C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- PROGRAM ZFFT2_PROGRAM !This example program performs a forward and reverse base-2 fast !Fourier transform (FFT) of a time series with uniform data !sampling and no missing data points. The subroutine invoked !(see code below) is "ZFFT2.F" and requires that the time series !length must be a positive integer power of 2 - i.e., 2, 4, 8, !16, 32, 64, ..., 1024, 2048, ..., 8192, ... ! !Time series lengths in data are seldom a power of 2, so a method !called "zero padding" is often used to extend the time series up !to the next power of 2 length by adding zero values after the !final data value if a base-2 FFT is used. One advantage of zero !padding is that frequency resolution is deemed better since !frequency resolution is inversely proportional to the number of !time series values. A disadvantage of zero padding is that you !are fitting your time series with a bunch of fictitious zeros to !sine and cosine harmonics, and if you add enough zeros you may !erroneously end up with precisely the power spectrum that you !wish for. The base-2 FFT is very fast and accurate as long as !the number of zeros padded do not exceed (conservatively) about !10-20 percent of the original time series length. ! !More extensive FFT algorithms invoke the fact that time series !length (an arbitrary natural number N) can be written as a !product of prime number factors with each prime number raised to !a positive integer value. In other words, every natural number !N greater than one can be written ! ! N = 2^k2 * 3^k3 * 5^k5 * 7^k7 * ... ! !These more advanced FFT algorithms have difficulty when even one !prime number factor is large. This can cause computation time !to increase by several orders of magnitude unless these routines !also use an internal "zero padding" or similar scheme. ! !The first step of FFT analysis is always to determine the prime !number factorization of the time series length - see the program !"PRIMEFAC.F" which provides prime number factorization of !arbitrary natural numbers N. For arbitrary N with more advanced !FFT routines, one can almost always repeat or delete the last !few time series values (creating a new time series length N) to !optimize computation speed without seriously jeopardizing the !derived Fourier coefficients. ! !Author: Dr. Jerry R. Ziemke PARAMETER(N=16) REAL X1(N),X2(N) !X1 below are the example time series values: DATA X1 /1., 7., 4., 3., 2., 9., 0., 2., 5., 6., 9., 2., & 5., 8., 0., 1./ ! WRITE(*,*)' ' WRITE(*,*)'Beginning time series values:' WRITE(*,*) X1 WRITE(*,*)' ' WRITE(*,*)'Forward FFT (Freq. index, cosine and sine coeff''s):' WRITE(*,*)' P C(P) S(P)' WRITE(*,*)'--------------------------' ISIGN=1 CALL ZFFT2(ISIGN,N,X1,X2) DO I=1,N/2+1 WRITE(*,'(I5, 2F10.5)') I, X1(I), X2(I) ENDDO WRITE(*,*)' ' WRITE(*,*)'Reverse FFT (should be original time series values):' ISIGN=-1 CALL ZFFT2(ISIGN,N,X1,X2) WRITE(*,*) X1 STOP END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- SUBROUTINE ZFFT2(ISIGN,N,X1,X2) C This subroutine is a base-2 Fast Fourier Transform (FFT). C C Calling arguments (all four are required): C C (1) ISIGN = 1 (Derive Fourier cosine and sine coefficients from C the original timeseries of length N): C C When making a call to ZFFT2 with ISIGN = 1, X1(1, 2, ..., N) C is the original time series, and X2 is a work vector of the C same size as X1 and need not have any assigned values. The value C of N MUST be a power of 2 - i.e., 2, 4, 8, 16, 32, 64, ... After C calling ZFFT2 with ISIGN = 1, X1 and X2 will be the Fourier cosine C and sine coefficients (indices one through N/2+1 respectively in C Fortran language). The mathematics is: C C n/2 C x(t;t=0,1,..,n-1) = SUM [c(p)*cos(w*t)+s(p)*sin(w*t)] C p=0 C C where x(t) is the original time series of length n, c and s are C the Fourier cosine and sine coefficients, and w is the circular C frequency = 2*PI*p/n (PI = 3.141592653589793...) C C The Fourier cosine and sine coefficients are identified from C either least-squares analysis or orthogonal relationships C between harmonic sine and cosine functions, and are: C C 1 n-1 1 n-1 t+1 C c(0) = - * SUM x(t) c(n/2) = - * SUM (-1) * x(t) C n t=0 n t=0 C C s(0) = s(n/2) = 0 C C 2 n-1 C c(p;p=1,2,..,n/2-1) = - * SUM x(t)*cos(w*t) C n t=0 C s(p; " ) = " " " sin( " ) C C (1) ISIGN = -1 (Derive time series from Fourier cosine and sine C coefficients): C C When calling ZFFT2 with ISIGN = -1, X1(1, 2, ..., N/2+1) and C X2(1, 2, ..., N/2+1) are the Fourier cosine and sine coefficients, C respectively. After the call, X1(1, 2, ..., N) will be the C reconstructed time series. C C (2) N: Time series length, as a power of 2. C C (3) and (4) - X1(N) and X2(N): C C With ISIGN = 1, X1 is the original time series, and values for X2 C need not be defined; X1 and X2 are the respective Fourier cosine C and sine coefficients following this forward FFT call. With C ISIGN = -1 (reverse FFT), X1 amd X2 are the Fourier cosine and C sine coefficients, respectively; X1 after this reverse FFT call C is the reconstructed timeseries. C C Author: Dr. Jerry R. Ziemke REAL X1(N),X2(N),MEAN IF (ISIGN.EQ.1) THEN MEAN=0. DO I=1,N MEAN=MEAN+X1(I) ENDDO MEAN=MEAN/N DO I=1,N X1(I)=X1(I)-MEAN X2(I)=0. ENDDO CALL UTF(ISIGN,N,X1,X2) DO I=2,N/2 X1(I)=2.0*X1(I) X2(I)=-2.0*X2(I) ENDDO X1(1)=X1(1)+MEAN ENDIF IF (ISIGN.EQ.-1) THEN MEAN=X1(1) X1(1)=0. DO I=2,N/2 X1(I)=0.5*X1(I) X2(I)=-0.5*X2(I) ENDDO DO I=2,N/2 X1(N/2+I)=X1(N/2+2-I) X2(N/2+I)=-X2(N/2+2-I) ENDDO X2(1)=0. X2(N/2+1)=0. CALL UTF(ISIGN,N,X1,X2) DO I=1,N X1(I)=X1(I)+MEAN ENDDO ENDIF RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- SUBROUTINE UTF(ISIGN,N,X1,X2) PARAMETER(PI=3.141592653589793) REAL X1(N),X2(N),CS(2),XA(2) INTEGER MS(20) NPOW=IFIX(ALOG(FLOAT(N))/ALOG(2.)+0.5) DT=1./N ZZ=(2.*PI*ISIGN)/N IF (ISIGN.GT.0) DELTA=DT IF (ISIGN.LT.0) DELTA=1. MS(1)=N/2 DO I=2,NPOW MS(I)=MS(I-1)/2 ENDDO NN=N MM=2 DO 60 LA=1,NPOW NN=NN/2 NW=0 IMAX=MM/2 DO 50 I0=1,IMAX I=1+2*(I0-1) II=NN*I W=NW*ZZ CS(1)=COS(W) CS(2)=SIN(W) DO 30 J=1,NN II=II+1 IJ=II-NN XA(1)=CS(1)*X1(II)-CS(2)*X2(II) XA(2)=CS(1)*X2(II)+CS(2)*X1(II) X1(II)=X1(IJ)-XA(1) X2(II)=X2(IJ)-XA(2) X1(IJ)=X1(IJ)+XA(1) X2(IJ)=X2(IJ)+XA(2) 30 CONTINUE DO 40 LOC=2,NPOW LL=NW-MS(LOC) IF (LL.LT.0) THEN NW=MS(LOC)+NW GOTO 50 ENDIF IF (LL.EQ.0) THEN NW=MS(LOC+1) GOTO 50 ENDIF NW=LL 40 CONTINUE NW=MS(LOC)+NW 50 CONTINUE MM=MM*2 60 CONTINUE NW=0 DO 80 I=1,N N1=NW+1 H1=X1(N1) H2=X2(N1) IF ((N1-I).GE.0) THEN IF ((N1-I).GT.0) THEN X1(N1)=X1(I)*DELTA X2(N1)=X2(I)*DELTA X1(I)=H1*DELTA X2(I)=H2*DELTA ENDIF IF ((N1-I).EQ.0) THEN X1(I)=H1*DELTA X2(I)=H2*DELTA ENDIF ENDIF DO 70 LOC=1,NPOW LL=NW-MS(LOC) IF (LL.LT.0) THEN NW=MS(LOC)+NW GOTO 80 ENDIF IF (LL.EQ.0) THEN NW=MS(LOC+1) GOTO 80 ENDIF IF (LL.GT.0) THEN NW=LL ENDIF 70 CONTINUE NW=MS(LOC+1) 80 CONTINUE RETURN END