| << Chapter < Page | Chapter >> Page > |
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
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 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
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
Notification Switch
Would you like to follow the 'Fast fourier transforms' conversation and receive update notifications?
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|