This page is optimized for mobile devices, if you would prefer the desktop version just click here

0.17 Appendix 3: fft computer programs

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

<< Chapter < Page Page > Chapter >>

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.