2012-05-22 13:09:48 -04:00
|
|
|
subroutine ffft(d,npts,isign,ireal)
|
|
|
|
|
|
|
|
|
|
C Fourier transform of length npts=2**k, performed in place.
|
|
|
|
|
C Input data in array d, treated as complex if ireal=0, and as real if ireal=1.
|
|
|
|
|
C In either case the transform values are returned in array d, treated as
|
|
|
|
|
C complex. The DC term is d(1), and d(npts/2+1) is the term at the Nyquist
|
|
|
|
|
C frequency. The basic algorithm is the same as Norm Brenner's FOUR1, and
|
|
|
|
|
C uses radix-2 transforms.
|
|
|
|
|
|
|
|
|
|
C J. H. Taylor, Princeton University.
|
|
|
|
|
|
|
|
|
|
complex d(npts),t,w,wstep,tt,uu
|
|
|
|
|
data pi/3.14159265359/
|
|
|
|
|
|
|
|
|
|
C Shuffle the data to bit-reversed order.
|
|
|
|
|
|
|
|
|
|
imax=npts/(ireal+1)
|
|
|
|
|
irev=1
|
|
|
|
|
do 5 i=1,imax
|
|
|
|
|
if(i.ge.irev) go to 2
|
|
|
|
|
t=d(i)
|
|
|
|
|
d(i)=d(irev)
|
|
|
|
|
d(irev)=t
|
|
|
|
|
2 mmax=imax/2
|
|
|
|
|
3 if(irev.le.mmax) go to 5
|
|
|
|
|
irev=irev-mmax
|
|
|
|
|
mmax=mmax/2
|
|
|
|
|
if(mmax.ge.1) go to 3
|
|
|
|
|
5 irev=irev+mmax
|
|
|
|
|
|
|
|
|
|
C The radix-2 transform begins here.
|
|
|
|
|
|
|
|
|
|
api=isign*pi/2.
|
|
|
|
|
mmax=1
|
|
|
|
|
6 istep=2*mmax
|
|
|
|
|
wstep=cmplx(-2.*sin(api/mmax)**2,sin(2.*api/mmax))
|
|
|
|
|
w=1.
|
|
|
|
|
do 9 m=1,mmax
|
|
|
|
|
|
|
|
|
|
C This in the inner-most loop -- optimization here is important!
|
|
|
|
|
do 8 i=m,imax,istep
|
|
|
|
|
t=w*d(i+mmax)
|
|
|
|
|
d(i+mmax)=d(i)-t
|
|
|
|
|
8 d(i)=d(i)+t
|
|
|
|
|
|
|
|
|
|
9 w=w*(1.+wstep)
|
|
|
|
|
mmax=istep
|
|
|
|
|
if(mmax.lt.imax) go to 6
|
|
|
|
|
|
|
|
|
|
if(ireal.eq.0) return
|
|
|
|
|
|
|
|
|
|
C Now complete the last stage of a doubled-up real transform.
|
|
|
|
|
|
|
|
|
|
jmax=imax/2 + 1
|
|
|
|
|
wstep=cmplx(-2.*sin(isign*pi/npts)**2,sin(isign*pi/imax))
|
|
|
|
|
w=1.0
|
|
|
|
|
d(imax+1)=d(1)
|
|
|
|
|
|
|
|
|
|
do 10 j=1,jmax
|
|
|
|
|
uu=cmplx(real(d(j))+real(d(2+imax-j)),aimag(d(j)) -
|
|
|
|
|
+ aimag(d(2+imax-j)))
|
|
|
|
|
tt=w*cmplx(aimag(d(j))+aimag(d(2+imax-j)),-real(d(j)) +
|
|
|
|
|
+ real(d(2+imax-j)))
|
|
|
|
|
d(j)=uu+tt
|
|
|
|
|
d(2+imax-j)=conjg(uu-tt)
|
|
|
|
|
10 w=w*(1.+wstep)
|
|
|
|
|
|
|
|
|
|
return
|
|
|
|
|
end
|
