Changeset 20936 for trunk/Ohana/src/opihi/lib.data/fft.c

Timestamp:

Dec 7, 2008, 3:31:01 PM (18 years ago)

Author:

eugene

Message:

big update from eam_branch_20081124 with updates to Opihi math

File:

• trunk/Ohana/src/opihi/lib.data/fft.c (modified) (1 diff)

trunk/Ohana/src/opihi/lib.data/fft.c

@@ -133 +133 @@
 }
 
+// fft based on code by Douglas L. Jones (see note at EOF). modified for Ohana C style
+void dfft1D (double *x, double *y, int n, int Nbit, int forward) {
+
+  int i,j,k,n1,n2;
+  double c,s,e,a,t1,t2;        
+  double factor;
+         
+  // bit-reverse
+  j = 0; 
+  n2 = n/2;
+  for (i = 1; i < n - 1; i++) {
+    n1 = n2;
+    while ( j >= n1 ) {
+      j -= n1;
+      n1 /= 2;
+    }
+    j += n1;
+               
+    if (i < j) {
+      t1 = x[i];
+      x[i] = x[j];
+      x[j] = t1;
+      t1 = y[i];
+      y[i] = y[j];
+      y[j] = t1;
+    }
+  }
+                                          
+  n1 = 0; /* FFT */
+  n2 = 1;
+                                             
+  if (forward) {
+    factor = +2.0*M_PI;
+  } else {
+    factor = -2.0*M_PI;
+  }
+
+  for (i=0; i < Nbit; i++) {
+    n1 = n2;
+    n2 = n2 + n2;
+    e = factor/n2;
+    a = 0.0;
+                                             
+    for (j=0; j < n1; j++) {
+      c = cos(a);
+      s = sin(a);
+      a = a + e;
+                                            
+      for (k=j; k < n; k=k+n2) {
+        t1 = c*x[k+n1] - s*y[k+n1];
+        t2 = s*x[k+n1] + c*y[k+n1];
+        x[k+n1] = x[k] - t1;
+        y[k+n1] = y[k] - t2;
+        x[k] = x[k] + t1;
+        y[k] = y[k] + t2;
+      }
+    }
+  }
+                                      
+  // re-normalize
+  for (i = 0; i < n; i++) {
+    x[i] /= n;
+    y[i] /= n;
+  }
+
+  return;
+}                          
+
+// This implementation uses the 1-D fft above for each of the vectors in each dimension.
+// This requires 2(Nx*Ny*...) mem copies, but the fft operations are likely to happen in
+// cache.
+int dfftND (double *x, double *y, int Ndim, int *Nsize, int forward) {
+
+  int i, nIndex, minor, major, iDim;
+  int step, Nmajor, Nminor, Nmax, Ntotal;
+  int *Nbit;
+  double *tmpX, *tmpY;
+
+  ALLOCATE (Nbit, int, Ndim);
+
+  // find the longest axis and allocate storage for that length
+  Nmax = 0;
+  Ntotal = 1;
+  for (i = 0; i < Ndim; i++) {
+    Nmax = MAX(Nmax, Nsize[i]);
+    Ntotal *= Nsize[i];
+    if (!IsBinary (Nsize[i], &Nbit[i])) {
+      free (Nbit);
+      return (FALSE);
+    }
+  }
+  ALLOCATE (tmpX, double, Nmax);
+  ALLOCATE (tmpY, double, Nmax);
+  
+  step = 1;
+  Nminor = 1;
+  Nmajor = Ntotal;
+  for (iDim = 0; iDim < Ndim; iDim++) {
+    step *= Nsize[iDim];
+    Nmajor /= Nsize[iDim];
+
+    // we perform the FFT along all other dimensions 
+    for (major = 0; major < Nmajor; major++) {
+      for (minor = 0; minor < Nminor; minor++) {
+        // nIndex = minor + i*Nminor + major*step;
+        // extract the data values to the temp vector
+        nIndex = minor + major*step;
+        for (i = 0; i < Nsize[iDim]; i++) {
+          tmpX[i] = x[nIndex];
+          tmpY[i] = y[nIndex];
+          nIndex += Nminor;
+        }
+
+        dfft1D (tmpX, tmpY, Nsize[iDim], Nbit[iDim], forward);
+
+        // replace the result vectors
+        nIndex = minor + major*step;
+        for (i = 0; i < Nsize[iDim]; i++) {
+          x[nIndex] = tmpX[i];
+          y[nIndex] = tmpY[i];
+          nIndex += Nminor;
+        }
+      }
+    }
+    Nminor *= Nsize[iDim];
+  }
+  free (Nbit);
+  free (tmpX);
+  free (tmpY);
+  return (TRUE);
+}
+
 // check that a number is binary (2^Nbit).  returns int(log_2(N)) in Nbit
 int IsBinary (int N, int *Nbit) {