1 files changed, 98 insertions, 0 deletions
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