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