<< Chapter < Page Chapter >> Page >
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

Get Jobilize Job Search Mobile App in your pocket Now!

Get it on Google Play Download on the App Store Now




Source:  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.

Notification Switch

Would you like to follow the 'Fast fourier transforms' conversation and receive update notifications?

Ask