This page is optimized for mobile devices, if you would prefer the desktop version just
click here
Fortran programs for efficient DFT, Cooley-Tukey, and Prime Factor Algorithm FFTs.
Goertzel algorithm
A FORTRAN implementation of the first-order Goertzel algorithm with in-order input as given in ( [link] ) and [link] is given below.
C----------------------------------------------
C GOERTZEL'S DFT ALGORITHMC First order, input inorder
C C. S. BURRUS, SEPT 1983C---------------------------------------------
SUBROUTINE DFT(X,Y,A,B,N)REAL X(260), Y(260), A(260), B(260)
Q = 6.283185307179586/NDO 20 J=1, N
C = COS(Q*(J-1))S = SIN(Q*(J-1))
AT = X(1)BT = Y(1)
DO 30 I = 2, NT = C*AT - S*BT + X(I)
BT = C*BT + S*AT + Y(I)AT = T
30 CONTINUEA(J) = C*AT - S*BT
B(J) = C*BT + S*AT20 CONTINUE
RETURNEND
First Order Goertzel Algorithm
Second order goertzel algorithm
Below is the program for a second order Goertzel algorithm.
C----------------------------------------------
C GOERTZEL'S DFT ALGORITHMC Second order, input inorder
C C. S. BURRUS, SEPT 1983C---------------------------------------------
SUBROUTINE DFT(X,Y,A,B,N)REAL X(260), Y(260), A(260), B(260)
CQ = 6.283185307179586/N
DO 20 J = 1, NC = COS(Q*(J-1))
S = SIN(Q*(J-1))CC = 2*C
A2 = 0B2 = 0
A1 = X(1)B1 = Y(1)
DO 30 I = 2, NT = A1
A1 = CC*A1 - A2 + X(I)A2 = T
T = B1B1 = CC*B1 - B2 + Y(I)
B2 = T30 CONTINUE
A(J) = C*A1 - A2 - S*B1B(J) = C*B1 - B2 + S*A1
20 CONTINUEC
RETURNEND
Second Order Goertzel Algorithm
Second order goertzel algorithm 2
Second order Goertzel algorithm that calculates two outputs at a time.
C-------------------------------------------------------
C GOERTZEL'S DFT ALGORITHM, Second orderC Input inorder, output by twos; C.S. Burrus, SEPT 1991
C-------------------------------------------------------SUBROUTINE DFT(X,Y,A,B,N)
REAL X(260), Y(260), A(260), B(260)Q = 6.283185307179586/N
DO 20 J = 1, N/2 + 1C = COS(Q*(J-1))
S = SIN(Q*(J-1))CC = 2*C
A2 = 0B2 = 0
A1 = X(1)B1 = Y(1)
DO 30 I = 2, NT = A1
A1 = CC*A1 - A2 + X(I)A2 = T
T = B1B1 = CC*B1 - B2 + Y(I)
B2 = T30 CONTINUE
A2 = C*A1 - A2T = S*B1
A(J) = A2 - TA(N-J+2) = A2 + T
B2 = C*B1 - B2T = S*A1
B(J) = B2 + TB(N-J+2) = B2 - T
20 CONTINUERETURN
ENDFigure. Second Order Goertzel Calculating Two Outputs at a Time
Basic qft algorithm
A FORTRAN implementation of the basic QFT algorithm is given below to show how the theory is implemented. The program is written for clarity, not tominimize the number of floating point operations.
C
SUBROUTINE QDFT(X,Y,XX,YY,NN)REAL X(0:260),Y(0:260),XX(0:260),YY(0:260)
CN1 = NN - 1
N2 = N1/2N21 = NN/2
Q = 6.283185308/NNDO 2 K = 0, N21
SSX = X(0)SSY = Y(0)
SDX = 0SDY = 0
IF (MOD(NN,2).EQ.0) THENSSX = SSX + COS(3.1426*K)*X(N21)
SSY = SSY + COS(3.1426*K)*Y(N21)ENDIF
DO 3 N = 1, N2SSX = SSX + (X(N) + X(NN-N))*COS(Q*N*K)
SSY = SSY + (Y(N) + Y(NN-N))*COS(Q*N*K)SDX = SDX + (X(N) - X(NN-N))*SIN(Q*N*K)
SDY = SDY + (Y(N) - Y(NN-N))*SIN(Q*N*K)3 CONTINUE
XX(K) = SSX + SDYYY(K) = SSY - SDX
XX(NN-K) = SSX - SDYYY(NN-K) = SSY + SDX
2 CONTINUERETURN
END
Simple QFT Fortran Program
Read also:
OpenStax, Fast fourier transforms. OpenStax CNX. Nov 18, 2012 Download for free at http://cnx.org/content/col10550/1.22
Google Play and the Google Play logo are trademarks of Google Inc.