I hope this isn't blocking up the net.
Hello, can anyone help me with this fft code. Original by Bergland & Dolan
published in Programs for Digital Signal Processing IEEE Press 1979
The main program tests fft routines, the problem is with the subroutine code
fft842 which complies but doesnot produce the right answers.
I am using a xlf complier on a RS6000. The original code was written to comply
with ANSI
fortran but I didn't have to translate too much. The subroutine call
statements have changed slightly from the original
--------------------------code begins below this line-------------------------
c
c main program foursubt
c
complex w,c(32),d,e,b(32),qb(32),a
open(unit=10,file='spectra.dat',status='unknown')
c
c sequence length n=2**mu
c method is to compute dft of known function; a**i, i=0,1,...,n-1
c output of program is largest difference between dft computed two ways
c a max difference large compared to the accuracy of the computer
c would probably indicate a program error
c
mu=5
c
ioutd=10
nn=2**mu
tpi=8.0*atan(1.0)
tpion=tpi/float(nn)
w=cmplx(cos(tpion),-sin(tpion))
c
c generate a**k as test function
c
a=(0.9,0.3)
b(1)=(1.0,0.0)
qb(1)=b(1)
do 10 k=2,nn
b(k)=a**(k-1)
qb(k)=b(k)
10 continue
c
c print complex input sequence
c
write(ioutd,9999)
9999 format(/24h complex input sequence )
write(ioutd,9998) (i,qb(i),i=1,nn)
9998 format(2(2x,1h(, i3,1h),2e14.6))
c
c b(1) contains a**0; b(k) contains a**k-1; etc
c
c compute dft of b in closed form
c
d=(1.0,0.0)-a**nn
do 20 k=1,nn
e=(1.0,0.0)-a*w**(k-1)
c(k)=d/e
20 continue
c
c dft of b is (1-a**nn)/(1-aw**k)
c
c now compute dft of b using fft842
c
call fft842(0,nn,x,y,b)
c
c print fft842 dft and theoretical dft
c
write(ioutd,9997)
9997 format(/12h fft842 dft )
write(ioutd,9998) (i,b(i),i=1,nn)
write(ioutd,9996)
9996 format(/17h theoretical dft )
write(ioutd,9998) (i,c(i),i=1,nn)
c
c find max difference between b and c
c
dd=0.0
do 30 i=1,nn
de=cabs(c(i)-b(i))
if (dd.gt.de) go to 30
dd=de
30 continue
c
c compute inverse dft of c
c
call fft842(1,nn,x,y,c)
c
c print inverse dft
c
write(ioutd,9995)
9995 format(/20h fft842 inverse dft )
write(ioutd,9998) (i,c(i),i=1,nn)
c
c compute max diff between input and inverse of theor dft
c
gg=0.0
do 40 i=1,nn
gg1=cabs(qb(i)-c(i))
if (gg1.le.gg) go to 40
gg=gg1
40 continue
c
c print maximum difference
c
write(ioutd,9994) dd
write(ioutd,9993) gg
write(6,9994) dd
write(6,9993) gg
9994 format(/42h max diff between theor and fft842 dft is , e12.3)
9993 format(/51h max diff between original data and inverse dft is ,
+ e12.3)
stop
end
cc
c subroutine fft842
c fft for n=2**m
c complex input
c
subroutine fft842(in,n,x,y,z)
c
c this subroutine replaces the vector z=x+iy by it's finite discrete
c complex fourier transform if in=0. the inverse transformation is
c calculated for in=1. it performs as many base 8 iterations as
c possible and then finishes with a base 4 iteration or
c base 2 iteration if needed
c
c the subroutine is called as subroutine fft842(in,n,x,y)
c the integer n (a power of 2), the n real location array x and
c the n real location y must be supplied to the subroutine
c
complex z(n)
dimension x(n),y(n),l(15)
common /con2/pi2,p7
equivalence (l15,l(1)),(l14,l(2)),(l13,l(3)),(l12,l(4)),
& (l11,l(5)),(l10,l(6)),(l9,l(7)),(l8,l(8)),(l7,l(9)),
& (l6,l(10)),(l5,l(11)),(l4,l(12)),(l3,l(13)),(l2,l(14)),
& (l1,l(15))
c
c iw is a machine dependent write device number
c
c iw=i1mach(2)
iw=6
c
pi2=8.0*atan(1.0)
p7=1.0/sqrt(2.0)
do 10 i=1,15
m=i
nt=2**i
if(n.eq.nt)go to 20
10 continue
write(iw,9999)
9999 format(35h n is not a power of two for fft842)
stop
20 n2pow=m
nthpo=n
fn=nthpo
do 5 i=1,n
x(i)=real(z(i))
y(i)=imag(z(i))
5 continue
if(in.eq.1)go to 40
do 30 i=1,nthpo
y(i)=-y(i)
30 continue
40 n8pow=n2pow/3
if(n8pow.eq.0)go to 60
c
c radix 8 passes, if any
c
do 50 ipass=1,n8pow
nxtlt=2**(n2pow-3*ipass)
lengt=8*nxtlt
call r8tx(nxtlt,nthpo,lengt,x,y)
50 continue
c
c is there a four factor left
c
60 if(n2pow-3*n8pow-1)90,70,80
c
c go through the base 2 iteration
c
70 call r2tx(nthpo,x,y)
go to 90
c
c go through the base 4 iteration
c
80 call r4tx(nthpo,x,y)
c
90 do 110 j=1,15
l(j)=1
if(j-n2pow)100,100,110
100 l(j)=2**(n2pow+1-j)
110 continue
ij=1
do 130 j1=1,l1
do 130 j2=j1,l2,l1
do 130 j3=j2,l3,l2
do 130 j4=j3,l4,l3
do 130 j5=j4,l5,l4
do 130 j6=j5,l6,l5
do 130 j7=j6,l7,l6
do 130 j8=j7,l8,l7
do 130 j9=j8,l9,l8
do 130 j10=j9,l10,l9
do 130 j11=j10,l11,l10
do 130 j12=j11,l12,l11
do 130 j13=j12,l13,l12
do 130 j14=j13,l14,l13
do 130 ji=j14,l15,l14
if(ij-ji)120,130,130
120 r=x(ij)
x(ij)=x(ji)
x(ji)=r
fi=y(ij)
y(ij)=y(ji)
y(ji)=fi
130 ij=ij+1
if(in.eq.1)go to 150
do 140 i=1,nthpo
y(i)=-y(i)
140 continue
go to 170
150 do 160 i=1,nthpo
x(i)=x(i)/fn
y(i)=y(i)/fn
160 continue
170 do 15 i=1,nthpo
z(i)=(x(i),y(i))
15 continue
return
end
c
c subroutine r2tx
c radix 2 iteration subroutine
c
subroutine r2tx(nthpo,x,y)
dimension x(nthpo),y(nthpo)
do 10 k0=1,nthpo,2
k1=k0+1
r1=x(k0)+x(k1)
x(k1)=x(k0)-x(k1)
x(k0)=r1
fi1=y(k0)+y(k1)
y(k1)=y(k0)-y(k1)
y(k0)=fi1
10 continue
return
end
c
c subroutine r4tx
c radix 4 iteration
c
subroutine r4tx(nthpo,x,y)
dimension x(nthpo),y(nthpo)
do 10 k0=1,nthpo,4
k1=k0+1
k2=k1+1
k3=k2+1
r1=x(k0)+x(k2)
r2=x(k0)-x(k2)
r3=x(k1)+x(k3)
r4=x(k1)-x(k3)
fi1=y(k0)+y(k2)
fi2=y(k0)-y(k2)
fi3=y(k1)+y(k3)
fi4=y(k1)-y(k3)
x(k0)=r1+r3
y(k0)=fi1+fi3
x(k1)=r1-r3
y(k1)=fi1-fi3
x(k2)=r2-fi4
y(k2)=fi2+r4
x(k3)=r2+fi4
y(k3)=fi2-r4
10 continue
return
end
c
c subroutine r8tx
c radix 8 iteration subroutine
c
subroutine r8tx(nxtlt,nthpo,lengt,x,y)
dimension x(nthpo),y(nthpo)
common /con2/pi2,p7
c
scale=pi2/float(lengt)
do 30 j=1,nxtlt
arg=float(j-1)*scale
c1=cos(arg)
s1=sin(arg)
c2=c1**2-s1**2
s2=c1*s1+c1*s1
c3=c1*c2-s1*s2
s3=c2*s1+s2*c1
c4=c2**2-s2**2
s4=c2*s2+c2*s2
c5=c2*c3-s2*s3
s5=c3*s2+s3*c2
c6=c3**2-s3**2
s6=c3*s3+c3*s3
c7=c3*c4-s3*s4
s7=c4*s3+s4*c3
do 20 k0=j,nthpo,lengt
k1=k0+nxtlt
k2=k1+nxtlt
k3=k2+nxtlt
k4=k3+nxtlt
k5=k4+nxtlt
k6=k5+nxtlt
k7=k6+nxtlt
ar0=x(k0)+x(k4)
ar1=x(k1)+x(k5)
ar2=x(k2)+x(k6)
ar3=x(k3)+x(k7)
ar4=x(k0)-x(k4)
ar5=x(k1)-x(k5)
ar6=x(k2)-x(k6)
ar7=x(k3)-x(k7)
ai0=y(k0)+y(k4)
ai1=y(k1)+y(k5)
ai2=y(k2)+y(k6)
ai3=y(k3)+y(k7)
ai4=y(k0)-y(k4)
ai5=y(k1)-y(k5)
ai6=y(k2)-y(k6)
ai7=y(k3)-y(k7)
br0=ar0+ar2
br1=ar1+ar3
br2=ar0-ar2
br3=ar1-ar3
br4=ar4-ai6
br5=ar5-ai7
br6=ar4+ai6
br7=ar5+ai7
bi0=ai0+ai2
bi1=ai1+ai3
bi2=ai0-ai2
bi3=ai1-ai3
bi4=ai4+ar6
bi5=ai5+ar7
bi6=ai4-ar6
bi7=ai5-ar7
x(k0)=br0+br1
y(k0)=bi0+bi1
if(j.le.1)go to 10
x(k1)=c4*(br0-br1)-s4*(bi0-bi1)
y(k1)=c4*(bi0-bi1)+s4*(br0-br1)
x(k2)=c2*(br2-bi3)-s2*(bi2+br3)
y(k2)=c2*(bi2+br3)+s2*(br2-bi3)
x(k3)=c6*(br2+bi3)-s6*(bi2-br3)
y(k3)=c6*(bi2-br3)+s6*(br2+bi3)
tr=p7*(br5-bi5)
ti=p7*(br5+bi5)
x(k4)=c1*(br4+tr)-s1*(bi4+ti)
y(k4)=c1*(bi4+ti)+s1*(br4+tr)
x(k5)=c5*(br4-tr)-s5*(bi4-ti)
y(k5)=c5*(bi4-ti)+s5*(br4-tr)
tr=-p7*(br7+bi7)
ti=p7*(br7-bi7)
x(k6)=c3*(br6+tr)-s3*(bi6-ti)
y(k6)=c3*(bi6-ti)+s3*(br6+tr)
x(k7)=c7*(br6-tr)-s7*(bi6-ti)
y(k7)=c7*(bi6-ti)+s7*(br6-tr)
go to 20
10 x(k1)=br0-br1
y(k1)=bi0-bi1
x(k2)=br2-bi3
y(k2)=bi2+br3
x(k3)=br2+bi3
y(k3)=bi2-br3
tr=p7*(br5-bi5)
ti=p7*(br5+bi5)
x(k4)=br4+tr
y(k4)=bi4+ti
x(k5)=br4-tr
y(k5)=bi4-ti
tr=-p7*(br7+bi7)
ti=p7*(br7-bi7)
x(k6)=br6+tr
y(k6)=bi6+ti
x(k7)=br6-tr
y(k7)=bi6-ti
20 continue
30 continue
return
end
fft842 subroutine by Bergland & Dolan IEEE 1979