diff options
Diffstat (limited to 'math/ieee/chap1')
46 files changed, 4833 insertions, 0 deletions
diff --git a/math/ieee/chap1/README b/math/ieee/chap1/README new file mode 100644 index 00000000..bf86b45f --- /dev/null +++ b/math/ieee/chap1/README @@ -0,0 +1,14 @@ + This directory contains the IEEE routines from Chapter 1 - +Discrete Fourier Transform Programs. Some of the routines are +provided here as separete files (fourea.f, wfta.f, ...) but are +purposely not put into the library, as these routines are not for +general use (fourea.f is an inefficient, demonstration version; +wfta.f is the Winograd DFT which is slower than the regular DFT +and uses far too much memory ona 16-bit machine to be of +practical use; ...). + The directory "test" contains the test routines from the chapter. +The directory "time" contains some routines to time various of the +routines. + The file "compall" compiles the appropiate routines and then +"mklib" will make the library, which one will probably want to +move to "/usr/lib/libieee.a". diff --git a/math/ieee/chap1/const.f b/math/ieee/chap1/const.f new file mode 100644 index 00000000..54eb4e41 --- /dev/null +++ b/math/ieee/chap1/const.f @@ -0,0 +1,114 @@ +c +c----------------------------------------------------------------------- +c subroutine: const +c computes the multipliers for the various modules +c----------------------------------------------------------------------- +c + subroutine const(co3,co8,co16,co9,cdc,cdd) + double precision dtheta,dtwopi,dsq32,dsq2 + double precision dcos1,dcos2,dcos3,dcos4 + double precision dsin1,dsin2,dsin3,dsin4 + dimension co3(3),co8(8),co16(18),co9(11),cdc(9),cdd(6) + dtwopi=8.0d0*datan(1.0d0) + dsq32=dsqrt(0.75d0) + dsq2=dsqrt(0.5d0) +c +c multipliers for the three point module +c + co3(1)=1.0 + co3(2)=-1.5 + co3(3)=-dsq32 +c +c multipliers for the five point module +c + dtheta=dtwopi/5.0d0 + dcos1=dcos(dtheta) + dcos2=dcos(2.0d0*dtheta) + dsin1=dsin(dtheta) + dsin2=dsin(2.0d0*dtheta) + cdd(1)=1.0 + cdd(2)=-1.25 + cdd(3)=-dsin1-dsin2 + cdd(4)=0.5*(dcos1-dcos2) + cdd(5)=dsin1-dsin2 + cdd(6)=dsin2 +c +c +c multipliers for the seven point module +c + dtheta=dtwopi/7.0d0 + dcos1=dcos(dtheta) + dcos2=dcos(2.0d0*dtheta) + dcos3=dcos(3.0d0*dtheta) + dsin1=dsin(dtheta) + dsin2=dsin(2.0d0*dtheta) + dsin3=dsin(3.0d0*dtheta) + cdc(1)=1.0 + cdc(2)=-7.0d0/6.0d0 + cdc(3)=-(dsin1+dsin2-dsin3)/3.0d0 + cdc(4)=(dcos1+dcos2-2.0d0*dcos3)/3.0d0 + cdc(5)=(2.0d0*dcos1-dcos2-dcos3)/3.0d0 + cdc(6)=-(2.0d0*dsin1-dsin2+dsin3)/3.0d0 + cdc(7)=-(dsin1+dsin2+2.0d0*dsin3)/3.0d0 + cdc(8)=(dcos1-2.0d0*dcos2+dcos3)/3.0d0 + cdc(9)=-(dsin1-2.0d0*dsin2-dsin3)/3.0d0 +c +c multipliers for the eight point module +c + co8(1)=1.0 + co8(2)=1.0 + co8(3)=1.0 + co8(4)=-1.0 + co8(5)=1.0 + co8(6)=-dsq2 + co8(7)=-1.0 + co8(8)=dsq2 +c +c multipliers for the nine point module +c + dtheta=dtwopi/9.0d0 + dcos1=dcos(dtheta) + dcos2=dcos(2.0d0*dtheta) + dcos4=dcos(4.0d0*dtheta) + dsin1=dsin(dtheta) + dsin2=dsin(2.0d0*dtheta) + dsin4=dsin(4.0d0*dtheta) + co9(1)=1.0 + co9(2)=-1.5 + co9(3)=-dsq32 + co9(4)=0.5 + co9(5)=(2.0d0*dcos1-dcos2-dcos4)/3.0d0 + co9(6)=(dcos1-2.0d0*dcos2+dcos4)/3.0d0 + co9(7)=(dcos1+dcos2-2.0d0*dcos4)/3.0d0 + co9(8)=-(2.0d0*dsin1+dsin2-dsin4)/3.0d0 + co9(9)=-(dsin1+2.0d0*dsin2+dsin4)/3.0d0 + co9(10)=-(dsin1-dsin2-2.0d0*dsin4)/3.0d0 + co9(11)=-dsq32 +c +c multipliers for the sixteen point module +c + dtheta=dtwopi/16.0d0 + dcos1=dcos(dtheta) + dcos3=dcos(3.0d0*dtheta) + dsin1=dsin(dtheta) + dsin3=dsin(3.0d0*dtheta) + co16(1)=1.0 + co16(2)=1.0 + co16(3)=1.0 + co16(4)=-1.0 + co16(5)=1.0 + co16(6)=-dsq2 + co16(7)=-1.0 + co16(8)=dsq2 + co16(9)=1.0 + co16(10)=-(dsin1-dsin3) + co16(11)=-dsq2 + co16(12)=-co16(10) + co16(13)=-1.0 + co16(14)=-(dsin1+dsin3) + co16(15)=dsq2 + co16(16)=-co16(14) + co16(17)=-dsin3 + co16(18)=dcos3 + return + end diff --git a/math/ieee/chap1/fast.f b/math/ieee/chap1/fast.f new file mode 100644 index 00000000..9a0ad713 --- /dev/null +++ b/math/ieee/chap1/fast.f @@ -0,0 +1,73 @@ +c +c----------------------------------------------------------------------- +c subroutine: fast +c replaces the real vector b(k), for k=1,2,...,n, +c with its finite discrete fourier transform +c----------------------------------------------------------------------- +c + subroutine fast(b, n) +c +c the dc term is returned in location b(1) with b(2) set to 0. +c thereafter the jth harmonic is returned as a complex +c number stored as b(2*j+1) + i b(2*j+2). +c the n/2 harmonic is returned in b(n+1) with b(n+2) set to 0. +c hence, b must be dimensioned to size n+2. +c the subroutine is called as fast(b,n) where n=2**m and +c b is the real array described above. +c + dimension b(2) + common /cons/ pii, p7, p7two, c22, s22, pi2 +c +c iw is a machine dependent write device number +c + iw = i1mach(2) +c + pii = 4.*atan(1.) + pi8 = pii/8. + p7 = 1./sqrt(2.) + p7two = 2.*p7 + c22 = cos(pi8) + s22 = sin(pi8) + pi2 = 2.*pii + do 10 i=1,15 + m = i + nt = 2**i + if (n.eq.nt) go to 20 + 10 continue + write (iw,9999) +9999 format (33h n is not a power of two for fast) + stop + 20 n4pow = m/2 +c +c do a radix 2 iteration first if one is required. +c + if (m-n4pow*2) 40, 40, 30 + 30 nn = 2 + int = n/nn + call fr2tr(int, b(1), b(int+1)) + go to 50 + 40 nn = 1 +c +c perform radix 4 iterations. +c + 50 if (n4pow.eq.0) go to 70 + do 60 it=1,n4pow + nn = nn*4 + int = n/nn + call fr4tr(int, nn, b(1), b(int+1), b(2*int+1), b(3*int+1), + * b(1), b(int+1), b(2*int+1), b(3*int+1)) + 60 continue +c +c perform in-place reordering. +c + 70 call ford1(m, b) + call ford2(m, b) + t = b(2) + b(2) = 0. + b(n+1) = t + b(n+2) = 0. + do 80 it=4,n,2 + b(it) = -b(it) + 80 continue + return + end diff --git a/math/ieee/chap1/ffa.f b/math/ieee/chap1/ffa.f new file mode 100644 index 00000000..dd62ea01 --- /dev/null +++ b/math/ieee/chap1/ffa.f @@ -0,0 +1,84 @@ +c +c----------------------------------------------------------------------- +c subroutine: ffa +c fast fourier analysis subroutine +c----------------------------------------------------------------------- +c + subroutine ffa(b, nfft) +c +c this subroutine replaces the real vector b(k), (k=1,2,...,n), +c with its finite discrete fourier transform. the dc term is +c returned in location b(1) with b(2) set to 0. thereafter, the +c jth harmonic is returned as a complex number stored as +c b(2*j+1) + i b(2*j+2). note that the n/2 harmonic is returned +c in b(n+1) with b(n+2) set to 0. hence, b must be dimensioned +c to size n+2. +c subroutine is called as ffa (b,n) where n=2**m and b is an +c n term real array. a real-valued, radix 8 algorithm is used +c with in-place reordering and the trig functions are computed as +c needed. +c + dimension b(2) + common /con/ pii, p7, p7two, c22, s22, pi2 +c +c iw is a machine dependent write device number +c + iw = i1mach(2) +c + pii = 4.*atan(1.) + pi8 = pii/8. + p7 = 1./sqrt(2.) + p7two = 2.*p7 + c22 = cos(pi8) + s22 = sin(pi8) + pi2 = 2.*pii + n = 1 + do 10 i=1,15 + m = i + n = n*2 + if (n.eq.nfft) go to 20 + 10 continue + write (iw,9999) +9999 format (30h nfft not a power of 2 for ffa) + stop + 20 continue + n8pow = m/3 +c +c do a radix 2 or radix 4 iteration first if one is required +c + if (m-n8pow*3-1) 50, 40, 30 + 30 nn = 4 + int = n/nn + call r4tr(int, b(1), b(int+1), b(2*int+1), b(3*int+1)) + go to 60 + 40 nn = 2 + int = n/nn + call r2tr(int, b(1), b(int+1)) + go to 60 + 50 nn = 1 +c +c perform radix 8 iterations +c + 60 if (n8pow) 90, 90, 70 + 70 do 80 it=1,n8pow + nn = nn*8 + int = n/nn + call r8tr(int, nn, b(1), b(int+1), b(2*int+1), b(3*int+1), + * b(4*int+1), b(5*int+1), b(6*int+1), b(7*int+1), b(1), + * b(int+1), b(2*int+1), b(3*int+1), b(4*int+1), b(5*int+1), + * b(6*int+1), b(7*int+1)) + 80 continue +c +c perform in-place reordering +c + 90 call ord1(m, b) + call ord2(m, b) + t = b(2) + b(2) = 0. + b(nfft+1) = t + b(nfft+2) = 0. + do 100 i=4,nfft,2 + b(i) = -b(i) + 100 continue + return + end diff --git a/math/ieee/chap1/ffs.f b/math/ieee/chap1/ffs.f new file mode 100644 index 00000000..2858fdb0 --- /dev/null +++ b/math/ieee/chap1/ffs.f @@ -0,0 +1,80 @@ +c +c----------------------------------------------------------------------- +c subroutine: ffs +c fast fourier synthesis subroutine +c radix 8-4-2 +c----------------------------------------------------------------------- +c + subroutine ffs(b, nfft) +c +c this subroutine synthesizes the real vector b(k), where +c k=1,2,...,n. the initial fourier coefficients are placed in +c the b array of size n+2. the dc term is in b(1) with +c b(2) equal to 0. +c the jth harmonic is stored as b(2*j+1) + i b(2*j+2). +c the n/2 harmonic is in b(n+1) with b(n+2) equal to 0. +c the subroutine is called as ffs(b,n) where n=2**m and +c b is the n term real array discussed above. +c + dimension b(2) + common /con1/ pii, p7, p7two, c22, s22, pi2 +c +c iw is a machine dependent write device number +c + iw = i1mach(2) +c + pii = 4.*atan(1.) + pi8 = pii/8. + p7 = 1./sqrt(2.) + p7two = 2.*p7 + c22 = cos(pi8) + s22 = sin(pi8) + pi2 = 2.*pii + n = 1 + do 10 i=1,15 + m = i + n = n*2 + if (n.eq.nfft) go to 20 + 10 continue + write (iw,9999) +9999 format (30h nfft not a power of 2 for ffs) + stop + 20 continue + b(2) = b(nfft+1) + do 30 i=1,nfft + b(i) = b(i)/float(nfft) + 30 continue + do 40 i=4,nfft,2 + b(i) = -b(i) + 40 continue + n8pow = m/3 +c +c reorder the input fourier coefficients +c + call ord2(m, b) + call ord1(m, b) +c + if (n8pow.eq.0) go to 60 +c +c perform the radix 8 iterations +c + nn = n + do 50 it=1,n8pow + int = n/nn + call r8syn(int, nn, b, b(int+1), b(2*int+1), b(3*int+1), + * b(4*int+1), b(5*int+1), b(6*int+1), b(7*int+1), b(1), + * b(int+1), b(2*int+1), b(3*int+1), b(4*int+1), b(5*int+1), + * b(6*int+1), b(7*int+1)) + nn = nn/8 + 50 continue +c +c do a radix 2 or radix 4 iteration if one is required +c + 60 if (m-n8pow*3-1) 90, 80, 70 + 70 int = n/4 + call r4syn(int, b(1), b(int+1), b(2*int+1), b(3*int+1)) + go to 90 + 80 int = n/2 + call r2tr(int, b(1), b(int+1)) + 90 return + end diff --git a/math/ieee/chap1/fft842.f b/math/ieee/chap1/fft842.f new file mode 100644 index 00000000..37a3a651 --- /dev/null +++ b/math/ieee/chap1/fft842.f @@ -0,0 +1,116 @@ +c +c----------------------------------------------------------------------- +c subroutine: fft842 +c fast fourier transform for n=2**m +c complex input +c----------------------------------------------------------------------- +c + subroutine fft842(in, n, x, y) +c +c this program replaces the vector z=x+iy by its finite +c discrete, complex fourier transform if in=0. the inverse transform +c is calculated for in=1. it performs as many base +c 8 iterations as possible and then finishes with a base 4 iteration +c or a 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 array y must be supplied to the subroutine. +c + dimension x(2), y(2), 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 +c iw is a machine dependent write device number +c + iw = i1mach(2) +c + pi2 = 8.*atan(1.) + p7 = 1./sqrt(2.) + 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 + 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(1), x(nxtlt+1), x(2*nxtlt+1), + * x(3*nxtlt+1), x(4*nxtlt+1), x(5*nxtlt+1), x(6*nxtlt+1), + * x(7*nxtlt+1), y(1), y(nxtlt+1), y(2*nxtlt+1), y(3*nxtlt+1), + * y(4*nxtlt+1), y(5*nxtlt+1), y(6*nxtlt+1), y(7*nxtlt+1)) + 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 +c + 70 call r2tx(nthpo, x(1), x(2), y(1), y(2)) + go to 90 +c +c go through the base 4 iteration +c + 80 call r4tx(nthpo, x(1), x(2), x(3), x(4), y(1), y(2), y(3), y(4)) +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 return + end diff --git a/math/ieee/chap1/fftaoh.f b/math/ieee/chap1/fftaoh.f new file mode 100644 index 00000000..f703c98b --- /dev/null +++ b/math/ieee/chap1/fftaoh.f @@ -0,0 +1,82 @@ +c +c----------------------------------------------------------------------- +c subroutine: fftaoh +c compute dft for real, antisymmetric, odd harmonic, n-point sequence +c using n/4-point fft +c antisymmetric sequence means x(m)=-x(n-m), m=1,...,n/2-1 +c odd harmonic means x(2*k)=0, all k, where x(k) is the dft of x(m) +c x(m) has the property x(m)=x(n/2-m), m=0,1,...,n/4-1, x(0)=0 +c note: index m is sequence index--not fortran index +c----------------------------------------------------------------------- +c + subroutine fftaoh(x, n, y) + dimension x(1), y(1) +c +c x = real array which on input contains the (n/4+1) points of the +c input sequence (antisymmetrical) +c on output x contains the n/4 imaginary points of the odd +c harmonics of the transform of the input--i.e. the zero +c valued real parts are not given nor are the zero-valued +c even harmonics +c n = true size of input +c y = scratch array of size n/4+2 +c +c +c handle n = 2 and n = 4 cases separately +c + if (n.gt.4) go to 20 + if (n.eq.4) go to 10 +c +c for n=2, assume x(1)=0, x(2)=0, compute dft directly +c + x(1) = 0. + return +c +c n = 4 case, assume x(1)=x(3)=0, x(2)=-x(4)=x0, compute dft directly +c + 10 x(1) = -2.*x(2) + return + 20 twopi = 8.*atan(1.0) +c +c form new sequence, y(m)=x(2*m)+(x(2*m+1)-x(2*m-1)) +c + no2 = n/2 + no4 = n/4 + no8 = n/8 + if (no8.eq.1) go to 40 + do 30 i=2,no8 + ind = 2*i + t1 = x(ind) - x(ind-2) + y(i) = x(ind-1) + t1 + ind1 = n/4 + 2 - i + y(ind1) = x(ind-1) - t1 + 30 continue + 40 y(1) = 2.*x(2) + y(no8+1) = x(no4+1) +c +c the sequence y (n/4 points) has only odd harmonics +c call subroutine fftohm to exploit odd harmonics +c + call fftohm(y, no2) +c +c form original dft from complex odd harmonics of y(k) +c by unscrambling y(k) +c + tpn = twopi/float(n) + cosi = 2.*cos(tpn) + sini = 2.*sin(tpn) + cosd = cos(tpn*2.) + sind = sin(tpn*2.) + do 50 i=1,no8 + ind = 2*i + bk = y(ind-1)/sini + temp = cosi*cosd - sini*sind + sini = cosi*sind + sini*cosd + cosi = temp + ak = y(ind) + x(i) = ak - bk + ind1 = n/4 + 1 - i + x(ind1) = -ak - bk + 50 continue + return + end diff --git a/math/ieee/chap1/fftasm.f b/math/ieee/chap1/fftasm.f new file mode 100644 index 00000000..30c84b79 --- /dev/null +++ b/math/ieee/chap1/fftasm.f @@ -0,0 +1,67 @@ +c +c----------------------------------------------------------------------- +c subroutine: fftasm +c compute dft for real, antisymmetric, n-point sequence x(m) using +c n/2-point fft +c antisymmetric sequence means x(m)=-x(n-m), m=1,...,n/2-1 +c note: index m is sequence index--not fortran index +c----------------------------------------------------------------------- +c + subroutine fftasm(x, n, y) + dimension x(1), y(1) +c +c x = real array which on input contains the n/2 points of the +c input sequence (asymmetrical) +c on output x contains the n/2+1 imaginary points of the transform +c of the input--i.e. the zero valued real parts are not returned +c n = true size of input +c y = scratch array of size n/2+2 +c +c +c for n = 2, assume x(1)=0, x(2)=0, compute dft directly +c + if (n.eq.2) go to 30 + twopi = 8.*atan(1.0) +c +c form new sequence, y(m)=x(2*m)+(x(2*m+1)-x(2*m-1)) +c + no2 = n/2 + no4 = n/4 + do 10 i=2,no4 + ind = 2*i + t1 = x(ind) - x(ind-2) + y(i) = x(ind-1) + t1 + ind1 = no2 + 2 - i + y(ind1) = -x(ind-1) + t1 + 10 continue + y(1) = 2.*x(2) + y(no4+1) = -2.*x(no2) +c +c take n/2 point (real) fft of y +c + call fast(y, no2) +c +c form original dft by unscrambling y(k) +c use recursion relation to generate sin(tpn*i) multiplier +c + tpn = twopi/float(n) + cosi = 2.*cos(tpn) + sini = 2.*sin(tpn) + cosd = cosi/2. + sind = sini/2. + nind = no4 + 1 + do 20 i=2,nind + ind = 2*i + bk = y(ind-1)/sini + ak = y(ind) + x(i) = ak - bk + ind1 = no2 + 2 - i + x(ind1) = -ak - bk + temp = cosi*cosd - sini*sind + sini = cosi*sind + sini*cosd + cosi = temp + 20 continue + 30 x(1) = 0. + x(no2+1) = 0. + return + end diff --git a/math/ieee/chap1/fftohm.f b/math/ieee/chap1/fftohm.f new file mode 100644 index 00000000..6e3ecc9d --- /dev/null +++ b/math/ieee/chap1/fftohm.f @@ -0,0 +1,101 @@ +c +c----------------------------------------------------------------------- +c subroutine: fftohm +c compute dft for real, n-point, odd harmonic sequences using an +c n/2 point fft +c odd harmonic means x(2*k)=0, all k where x(k) is the dft of x(m) +c note: index m is sequence index--not fortran index +c----------------------------------------------------------------------- +c + subroutine fftohm(x, n) + dimension x(1) +c +c x = real array which on input contains the first n/2 points of the +c input +c on output x contains the n/4 complex values of the odd +c harmonics of the input--stored in the sequence re(x(1)),im(x(1)), +c re(x(2)),im(x(2)),... +c ****note: x must be dimensioned to size n/2+2 for fft routine +c n = true size of x sequence +c +c first compute real(x(1)) and real(x(n/2-1)) separately +c also simultaneously multiply original sequence by sin(twopi*(m-1)/n) +c sin and cos are computed recursively +c +c +c for n = 2, assume x(1)=x0, x(2)=-x0, compute dft directly +c + if (n.gt.2) go to 10 + x(1) = 2.*x(1) + x(2) = 0. + return + 10 twopi = 8.*atan(1.0) + tpn = twopi/float(n) +c +c compute x1=real(x(1)) and x2=imaginary(x(n/2-1)) +c x(n) = x(n)*4.*sin(twopi*(i-1)/n) +c + t1 = 0. +c +c cosd and sind are multipliers for recursion for sin and cos +c cosi and sini are initial conditions for recursion for sin and cos +c + cosd = cos(tpn*2.) + sind = sin(tpn*2.) + cosi = 1. + sini = 0. + no2 = n/2 + do 20 i=1,no2,2 + t = x(i)*cosi + x(i) = x(i)*4.*sini + temp = cosi*cosd - sini*sind + sini = cosi*sind + sini*cosd + cosi = temp + t1 = t1 + t + 20 continue +c +c reset initial conditions (cosi,sini) for new recursion +c + cosi = cos(tpn) + sini = sin(tpn) + t2 = 0. + do 30 i=2,no2,2 + t = x(i)*cosi + x(i) = x(i)*4.*sini + temp = cosi*cosd - sini*sind + sini = cosi*sind + sini*cosd + cosi = temp + t2 = t2 + t + 30 continue + x1 = 2.*(t1+t2) + x2 = 2.*(t1-t2) +c +c take n/2 point (real) fft of preprocessed sequence x +c + call fast(x, no2) +c +c for n = 4--skip recursion and initial conditions +c + if (n.eq.4) go to 50 +c +c initial conditions for recursion +c + x(2) = -x(1)/2. + x(1) = x1 +c +c for n = 8, skip recursion +c + if (n.eq.8) go to 50 +c +c unscramble y(k) using recursion formula +c + nind = no2 - 2 + do 40 i=3,nind,2 + t = x(i) + x(i) = x(i-2) + x(i+1) + x(i+1) = x(i-1) - t + 40 continue + 50 x(no2) = x(no2+1)/2. + x(no2-1) = x2 + return + end diff --git a/math/ieee/chap1/fftsoh.f b/math/ieee/chap1/fftsoh.f new file mode 100644 index 00000000..afc18c17 --- /dev/null +++ b/math/ieee/chap1/fftsoh.f @@ -0,0 +1,81 @@ +c +c----------------------------------------------------------------------- +c subroutine: fftsoh +c compute dft for real, symmetric, odd harmonic, n-point sequence +c using n/4-point fft +c symmetric sequence means x(m)=x(n-m), m=1,...,n/2-1 +c odd harmonic means x(2*k)=0, all k, where x(k) is the dft of x(m) +c x(m) has the property x(m)=-x(n/2-m), m=0,1,...,n/4-1, x(n/4)=0 +c note: index m is sequence index--not fortran index +c----------------------------------------------------------------------- +c + subroutine fftsoh(x, n, y) + dimension x(1), y(1) +c +c x = real array which on input contains the n/4 points of the +c input sequence (symmetrical) +c on output x contains the n/4 real points of the odd harmonics +c of the transform of the input--i.e. the zero valued imaginary +c parts are not given nor are the zero-valued even harmonics +c n = true size of input +c y = scratch array of size n/4+2 +c +c +c handle n = 2 and n = 4 cases separately +c + if (n.gt.4) go to 20 + if (n.eq.4) go to 10 +c +c for n=2, assume x(1)=x0, x(2)=-x0, compute dft directly +c + x(1) = 2.*x(1) + return +c +c n = 4 case, compute dft directly +c + 10 x(1) = 2.*x(1) + return + 20 twopi = 8.*atan(1.0) +c +c form new sequence, y(m)=x(2*m)+(x(2*m+1)-x(2*m-1)) +c + no2 = n/2 + no4 = n/4 + no8 = n/8 + if (no8.eq.1) go to 40 + do 30 i=2,no8 + ind = 2*i + t1 = x(ind) - x(ind-2) + y(i) = x(ind-1) + t1 + ind1 = n/4 + 2 - i + y(ind1) = -x(ind-1) + t1 + 30 continue + 40 y(1) = x(1) + y(no8+1) = -2.*x(no4) +c +c the sequence y (n/4 points) has only odd harmonics +c call subroutine fftohm to exploit odd harmonics +c + call fftohm(y, no2) +c +c form original dft from complex odd harmonics of y(k) +c by unscrambling y(k) +c + tpn = twopi/float(n) + cosi = 2.*cos(tpn) + sini = 2.*sin(tpn) + cosd = cos(tpn*2.) + sind = sin(tpn*2.) + do 50 i=1,no8 + ind = 2*i + bk = y(ind)/sini + temp = cosi*cosd - sini*sind + sini = cosi*sind + sini*cosd + cosi = temp + ak = y(ind-1) + x(i) = ak + bk + ind1 = n/4 + 1 - i + x(ind1) = ak - bk + 50 continue + return + end diff --git a/math/ieee/chap1/fftsym.f b/math/ieee/chap1/fftsym.f new file mode 100644 index 00000000..b1e795bb --- /dev/null +++ b/math/ieee/chap1/fftsym.f @@ -0,0 +1,84 @@ +c +c----------------------------------------------------------------------- +c subroutine: fftsym +c compute dft for real, symmetric, n-point sequence x(m) using +c n/2-point fft +c symmetric sequence means x(m)=x(n-m), m=1,...,n/2-1 +c note: index m is sequence index--not fortran index +c----------------------------------------------------------------------- +c + subroutine fftsym(x, n, y) + dimension x(1), y(1) +c +c x = real array which on input contains the n/2+1 points of the +c input sequence (symmetrical) +c on output x contains the n/2+1 real points of the transform of +c the input--i.e. the zero valued imaginary parts are not returned +c n = true size of input +c y = scratch array of size n/2+2 +c +c +c for n = 2, compute dft directly +c + if (n.gt.2) go to 10 + t = x(1) + x(2) + x(2) = x(1) - x(2) + x(1) = t + return + 10 twopi = 8.*atan(1.0) +c +c first compute b0 term, where b0=sum of odd values of x(m) +c + no2 = n/2 + no4 = n/4 + nind = no2 + 1 + b0 = 0. + do 20 i=2,nind,2 + b0 = b0 + x(i) + 20 continue + b0 = b0*2. +c +c for n = 4 skip recursion loop +c + if (n.eq.4) go to 40 +c +c form new sequence, y(m)=x(2*m)+(x(2*m+1)-x(2*m-1)) +c + do 30 i=2,no4 + ind = 2*i + t1 = x(ind) - x(ind-2) + y(i) = x(ind-1) + t1 + ind1 = no2 + 2 - i + y(ind1) = x(ind-1) - t1 + 30 continue + 40 y(1) = x(1) + y(no4+1) = x(no2+1) +c +c take n/2 point (real) fft of y +c + call fast(y, no2) +c +c form original dft by unscrambling y(k) +c use recursion to give sin(tpn*i) multiplier +c + tpn = twopi/float(n) + cosi = 2.*cos(tpn) + sini = 2.*sin(tpn) + cosd = cosi/2. + sind = sini/2. + nind = no4 + 1 + do 50 i=2,nind + ind = 2*i + bk = y(ind)/sini + ak = y(ind-1) + x(i) = ak + bk + nind1 = n/2 + 2 - i + x(nind1) = ak - bk + temp = cosi*cosd - sini*sind + sini = cosi*sind + sini*cosd + cosi = temp + 50 continue + x(1) = b0 + y(1) + x(no2+1) = y(1) - b0 + return + end diff --git a/math/ieee/chap1/ford1.f b/math/ieee/chap1/ford1.f new file mode 100644 index 00000000..4982246f --- /dev/null +++ b/math/ieee/chap1/ford1.f @@ -0,0 +1,24 @@ +c +c----------------------------------------------------------------------- +c subroutine: ford1 +c in-place reordering subroutine +c----------------------------------------------------------------------- +c + subroutine ford1(m, b) + dimension b(2) +c + k = 4 + kl = 2 + n = 2**m + do 40 j=4,n,2 + if (k-j) 20, 20, 10 + 10 t = b(j) + b(j) = b(k) + b(k) = t + 20 k = k - 2 + if (k-kl) 30, 30, 40 + 30 k = 2*j + kl = j + 40 continue + return + end diff --git a/math/ieee/chap1/ford2.f b/math/ieee/chap1/ford2.f new file mode 100644 index 00000000..8ee26d69 --- /dev/null +++ b/math/ieee/chap1/ford2.f @@ -0,0 +1,46 @@ +c +c----------------------------------------------------------------------- +c subroutine: ford2 +c in-place reordering subroutine +c----------------------------------------------------------------------- +c + subroutine ford2(m, b) + dimension l(15), b(2) + 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)) + n = 2**m + l(1) = n + do 10 k=2,m + l(k) = l(k-1)/2 + 10 continue + do 20 k=m,14 + l(k+1) = 2 + 20 continue + ij = 2 + do 40 j1=2,l1,2 + do 40 j2=j1,l2,l1 + do 40 j3=j2,l3,l2 + do 40 j4=j3,l4,l3 + do 40 j5=j4,l5,l4 + do 40 j6=j5,l6,l5 + do 40 j7=j6,l7,l6 + do 40 j8=j7,l8,l7 + do 40 j9=j8,l9,l8 + do 40 j10=j9,l10,l9 + do 40 j11=j10,l11,l10 + do 40 j12=j11,l12,l11 + do 40 j13=j12,l13,l12 + do 40 j14=j13,l14,l13 + do 40 ji=j14,l15,l14 + if (ij-ji) 30, 40, 40 + 30 t = b(ij-1) + b(ij-1) = b(ji-1) + b(ji-1) = t + t = b(ij) + b(ij) = b(ji) + b(ji) = t + 40 ij = ij + 2 + return + end diff --git a/math/ieee/chap1/fourea.f b/math/ieee/chap1/fourea.f new file mode 100644 index 00000000..d412c0e2 --- /dev/null +++ b/math/ieee/chap1/fourea.f @@ -0,0 +1,98 @@ +c +c----------------------------------------------------------------------- +c subroutine: fourea +c performs cooley-tukey fast fourier transform +c----------------------------------------------------------------------- +c + subroutine fourea(data, n, isi) +c +c the cooley-tukey fast fourier transform in ansi fortran +c +c data is a one-dimensional complex array whose length, n, is a +c power of two. isi is +1 for an inverse transform and -1 for a +c forward transform. transform values are returned in the input +c array, replacing the input. +c transform(j)=sum(data(i)*w**((i-1)*(j-1))), where i and j run +c from 1 to n and w = exp (isi*2*pi*sqrt(-1)/n). program also +c computes inverse transform, for which the defining expression +c is invtr(j)=(1/n)*sum(data(i)*w**((i-1)*(j-1))). +c running time is proportional to n*log2(n), rather than to the +c classical n**2. +c after program by brenner, june 1967. this is a very short version +c of the fft and is intended mainly for demonstration. programs +c are available in this collection which run faster and are not +c restricted to powers of 2 or to one-dimensional arrays. +c see -- ieee trans audio (june 1967), special issue on fft. +c + complex data(1) + complex temp, w + ioutd = i1mach(2) +c +c check for power of two up to 15 +c + nn = 1 + do 10 i=1,15 + m = i + nn = nn*2 + if (nn.eq.n) go to 20 + 10 continue + write (ioutd,9999) +9999 format (30h n not a power of 2 for fourea) + stop + 20 continue +c + pi = 4.*atan(1.) + fn = n +c +c this section puts data in bit-reversed order +c + j = 1 + do 80 i=1,n +c +c at this point, i and j are a bit reversed pair (except for the +c displacement of +1) +c + if (i-j) 30, 40, 40 +c +c exchange data(i) with data(j) if i.lt.j. +c + 30 temp = data(j) + data(j) = data(i) + data(i) = temp +c +c implement j=j+1, bit-reversed counter +c + 40 m = n/2 + 50 if (j-m) 70, 70, 60 + 60 j = j - m + m = (m+1)/2 + go to 50 + 70 j = j + m + 80 continue +c +c now compute the butterflies +c + mmax = 1 + 90 if (mmax-n) 100, 130, 130 + 100 istep = 2*mmax + do 120 m=1,mmax + theta = pi*float(isi*(m-1))/float(mmax) + w = cmplx(cos(theta),sin(theta)) + do 110 i=m,n,istep + j = i + mmax + temp = w*data(j) + data(j) = data(i) - temp + data(i) = data(i) + temp + 110 continue + 120 continue + mmax = istep + go to 90 + 130 if (isi) 160, 140, 140 +c +c for inv trans -- isi=1 -- multiply output by 1/n +c + 140 do 150 i=1,n + data(i) = data(i)/fn + 150 continue + 160 return + end diff --git a/math/ieee/chap1/fr2tr.f b/math/ieee/chap1/fr2tr.f new file mode 100644 index 00000000..24a103ff --- /dev/null +++ b/math/ieee/chap1/fr2tr.f @@ -0,0 +1,15 @@ +c +c----------------------------------------------------------------------- +c subroutine: fr2tr +c radix 2 iteration subroutine +c----------------------------------------------------------------------- +c + subroutine fr2tr(int, b0, b1) + dimension b0(2), b1(2) + do 10 k=1,int + t = b0(k) + b1(k) + b1(k) = b0(k) - b1(k) + b0(k) = t + 10 continue + return + end diff --git a/math/ieee/chap1/fr4syn.f b/math/ieee/chap1/fr4syn.f new file mode 100644 index 00000000..5ce978f0 --- /dev/null +++ b/math/ieee/chap1/fr4syn.f @@ -0,0 +1,109 @@ +c +c----------------------------------------------------------------------- +c subroutine: fr4syn +c radix 4 synthesis +c----------------------------------------------------------------------- +c +c + subroutine fr4syn(int, nn, b0, b1, b2, b3, b4, b5, b6, b7) + dimension l(15), b0(2), b1(2), b2(2), b3(2), b4(2), b5(2), b6(2), + * b7(2) + common /const/ pii, p7, p7two, c22, s22, pi2 + 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 + l(1) = nn/4 + do 40 k=2,15 + if (l(k-1)-2) 10, 20, 30 + 10 l(k-1) = 2 + 20 l(k) = 2 + go to 40 + 30 l(k) = l(k-1)/2 + 40 continue +c + piovn = pii/float(nn) + ji = 3 + jl = 2 + jr = 2 +c + do 120 j1=2,l1,2 + do 120 j2=j1,l2,l1 + do 120 j3=j2,l3,l2 + do 120 j4=j3,l4,l3 + do 120 j5=j4,l5,l4 + do 120 j6=j5,l6,l5 + do 120 j7=j6,l7,l6 + do 120 j8=j7,l8,l7 + do 120 j9=j8,l9,l8 + do 120 j10=j9,l10,l9 + do 120 j11=j10,l11,l10 + do 120 j12=j11,l12,l11 + do 120 j13=j12,l13,l12 + do 120 j14=j13,l14,l13 + do 120 jthet=j14,l15,l14 + th2 = jthet - 2 + if (th2) 50, 50, 90 + 50 do 60 k=1,int + t0 = b0(k) + b1(k) + t1 = b0(k) - b1(k) + t2 = b2(k)*2.0 + t3 = b3(k)*2.0 + b0(k) = t0 + t2 + b2(k) = t0 - t2 + b1(k) = t1 + t3 + b3(k) = t1 - t3 + 60 continue +c + if (nn-4) 120, 120, 70 + 70 k0 = int*4 + 1 + kl = k0 + int - 1 + do 80 k=k0,kl + t2 = b0(k) - b2(k) + t3 = b1(k) + b3(k) + b0(k) = (b0(k)+b2(k))*2.0 + b2(k) = (b3(k)-b1(k))*2.0 + b1(k) = (t2+t3)*p7two + b3(k) = (t3-t2)*p7two + 80 continue + go to 120 + 90 arg = th2*piovn + 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 +c + int4 = int*4 + j0 = jr*int4 + 1 + k0 = ji*int4 + 1 + jlast = j0 + int - 1 + do 100 j=j0,jlast + k = k0 + j - j0 + t0 = b0(j) + b6(k) + t1 = b7(k) - b1(j) + t2 = b0(j) - b6(k) + t3 = b7(k) + b1(j) + t4 = b2(j) + b4(k) + t5 = b5(k) - b3(j) + t6 = b5(k) + b3(j) + t7 = b4(k) - b2(j) + b0(j) = t0 + t4 + b4(k) = t1 + t5 + b1(j) = (t2+t6)*c1 - (t3+t7)*s1 + b5(k) = (t2+t6)*s1 + (t3+t7)*c1 + b2(j) = (t0-t4)*c2 - (t1-t5)*s2 + b6(k) = (t0-t4)*s2 + (t1-t5)*c2 + b3(j) = (t2-t6)*c3 - (t3-t7)*s3 + b7(k) = (t2-t6)*s3 + (t3-t7)*c3 + 100 continue + jr = jr + 2 + ji = ji - 2 + if (ji-jl) 110, 110, 120 + 110 ji = 2*jr - 1 + jl = jr + 120 continue + return + end diff --git a/math/ieee/chap1/fr4tr.f b/math/ieee/chap1/fr4tr.f new file mode 100644 index 00000000..16ba3fe6 --- /dev/null +++ b/math/ieee/chap1/fr4tr.f @@ -0,0 +1,118 @@ +c +c----------------------------------------------------------------------- +c subroutine: fr4tr +c radix 4 iteration subroutine +c----------------------------------------------------------------------- +c + subroutine fr4tr(int, nn, b0, b1, b2, b3, b4, b5, b6, b7) + dimension l(15), b0(2), b1(2), b2(2), b3(2), b4(2), b5(2), b6(2), + * b7(2) + common /cons/ pii, p7, p7two, c22, s22, pi2 + 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 jthet is a reversed binary counter, jr steps two at a time to +c locate the real parts of intermediate results, and ji locates +c the imaginary part corresponding to jr. +c + l(1) = nn/4 + do 40 k=2,15 + if (l(k-1)-2) 10, 20, 30 + 10 l(k-1) = 2 + 20 l(k) = 2 + go to 40 + 30 l(k) = l(k-1)/2 + 40 continue +c + piovn = pii/float(nn) + ji = 3 + jl = 2 + jr = 2 +c + do 120 j1=2,l1,2 + do 120 j2=j1,l2,l1 + do 120 j3=j2,l3,l2 + do 120 j4=j3,l4,l3 + do 120 j5=j4,l5,l4 + do 120 j6=j5,l6,l5 + do 120 j7=j6,l7,l6 + do 120 j8=j7,l8,l7 + do 120 j9=j8,l9,l8 + do 120 j10=j9,l10,l9 + do 120 j11=j10,l11,l10 + do 120 j12=j11,l12,l11 + do 120 j13=j12,l13,l12 + do 120 j14=j13,l14,l13 + do 120 jthet=j14,l15,l14 + th2 = jthet - 2 + if (th2) 50, 50, 90 + 50 do 60 k=1,int + t0 = b0(k) + b2(k) + t1 = b1(k) + b3(k) + b2(k) = b0(k) - b2(k) + b3(k) = b1(k) - b3(k) + b0(k) = t0 + t1 + b1(k) = t0 - t1 + 60 continue +c + if (nn-4) 120, 120, 70 + 70 k0 = int*4 + 1 + kl = k0 + int - 1 + do 80 k=k0,kl + pr = p7*(b1(k)-b3(k)) + pi = p7*(b1(k)+b3(k)) + b3(k) = b2(k) + pi + b1(k) = pi - b2(k) + b2(k) = b0(k) - pr + b0(k) = b0(k) + pr + 80 continue + go to 120 +c + 90 arg = th2*piovn + 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 +c + int4 = int*4 + j0 = jr*int4 + 1 + k0 = ji*int4 + 1 + jlast = j0 + int - 1 + do 100 j=j0,jlast + k = k0 + j - j0 + r1 = b1(j)*c1 - b5(k)*s1 + r5 = b1(j)*s1 + b5(k)*c1 + t2 = b2(j)*c2 - b6(k)*s2 + t6 = b2(j)*s2 + b6(k)*c2 + t3 = b3(j)*c3 - b7(k)*s3 + t7 = b3(j)*s3 + b7(k)*c3 + t0 = b0(j) + t2 + t4 = b4(k) + t6 + t2 = b0(j) - t2 + t6 = b4(k) - t6 + t1 = r1 + t3 + t5 = r5 + t7 + t3 = r1 - t3 + t7 = r5 - t7 + b0(j) = t0 + t1 + b7(k) = t4 + t5 + b6(k) = t0 - t1 + b1(j) = t5 - t4 + b2(j) = t2 - t7 + b5(k) = t6 + t3 + b4(k) = t2 + t7 + b3(j) = t3 - t6 + 100 continue +c + jr = jr + 2 + ji = ji - 2 + if (ji-jl) 110, 110, 120 + 110 ji = 2*jr - 1 + jl = jr + 120 continue + return + end diff --git a/math/ieee/chap1/fsst.f b/math/ieee/chap1/fsst.f new file mode 100644 index 00000000..75ce56cb --- /dev/null +++ b/math/ieee/chap1/fsst.f @@ -0,0 +1,71 @@ +c +c----------------------------------------------------------------------- +c subroutine: fsst +c fourier synthesis subroutine +c----------------------------------------------------------------------- +c + subroutine fsst(b, n) +c +c this subroutine synthesizes the real vector b(k), for +c k=1,2,...,n, from the fourier coefficients stored in the +c b array of size n+2. the dc term is in b(1) with b(2) equal +c to 0. the jth harmonic is stored as b(2*j+1) + i b(2*j+2). +c the n/2 harmonic is in b(n+1) with b(n+2) equal to 0. +c the subroutine is called as fsst(b,n) where n=2**m and +c b is the real array discussed above. +c + dimension b(2) + common /const/ pii, p7, p7two, c22, s22, pi2 +c +c iw is a machine dependent write device number +c + iw = i1mach(2) +c + pii = 4.*atan(1.) + pi8 = pii/8. + p7 = 1./sqrt(2.) + p7two = 2.*p7 + c22 = cos(pi8) + s22 = sin(pi8) + pi2 = 2.*pii + do 10 i=1,15 + m = i + nt = 2**i + if (n.eq.nt) go to 20 + 10 continue + write (iw,9999) +9999 format (33h n is not a power of two for fsst) + stop + 20 b(2) = b(n+1) + do 30 i=4,n,2 + b(i) = -b(i) + 30 continue +c +c scale the input by n +c + do 40 i=1,n + b(i) = b(i)/float(n) + 40 continue + n4pow = m/2 +c +c scramble the inputs +c + call ford2(m, b) + call ford1(m, b) +c + if (n4pow.eq.0) go to 60 + nn = 4*n + do 50 it=1,n4pow + nn = nn/4 + int = n/nn + call fr4syn(int, nn, b(1), b(int+1), b(2*int+1), b(3*int+1), + * b(1), b(int+1), b(2*int+1), b(3*int+1)) + 50 continue +c +c do a radix 2 iteration if one is required +c + 60 if (m-n4pow*2) 80, 80, 70 + 70 int = n/2 + call fr2tr(int, b(1), b(int+1)) + 80 return + end diff --git a/math/ieee/chap1/iftaoh.f b/math/ieee/chap1/iftaoh.f new file mode 100644 index 00000000..8cb0a908 --- /dev/null +++ b/math/ieee/chap1/iftaoh.f @@ -0,0 +1,87 @@ +c +c----------------------------------------------------------------------- +c subroutine: iftaoh +c compute idft for real, antisymmetric, odd harmonic, n-point sequence +c using n/4-point fft +c antisymmetric sequence means x(m)=-x(n-m), m=1,...,n/2-1 +c odd harmonic means x(2*k)=0, all k, where x(k) is the dft of x(m) +c x(m)has the property x(m)=x(n/2-m), m=0,1,...,n/4-1, x(0)=0 +c note: index m is sequence index--not fortran index +c----------------------------------------------------------------------- +c + subroutine iftaoh(x, n, y) + dimension x(1), y(1) +c +c x = real array which on input contains the n/4 imaginary points +c of the odd harmonics of the transform of the original time +c sequence--i.e. the zero valued real parts are not input nor +c are the zero-valued even harmonics +c on output x contains the first (n/4+1) points of the original +c time sequence (antisymmetrical) +c n = true size of input +c y = scratch array of size n/4+2 +c +c +c handle n = 2 and n = 4 cases separately +c + if (n.gt.4) go to 20 + if (n.eq.4) go to 10 +c +c for n=2 assume x(1)=0, x(2)=0, compute idft directly +c + x(1) = 0. + return +c +c for n=4, assume x(1)=x(3)=0, x(2)=-x(4)=x0, compute idft directly +c + 10 x(2) = -x(1)/2. + x(1) = 0. + return +c +c code for values of n which are multiples of 8 +c + 20 twopi = 8.*atan(1.0) + no2 = n/2 + no4 = n/4 + no8 = n/8 + tpn = twopi/float(n) +c +c scramble original dft (x(k)) to give y(k) +c use recursion to give sin multipliers +c + cosi = cos(tpn) + sini = sin(tpn) + cosd = cos(tpn*2.) + sind = sin(tpn*2.) + do 30 i=1,no8 + ind = 2*i + ind1 = no4 + 1 - i + ak = (x(i)-x(ind1))/2. + bk = -(x(i)+x(ind1)) + y(ind) = ak + y(ind-1) = bk*sini + temp = cosi*cosd - sini*sind + sini = cosi*sind + sini*cosd + cosi = temp + 30 continue +c +c the sequence y(k) is an odd harmonic sequence +c use subroutine iftohm to give y(m) +c + call iftohm(y, no2) +c +c form x sequence from y sequence +c + x(2) = y(1)/2. + x(1) = 0. + if (n.eq.8) return + do 40 i=2,no8 + ind = 2*i + ind1 = no4 + 2 - i + x(ind-1) = (y(i)+y(ind1))/2. + t1 = (y(i)-y(ind1))/2. + x(ind) = t1 + x(ind-2) + 40 continue + x(no4+1) = y(no8+1) + return + end diff --git a/math/ieee/chap1/iftasm.f b/math/ieee/chap1/iftasm.f new file mode 100644 index 00000000..6ec12c7d --- /dev/null +++ b/math/ieee/chap1/iftasm.f @@ -0,0 +1,77 @@ +c +c----------------------------------------------------------------------- +c subroutine: iftasm +c compute idft for real, antisymmetric, n-point sequence x(m) using +c n/2-point fft +c antisymmetric sequence means x(m)=-x(n-m), m=1,...,n/2-1 +c note: index m is sequence index--not fortran index +c----------------------------------------------------------------------- +c + subroutine iftasm(x, n, y) + dimension x(1), y(1) +c +c x = imaginary array which on input contains the n/2+1 real points of +c the transform of the input--i.e. the zero valued real parts +c are not given as input +c on output x contains the n/2 points of the time sequence +c (antisymmetrical) +c n = true size of input +c y = scratch array of size n/2+2 +c +c +c for n = 2, assume x(1)=0, x(2)=0 +c + if (n.gt.2) go to 10 + x(1) = 0 + x(2) = 0 + return + 10 twopi = 8.*atan(1.0) +c +c first compute x1=x(1) term directly +c use recursion on the sine cosine terms +c + no2 = n/2 + no4 = n/4 + tpn = twopi/float(n) +c +c scramble original dft (x(k)) to give y(k) +c use recursion relation to give sin(tpn*i) multiplier +c + cosi = cos(tpn) + sini = sin(tpn) + cosd = cosi + sind = sini + nind = no4 + 1 + do 20 i=2,nind + ind = 2*i + ind1 = no2 + 2 - i + ak = (x(i)-x(ind1))/2. + bk = -(x(i)+x(ind1)) + y(ind) = ak + y(ind-1) = bk*sini + temp = cosi*cosd - sini*sind + sini = cosi*sind + sini*cosd + cosi = temp + 20 continue + y(1) = 0. + y(2) = 0. +c +c take n/2 point idft of y +c + call fsst(y, no2) +c +c form x sequence from y sequence +c + x(2) = y(1)/2. + x(1) = 0. + if (n.eq.4) go to 40 + do 30 i=2,no4 + ind = 2*i + ind1 = no2 + 2 - i + x(ind-1) = (y(i)-y(ind1))/2. + t1 = (y(i)+y(ind1))/2. + x(ind) = t1 + x(ind-2) + 30 continue + 40 x(no2) = -y(no4+1)/2. + return + end diff --git a/math/ieee/chap1/iftohm.f b/math/ieee/chap1/iftohm.f new file mode 100644 index 00000000..df084e87 --- /dev/null +++ b/math/ieee/chap1/iftohm.f @@ -0,0 +1,83 @@ +c +c----------------------------------------------------------------------- +c subroutine: iftohm +c compute idft for real, n-point, odd harmonic sequences using an +c n/2 point fft +c odd harmonic means x(2*k)=0, all k where x(k) is the dft of x(m) +c note: index m is sequence index--not fortran index +c----------------------------------------------------------------------- +c + subroutine iftohm(x, n) + dimension x(1) +c +c x = real array which on input contains the n/4 complex values of the +c odd harmonics of the input--stored in the sequence re(x(1)), +c im(x(1)),re(x(2)),im(x(2)),... +c on output x contains the first n/2 points of the input +c ****note: x must be dimensioned to size n/2+2 for fft routine +c n = true size of x sequence +c +c first compute real(x(1)) and real(x(n/2-1)) separately +c also simultaneously multiply original sequence by sin(twopi*(m-1)/n) +c sin and cos are computed recursively +c +c +c for n = 2, assume x(1)=x0, x(2)=-x0, compute idft directly +c + if (n.gt.2) go to 10 + x(1) = 0.5*x(1) + x(2) = -x(1) + return + 10 twopi = 8.*atan(1.0) + tpn = twopi/float(n) + no2 = n/2 + no4 = n/4 + nind = no2 +c +c solve for x(0)=x0 directly +c + x0 = 0. + do 20 i=1,no2,2 + x0 = x0 + 2.*x(i) + 20 continue + x0 = x0/float(n) +c +c form y(k)=j*(x(2k+1)-x(2k-1)) +c overwrite x array with y sequence +c + xpr = x(1) + xpi = x(2) + x(1) = -2.*x(2) + x(2) = 0. + if (no4.eq.1) go to 40 + do 30 i=3,nind,2 + ti = x(i) - xpr + tr = -x(i+1) + xpi + xpr = x(i) + xpi = x(i+1) + x(i) = tr + x(i+1) = ti + 30 continue + 40 x(no2+1) = 2.*xpi + x(no2+2) = 0. +c +c take n/2 point (real) ifft of preprocessed sequence x +c + call fsst(x, no2) +c +c solve for x(m) by dividing by 4*sin(twopi*m/n) for m=1,2,...,n/2-1 +c for m=0 substitute precomputed value x0 +c + cosi = 4. + sini = 0. + cosd = cos(tpn) + sind = sin(tpn) + do 50 i=2,no2 + temp = cosi*cosd - sini*sind + sini = cosi*sind + sini*cosd + cosi = temp + x(i) = x(i)/sini + 50 continue + x(1) = x0 + return + end diff --git a/math/ieee/chap1/iftsoh.f b/math/ieee/chap1/iftsoh.f new file mode 100644 index 00000000..6098e6e4 --- /dev/null +++ b/math/ieee/chap1/iftsoh.f @@ -0,0 +1,94 @@ +c +c----------------------------------------------------------------------- +c subroutine: iftsoh +c compute idft for real, symmetric, odd harmonic, n-point sequence +c using n/4-point fft +c symmetric sequence means x(m)=x(n-m), m=1,...,n/2-1 +c odd harmonic means x(2*k)=0, all k, where x(k) is the dft of x(m) +c x(m) has the property x(m)=-x(n/2-m), m=0,1,...,n/4-1, x(n/4)=0 +c note: index m is sequence index--not fortran index +c----------------------------------------------------------------------- +c + subroutine iftsoh(x, n, y) + dimension x(1), y(1) +c +c x = real array which on input contains the n/4 real points of +c the odd harmonics of the transform of the original time sequence +c i.e. the zero valued imaginary parts are not given nor are the +c zero valued even harmonics +c on output x contains the first n/4 points of the original input +c sequence (symmetrical) +c n = true size of input +c y = scratch array of size n/4+2 +c +c +c handle n = 2 and n = 4 cases separately +c + if (n.gt.4) go to 10 +c +c for n=2, 4 assume x(1)=x0, x(2)=-x0, compute idft directly +c + x(1) = x(1)/2. + return +c +c code for values of n which are multiples of 8 +c + 10 twopi = 8.*atan(1.0) + no2 = n/2 + no4 = n/4 + no8 = n/8 + tpn = twopi/float(n) +c +c first compute x1=x(1) term directly +c use recursion on the sine cosine terms +c + cosd = cos(tpn*2.) + sind = sin(tpn*2.) + cosi = 2.*cos(tpn) + sini = 2.*sin(tpn) + x1 = 0. + do 20 i=1,no4 + x1 = x1 + x(i)*cosi + temp = cosi*cosd - sini*sind + sini = cosi*sind + sini*cosd + cosi = temp + 20 continue + x1 = x1/float(n) +c +c scramble original dft (x(k)) to give y(k) +c use recursion relation to give sin multipliers +c + cosi = cos(tpn) + sini = sin(tpn) + do 30 i=1,no8 + ind = 2*i + ind1 = no4 + 1 - i + ak = (x(i)+x(ind1))/2. + bk = (x(i)-x(ind1)) + y(ind-1) = ak + y(ind) = bk*sini + temp = cosi*cosd - sini*sind + sini = cosi*sind + sini*cosd + cosi = temp + 30 continue +c +c the sequence y(k) is the odd harmonics dft output +c use subroutine iftohm to obtain y(m), the inverse transform +c + call iftohm(y, no2) +c +c form x(m) sequence from y(m) sequence +c use x1 initial condition on the recursion +c + x(1) = y(1) + x(2) = x1 + if (no8.eq.1) return + do 40 i=2,no8 + ind = 2*i + ind1 = no4 + 2 - i + t1 = (y(i)+y(ind1))/2. + x(ind-1) = (y(i)-y(ind1))/2. + x(ind) = t1 + x(ind-2) + 40 continue + return + end diff --git a/math/ieee/chap1/iftsym.f b/math/ieee/chap1/iftsym.f new file mode 100644 index 00000000..e45c9a37 --- /dev/null +++ b/math/ieee/chap1/iftsym.f @@ -0,0 +1,90 @@ +c +c----------------------------------------------------------------------- +c subroutine: iftsym +c compute idft for real, symmetric, n-point sequence x(m) using +c n/2-point fft +c symmetric sequence means x(m)=x(n-m), m=1,...,n/2-1 +c note: index m is sequence index--not fortran index +c----------------------------------------------------------------------- +c + subroutine iftsym(x, n, y) + dimension x(1), y(1) +c +c x = real array which on input contains the n/2+1 real points of the +c transform of the input--i.e. the zero valued imaginary parts +c are not given as input +c on output x contains the n/2+1 points of the time sequence +c (symmetrical) +c n = true size of input +c y = scratch array of size n/2+2 +c +c +c for n = 2, compute idft directly +c + if (n.gt.2) go to 10 + t = (x(1)+x(2))/2. + x(2) = (x(1)-x(2))/2. + x(1) = t + return + 10 twopi = 8.*atan(1.0) +c +c first compute x1=x(1) term directly +c use recursion on the sine cosine terms +c + no2 = n/2 + no4 = n/4 + tpn = twopi/float(n) + cosd = cos(tpn) + sind = sin(tpn) + cosi = 2. + sini = 0. + x1 = x(1) - x(no2+1) + do 20 i=2,no2 + temp = cosi*cosd - sini*sind + sini = cosi*sind + sini*cosd + cosi = temp + x1 = x1 + x(i)*cosi + 20 continue + x1 = x1/float(n) +c +c scramble original dft (x(k)) to give y(k) +c use recursion relation to generate sin(tpn*i) multiplier +c + cosi = cos(tpn) + sini = sin(tpn) + cosd = cosi + sind = sini + y(1) = (x(1)+x(no2+1))/2. + y(2) = 0. + nind = no4 + 1 + do 30 i=2,nind + ind = 2*i + nind1 = no2 + 2 - i + ak = (x(i)+x(nind1))/2. + bk = (x(i)-x(nind1)) + y(ind-1) = ak + y(ind) = bk*sini + temp = cosi*cosd - sini*sind + sini = cosi*sind + sini*cosd + cosi = temp + 30 continue +c +c take n/2 point idft of y +c + call fsst(y, no2) +c +c form x sequence from y sequence +c + x(1) = y(1) + x(2) = x1 + if (n.eq.4) go to 50 + do 40 i=2,no4 + ind = 2*i + ind1 = no2 + 2 - i + x(ind-1) = (y(i)+y(ind1))/2. + t1 = (y(i)-y(ind1))/2. + x(ind) = t1 + x(ind-2) + 40 continue + 50 x(no2+1) = y(no4+1) + return + end diff --git a/math/ieee/chap1/inishl.f b/math/ieee/chap1/inishl.f new file mode 100644 index 00000000..925cd657 --- /dev/null +++ b/math/ieee/chap1/inishl.f @@ -0,0 +1,179 @@ +c +c----------------------------------------------------------------------- +c subroutine: inishl +c this subroutine initializes the wfta routine for a given +c value of the transform length n. the factors of n are +c determined, the multiplication coefficients are calculated +c and stored in the array coef(.), the input and output +c permutation vectors are computed and stored in the arrays +c indx1(.) and indx2(.) +c +c----------------------------------------------------------------------- +c + subroutine inishl(n,coef,xr,xi,indx1,indx2,ierr) + dimension coef(1),xr(1),xi(1) + integer s1,s2,s3,s4,indx1(1),indx2(1),p1 + dimension co3(3),co4(4),co8(8),co9(11),co16(18),cda(18),cdb(11), + 1cdc(9),cdd(6) + common na,nb,nc,nd,nd1,nd2,nd3,nd4 +c +c data statements assign short dft coefficients. +c + data co4(1),co4(2),co4(3),co4(4)/4*1.0/ +c + data cda(1),cda(2),cda(3),cda(4),cda(5),cda(6),cda(7), + 1 cda(8),cda(9),cda(10),cda(11),cda(12),cda(13),cda(14), + 2 cda(15),cda(16),cda(17),cda(18)/18*1.0/ +c + data cdb(1),cdb(2),cdb(3),cdb(4),cdb(5),cdb(6),cdb(7),cdb(8), + 1 cdb(9),cdb(10),cdb(11)/11*1.0/ +c + data ionce/1/ +c +c get multiplier constants +c + if(ionce.ne.1) go to 20 + call const(co3,co8,co16,co9,cdc,cdd) +20 ionce=-1 +c +c following segment determines factors of n and chooses +c the appropriate short dft coefficients. +c + iout=i1mach(2) + ierr=0 + nd1=1 + na=1 + nb=1 + nd2=1 + nc=1 + nd3=1 + nd=1 + nd4=1 + if(n.le.0 .or. n.gt.5040) go to 190 + if(16*(n/16).eq.n) go to 30 + if(8*(n/8).eq.n) go to 40 + if(4*(n/4).eq.n) go to 50 + if(2*(n/2).ne.n) go to 70 + nd1=2 + na=2 + cda(2)=1.0 + go to 70 +30 nd1=18 + na=16 + do 31 j=1,18 +31 cda(j)=co16(j) + go to 70 +40 nd1=8 + na=8 + do 41 j=1,8 +41 cda(j)=co8(j) + go to 70 +50 nd1=4 + na=4 + do 51 j=1,4 +51 cda(j)=co4(j) +70 if(3*(n/3).ne.n) go to 120 + if(9*(n/9).eq.n) go to 100 + nd2=3 + nb=3 + do 71 j=1,3 +71 cdb(j)=co3(j) + go to 120 +100 nd2=11 + nb=9 + do 110 j=1,11 +110 cdb(j)=co9(j) +120 if(7*(n/7).ne.n) go to 160 + nd3=9 + nc=7 +160 if(5*(n/5).ne.n) go to 190 + nd4=6 + nd=5 +190 m=na*nb*nc*nd + if(m.eq.n) go to 250 + write(iout,210) +210 format(21h this n does not work) + ierr=-1 + return +c +c next segment generates the dft coefficients by +c multiplying together the short dft coefficients +c +250 j=1 + do 300 n4=1,nd4 + do 300 n3=1,nd3 + do 300 n2=1,nd2 + do 300 n1=1,nd1 + coef(j)=cda(n1)*cdb(n2)*cdc(n3)*cdd(n4) + j=j+1 +300 continue +c +c following segment forms the input indexing vector +c + j=1 + nu=nb*nc*nd + nv=na*nc*nd + nw=na*nb*nd + ny=na*nb*nc + k=1 + do 440 n4=1,nd + do 430 n3=1,nc + do 420 n2=1,nb + do 410 n1=1,na +405 if(k.le.n) go to 408 + k=k-n + go to 405 +408 indx1(j)=k + j=j+1 +410 k=k+nu +420 k=k+nv +430 k=k+nw +440 k=k+ny +c +c following segment forms the output indexing vector +c + m=1 + s1=0 + s2=0 + s3=0 + s4=0 + if(na.eq.1) go to 530 +520 p1=m*nu-1 + if((p1/na)*na.eq.p1) go to 510 + m=m+1 + go to 520 +510 s1=p1+1 +530 if(nb.eq.1) go to 540 + m=1 +550 p1=m*nv-1 + if((p1/nb)*nb.eq.p1) go to 560 + m=m+1 + go to 550 +560 s2=p1+1 +540 if(nc.eq.1) go to 630 + m=1 +620 p1=m*nw-1 + if((p1/nc)*nc.eq.p1) go to 610 + m=m+1 + go to 620 +610 s3=p1+1 +630 if(nd.eq.1) go to 660 + m=1 +640 p1=m*ny-1 + if((p1/nd)*nd.eq.p1) go to 650 + m=m+1 + go to 640 +650 s4=p1+1 +660 j=1 + do 810 n4=1,nd + do 810 n3=1,nc + do 810 n2=1,nb + do 810 n1=1,na + indx2(j)=s1*(n1-1)+s2*(n2-1)+s3*(n3-1)+s4*(n4-1)+1 +900 if(indx2(j).le.n) go to 910 + indx2(j)=indx2(j)-n + go to 900 +910 j=j+1 +810 continue + return + end diff --git a/math/ieee/chap1/ord1.f b/math/ieee/chap1/ord1.f new file mode 100644 index 00000000..21e4494e --- /dev/null +++ b/math/ieee/chap1/ord1.f @@ -0,0 +1,24 @@ +c +c----------------------------------------------------------------------- +c subroutine: ord1 +c in-place reordering subroutine +c----------------------------------------------------------------------- +c + subroutine ord1(m, b) + dimension b(2) +c + k = 4 + kl = 2 + n = 2**m + do 40 j=4,n,2 + if (k-j) 20, 20, 10 + 10 t = b(j) + b(j) = b(k) + b(k) = t + 20 k = k - 2 + if (k-kl) 30, 30, 40 + 30 k = 2*j + kl = j + 40 continue + return + end diff --git a/math/ieee/chap1/ord2.f b/math/ieee/chap1/ord2.f new file mode 100644 index 00000000..10f03ec7 --- /dev/null +++ b/math/ieee/chap1/ord2.f @@ -0,0 +1,46 @@ +c +c----------------------------------------------------------------------- +c subroutine: ord2 +c in-place reordering subroutine +c----------------------------------------------------------------------- +c + subroutine ord2(m, b) + dimension l(15), b(2) + 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)) + n = 2**m + l(1) = n + do 10 k=2,m + l(k) = l(k-1)/2 + 10 continue + do 20 k=m,14 + l(k+1) = 2 + 20 continue + ij = 2 + do 40 j1=2,l1,2 + do 40 j2=j1,l2,l1 + do 40 j3=j2,l3,l2 + do 40 j4=j3,l4,l3 + do 40 j5=j4,l5,l4 + do 40 j6=j5,l6,l5 + do 40 j7=j6,l7,l6 + do 40 j8=j7,l8,l7 + do 40 j9=j8,l9,l8 + do 40 j10=j9,l10,l9 + do 40 j11=j10,l11,l10 + do 40 j12=j11,l12,l11 + do 40 j13=j12,l13,l12 + do 40 j14=j13,l14,l13 + do 40 ji=j14,l15,l14 + if (ij-ji) 30, 40, 40 + 30 t = b(ij-1) + b(ij-1) = b(ji-1) + b(ji-1) = t + t = b(ij) + b(ij) = b(ji) + b(ji) = t + 40 ij = ij + 2 + return + end diff --git a/math/ieee/chap1/r2tr.f b/math/ieee/chap1/r2tr.f new file mode 100644 index 00000000..510763fe --- /dev/null +++ b/math/ieee/chap1/r2tr.f @@ -0,0 +1,16 @@ +c +c----------------------------------------------------------------------- +c subroutine: r2tr +c radix 2 iteration subroutine +c----------------------------------------------------------------------- +c +c + subroutine r2tr(int, b0, b1) + dimension b0(2), b1(2) + do 10 k=1,int + t = b0(k) + b1(k) + b1(k) = b0(k) - b1(k) + b0(k) = t + 10 continue + return + end diff --git a/math/ieee/chap1/r2tx.f b/math/ieee/chap1/r2tx.f new file mode 100644 index 00000000..425c10f8 --- /dev/null +++ b/math/ieee/chap1/r2tx.f @@ -0,0 +1,18 @@ +c +c----------------------------------------------------------------------- +c subroutine: r2tx +c radix 2 iteration subroutine +c----------------------------------------------------------------------- +c + subroutine r2tx(nthpo, cr0, cr1, ci0, ci1) + dimension cr0(2), cr1(2), ci0(2), ci1(2) + do 10 k=1,nthpo,2 + r1 = cr0(k) + cr1(k) + cr1(k) = cr0(k) - cr1(k) + cr0(k) = r1 + fi1 = ci0(k) + ci1(k) + ci1(k) = ci0(k) - ci1(k) + ci0(k) = fi1 + 10 continue + return + end diff --git a/math/ieee/chap1/r4syn.f b/math/ieee/chap1/r4syn.f new file mode 100644 index 00000000..77c97bba --- /dev/null +++ b/math/ieee/chap1/r4syn.f @@ -0,0 +1,20 @@ +c +c----------------------------------------------------------------------- +c subroutine: r4syn +c radix 4 synthesis +c----------------------------------------------------------------------- +c + subroutine r4syn(int, b0, b1, b2, b3) + dimension b0(2), b1(2), b2(2), b3(2) + do 10 k=1,int + t0 = b0(k) + b1(k) + t1 = b0(k) - b1(k) + t2 = b2(k) + b2(k) + t3 = b3(k) + b3(k) + b0(k) = t0 + t2 + b2(k) = t0 - t2 + b1(k) = t1 + t3 + b3(k) = t1 - t3 + 10 continue + return + end diff --git a/math/ieee/chap1/r4tr.f b/math/ieee/chap1/r4tr.f new file mode 100644 index 00000000..5fc46eab --- /dev/null +++ b/math/ieee/chap1/r4tr.f @@ -0,0 +1,18 @@ +c +c----------------------------------------------------------------------- +c subroutine: r4tr +c radix 4 iteration subroutine +c----------------------------------------------------------------------- +c + subroutine r4tr(int, b0, b1, b2, b3) + dimension b0(2), b1(2), b2(2), b3(2) + do 10 k=1,int + r0 = b0(k) + b2(k) + r1 = b1(k) + b3(k) + b2(k) = b0(k) - b2(k) + b3(k) = b1(k) - b3(k) + b0(k) = r0 + r1 + b1(k) = r0 - r1 + 10 continue + return + end diff --git a/math/ieee/chap1/r4tx.f b/math/ieee/chap1/r4tx.f new file mode 100644 index 00000000..4e5649e8 --- /dev/null +++ b/math/ieee/chap1/r4tx.f @@ -0,0 +1,29 @@ +c +c----------------------------------------------------------------------- +c subroutine: r4tx +c radix 4 iteration subroutine +c----------------------------------------------------------------------- +c + subroutine r4tx(nthpo, cr0, cr1, cr2, cr3, ci0, ci1, ci2, ci3) + dimension cr0(2), cr1(2), cr2(2), cr3(2), ci0(2), ci1(2), ci2(2), + * ci3(2) + do 10 k=1,nthpo,4 + r1 = cr0(k) + cr2(k) + r2 = cr0(k) - cr2(k) + r3 = cr1(k) + cr3(k) + r4 = cr1(k) - cr3(k) + fi1 = ci0(k) + ci2(k) + fi2 = ci0(k) - ci2(k) + fi3 = ci1(k) + ci3(k) + fi4 = ci1(k) - ci3(k) + cr0(k) = r1 + r3 + ci0(k) = fi1 + fi3 + cr1(k) = r1 - r3 + ci1(k) = fi1 - fi3 + cr2(k) = r2 - fi4 + ci2(k) = fi2 + r4 + cr3(k) = r2 + fi4 + ci3(k) = fi2 - r4 + 10 continue + return + end diff --git a/math/ieee/chap1/r8syn.f b/math/ieee/chap1/r8syn.f new file mode 100644 index 00000000..054413d7 --- /dev/null +++ b/math/ieee/chap1/r8syn.f @@ -0,0 +1,186 @@ +c +c----------------------------------------------------------------------- +c subroutine: r8syn +c radix 8 synthesis subroutine +c----------------------------------------------------------------------- +c + subroutine r8syn(int, nn, br0, br1, br2, br3, br4, br5, br6, br7, + * bi0, bi1, bi2, bi3, bi4, bi5, bi6, bi7) + dimension l(15), br0(2), br1(2), br2(2), br3(2), br4(2), br5(2), + * br6(2), br7(2), bi0(2), bi1(2), bi2(2), bi3(2), bi4(2), + * bi5(2), bi6(2), bi7(2) + common /con1/ pii, p7, p7two, c22, s22, pi2 + 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)) + l(1) = nn/8 + do 40 k=2,15 + if (l(k-1)-2) 10, 20, 30 + 10 l(k-1) = 2 + 20 l(k) = 2 + go to 40 + 30 l(k) = l(k-1)/2 + 40 continue + piovn = pii/float(nn) + ji = 3 + jl = 2 + jr = 2 +c + do 120 j1=2,l1,2 + do 120 j2=j1,l2,l1 + do 120 j3=j2,l3,l2 + do 120 j4=j3,l4,l3 + do 120 j5=j4,l5,l4 + do 120 j6=j5,l6,l5 + do 120 j7=j6,l7,l6 + do 120 j8=j7,l8,l7 + do 120 j9=j8,l9,l8 + do 120 j10=j9,l10,l9 + do 120 j11=j10,l11,l10 + do 120 j12=j11,l12,l11 + do 120 j13=j12,l13,l12 + do 120 j14=j13,l14,l13 + do 120 jthet=j14,l15,l14 + th2 = jthet - 2 + if (th2) 50, 50, 90 + 50 do 60 k=1,int + t0 = br0(k) + br1(k) + t1 = br0(k) - br1(k) + t2 = br2(k) + br2(k) + t3 = br3(k) + br3(k) + t4 = br4(k) + br6(k) + t6 = br7(k) - br5(k) + t5 = br4(k) - br6(k) + t7 = br7(k) + br5(k) + pr = p7*(t7+t5) + pi = p7*(t7-t5) + tt0 = t0 + t2 + tt1 = t1 + t3 + t2 = t0 - t2 + t3 = t1 - t3 + t4 = t4 + t4 + t5 = pr + pr + t6 = t6 + t6 + t7 = pi + pi + br0(k) = tt0 + t4 + br1(k) = tt1 + t5 + br2(k) = t2 + t6 + br3(k) = t3 + t7 + br4(k) = tt0 - t4 + br5(k) = tt1 - t5 + br6(k) = t2 - t6 + br7(k) = t3 - t7 + 60 continue + if (nn-8) 120, 120, 70 + 70 k0 = int*8 + 1 + kl = k0 + int - 1 + do 80 k=k0,kl + t1 = bi0(k) + bi6(k) + t2 = bi7(k) - bi1(k) + t3 = bi0(k) - bi6(k) + t4 = bi7(k) + bi1(k) + pr = t3*c22 + t4*s22 + pi = t4*c22 - t3*s22 + t5 = bi2(k) + bi4(k) + t6 = bi5(k) - bi3(k) + t7 = bi2(k) - bi4(k) + t8 = bi5(k) + bi3(k) + rr = t8*c22 - t7*s22 + ri = -t8*s22 - t7*c22 + bi0(k) = (t1+t5) + (t1+t5) + bi4(k) = (t2+t6) + (t2+t6) + bi1(k) = (pr+rr) + (pr+rr) + bi5(k) = (pi+ri) + (pi+ri) + t5 = t1 - t5 + t6 = t2 - t6 + bi2(k) = p7two*(t6+t5) + bi6(k) = p7two*(t6-t5) + rr = pr - rr + ri = pi - ri + bi3(k) = p7two*(ri+rr) + bi7(k) = p7two*(ri-rr) + 80 continue + go to 120 + 90 arg = th2*piovn + 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 + int8 = int*8 + j0 = jr*int8 + 1 + k0 = ji*int8 + 1 + jlast = j0 + int - 1 + do 100 j=j0,jlast + k = k0 + j - j0 + tr0 = br0(j) + bi6(k) + ti0 = bi7(k) - br1(j) + tr1 = br0(j) - bi6(k) + ti1 = bi7(k) + br1(j) + tr2 = br2(j) + bi4(k) + ti2 = bi5(k) - br3(j) + tr3 = bi5(k) + br3(j) + ti3 = bi4(k) - br2(j) + tr4 = br4(j) + bi2(k) + ti4 = bi3(k) - br5(j) + t0 = br4(j) - bi2(k) + t1 = bi3(k) + br5(j) + tr5 = p7*(t0+t1) + ti5 = p7*(t1-t0) + tr6 = br6(j) + bi0(k) + ti6 = bi1(k) - br7(j) + t0 = br6(j) - bi0(k) + t1 = bi1(k) + br7(j) + tr7 = -p7*(t0-t1) + ti7 = -p7*(t1+t0) + t0 = tr0 + tr2 + t1 = ti0 + ti2 + t2 = tr1 + tr3 + t3 = ti1 + ti3 + tr2 = tr0 - tr2 + ti2 = ti0 - ti2 + tr3 = tr1 - tr3 + ti3 = ti1 - ti3 + t4 = tr4 + tr6 + t5 = ti4 + ti6 + t6 = tr5 + tr7 + t7 = ti5 + ti7 + ttr6 = ti4 - ti6 + ti6 = tr6 - tr4 + ttr7 = ti5 - ti7 + ti7 = tr7 - tr5 + br0(j) = t0 + t4 + bi0(k) = t1 + t5 + br1(j) = c1*(t2+t6) - s1*(t3+t7) + bi1(k) = c1*(t3+t7) + s1*(t2+t6) + br2(j) = c2*(tr2+ttr6) - s2*(ti2+ti6) + bi2(k) = c2*(ti2+ti6) + s2*(tr2+ttr6) + br3(j) = c3*(tr3+ttr7) - s3*(ti3+ti7) + bi3(k) = c3*(ti3+ti7) + s3*(tr3+ttr7) + br4(j) = c4*(t0-t4) - s4*(t1-t5) + bi4(k) = c4*(t1-t5) + s4*(t0-t4) + br5(j) = c5*(t2-t6) - s5*(t3-t7) + bi5(k) = c5*(t3-t7) + s5*(t2-t6) + br6(j) = c6*(tr2-ttr6) - s6*(ti2-ti6) + bi6(k) = c6*(ti2-ti6) + s6*(tr2-ttr6) + br7(j) = c7*(tr3-ttr7) - s7*(ti3-ti7) + bi7(k) = c7*(ti3-ti7) + s7*(tr3-ttr7) + 100 continue + jr = jr + 2 + ji = ji - 2 + if (ji-jl) 110, 110, 120 + 110 ji = 2*jr - 1 + jl = jr + 120 continue + return + end diff --git a/math/ieee/chap1/r8tr.f b/math/ieee/chap1/r8tr.f new file mode 100644 index 00000000..d49ef58c --- /dev/null +++ b/math/ieee/chap1/r8tr.f @@ -0,0 +1,201 @@ +c +c----------------------------------------------------------------------- +c subroutine: r8tr +c radix 8 iteration subroutine +c----------------------------------------------------------------------- +c + subroutine r8tr(int, nn, br0, br1, br2, br3, br4, br5, br6, br7, + * bi0, bi1, bi2, bi3, bi4, bi5, bi6, bi7) + dimension l(15), br0(2), br1(2), br2(2), br3(2), br4(2), br5(2), + * br6(2), br7(2), bi0(2), bi1(2), bi2(2), bi3(2), bi4(2), + * bi5(2), bi6(2), bi7(2) + common /con/ pii, p7, p7two, c22, s22, pi2 + 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 set up counters such that jthet steps through the arguments +c of w, jr steps through starting locations for the real part of the +c intermediate results and ji steps through starting locations +c of the imaginary part of the intermediate results. +c + l(1) = nn/8 + do 40 k=2,15 + if (l(k-1)-2) 10, 20, 30 + 10 l(k-1) = 2 + 20 l(k) = 2 + go to 40 + 30 l(k) = l(k-1)/2 + 40 continue + piovn = pii/float(nn) + ji = 3 + jl = 2 + jr = 2 + do 120 j1=2,l1,2 + do 120 j2=j1,l2,l1 + do 120 j3=j2,l3,l2 + do 120 j4=j3,l4,l3 + do 120 j5=j4,l5,l4 + do 120 j6=j5,l6,l5 + do 120 j7=j6,l7,l6 + do 120 j8=j7,l8,l7 + do 120 j9=j8,l9,l8 + do 120 j10=j9,l10,l9 + do 120 j11=j10,l11,l10 + do 120 j12=j11,l12,l11 + do 120 j13=j12,l13,l12 + do 120 j14=j13,l14,l13 + do 120 jthet=j14,l15,l14 + th2 = jthet - 2 + if (th2) 50, 50, 90 + 50 do 60 k=1,int + t0 = br0(k) + br4(k) + t1 = br1(k) + br5(k) + t2 = br2(k) + br6(k) + t3 = br3(k) + br7(k) + t4 = br0(k) - br4(k) + t5 = br1(k) - br5(k) + t6 = br2(k) - br6(k) + t7 = br3(k) - br7(k) + br2(k) = t0 - t2 + br3(k) = t1 - t3 + t0 = t0 + t2 + t1 = t1 + t3 + br0(k) = t0 + t1 + br1(k) = t0 - t1 + pr = p7*(t5-t7) + pi = p7*(t5+t7) + br4(k) = t4 + pr + br7(k) = t6 + pi + br6(k) = t4 - pr + br5(k) = pi - t6 + 60 continue + if (nn-8) 120, 120, 70 + 70 k0 = int*8 + 1 + kl = k0 + int - 1 + do 80 k=k0,kl + pr = p7*(bi2(k)-bi6(k)) + pi = p7*(bi2(k)+bi6(k)) + tr0 = bi0(k) + pr + ti0 = bi4(k) + pi + tr2 = bi0(k) - pr + ti2 = bi4(k) - pi + pr = p7*(bi3(k)-bi7(k)) + pi = p7*(bi3(k)+bi7(k)) + tr1 = bi1(k) + pr + ti1 = bi5(k) + pi + tr3 = bi1(k) - pr + ti3 = bi5(k) - pi + pr = tr1*c22 - ti1*s22 + pi = ti1*c22 + tr1*s22 + bi0(k) = tr0 + pr + bi6(k) = tr0 - pr + bi7(k) = ti0 + pi + bi1(k) = pi - ti0 + pr = -tr3*s22 - ti3*c22 + pi = tr3*c22 - ti3*s22 + bi2(k) = tr2 + pr + bi4(k) = tr2 - pr + bi5(k) = ti2 + pi + bi3(k) = pi - ti2 + 80 continue + go to 120 + 90 arg = th2*piovn + 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 + int8 = int*8 + j0 = jr*int8 + 1 + k0 = ji*int8 + 1 + jlast = j0 + int - 1 + do 100 j=j0,jlast + k = k0 + j - j0 + tr1 = br1(j)*c1 - bi1(k)*s1 + ti1 = br1(j)*s1 + bi1(k)*c1 + tr2 = br2(j)*c2 - bi2(k)*s2 + ti2 = br2(j)*s2 + bi2(k)*c2 + tr3 = br3(j)*c3 - bi3(k)*s3 + ti3 = br3(j)*s3 + bi3(k)*c3 + tr4 = br4(j)*c4 - bi4(k)*s4 + ti4 = br4(j)*s4 + bi4(k)*c4 + tr5 = br5(j)*c5 - bi5(k)*s5 + ti5 = br5(j)*s5 + bi5(k)*c5 + tr6 = br6(j)*c6 - bi6(k)*s6 + ti6 = br6(j)*s6 + bi6(k)*c6 + tr7 = br7(j)*c7 - bi7(k)*s7 + ti7 = br7(j)*s7 + bi7(k)*c7 +c + t0 = br0(j) + tr4 + t1 = bi0(k) + ti4 + tr4 = br0(j) - tr4 + ti4 = bi0(k) - ti4 + t2 = tr1 + tr5 + t3 = ti1 + ti5 + tr5 = tr1 - tr5 + ti5 = ti1 - ti5 + t4 = tr2 + tr6 + t5 = ti2 + ti6 + tr6 = tr2 - tr6 + ti6 = ti2 - ti6 + t6 = tr3 + tr7 + t7 = ti3 + ti7 + tr7 = tr3 - tr7 + ti7 = ti3 - ti7 +c + tr0 = t0 + t4 + ti0 = t1 + t5 + tr2 = t0 - t4 + ti2 = t1 - t5 + tr1 = t2 + t6 + ti1 = t3 + t7 + tr3 = t2 - t6 + ti3 = t3 - t7 + t0 = tr4 - ti6 + t1 = ti4 + tr6 + t4 = tr4 + ti6 + t5 = ti4 - tr6 + t2 = tr5 - ti7 + t3 = ti5 + tr7 + t6 = tr5 + ti7 + t7 = ti5 - tr7 + br0(j) = tr0 + tr1 + bi7(k) = ti0 + ti1 + bi6(k) = tr0 - tr1 + br1(j) = ti1 - ti0 + br2(j) = tr2 - ti3 + bi5(k) = ti2 + tr3 + bi4(k) = tr2 + ti3 + br3(j) = tr3 - ti2 + pr = p7*(t2-t3) + pi = p7*(t2+t3) + br4(j) = t0 + pr + bi3(k) = t1 + pi + bi2(k) = t0 - pr + br5(j) = pi - t1 + pr = -p7*(t6+t7) + pi = p7*(t6-t7) + br6(j) = t4 + pr + bi1(k) = t5 + pi + bi0(k) = t4 - pr + br7(j) = pi - t5 + 100 continue + jr = jr + 2 + ji = ji - 2 + if (ji-jl) 110, 110, 120 + 110 ji = 2*jr - 1 + jl = jr + 120 continue + return + end diff --git a/math/ieee/chap1/r8tx.f b/math/ieee/chap1/r8tx.f new file mode 100644 index 00000000..9cc4f591 --- /dev/null +++ b/math/ieee/chap1/r8tx.f @@ -0,0 +1,107 @@ +c +c----------------------------------------------------------------------- +c subroutine: r8tx +c radix 8 iteration subroutine +c----------------------------------------------------------------------- +c + subroutine r8tx(nxtlt, nthpo, lengt, cr0, cr1, cr2, cr3, cr4, + * cr5, cr6, cr7, ci0, ci1, ci2, ci3, ci4, ci5, ci6, ci7) + dimension cr0(2), cr1(2), cr2(2), cr3(2), cr4(2), cr5(2), cr6(2), + * cr7(2), ci1(2), ci2(2), ci3(2), ci4(2), ci5(2), ci6(2), + * ci7(2), ci0(2) + 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 k=j,nthpo,lengt + ar0 = cr0(k) + cr4(k) + ar1 = cr1(k) + cr5(k) + ar2 = cr2(k) + cr6(k) + ar3 = cr3(k) + cr7(k) + ar4 = cr0(k) - cr4(k) + ar5 = cr1(k) - cr5(k) + ar6 = cr2(k) - cr6(k) + ar7 = cr3(k) - cr7(k) + ai0 = ci0(k) + ci4(k) + ai1 = ci1(k) + ci5(k) + ai2 = ci2(k) + ci6(k) + ai3 = ci3(k) + ci7(k) + ai4 = ci0(k) - ci4(k) + ai5 = ci1(k) - ci5(k) + ai6 = ci2(k) - ci6(k) + ai7 = ci3(k) - ci7(k) + 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 + cr0(k) = br0 + br1 + ci0(k) = bi0 + bi1 + if (j.le.1) go to 10 + cr1(k) = c4*(br0-br1) - s4*(bi0-bi1) + ci1(k) = c4*(bi0-bi1) + s4*(br0-br1) + cr2(k) = c2*(br2-bi3) - s2*(bi2+br3) + ci2(k) = c2*(bi2+br3) + s2*(br2-bi3) + cr3(k) = c6*(br2+bi3) - s6*(bi2-br3) + ci3(k) = c6*(bi2-br3) + s6*(br2+bi3) + tr = p7*(br5-bi5) + ti = p7*(br5+bi5) + cr4(k) = c1*(br4+tr) - s1*(bi4+ti) + ci4(k) = c1*(bi4+ti) + s1*(br4+tr) + cr5(k) = c5*(br4-tr) - s5*(bi4-ti) + ci5(k) = c5*(bi4-ti) + s5*(br4-tr) + tr = -p7*(br7+bi7) + ti = p7*(br7-bi7) + cr6(k) = c3*(br6+tr) - s3*(bi6+ti) + ci6(k) = c3*(bi6+ti) + s3*(br6+tr) + cr7(k) = c7*(br6-tr) - s7*(bi6-ti) + ci7(k) = c7*(bi6-ti) + s7*(br6-tr) + go to 20 + 10 cr1(k) = br0 - br1 + ci1(k) = bi0 - bi1 + cr2(k) = br2 - bi3 + ci2(k) = bi2 + br3 + cr3(k) = br2 + bi3 + ci3(k) = bi2 - br3 + tr = p7*(br5-bi5) + ti = p7*(br5+bi5) + cr4(k) = br4 + tr + ci4(k) = bi4 + ti + cr5(k) = br4 - tr + ci5(k) = bi4 - ti + tr = -p7*(br7+bi7) + ti = p7*(br7-bi7) + cr6(k) = br6 + tr + ci6(k) = bi6 + ti + cr7(k) = br6 - tr + ci7(k) = bi6 - ti + 20 continue + 30 continue + return + end diff --git a/math/ieee/chap1/rad4sb.f b/math/ieee/chap1/rad4sb.f new file mode 100644 index 00000000..dfebf0cc --- /dev/null +++ b/math/ieee/chap1/rad4sb.f @@ -0,0 +1,38 @@ +c +c----------------------------------------------------------------------- +c subroutine: rad4sb +c used by subroutine radix4. never directly accessed by user. +c----------------------------------------------------------------------- +c + subroutine rad4sb(ntype) +c +c input: ntype = type of butterfly invoked +c output: parameters used by subroutine radix4 +c + dimension ix(2996) + common /xx/ix + common ntypl,kkp,index,ixc + if(ntype.eq.ntypl) go to 7 + ix(ixc)=0 + ix(ixc+1)=ntype + ixc=ixc+2 + if(ntype.ne.4) go to 4 + indexp=(index-1)*9 + ix(ixc)=kkp+1 + ix(ixc+1)=indexp+1 + ixc=ixc+2 + go to 6 +4 ix(ixc)=kkp+1 + ixc=ixc+1 +6 ntypl=ntype + return +7 if(ntype.ne.4) go to 8 + indexp=(index-1)*9 + ix(ixc)=kkp+1 + ix(ixc+1)=indexp+1 + ixc=ixc+2 + return +8 ix(ixc)=kkp+1 + ixc=ixc+1 + return + end diff --git a/math/ieee/chap1/radix4.f b/math/ieee/chap1/radix4.f new file mode 100644 index 00000000..cd703305 --- /dev/null +++ b/math/ieee/chap1/radix4.f @@ -0,0 +1,488 @@ +c +c----------------------------------------------------------------------- +c subroutine: radix4 +c computes forward or inverse complex dft via radix-4 fft. +c uses autogen technique to yield time efficient program. +c----------------------------------------------------------------------- +c + subroutine radix4(mm,iflag,jflag) +c +c input: +c mm = power of 4 (i.e., n = 4**mm complex point transform) +c (mm.ge.2 and mm.le.5) +c +c iflag = 1 on first pass for given n +c = 0 on subsequent passes for given n +c +c jflag = -1 for forward transform +c = +1 for inverse transform +c +c input/output: +c a = array of dimensions 2*n with real and imaginary parts +c of dft input/output in odd, even array components. +c +c for optimal time efficiency, common is used to pass arrays. +c this means that dimensions of arrays a, ix, and t can be +c modified to reflect maximum value of n = 4**mm to be used. note +c that array "ix" is also dimensioned in subroutine "rad4sb". +c +c i.e., a( ) ix( ) t( ) +c +c m =2 32 38 27 +c m<=3 128 144 135 +c m<=4 512 658 567 +c m<=5 2048 2996 2295 +c + dimension a(2048),ix(2996),t(2295) + dimension nfac(11),np(209) + common ntypl,kkp,index,ixc + common /aa/a + common /xx/ix +c +c check for mm<2 or mm>5 +c + if(mm.lt.2.or.mm.gt.5)stop +c +c initialize on first pass """""""""""""""""""""""""""""""""""""""" +c + if(iflag.eq.1) go to 9999 +c +c fast fourier transform start #################################### +c +8885 kspan=2*4**mm + if(jflag.eq.1) go to 8887 +c +c conjugate data for forward transform +c + do 8886 j=2,n2,2 +8886 a(j)=-a(j) + go to 8889 +c +c multiply data by n**(-1) if inverse transform +c +8887 do 8888 j=1,n2,2 + a(j)=a(j)*xp +8888 a(j+1)=a(j+1)*xp +8889 i=3 + it=ix(i-1) + go to (1,2,3,4,5,6,7,8),it +c*********************************************************************** +c +c 8 multiply butterfly +c +1 kk=ix(i) +c +11 k1=kk+kspan + k2=k1+kspan + k3=k2+kspan +c + akp=a(kk)+a(k2) + akm=a(kk)-a(k2) + ajp=a(k1)+a(k3) + ajm=a(k1)-a(k3) + a(kk)=akp+ajp +c + bkp=a(kk+1)+a(k2+1) + bkm=a(kk+1)-a(k2+1) + bjp=a(k1+1)+a(k3+1) + bjm=a(k1+1)-a(k3+1) + a(kk+1)=bkp+bjp +c + bjp=bkp-bjp +c + a(k2+1)=(akp+bjp-ajp)*c707 + a(k2)=a(k2+1)+bjp*cm141 +c + bkp=bkm+ajm + akp=akm-bjm +c + ac0=(akp+bkp)*c924 + a(k1+1)=ac0+akp*cm541 + a(k1) =ac0+bkp*cm131 +c + bkm=bkm-ajm + akm=akm+bjm +c + ac0=(akm+bkm)*c383 + a(k3+1)=ac0+akm*c541 + a(k3) =ac0+bkm*cm131 +c + i=i+1 + kk=ix(i) + if (kk) 111,111,11 +111 i=i+2 + it=ix(i-1) + go to (1,2,3,4,5,6,7,8), it +c*********************************************************************** +c +c 4 multiply butterfly +c +2 kk=ix(i) +c +22 k1=kk+kspan + k2=k1+kspan + k3=k2+kspan +c + akp=a(kk)+a(k2) + akm=a(kk)-a(k2) + ajp=a(k1)+a(k3) + ajm=a(k1)-a(k3) + a(kk)=akp+ajp +c + bkp=a(kk+1)+a(k2+1) + bkm=a(kk+1)-a(k2+1) + bjp=a(k1+1)+a(k3+1) + bjm=a(k1+1)-a(k3+1) + a(kk+1)=bkp+bjp + a(k2)=-bkp+bjp + a(k2+1)=akp-ajp +c + bkp=bkm+ajm +c + a(k1+1)=(bkp+akm-bjm)*c707 + a(k1)=a(k1+1)+bkp*cm141 +c + akm=akm+bjm +c + a(k3+1)=(akm+ajm-bkm)*c707 + a(k3)=a(k3+1)+akm*cm141 +c + i=i+1 + kk=ix(i) + if (kk) 222,222,22 +222 i=i+2 + it=ix(i-1) + go to (1,2,3,4,5,6,7,8), it +c*********************************************************************** +c +c 8 multiply butterfly +c +3 kk=ix(i) +c +33 k1=kk+kspan + k2=k1+kspan + k3=k2+kspan +c + akp=a(kk)+a(k2) + akm=a(kk)-a(k2) + ajp=a(k1)+a(k3) + ajm=a(k1)-a(k3) + a(kk)=akp+ajp +c + bkp=a(kk+1)+a(k2+1) + bkm=a(kk+1)-a(k2+1) + bjp=a(k1+1)+a(k3+1) + bjm=a(k1+1)-a(k3+1) + a(kk+1)=bkp+bjp +c + ajp=akp-ajp +c + a(k2+1)=(ajp+bjp-bkp)*c707 + a(k2)=a(k2+1)+ajp*cm141 +c + bkp=bkm+ajm + akp=akm-bjm +c + ac0=(akp+bkp)*c383 + a(k1+1)=ac0+akp*c541 + a(k1) =ac0+bkp*cm131 +c + bkm=bkm-ajm + akm=akm+bjm +c + ac0=(akm+bkm)*cm924 + a(k3+1)=ac0+akm*c541 + a(k3) =ac0+bkm*c131 +c + i=i+1 + kk=ix(i) + if (kk) 333,333,33 +333 i=i+2 + it=ix(i-1) + go to (1,2,3,4,5,6,7,8), it +c*********************************************************************** +c +c general 9 multiply butterfly +c +4 kk=ix(i) +c +44 k1=kk+kspan + k2=k1+kspan + k3=k2+kspan +c + akp=a(kk)+a(k2) + akm=a(kk)-a(k2) + ajp=a(k1)+a(k3) + ajm=a(k1)-a(k3) + a(kk)=akp+ajp +c + bkp=a(kk+1)+a(k2+1) + bkm=a(kk+1)-a(k2+1) + bjp=a(k1+1)+a(k3+1) + bjm=a(k1+1)-a(k3+1) + a(kk+1)=bkp+bjp +c + ajp=akp-ajp + bjp=bkp-bjp +c + j=ix(i+1) +c + ac0=(ajp+bjp)*t(j+8) + a(k2+1)=ac0+ajp*t(j+6) + a(k2) =ac0+bjp*t(j+7) +c + bkp=bkm+ajm + akp=akm-bjm +c + ac0=(akp+bkp)*t(j+5) + a(k1+1)=ac0+akp*t(j+3) + a(k1) =ac0+bkp*t(j+4) +c + bkm=bkm-ajm + akm=akm+bjm +c + ac0=(akm+bkm)*t(j+2) + a(k3+1)=ac0+akm*t(j) + a(k3) =ac0+bkm*t(j+1) +c + i=i+2 + kk=ix(i) + if (kk) 444,444,44 +444 i=i+2 + it=ix(i-1) + go to (1,2,3,4,5,6,7,8), it +c*********************************************************************** +c +c 0 multiply butterfly +c +5 kk=ix(i) +c +55 k1=kk+kspan + k2=k1+kspan + k3=k2+kspan +c + akp=a(kk)+a(k2) + akm=a(kk)-a(k2) + ajp=a(k1)+a(k3) + ajm=a(k1)-a(k3) + a(kk)=akp+ajp + a(k2)=akp-ajp +c + bkp=a(kk+1)+a(k2+1) + bkm=a(kk+1)-a(k2+1) + bjp=a(k1+1)+a(k3+1) + bjm=a(k1+1)-a(k3+1) + a(kk+1)=bkp+bjp + a(k2+1)=bkp-bjp +c + a(k3+1)=bkm-ajm + a(k1+1)=bkm+ajm + a(k3)=akm+bjm + a(k1)=akm-bjm +c + i=i+1 + kk=ix(i) + if (kk) 555,555,55 +555 i=i+2 + it=ix(i-1) + go to (1,2,3,4,5,6,7,8), it +c*********************************************************************** +c +c offset reduced +c +6 kspan=kspan/4 + i=i+2 + it=ix(i-1) + go to (1,2,3,4,5,6,7,8), it +c*********************************************************************** +c +c bit reversal (shuffling) +c +7 ip1=ix(i) +77 ip2=ix(i+1) + t1=a(ip2) + a(ip2)=a(ip1) + a(ip1)=t1 + t1=a(ip2+1) + a(ip2+1)=a(ip1+1) + a(ip1+1)=t1 + i=i+2 + ip1=ix(i) + if (ip1) 777,777,77 +777 i=i+2 + it=ix(i-1) + go to (1,2,3,4,5,6,7,8), it +c*********************************************************************** +8 if(jflag.eq.1) go to 888 +c +c conjugate output if forward transform +c + do 88 j=2,n2,2 +88 a(j)=-a(j) +888 return +c +c fast fourier transform ends ##################################### +c +c initialization phase starts. done only once +c +9999 ixc=1 + n=4**mm + xp=n + xp=1./xp + ntot=n + n2=n*2 + nspan=n + n1test=n/16 + n2test=n/8 + n3test=(3*n)/16 + nspan4=nspan/4 + ibase=0 + isn=1 + inc=isn + rad=8.0*atan(1.0) + pi=4.*atan(1.0) + c707=sin(pi/4.) + cm141=-2.*c707 + c383=sin(pi/8.) + c924=cos(pi/8.) + cm924=-c924 + c541=c924-c383 + cm541=-c541 + c131=c924+c383 + cm131=-c131 +10 nt=inc*ntot + ks=inc*nspan + kspan=ks + jc=ks/n + radf=rad*float(jc)*.5 + i=0 +c +c determine the factors of n +c all factors must be 4 for this version +c + m=0 + k=n +15 m=m+1 + nfac(m)=4 + k=k/4 +20 if(k-(k/4)*4.eq.0) go to 15 + kt=1 + if(n.ge.256) kt=2 + kspan0=kspan + ntypl=0 +c +100 ndelta=kspan0/kspan + index=0 + sd=radf/float(kspan) + cd=2.0*sin(sd)**2 + sd=sin(sd+sd) + kk=1 + i=i+1 +c +c transform for a factor of 4 +c + kspan=kspan/4 + ix(ixc)=0 + ix(ixc+1)=6 + ixc=ixc+2 +c +410 c1=1.0 + s1=0.0 +420 k1=kk+kspan + k2=k1+kspan + k3=k2+kspan + if(s1.eq.0.0) go to 460 +430 if(kspan.ne.nspan4) go to 431 + t(ibase+5)=-(s1+c1) + t(ibase+6)=c1 + t(ibase+4)=s1-c1 + t(ibase+8)=-(s2+c2) + t(ibase+9)=c2 + t(ibase+7)=s2-c2 + t(ibase+2)=-(s3+c3) + t(ibase+3)=c3 + t(ibase+1)=s3-c3 + ibase=ibase+9 +c +431 kkp=(kk-1)*2 + if(index.ne.n1test) go to 150 + call rad4sb(1) + go to 5035 +150 if(index.ne.n2test) go to 160 + call rad4sb(2) + go to 5035 +160 if(index.ne.n3test) go to 170 + call rad4sb(3) + go to 5035 +170 call rad4sb(4) +5035 kk=k3+kspan + if(kk.le.nt) go to 420 +440 index=index+ndelta + c2=c1-(cd*c1+sd*s1) + s1=(sd*c1-cd*s1)+s1 + c1=c2 + c2=c1*c1-s1*s1 + s2=c1*s1+c1*s1 + c3=c2*c1-s2*s1 + s3=c2*s1+s2*c1 + kk=kk-nt+jc + if(kk.le.kspan) go to 420 + kk=kk-kspan+inc + if(kk.le.jc) go to 410 + if(kspan.eq.jc) go to 800 + go to 100 +460 kkp=(kk-1)*2 + call rad4sb(5) +5050 kk=k3+kspan + if(kk.le.nt) go to 420 + go to 440 +c +800 ix(ixc)=0 + ix(ixc+1)=7 + ixc=ixc+2 +c +c compute parameters to permute the results to normal order +c done in two steps +c permutation for square factors of n +c + np(1)=ks + k=kt+kt+1 + if(m.lt.k) k=k-1 + j=1 + np(k+1)=jc +810 np(j+1)=np(j)/nfac(j) + np(k)=np(k+1)*nfac(j) + j=j+1 + k=k-1 + if(j.lt.k) go to 810 + k3=np(k+1) + kspan=np(2) + kk=jc+1 + k2=kspan+1 + j=1 +c +c permutation for single variate transform +c +820 kkp=(kk-1)*2 + k2p=(k2-1)*2 + ix(ixc)=kkp+1 + ix(ixc+1)=k2p+1 + ixc=ixc+2 + kk=kk+inc + k2=kspan+k2 + if(k2.lt.ks) go to 820 +830 k2=k2-np(j) + j=j+1 + k2=np(j+1)+k2 + if(k2.gt.np(j)) go to 830 + j=1 +840 if(kk.lt.k2) go to 820 + kk=kk+inc + k2=kspan+k2 + if(k2.lt.ks) go to 840 + if(kk.lt.ks) go to 830 + jc=k3 + ix(ixc)=0 + ix(ixc+1)=8 + go to 8885 + end 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 diff --git a/math/ieee/chap1/time/time12.f b/math/ieee/chap1/time/time12.f new file mode 100644 index 00000000..774c2cc9 --- /dev/null +++ b/math/ieee/chap1/time/time12.f @@ -0,0 +1,53 @@ +c----------------------------------------------------------------------- +c +c----------------------------------------------------------------------- +c + parameter (SIZE = 1024, ILOOP = 100) + complex a, w + real breal(SIZE), bimag(SIZE), qbreal(SIZE), qbimag(SIZE) +c + ioutd = i1mach(2) + nn = SIZE + tpi = 8.*atan(1.) + tpion = tpi/float(nn) + w = cmplx(cos(tpion),-sin(tpion)) +c +c generate a**k as test function +c result to b(i) for modification by dft and idft subroutines and +c a copy to qb(i) to compare final result with for error difference. +c + a = (.9,.3) + breal(1) = 1.0 + bimag(1) = 0.0 + qbreal(1) = 1.0 + qbimag(1) = 0.0 + do 10 k=2,nn + w = a**(k-1) + breal(k) = real(w) + bimag(k) = aimag(w) + qbreal(k) = breal(k) + qbimag(k) = bimag(k) + 10 continue +c +c now compute dft, idft, dft, idft, ... +c first dft is computed specially, in case subroutine needs to be started. +c + call fft842(0, SIZE, breal, bimag) + do 25 icount = 1, ILOOP + call fft842(1, SIZE, breal, bimag) + call fft842(0, SIZE, breal, bimag) + 25 continue + call fft842(1, SIZE, breal, bimag) +c +c calculate rms error between b(i) and qb(i). +c + err = 0. + do 30 i=1,nn + err = err + (breal(i)-qbreal(i))**2 + * + (bimag(i)-qbimag(i))**2 + 30 continue + err = sqrt(err / float(nn)) + write (ioutd,9994) ILOOP, err + 9994 format(' rms error, after ', i5, ' loops = ', e15.8) + stop + end diff --git a/math/ieee/chap1/time/time17.f b/math/ieee/chap1/time/time17.f new file mode 100644 index 00000000..adcce4c2 --- /dev/null +++ b/math/ieee/chap1/time/time17.f @@ -0,0 +1,53 @@ +c----------------------------------------------------------------------- +c +c----------------------------------------------------------------------- +c + parameter (SIZE = 1008, ILOOP = 100) + complex a, w + real breal(SIZE), bimag(SIZE), qbreal(SIZE), qbimag(SIZE) +c + ioutd = i1mach(2) + nn = SIZE + tpi = 8.*atan(1.) + tpion = tpi/float(nn) + w = cmplx(cos(tpion),-sin(tpion)) +c +c generate a**k as test function +c result to b(i) for modification by dft and idft subroutines and +c a copy to qb(i) to compare final result with for error difference. +c + a = (.9,.3) + breal(1) = 1.0 + bimag(1) = 0.0 + qbreal(1) = 1.0 + qbimag(1) = 0.0 + do 10 k=2,nn + w = a**(k-1) + breal(k) = real(w) + bimag(k) = aimag(w) + qbreal(k) = breal(k) + qbimag(k) = bimag(k) + 10 continue +c +c now compute dft, idft, dft, idft, ... +c first dft is computed specially, in case subroutine needs to be started. +c + call wfta(breal, bimag, SIZE, 0, 0, ierr) + do 25 icount = 1, ILOOP + call wfta(breal, bimag, SIZE, 1, 1, ierr) + call wfta(breal, bimag, SIZE, 0, 1, ierr) + 25 continue + call wfta(breal, bimag, SIZE, 1, 1, ierr) +c +c calculate rms error between b(i) and qb(i). +c + err = 0. + do 30 i=1,nn + err = err + (breal(i)-qbreal(i))**2 + * + (bimag(i)-qbimag(i))**2 + 30 continue + err = sqrt(err / float(nn)) + write (ioutd,9994) ILOOP, err + 9994 format(' rms error, after ', i5, ' loops = ', e15.8) + stop + end diff --git a/math/ieee/chap1/time/time18.f b/math/ieee/chap1/time/time18.f new file mode 100644 index 00000000..87680554 --- /dev/null +++ b/math/ieee/chap1/time/time18.f @@ -0,0 +1,48 @@ +c----------------------------------------------------------------------- +c +c----------------------------------------------------------------------- +c + parameter (SIZE = 1024, ILOOP = 100) + complex w, b(SIZE), qb(SIZE), a + common /aa/ b +c + ioutd = i1mach(2) + nn = SIZE + tpi = 8.*atan(1.) + tpion = tpi/float(nn) + w = cmplx(cos(tpion),-sin(tpion)) +c +c generate a**k as test function +c result to b(i) for modification by dft and idft subroutines and +c a copy to qb(i) to compare final result with for error difference. +c + a = (.9,.3) + b(1) = (1.,0.) + qb(1) = b(1) + do 10 k=2,nn + b(k) = a**(k-1) + qb(k) = b(k) + 10 continue +c +c now compute dft, idft, dft, idft, ... +c first dft is computed specially, in case subroutine needs to be started. +c + call radix4(5, 1, -1) + do 25 icount = 1, ILOOP + call radix4(5, 0, 1) + call radix4(5, 0, -1) + 25 continue + call radix4(5, 0, 1) +c +c calculate rms error between b(i) and qb(i). +c + err = 0. + do 30 i=1,nn + de = cabs(qb(i)-b(i)) + err = err + de**2 + 30 continue + err = sqrt(err / float(nn)) + write (ioutd,9994) ILOOP, err + 9994 format(' rms error, after ', i5, ' loops = ', e15.8) + stop + end diff --git a/math/ieee/chap1/weave1.f b/math/ieee/chap1/weave1.f new file mode 100644 index 00000000..9c83d0ea --- /dev/null +++ b/math/ieee/chap1/weave1.f @@ -0,0 +1,371 @@ +c +c----------------------------------------------------------------------- +c subroutine: weave1 +c this subroutine implements the different pre-weave +c modules of the wfta. the working arrays are sr and si. +c the routine checks to see which factors are present +c in the transform length n = na*nb*nc*nd and executes +c the pre-weave code for these factors. +c +c----------------------------------------------------------------------- +c + subroutine weave1(sr,si) + common na,nb,nc,nd,nd1,nd2,nd3,nd4 + dimension q(8),t(16) + dimension sr(1),si(1) + if(na.eq.1) go to 300 + if(na.ne.2) go to 800 +c +c ********************************************************************** +c +c the following code implements the 2 point pre-weave module +c +c ********************************************************************** +c + nlup2=2*(nd2-nb) + nlup23=2*nd2*(nd3-nc) + nbase=1 + do 240 n4=1,nd + do 230 n3=1,nc + do 220 n2=1,nb + nr1=nbase+1 + t0=sr(nbase)+sr(nr1) + sr(nr1)=sr(nbase)-sr(nr1) + sr(nbase)=t0 + t0=si(nbase)+si(nr1) + si(nr1)=si(nbase)-si(nr1) + si(nbase)=t0 +220 nbase=nbase+2 +230 nbase=nbase+nlup2 +240 nbase=nbase+nlup23 +800 if(na.ne.8) go to 1600 +c +c ********************************************************************** +c +c the following code implements the 8 point pre-weave module +c +c ********************************************************************** +c + nlup2=8*(nd2-nb) + nlup23=8*nd2*(nd3-nc) + nbase=1 + do 840 n4=1,nd + do 830 n3=1,nc + do 820 n2=1,nb + nr1=nbase+1 + nr2=nr1+1 + nr3=nr2+1 + nr4=nr3+1 + nr5=nr4+1 + nr6=nr5+1 + nr7=nr6+1 + t3=sr(nr3)+sr(nr7) + t7=sr(nr3)-sr(nr7) + t0=sr(nbase)+sr(nr4) + sr(nr4)=sr(nbase)-sr(nr4) + t1=sr(nr1)+sr(nr5) + t5=sr(nr1)-sr(nr5) + t2=sr(nr2)+sr(nr6) + sr(nr6)=sr(nr2)-sr(nr6) + sr(nbase)=t0+t2 + sr(nr2)=t0-t2 + sr(nr1)=t1+t3 + sr(nr3)=t1-t3 + sr(nr5)=t5+t7 + sr(nr7)=t5-t7 + t3=si(nr3)+si(nr7) + t7=si(nr3)-si(nr7) + t0=si(nbase)+si(nr4) + si(nr4)=si(nbase)-si(nr4) + t1=si(nr1)+si(nr5) + t5=si(nr1)-si(nr5) + t2=si(nr2)+si(nr6) + si(nr6)=si(nr2)-si(nr6) + si(nbase)=t0+t2 + si(nr2)=t0-t2 + si(nr1)=t1+t3 + si(nr3)=t1-t3 + si(nr5)=t5+t7 + si(nr7)=t5-t7 +820 nbase=nbase+8 +830 nbase=nbase+nlup2 +840 nbase=nbase+nlup23 +1600 if(na.ne.16) go to 300 +c +c ********************************************************************** +c +c the following code implements the 16 point pre-weave module +c +c ********************************************************************** +c + nlup2=18*(nd2-nb) + nlup23=18*nd2*(nd3-nc) + nbase=1 + do 1640 n4=1,nd + do 1630 n3=1,nc + do 1620 n2=1,nb + nr1=nbase+1 + nr2=nr1+1 + nr3=nr2+1 + nr4=nr3+1 + nr5=nr4+1 + nr6=nr5+1 + nr7=nr6+1 + nr8=nr7+1 + nr9=nr8+1 + nr10=nr9+1 + nr11=nr10+1 + nr12=nr11+1 + nr13=nr12+1 + nr14=nr13+1 + nr15=nr14+1 + nr16=nr15+1 + nr17=nr16+1 + jbase=nbase + do 1645 j=1,8 + t(j)=sr(jbase)+sr(jbase+8) + t(j+8)=sr(jbase)-sr(jbase+8) + jbase=jbase+1 +1645 continue + do 1650 j=1,4 + q(j)=t(j)+t(j+4) + q(j+4)=t(j)-t(j+4) +1650 continue + sr(nbase)=q(1)+q(3) + sr(nr2)=q(1)-q(3) + sr(nr1)=q(2)+q(4) + sr(nr3)=q(2)-q(4) + sr(nr5)=q(6)+q(8) + sr(nr7)=q(6)-q(8) + sr(nr4)=q(5) + sr(nr6)=q(7) + sr(nr8)=t(9) + sr(nr9)=t(10)+t(16) + sr(nr15)=t(10)-t(16) + sr(nr13)=t(14)+t(12) + sr(nr11)=t(14)-t(12) + sr(nr17)=sr(nr11)+sr(nr15) + sr(nr16)=sr(nr9)+sr(nr13) + sr(nr10)=t(11)+t(15) + sr(nr14)=t(11)-t(15) + sr(nr12)=t(13) + jbase=nbase + do 1745 j=1,8 + t(j)=si(jbase)+si(jbase+8) + t(j+8)=si(jbase)-si(jbase+8) + jbase=jbase+1 +1745 continue + do 1750 j=1,4 + q(j)=t(j)+t(j+4) + q(j+4)=t(j)-t(j+4) +1750 continue + si(nbase)=q(1)+q(3) + si(nr2)=q(1)-q(3) + si(nr1)=q(2)+q(4) + si(nr3)=q(2)-q(4) + si(nr5)=q(6)+q(8) + si(nr7)=q(6)-q(8) + si(nr4)=q(5) + si(nr6)=q(7) + si(nr8)=t(9) + si(nr9)=t(10)+t(16) + si(nr15)=t(10)-t(16) + si(nr13)=t(14)+t(12) + si(nr11)=t(14)-t(12) + si(nr17)=si(nr11)+si(nr15) + si(nr16)=si(nr9)+si(nr13) + si(nr10)=t(11)+t(15) + si(nr14)=t(11)-t(15) + si(nr12)=t(13) +1620 nbase=nbase+18 +1630 nbase=nbase+nlup2 +1640 nbase=nbase+nlup23 +300 if(nb.eq.1) go to 700 + if(nb.ne.3) go to 900 +c +c ********************************************************************** +c +c the following code implements the 3 point pre-weave module +c +c ********************************************************************** +c + nlup2=2*nd1 + nlup23=3*nd1*(nd3-nc) + nbase=1 + noff=nd1 + do 340 n4=1,nd + do 330 n3=1,nc + do 310 n2=1,nd1 + nr1=nbase+noff + nr2=nr1+noff + t1=sr(nr1)+sr(nr2) + sr(nbase)=sr(nbase)+t1 + sr(nr2)=sr(nr1)-sr(nr2) + sr(nr1)=t1 + t1=si(nr1)+si(nr2) + si(nbase)=si(nbase)+t1 + si(nr2)=si(nr1)-si(nr2) + si(nr1)=t1 +310 nbase=nbase+1 +330 nbase=nbase+nlup2 +340 nbase=nbase+nlup23 +900 if(nb.ne.9) go to 700 +c +c ********************************************************************** +c +c the following code implements the 9 point pre-weave module +c +c ********************************************************************** +c + nlup2=10*nd1 + nlup23=11*nd1*(nd3-nc) + nbase=1 + noff=nd1 + do 940 n4=1,nd + do 930 n3=1,nc + do 910 n2=1,nd1 + nr1=nbase+noff + nr2=nr1+noff + nr3=nr2+noff + nr4=nr3+noff + nr5=nr4+noff + nr6=nr5+noff + nr7=nr6+noff + nr8=nr7+noff + nr9=nr8+noff + nr10=nr9+noff + t3=sr(nr3)+sr(nr6) + t6=sr(nr3)-sr(nr6) + sr(nbase)=sr(nbase)+t3 + t7=sr(nr7)+sr(nr2) + t2=sr(nr7)-sr(nr2) + sr(nr2)=t6 + t1=sr(nr1)+sr(nr8) + t8=sr(nr1)-sr(nr8) + sr(nr1)=t3 + t4=sr(nr4)+sr(nr5) + t5=sr(nr4)-sr(nr5) + sr(nr3)=t1+t4+t7 + sr(nr4)=t1-t7 + sr(nr5)=t4-t1 + sr(nr6)=t7-t4 + sr(nr10)=t2+t5+t8 + sr(nr7)=t8-t2 + sr(nr8)=t5-t8 + sr(nr9)=t2-t5 + t3=si(nr3)+si(nr6) + t6=si(nr3)-si(nr6) + si(nbase)=si(nbase)+t3 + t7=si(nr7)+si(nr2) + t2=si(nr7)-si(nr2) + si(nr2)=t6 + t1=si(nr1)+si(nr8) + t8=si(nr1)-si(nr8) + si(nr1)=t3 + t4=si(nr4)+si(nr5) + t5=si(nr4)-si(nr5) + si(nr3)=t1+t4+t7 + si(nr4)=t1-t7 + si(nr5)=t4-t1 + si(nr6)=t7-t4 + si(nr10)=t2+t5+t8 + si(nr7)=t8-t2 + si(nr8)=t5-t8 + si(nr9)=t2-t5 +910 nbase=nbase+1 +930 nbase=nbase+nlup2 +940 nbase=nbase+nlup23 +700 if(nc.ne.7) go to 500 +c +c ********************************************************************** +c +c the following code implements the 7 point pre-weave module +c +c ********************************************************************** +c + noff=nd1*nd2 + nbase=1 + nlup2=8*noff + do 740 n4=1,nd + do 710 n1=1,noff + nr1=nbase+noff + nr2=nr1+noff + nr3=nr2+noff + nr4=nr3+noff + nr5=nr4+noff + nr6=nr5+noff + nr7=nr6+noff + nr8=nr7+noff + t1=sr(nr1)+sr(nr6) + t6=sr(nr1)-sr(nr6) + t4=sr(nr4)+sr(nr3) + t3=sr(nr4)-sr(nr3) + t2=sr(nr2)+sr(nr5) + t5=sr(nr2)-sr(nr5) + sr(nr5)=t6-t3 + sr(nr2)=t5+t3+t6 + sr(nr6)=t5-t6 + sr(nr8)=t3-t5 + sr(nr3)=t2-t1 + sr(nr4)=t1-t4 + sr(nr7)=t4-t2 + t1=t1+t4+t2 + sr(nbase)=sr(nbase)+t1 + sr(nr1)=t1 + t1=si(nr1)+si(nr6) + t6=si(nr1)-si(nr6) + t4=si(nr4)+si(nr3) + t3=si(nr4)-si(nr3) + t2=si(nr2)+si(nr5) + t5=si(nr2)-si(nr5) + si(nr5)=t6-t3 + si(nr2)=t5+t3+t6 + si(nr6)=t5-t6 + si(nr8)=t3-t5 + si(nr3)=t2-t1 + si(nr4)=t1-t4 + si(nr7)=t4-t2 + t1=t1+t4+t2 + si(nbase)=si(nbase)+t1 + si(nr1)=t1 +710 nbase=nbase+1 +740 nbase=nbase+nlup2 +500 if(nd.ne.5) return +c +c ********************************************************************** +c +c the following code implements the 5 point pre-weave module +c +c ********************************************************************** +c + noff=nd1*nd2*nd3 + nbase=1 + do 510 n1=1,noff + nr1=nbase+noff + nr2=nr1+noff + nr3=nr2+noff + nr4=nr3+noff + nr5=nr4+noff + t4=sr(nr1)-sr(nr4) + t1=sr(nr1)+sr(nr4) + t3=sr(nr3)+sr(nr2) + t2=sr(nr3)-sr(nr2) + sr(nr3)=t1-t3 + sr(nr1)=t1+t3 + sr(nbase)=sr(nbase)+sr(nr1) + sr(nr5)=t2+t4 + sr(nr2)=t4 + sr(nr4)=t2 + t4=si(nr1)-si(nr4) + t1=si(nr1)+si(nr4) + t3=si(nr3)+si(nr2) + t2=si(nr3)-si(nr2) + si(nr3)=t1-t3 + si(nr1)=t1+t3 + si(nbase)=si(nbase)+si(nr1) + si(nr5)=t2+t4 + si(nr2)=t4 + si(nr4)=t2 +510 nbase=nbase+1 + return + end diff --git a/math/ieee/chap1/weave2.f b/math/ieee/chap1/weave2.f new file mode 100644 index 00000000..5c04ff91 --- /dev/null +++ b/math/ieee/chap1/weave2.f @@ -0,0 +1,412 @@ +c +c----------------------------------------------------------------------- +c +c subroutine: weave2 +c this subroutine implements the post-weave modules +c of the wfta. the working arrays are sr and si. +c the routine checks to see which factors are present +c in the transform length n = na*nb*nc*nd and executes +c the post-weave code for these factors. +c +c----------------------------------------------------------------------- +c + subroutine weave2(sr,si) + common na,nb,nc,nd,nd1,nd2,nd3,nd4 + dimension sr(1),si(1) + dimension q(8),t(16) + if(nd.ne.5) go to 700 +c +c ********************************************************************** +c +c the following code implements the 5 point post-weave module +c +c ********************************************************************** +c + noff=nd1*nd2*nd3 + nbase=1 + do 510 n1=1,noff + nr1=nbase+noff + nr2=nr1+noff + nr3=nr2+noff + nr4=nr3+noff + nr5=nr4+noff + t1=sr(nbase)+sr(nr1) + t3=t1-sr(nr3) + t1=t1+sr(nr3) + t4=si(nr2)+si(nr5) + t2=si(nr4)+si(nr5) + sr(nr1)=t1-t4 + sr4=t1+t4 + sr2=t3+t2 + sr(nr3)=t3-t2 + t1=si(nbase)+si(nr1) + t3=t1-si(nr3) + t1=t1+si(nr3) + t4=sr(nr2)+sr(nr5) + t2=sr(nr4)+sr(nr5) + si(nr1)=t1+t4 + si(nr4)=t1-t4 + si(nr2)=t3-t2 + si(nr3)=t3+t2 + sr(nr2)=sr2 + sr(nr4)=sr4 +510 nbase=nbase+1 +700 if(nc.ne.7) go to 300 +c +c ********************************************************************** +c +c the following code implements the 7 point post-weave module +c +c ********************************************************************** +c + noff=nd1*nd2 + nbase=1 + nlup2=8*noff + do 740 n4=1,nd + do 710 n1=1,noff + nr1=nbase+noff + nr2=nr1+noff + nr3=nr2+noff + nr4=nr3+noff + nr5=nr4+noff + nr6=nr5+noff + nr7=nr6+noff + nr8=nr7+noff + t1=sr(nr1)+sr(nbase) + t2=t1-sr(nr3)-sr(nr4) + t4=t1+sr(nr3)-sr(nr7) + t1=t1+sr(nr4)+sr(nr7) + t6=si(nr2)+si(nr5)+si(nr8) + t5=si(nr2)-si(nr5)-si(nr6) + t3=si(nr2)+si(nr6)-si(nr8) + sr(nr1)=t1-t6 + sr6=t1+t6 + sr2=t2-t5 + sr5=t2+t5 + sr(nr4)=t4-t3 + sr(nr3)=t4+t3 + t1=si(nr1)+si(nbase) + t2=t1-si(nr3)-si(nr4) + t4=t1+si(nr3)-si(nr7) + t1=t1+si(nr4)+si(nr7) + t6=sr(nr2)+sr(nr5)+sr(nr8) + t5=sr(nr2)-sr(nr5)-sr(nr6) + t3=sr(nr2)+sr(nr6)-sr(nr8) + si(nr1)=t1+t6 + si(nr6)=t1-t6 + si(nr2)=t2+t5 + si(nr5)=t2-t5 + si(nr4)=t4+t3 + si(nr3)=t4-t3 + sr(nr2)=sr2 + sr(nr5)=sr5 + sr(nr6)=sr6 +710 nbase=nbase+1 +740 nbase=nbase+nlup2 +300 if(nb.eq.1) go to 400 + if(nb.ne.3) go to 900 +c +c ********************************************************************** +c +c the following code implements the 3 point post-weave module +c +c ********************************************************************** +c + nlup2=2*nd1 + nlup23=3*nd1*(nd3-nc) + nbase=1 + noff=nd1 + do 340 n5=1,nd + do 330 n4=1,nc + do 310 n2=1,nd1 + nr1=nbase+noff + nr2=nr1+noff + t1=sr(nbase)+sr(nr1) + sr(nr1)=t1-si(nr2) + sr2=t1+si(nr2) + t1=si(nbase)+si(nr1) + si(nr1)=t1+sr(nr2) + si(nr2)=t1-sr(nr2) + sr(nr2)=sr2 +310 nbase=nbase+1 +330 nbase=nbase+nlup2 +340 nbase=nbase+nlup23 +900 if(nb.ne.9) go to 400 +c +c ********************************************************************** +c +c the following code implements the 9 point post-weave module +c +c ********************************************************************** +c + nlup2=10*nd1 + nlup23=11*nd1*(nd3-nc) + nbase=1 + noff=nd1 + do 940 n4=1,nd + do 930 n3=1,nc + do 910 n2=1,nd1 + nr1=nbase+noff + nr2=nr1+noff + nr3=nr2+noff + nr4=nr3+noff + nr5=nr4+noff + nr6=nr5+noff + nr7=nr6+noff + nr8=nr7+noff + nr9=nr8+noff + nr10=nr9+noff + t3=sr(nbase)-sr(nr3) + t7=sr(nbase)+sr(nr1) + sr(nbase)=sr(nbase)+sr(nr3)+sr(nr3) + t6=t3+si(nr10) + sr(nr3)=t3-si(nr10) + t4=t7+sr(nr5)-sr(nr6) + t1=t7-sr(nr4)-sr(nr5) + t7=t7+sr(nr4)+sr(nr6) + sr(nr6)=t6 + t8=si(nr2)-si(nr7)-si(nr8) + t5=si(nr2)+si(nr8)-si(nr9) + t2=si(nr2)+si(nr7)+si(nr9) + sr(nr1)=t7-t2 + sr8=t7+t2 + sr(nr4)=t1-t8 + sr(nr5)=t1+t8 + sr7=t4-t5 + sr2=t4+t5 + t3=si(nbase)-si(nr3) + t7=si(nbase)+si(nr1) + si(nbase)=si(nbase)+si(nr3)+si(nr3) + t6=t3-sr(nr10) + si(nr3)=t3+sr(nr10) + t4=t7+si(nr5)-si(nr6) + t1=t7-si(nr4)-si(nr5) + t7=t7+si(nr4)+si(nr6) + si(nr6)=t6 + t8=sr(nr2)-sr(nr7)-sr(nr8) + t5=sr(nr2)+sr(nr8)-sr(nr9) + t2=sr(nr2)+sr(nr7)+sr(nr9) + si(nr1)=t7+t2 + si(nr8)=t7-t2 + si(nr4)=t1+t8 + si(nr5)=t1-t8 + si(nr7)=t4+t5 + si(nr2)=t4-t5 + sr(nr2)=sr2 + sr(nr7)=sr7 + sr(nr8)=sr8 +910 nbase=nbase+1 +930 nbase=nbase+nlup2 +940 nbase=nbase+nlup23 +400 if(na.eq.1) return + if(na.ne.4) go to 800 +c +c ********************************************************************** +c +c the following code implements the 4 point post-weave module +c +c ********************************************************************** +c + nlup2=4*(nd2-nb) + nlup23=4*nd2*(nd3-nc) + nbase=1 + do 440 n4=1,nd + do 430 n3=1,nc + do 420 n2=1,nb + nr1=nbase+1 + nr2=nr1+1 + nr3=nr2+1 + tr0=sr(nbase)+sr(nr2) + tr2=sr(nbase)-sr(nr2) + tr1=sr(nr1)+sr(nr3) + tr3=sr(nr1)-sr(nr3) + ti1=si(nr1)+si(nr3) + ti3=si(nr1)-si(nr3) + sr(nbase)=tr0+tr1 + sr(nr2)=tr0-tr1 + sr(nr1)=tr2+ti3 + sr(nr3)=tr2-ti3 + ti0=si(nbase)+si(nr2) + ti2=si(nbase)-si(nr2) + si(nbase)=ti0+ti1 + si(nr2)=ti0-ti1 + si(nr1)=ti2-tr3 + si(nr3)=ti2+tr3 +420 nbase=nbase+4 +430 nbase=nbase+nlup2 +440 nbase=nbase+nlup23 +800 if(na.ne.8) go to 1600 +c +c ********************************************************************** +c +c the following code implements the 8 point post-weave module +c +c ********************************************************************** +c + nlup2=8*(nd2-nb) + nlup23=8*nd2*(nd3-nc) + nbase=1 + do 840 n4=1,nd + do 830 n3=1,nc + do 820 n2=1,nb + nr1=nbase+1 + nr2=nr1+1 + nr3=nr2+1 + nr4=nr3+1 + nr5=nr4+1 + nr6=nr5+1 + nr7=nr6+1 + t1=sr(nbase)-sr(nr1) + sr(nbase)=sr(nbase)+sr(nr1) + sr6=sr(nr2)+si(nr3) + sr(nr2)=sr(nr2)-si(nr3) + t4=sr(nr4)-si(nr5) + t5=sr(nr4)+si(nr5) + t6=sr(nr7)-si(nr6) + t7=sr(nr7)+si(nr6) + sr(nr4)=t1 + sr(nr1)=t4+t6 + sr3=t4-t6 + sr5=t5-t7 + sr(nr7)=t5+t7 + t1=si(nbase)-si(nr1) + si(nbase)=si(nbase)+si(nr1) + t3=si(nr2)-sr(nr3) + si(nr2)=si(nr2)+sr(nr3) + t4=si(nr4)+sr(nr5) + t5=si(nr4)-sr(nr5) + si(nr6)=t3 + t6=sr(nr6)+si(nr7) + t7=sr(nr6)-si(nr7) + si(nr4)=t1 + si(nr1)=t4+t6 + si(nr3)=t4-t6 + si(nr5)=t5+t7 + si(nr7)=t5-t7 + sr(nr3)=sr3 + sr(nr5)=sr5 + sr(nr6)=sr6 +820 nbase=nbase+8 +830 nbase=nbase+nlup2 +840 nbase=nbase+nlup23 +1600 if(na.ne.16) return +c +c ********************************************************************** +c +c the following code implements the 16 point post-weave module +c +c ********************************************************************** +c + nlup2=18*(nd2-nb) + nlup23=18*nd2*(nd3-nc) + nbase=1 + do 1640 n4=1,nd + do 1630 n3=1,nc + do 1620 n2=1,nb + nr1=nbase+1 + nr2=nr1+1 + nr3=nr2+1 + nr4=nr3+1 + nr5=nr4+1 + nr6=nr5+1 + nr7=nr6+1 + nr8=nr7+1 + nr9=nr8+1 + nr10=nr9+1 + nr11=nr10+1 + nr12=nr11+1 + nr13=nr12+1 + nr14=nr13+1 + nr15=nr14+1 + nr16=nr15+1 + nr17=nr16+1 + t(2)=sr(nbase)-sr(nr1) + sr(nbase)=sr(nr1)+sr(nbase) + t(4)=sr(nr2)+si(nr3) + t(3)=sr(nr2)-si(nr3) + t(6)=sr(nr4)+si(nr5) + t(5)=sr(nr4)-si(nr5) + t(8)=-si(nr6)-sr(nr7) + t(7)=-si(nr6)+sr(nr7) + t(9)=sr(nr8)+sr(nr14) + t(15)=sr(nr8)-sr(nr14) + t(13)=-si(nr10)-si(nr12) + t(11)=si(nr10)-si(nr12) + t(16)=sr(nr15)-sr(nr17) + t(12)=sr(nr11)-sr(nr17) + t(10)=-si(nr9)-si(nr16) + t(14)=-si(nr16)+si(nr13) + sr(nr2)=t(5)+t(7) + sr6=t(5)-t(7) + sr10=t(6)+t(8) + sr(nr14)=t(6)-t(8) + q(7)=t(9)+t(10) + q(8)=t(9)-t(10) + q(1)=t(11)+t(12) + q(2)=t(11)-t(12) + q(4)=t(14)+t(15) + q(5)=t(15)-t(14) + q(3)=t(13)+t(16) + q(6)=t(13)-t(16) + sr(nr1)=q(3)+q(7) + sr(nr7)=q(7)-q(3) + sr9=q(8)+q(6) + sr(nr15)=q(8)-q(6) + sr5=q(1)+q(4) + sr3=q(4)-q(1) + sr13=q(2)+q(5) + sr11=q(5)-q(2) + sr(nr8)=t(2) + sr(nr4)=t(3) + sr12=t(4) + t(2)=si(nbase)-si(nr1) + si(nbase)=si(nr1)+si(nbase) + t(4)=si(nr2)-sr(nr3) + t(3)=si(nr2)+sr(nr3) + t(6)=si(nr4)-sr(nr5) + t(5)=si(nr4)+sr(nr5) + t(8)=sr(nr6)-si(nr7) + t(7)=sr(nr6)+si(nr7) + t(9)=si(nr8)+si(nr14) + t(15)=si(nr8)-si(nr14) + t(13)=sr(nr10)+sr(nr12) + t(11)=sr(nr12)-sr(nr10) + t(16)=si(nr15)-si(nr17) + t(12)=si(nr11)-si(nr17) + t(10)=sr(nr9)+sr(nr16) + si(nr2)=t(5)+t(7) + si(nr6)=t(5)-t(7) + si(nr10)=t(6)+t(8) + si(nr14)=t(6)-t(8) + q(7)=t(9)+t(10) + q(8)=t(9)-t(10) + q(1)=t(11)+t(12) + q(2)=t(11)-t(12) + q(4)=t(14)+t(15) + q(5)=t(15)-t(14) + q(3)=t(13)+t(16) + q(6)=t(13)-t(16) + si(nr1)=q(3)+q(7) + si(nr7)=q(7)-q(3) + si(nr9)=q(8)+q(6) + si(nr15)=q(8)-q(6) + si(nr5)=q(1)+q(4) + si(nr3)=q(4)-q(1) + si(nr13)=q(2)+q(5) + si(nr11)=q(5)-q(2) + si(nr8)=t(2) + si(nr4)=t(3) + si(nr12)=t(4) + sr(nr3)=sr3 + sr(nr5)=sr5 + sr(nr6)=sr6 + sr(nr9)=sr9 + sr(nr10)=sr10 + sr(nr11)=sr11 + sr(nr12)=sr12 + sr(nr13)=sr13 +1620 nbase=nbase+18 +1630 nbase=nbase+nlup2 +1640 nbase=nbase+nlup23 + return + end diff --git a/math/ieee/chap1/wfta.f b/math/ieee/chap1/wfta.f new file mode 100644 index 00000000..9e819941 --- /dev/null +++ b/math/ieee/chap1/wfta.f @@ -0,0 +1,150 @@ +c +c----------------------------------------------------------------------- +c subroutine: wfta +c winograd fourier transform algorithm +c----------------------------------------------------------------------- +c + subroutine wfta(xr,xi,n,invrs,init,ierr) + dimension xr(1),xi(1) +c +c inputs: +c n-- transform length. must be formed as the product of +c relatively prime integers from the set: +c 2,3,4,5,7,8,9,16 +c thus the largest possible value of n is 5040. +c xr(.)-- array that holds the real part of the data +c to be transformed. +c xi(.)-- array that holds the imaginary part of the +c data to be transformed. +c invrs-- parameter that flags whether or not the inverse +c transform is to be calculated. a division by n +c is included in the inverse. +c invrs = 1 yields inverse transform +c invrs .ne. 1 gives forward transform +c init-- parameter that flags whether or not the program +c is to be initialized for this value of n. the +c initialization is performed only once in order to +c to speed up the computation on succeeding calls +c to the wfta routine, when n is held fixed. +c init = 0 results in initialization. +c ierr-- error code that is negative when the wfta +c terminates incorrectly. +c 0 = successful completion +c -1 = this value of n does not factor properly +c -2 = an initialization has not been done for +c this value of n. +c +c +c the following two cards may be changed if the maximum +c desired transform length is less than 5040 +c +c ********************************************************************* + dimension sr(1782),si(1782),coef(1782) + integer indx1(1008),indx2(1008) +c ********************************************************************* +c + common na,nb,nc,nd,nd1,nd2,nd3,nd4 +c +c test for initial run +c + if(init.eq.0) call inishl(n,coef,xr,xi,indx1,indx2,ierr) +c + if(ierr.lt.0) return + m=na*nb*nc*nd + if(m.eq.n) go to 100 + ierr=-2 + return +c +c error(-2)-- program not initialized for this value of n +c +100 nmult=nd1*nd2*nd3*nd4 +c +c the following code maps the data arrays xr and xi to +c the working arrays sr and si via the mapping vector +c indx1(.). the permutation of the data follows the +c sino correspondence of the chinese remainder theorem. +c + j=1 + k=1 + inc1=nd1-na + inc2=nd1*(nd2-nb) + inc3=nd1*nd2*(nd3-nc) + do 140 n4=1,nd + do 130 n3=1,nc + do 120 n2=1,nb + do 110 n1=1,na + ind=indx1(k) + sr(j)=xr(ind) + si(j)=xi(ind) + k=k+1 +110 j=j+1 +120 j=j+inc1 +130 j=j+inc2 +140 j=j+inc3 +c +c do the pre-weave modules +c + call weave1(sr,si) +c +c the following loop performs all the multiplications of the +c winograd fourier transform algorithm. the multiplication +c coefficients are stored on the initialization pass in the +c array coef(.). +c + do 200 j=1,nmult + sr(j)=sr(j)*coef(j) + si(j)=si(j)*coef(j) + 200 continue +c +c do the post-weave modules +c + call weave2(sr,si) +c +c +c the following code maps the working arrays sr and si +c to the data arrays xr and xi via the mapping vector +c indx2(.). the permutation of the data follows the +c chinese remainder theorem. +c + j=1 + k=1 + inc1=nd1-na + inc2=nd1*(nd2-nb) + inc3=nd1*nd2*(nd3-nc) +c +c check for inverse +c + if(invrs.eq.1) go to 400 + do 340 n4=1,nd + do 330 n3=1,nc + do 320 n2=1,nb + do 310 n1=1,na + kndx=indx2(k) + xr(kndx)=sr(j) + xi(kndx)=si(j) + k=k+1 +310 j=j+1 +320 j=j+inc1 +330 j=j+inc2 +340 j=j+inc3 + return +c +c different permutation for the inverse +c +400 fn=float(n) + np2=n+2 + indx2(1)=n+1 + do 440 n4=1,nd + do 430 n3=1,nc + do 420 n2=1,nb + do 410 n1=1,na + kndx=np2-indx2(k) + xr(kndx)=sr(j)/fn + xi(kndx)=si(j)/fn + k=k+1 +410 j=j+1 +420 j=j+inc1 +430 j=j+inc2 +440 j=j+inc3 + return + end |