From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- math/ieee/chap1/test/test13.f | 260 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 260 insertions(+) create mode 100644 math/ieee/chap1/test/test13.f (limited to 'math/ieee/chap1/test/test13.f') 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 -- cgit