aboutsummaryrefslogtreecommitdiff
path: root/math/ieee/chap1/test
diff options
context:
space:
mode:
Diffstat (limited to 'math/ieee/chap1/test')
-rw-r--r--math/ieee/chap1/test/test12.f90
-rw-r--r--math/ieee/chap1/test/test13.f260
-rw-r--r--math/ieee/chap1/test/test17.f147
-rw-r--r--math/ieee/chap1/test/test18.f71
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
generated by cgit (git 2.34.1) at 2025-09-20 13:56:29 -0400