Initial commit

1 files changed, 84 insertions, 0 deletions

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