Initial commit

4 files changed, 568 insertions, 0 deletions

diff --git a/math/ieee/chap1/test/test12.f b/math/ieee/chap1/test/test12.f
new file mode 100644
index 00000000..a830ab2c
--- /dev/null
+++ b/math/ieee/chap1/test/test12.f
@@ -0,0 +1,90 @@
+c
+c-----------------------------------------------------------------------
+c main program: fastmain  -  fast fourier transforms
+c authors:      g. d. bergland and m. t. dolan
+c               bell laboratories, murray hill, new jersey 07974
+c
+c input:        the program calls on a random number
+c               generator for input and checks dft and
+c               idft with a 32-point sequence
+c-----------------------------------------------------------------------
+c
+      dimension x(32), y(32), b(34)
+c
+c generate random numbers and store array in b so
+c the same sequence can be used in all tests.
+c note that b is dimensioned to size n+2.
+c
+c iw is a machine dependent write device number
+c
+      iw = i1mach(2)
+c
+      do 10 i=1,32
+        x(i) = uni(0)
+        b(i) = x(i)
+  10  continue
+      m = 5
+      n = 2**m
+      np1 = n + 1
+      np2 = n + 2
+      knt = 1
+c
+c test fast-fsst then ffa-ffs
+c
+      write (iw,9999)
+  20  write (iw,9998) (b(i),i=1,n)
+      if (knt.eq.1) call fast(b, n)
+      if (knt.eq.2) call ffa(b, n)
+      write (iw,9997) (b(i),i=1,np1,2)
+      write (iw,9996) (b(i),i=2,np2,2)
+      if (knt.eq.1) call fsst(b, n)
+      if (knt.eq.2) call ffs(b, n)
+      write (iw,9995) (b(i),i=1,n)
+      knt = knt + 1
+      if (knt.eq.3) go to 40
+c
+      write (iw,9994)
+      do 30 i=1,n
+        b(i) = x(i)
+  30  continue
+      go to 20
+c
+c test fft842 with real input then complex
+c
+  40  write (iw,9993)
+      do 50 i=1,n
+        b(i) = x(i)
+        y(i) = 0.
+  50  continue
+  60  write (iw,9992) (b(i),i=1,n)
+      write (iw,9991) (y(i),i=1,n)
+      call fft842(0, n, b, y)
+      write (iw,9997) (b(i),i=1,n)
+      write (iw,9996) (y(i),i=1,n)
+      call fft842(1, n, b, y)
+      write (iw,9990) (b(i),i=1,n)
+      write (iw,9989) (y(i),i=1,n)
+      knt = knt + 1
+      if (knt.eq.5) go to 80
+c
+      write (iw,9988)
+      do 70 i=1,n
+        b(i) = x(i)
+        y(i) = uni(0)
+  70  continue
+      go to 60
+c
+9999  format (19h1test fast and fsst)
+9998  format (20h0real input sequence/(4e17.8))
+9997  format (29h0real components of transform/(4e17.8))
+9996  format (29h0imag components of transform/(4e17.8))
+9995  format (23h0real inverse transform/(4e17.8))
+9994  format (17h1test ffa and ffs)
+9993  format (37h1test fft842 with real input sequence/(4e17.8))
+9992  format (34h0real components of input sequence/(4e17.8))
+9991  format (34h0imag components of input sequence/(4e17.8))
+9990  format (37h0real components of inverse transform/(4e17.8))
+9989  format (37h0imag components of inverse transform/(4e17.8))
+9988  format (40h1test fft842 with complex input sequence)
+  80  stop
+      end
diff --git a/math/ieee/chap1/test/test13.f b/math/ieee/chap1/test/test13.f
new file mode 100644
index 00000000..d7157069
--- /dev/null
+++ b/math/ieee/chap1/test/test13.f
@@ -0,0 +1,260 @@
+c
+c-----------------------------------------------------------------------
+c main program: test program for fft subroutines
+c author:      l r rabiner
+c      bell laboratories, murray hill, new jersey, 07974
+c input:       randomly chosen sequences to test fft subroutines
+c      for sequences with special properties
+c      n is the fft length (n must be a power of 2)
+c          2<= n <= 4096
+c-----------------------------------------------------------------------
+c
+      dimension x(4098), y(4098)
+c
+c  define i/0 device codes
+c  input: input to this program is user-interactive
+c        that is - a question is written on the user
+c        terminal (iout1) ad the user types in the answer.
+c
+c output: all output is written on the standard
+c         output unit (iout2).
+c
+      ind = i1mach(1)
+      iout1 = i1mach(4)
+      iout2 = i1mach(2)
+c
+c read in analysis size for fft
+c
+  10  write (iout1,9999)
+9999  format (30h fft size(2.le.n.le.4096)(i4)=)
+      read (ind,9998) n
+9998  format (i4)
+      if (n.eq.0) stop
+      do 20 i=1,12
+        itest = 2**i
+        if (n.eq.itest) go to 30
+  20  continue
+      write (iout1,9997)
+9997  format (45h n is not a power of 2 in the range 2 to 4096)
+      go to 10
+  30  write (iout2,9996) n
+9996  format (11h testing n=, i5, 17h random sequences)
+      write (iout2,9992)
+      np2 = n + 2
+      no2 = n/2
+      no2p1 = no2 + 1
+      no4 = n/4
+      no4p1 = no4 + 1
+c
+c create symmetrical sequence of size n
+c
+      do 40 i=2,no2
+        x(i) = uni(0) - 0.5
+        ind1 = np2 - i
+        x(ind1) = x(i)
+  40  continue
+      x(1) = uni(0) - 0.5
+      x(no2p1) = uni(0) - 0.5
+      do 50 i=1,no2p1
+        y(i) = x(i)
+  50  continue
+      write (iout2,9995)
+9995  format (28h original symmetric sequence)
+      write (iout2,9993) (x(i),i=1,n)
+      write (iout2,9992)
+c
+c compute true fft of n point sequence
+c
+      call fast(x, n)
+      write (iout2,9994) n
+9994  format (1h , i4, 32h point fft of symmetric sequence)
+      write (iout2,9993) (x(i),i=1,np2)
+9993  format (1h , 5e13.5)
+      write (iout2,9992)
+9992  format (1h /1h )
+c
+c use subroutine fftsym to obtain dft from no2 point fft
+c
+      do 60 i=1,no2p1
+        x(i) = y(i)
+  60  continue
+      call fftsym(x, n, y)
+      write (iout2,9991)
+9991  format (17h output of fftsym)
+      write (iout2,9993) (x(i),i=1,no2p1)
+      write (iout2,9992)
+c
+c use subroutine iftsym to obtain original sequence from no2 point dft
+c
+      call iftsym(x, n, y)
+      write (iout2,9990)
+9990  format (17h output of iftsym)
+      write (iout2,9993) (x(i),i=1,no2p1)
+      write (iout2,9992)
+c
+c create antisymmetric n point sequence
+c
+      do 70 i=2,no2
+        x(i) = uni(0) - 0.5
+        ind1 = np2 - i
+        x(ind1) = -x(i)
+  70  continue
+      x(1) = 0.
+      x(no2p1) = 0.
+      do 80 i=1,no2p1
+        y(i) = x(i)
+  80  continue
+      write (iout2,9989)
+9989  format (32h original antisymmetric sequence)
+      write (iout2,9993) (x(i),i=1,n)
+      write (iout2,9992)
+c
+c obtain n point dft of antisymmetric sequence
+c
+      call fast(x, n)
+      write (iout2,9988) n
+9988  format (1h , i4, 36h point fft of antisymmetric sequence)
+      write (iout2,9993) (x(i),i=1,np2)
+      write (iout2,9992)
+c
+c use subroutine fftasm to obtain dft from no2 point fft
+c
+      do 90 i=1,no2
+        x(i) = y(i)
+  90  continue
+      call fftasm(x, n, y)
+      write (iout2,9987)
+9987  format (17h output of fftasm)
+      write (iout2,9993) (x(i),i=1,no2p1)
+      write (iout2,9992)
+c
+c use subroutine iftasm to obtain original sequence from no2 point dft
+c
+      call iftasm(x, n, y)
+      write (iout2,9986)
+9986  format (17h output of iftasm)
+      write (iout2,9993) (x(i),i=1,no2)
+      write (iout2,9992)
+c
+c create sequence with only odd harmonics--begin in frequency domain
+c
+      do 100 i=1,np2,2
+        x(i) = 0.
+        x(i+1) = 0.
+        if (mod(i,4).eq.1) go to 100
+        x(i) = uni(0) - 0.5
+        x(i+1) = uni(0) - 0.5
+        if (n.eq.2) x(i+1) = 0.
+ 100  continue
+      write (iout2,9985) n
+9985  format (1h , i4, 35h point fft of odd harmonic sequence)
+      write (iout2,9993) (x(i),i=1,np2)
+      write (iout2,9992)
+c
+c transform back to time sequence
+c
+      call fsst(x, n)
+      write (iout2,9984)
+9984  format (31h original odd harmonic sequence)
+      write (iout2,9993) (x(i),i=1,n)
+      write (iout2,9992)
+c
+c use subroutine fftohm to obtain dft from no2 point fft
+c
+      call fftohm(x, n)
+      write (iout2,9983)
+9983  format (17h output of fftohm)
+      write (iout2,9993) (x(i),i=1,no2)
+      write (iout2,9992)
+c
+c use subroutine iftohm to obtain original sequence from no2 point dft
+c
+      call iftohm(x, n)
+      write (iout2,9982)
+9982  format (17h output of iftohm)
+      write (iout2,9993) (x(i),i=1,no2)
+      write (iout2,9992)
+c
+c create sequence with only real valued odd harmonics
+c
+      do 110 i=1,np2,2
+        x(i) = 0.
+        x(i+1) = 0.
+        if (mod(i,4).eq.1) go to 110
+        x(i) = uni(0) - 0.5
+ 110  continue
+      write (iout2,9981) n
+9981  format (1h , i4, 45h point fft of odd harmonic, symmetric sequenc,
+     *    1he)
+      write (iout2,9993) (x(i),i=1,np2)
+      write (iout2,9992)
+c
+c transform back to time sequence
+c
+      call fsst(x, n)
+      write (iout2,9980)
+9980  format (42h original odd harmonic, symmetric sequence)
+      write (iout2,9993) (x(i),i=1,n)
+      write (iout2,9992)
+c
+c use subroutine fftsoh to obtain dft from no4 point fft
+c
+      call fftsoh(x, n, y)
+      write (iout2,9979)
+9979  format (17h output of fftsoh)
+      write (iout2,9993) (x(i),i=1,no4)
+      write (iout2,9992)
+c
+c use subroutine iftsoh to obtain original sequence from no4 point dft
+c
+      call iftsoh(x, n, y)
+      write (iout2,9978)
+9978  format (17h output of iftsoh)
+      write (iout2,9993) (x(i),i=1,no4)
+      write (iout2,9992)
+c
+c create sequence with only imaginary valued odd harmonics--begin
+c in frequency domain
+c
+      do 120 i=1,np2,2
+        x(i) = 0.
+        x(i+1) = 0.
+        if (mod(i,4).eq.1) go to 120
+        x(i+1) = uni(0) - 0.5
+ 120  continue
+      write (iout2,9977) n
+9977  format (1h , i4, 41h point fft of odd harmonic, antisymmetric,
+     *    9h sequence)
+      write (iout2,9993) (x(i),i=1,np2)
+      write (iout2,9992)
+c
+c transform back to time sequence
+c
+      call fsst(x, n)
+      write (iout2,9976)
+9976  format (46h original odd harmonic, antisymmetric sequence)
+      write (iout2,9993) (x(i),i=1,n)
+      write (iout2,9992)
+c
+c use subroutine fftaoh to obtain dft from no4 point fft
+c
+      call fftaoh(x, n, y)
+      write (iout2,9975)
+9975  format (17h output of fftaoh)
+      write (iout2,9993) (x(i),i=1,no4)
+      write (iout2,9992)
+c
+c use subroutine iftaoh to obtain original sequence from n/4 point dft
+c
+      call iftaoh(x, n, y)
+      write (iout2,9974)
+9974  format (17h output of iftaoh)
+      write (iout2,9993) (x(i),i=1,no4p1)
+      write (iout2,9992)
+c
+c begin a new page
+c
+      write (iout2,9973)
+9973  format (1h1)
+      go to 10
+      end
diff --git a/math/ieee/chap1/test/test17.f b/math/ieee/chap1/test/test17.f
new file mode 100644
index 00000000..15a34788
--- /dev/null
+++ b/math/ieee/chap1/test/test17.f
@@ -0,0 +1,147 @@
+c
+c-----------------------------------------------------------------------
+c main program: test program to exercise the wfta subroutine
+c   the test waveform is a complex exponential a**i whose
+c   transform is known analytically to be (1 - a**n)/(1 - a*w**k).
+c
+c authors:
+c   james h. mcclellan     and     hamid nawab
+c   department of electrical engineering and computer science
+c   massachusetts of technology
+c   cambridge, mass.  02139
+c
+c inputs:
+c   n-- transform length. it must be formed as the product of
+c       relatively prime integers from the set:
+c           2,3,4,5,7,8,9,16
+c   invrs is the flag for forward or inverse transform.
+c           invrs = 1 yields inverse transform
+c           invrs .ne. 1 gives forward transform
+c   rad and phi are the magnitude and angle (as a fraction of
+c       2*pi/n) of the complex exponential test signal.
+c          suggestion: rad = 0.98, phi = 0.5.
+c-----------------------------------------------------------------------
+c
+      double precision pi2,pin,xn,xj,xt
+      dimension xr(1260),xi(1260)
+      complex cone,ca,can,cnum,cden
+c
+c   output will be punched
+c
+      iout=i1mach(3)
+      input=i1mach(1)
+      cone=cmplx(1.0,0.0)
+      pi2=8.0d0*datan(1.0d0)
+50    continue
+      read(input,130)n
+130   format(i5)
+      write(iout,150) n
+150   format(10h length = ,i5)
+      if(n.le.0 .or. n.gt.1260) stop
+c
+c   enter a 1 to perform the inverse
+c
+c     read(input,130) invrs
+      invrs = 0
+c
+c   enter magnitude and angle (in fraction of 2*pi/n)
+c   avoid multiples of n for the angle if the radius is
+c   close to one.  suggestion: rad = 0.98, phi = 0.5.
+c
+c     read(input,160) rad,phi
+      rad = 0.98
+      phi = 0.5
+160   format(2f15.10)
+      xn=float(n)
+      pin=phi
+      pin=pin*pi2/xn
+c
+c   generate z**j
+c
+      init=0
+      do 200 j=1,n
+      an=rad**(j-1)
+      xj=j-1
+      xj=xj*pin
+      xt=dcos(xj)
+      xr(j)=xt
+      xr(j)=xr(j)*an
+      xt=dsin(xj)
+      xi(j)=xt
+      xi(j)=xi(j)*an
+200   continue
+      can=cmplx(xr(n),xi(n))
+      ca=cmplx(xr(2),xi(2))
+      can=can*ca
+c
+c   print first 50 values of input sequence
+c
+      max=50
+      if(n.lt.50)max=n
+      write(iout,300)(j,xr(j),xi(j),j=1,max)
+c
+c   call the winograd fourier transform algorithm
+c
+      call wfta(xr,xi,n,invrs,init,ierr)
+c
+c   check for error return
+c
+      if(ierr.lt.0) write(iout,250) ierr
+250   format(1x,5herror,i5)
+      if(ierr.lt.0) go to 50
+c
+c   print first 50 values of the transformed sequence
+c
+      write(iout,300)(j,xr(j),xi(j),j=1,max)
+300   format(1x,3hj =,i3,6hreal =,e20.12,6himag =,e20.12)
+c
+c   calculate absolute and relative deviations
+c
+      devabs=0.0
+      devrel=0.0
+      cnum=cone-can
+      pin=pi2/xn
+      do 350 j=1,n
+      xj=j-1
+      xj=-xj*pin
+      if(invrs.eq.1) xj=-xj
+      tr=dcos(xj)
+      ti=dsin(xj)
+      can=cmplx(tr,ti)
+      cden=cone-ca*can
+      cden=cnum/cden
+c
+c   true value of the transform (1. - a**n)/(1. - a*w**k),
+c   where a = rad*exp(j*phi*(2*pi/n)), w = exp(-j*2*pi/n).
+c   for the inverse transform the complex exponential w
+c   is conjugated.
+c
+      tr=real(cden)
+      ti=aimag(cden)
+      if(invrs.ne.1) go to 330
+c
+c   scale inverse transform by 1/n
+c
+      tr=tr/float(n)
+      ti=ti/float(n)
+330   tr=xr(j)-tr
+      ti=xi(j)-ti
+      devabs=sqrt(tr*tr+ti*ti)
+      xmag=sqrt(xr(j)*xr(j)+xi(j)*xi(j))
+      devrel=100.0*devabs/xmag
+      if(devabs.le.devmx1)go to 340
+      devmx1=devabs
+      labs=j-1
+340   if(devrel.le.devmx2)go to 350
+      devmx2=devrel
+      lrel=j-1
+350   continue
+c
+c   print the absolute and relative deviations together
+c   with their locations.
+c
+      write(iout,380) devabs,labs,devrel,lrel
+380   format(1x,21habsolute deviation = ,e20.12,9h at index,i5/
+     1 1x,21hrelative deviation = ,f11.7,8h percent,1x,9h at index,i5)
+      go to 50
+      end
diff --git a/math/ieee/chap1/test/test18.f b/math/ieee/chap1/test/test18.f
new file mode 100644
index 00000000..cf550e01
--- /dev/null
+++ b/math/ieee/chap1/test/test18.f
@@ -0,0 +1,71 @@
+c
+c-----------------------------------------------------------------------
+c main program: time-efficient radix-4 fast fourier transform
+c author:       l. robert morris
+c               department of systems engineering and computing science
+c               carleton university, ottawa, canada k1s 5b6
+c
+c input:        the array "a" contains the data to be transformed
+c-----------------------------------------------------------------------
+c
+c       test program for autogen radix-4 fft
+c
+        dimension a(2048),b(2048)
+        common /aa/a
+c
+        ioutd=i1mach(2)
+c
+c       compute dft and idft for n = 64, 256, and 1024 complex points
+c
+        do 1 mm=3,5
+        n=4**mm
+        do 2 j=1,n
+        a(2*j-1)=uni(0)
+        a(2*j  )=uni(0)
+        b(2*j-1)=a(2*j-1)
+2       b(2*j  )=a(2*j)
+c
+c       forward dft
+c
+        call radix4(mm,1,-1)
+c
+        if(mm.ne.3) go to 5
+c
+c       list dft input, output for n = 64 only
+c
+        write(ioutd,98)
+        write(ioutd,100)
+        do 3 j=1,n
+        write(ioutd,96) b(2*j-1),b(2*j),a(2*j-1),a(2*j)
+3       continue
+c
+c       inverse dft
+c
+5       call radix4(mm,0, 1)
+c
+c       list dft input, idft output for n = 64 only
+c
+        if(mm.ne.3) go to 7
+c
+        write(ioutd,99)
+        write(ioutd,100)
+        do 6 j=1,n
+        write(ioutd,96) b(2*j-1),b(2*j),a(2*j-1),a(2*j)
+6       continue
+c
+c       calculate rms error
+c
+7       err=0.0
+        do 8 j=1,n
+8       err=err+(a(2*j-1)-b(2*j-1))**2+(a(2*j)-b(2*j))**2
+        err=sqrt(err/float(n))
+        write(ioutd,97) mm,err
+1       continue
+c
+96      format(1x,4(f10.6,2x))
+97      format(1x,20h   rms error for m =,i2,4h is ,e14.6/)
+98      format(1x,43h       dft input                dft  output/)
+99      format(1x,43h       dft input                idft output/)
+100     format(1x,44h   real        imag       real          imag/)
+        stop
+        end