fft842 subroutine by Bergland & Dolan IEEE 1979 
Author Message
 fft842 subroutine by Bergland & Dolan IEEE 1979

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



Sat, 18 Jan 1997 13:01:44 GMT  
 fft842 subroutine by Bergland & Dolan IEEE 1979
The routine will work properly if the arrays x and y are dimensioned in main

eg DIMENSION X(32),Y(32)

I also spotted errors in the subroutine R8TX

        X(K6)=C3*(BR6+TR)-S3*(BI6+TI)
                                 ^this is a plus, not a minus
        Y(K6)=C3*(BI6+TI)+S3*(BR6+TR)
                     ^this is a plus, not a minus

These changes give a perfectly functional fft routine, it took about 4 seconds
to do a 32768 point fft and 25 seconds 131072 points. On a RS6000 though.



Mon, 20 Jan 1997 09:57:45 GMT  
 
 [ 3 post ] 

 Relevant Pages 

1. IEEE Standard Subroutines for CAMAC (ANSI/IEEE Std 758-1979)

2. array-conversion: $arr: [1979]=>1979, [2000]=>2000, [2002]=>2002 zu $arr2: [0]=>1979, [1]=>2000, [2]=>2002 ?

3. Charles Moore (1979) on the Future.

4. Some Q&A with Kahan on IEEE 754 (long)

5. CFP: Virtual Environments & IEEE Yuforic - Final CFP

6. IEEE P1481 Delay & Power Calculation Standard Update

7. DAC Tutorial: IEEE Delay & Power Calculation System

8. IEEE 1284 &Fast Ethernet MAC

9. IEEE 1284 & Fast Ethernet Mac

10. Ada/IEEE binding (IEEE 754)

11. DAC Tutorial: IEEE Delay & Power Calculation System

12. IEEE P1481 Delay & Power Calculation Standard Update

 

 
Powered by phpBB® Forum Software