|
C Here is an FFT routine I have used for some time
c***************************************************
C Fast Fourier Transform routine from Stearns p266
c "Digital Signal Analysis" Hayden Books - Samuel D. Stearns
c ISBN 0-8104-5828-4 1975
C Typed in by D. Drake Dec 26, 1986
c****************************************************
c fr and fi are the real and imaginary input arrays. N is the
c length of each array, ex: 1024. K is the log base 2 of N,
c ex: 10. fr and fi have the FFT on return, again in real and
c imaginary format.
c****************************************************
Subroutine FFT(fr,fi,k)
real fr(1),fi(1)
n=2**k
mr=0
nn=n-1
pi=3.1415926535
DO m=1,nn !bit reversal
l=n
1 l=l/2
if (mr+l.gt.nn) goto 1
mr=mod(mr,l) + l
if (mr.le.m) goto 2
tr=fr(m+1)
fr(m+1)=fr(mr+1)
fr(mr+1)=tr
ti=fi(m+1)
fi(m+1)=fi(mr+1)
fi(mr+1)=ti
2 END DO
l=1
3 if(l.ge.n) return !return to main routine
istep=2*l
el=float(l)
DO m=1,l
a=pi*float(1-m)/el
wr=cos(a)
wi=sin(a)
DO i=m,n,istep
j=i+l
tr=wr*fr(j)-wi*fi(j)
ti=wr*fi(j)+wi*fr(j)
fr(j)=fr(i)-tr
fi(j)=fi(i)-ti
fr(i)=fr(i)+tr
fi(i)=fi(i)+ti
END DO
END DO
l=istep
goto 3
END
|