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

Timestamp:

Jan 17, 2008, 8:41:08 AM (19 years ago)

Author:

eugene

Message:

converting fft and fftn to code from Jones, with Magnier version of multi-dim

File:

• trunk/Ohana/src/opihi/lib.data/fft.c (modified) (3 diffs)

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

@@ -2 +2 @@
 
 // fft based on code by Douglas L. Jones (see note at EOF). modified for Ohana C style
-void fftN (float *data, int *nn, int ndim, int isign) {
-  
-  fprintf (stderr, "broken\n");
-  return;
-}
-
-void fft (float *x, float *y, int n, int Nbit, int forward) {
+void fft1D (float *x, float *y, int n, int Nbit, int forward) {
 
   int i,j,k,n1,n2;
@@ -75 +69 @@
 }                          
 
+// 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 fftND (float *x, float *y, int Ndim, int *Nsize, int forward) {
+
+  int i, nIndex, minor, major, iDim;
+  int step, Nmajor, Nminor, Nmax, Ntotal;
+  int *Nbit;
+  float *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, float, Nmax);
+  ALLOCATE (tmpY, float, 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;
+        }
+
+        fft1D (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 (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) {
@@ -80 +137 @@
   int i, Nset;
 
-  *Nbit = 0;
+  if (Nbit != NULL) *Nbit = 0;
   Nset = 0;
   for (i = 0; i < 8*sizeof(int); i++) {
     if (N & 0x01) {
       Nset ++;
-      *Nbit = i;
+      if (Nbit != NULL) *Nbit = i;
     }
     N >>= 1;
   }
   if (Nset > 1) return (FALSE);
+  return (TRUE);  
 }