Changeset 16099 for trunk/Ohana/src/opihi/lib.data/fftold.c

Timestamp:

Jan 16, 2008, 1:10:16 PM (19 years ago)

Author:

eugene

Message:

working on the fft function

File:

• trunk/Ohana/src/opihi/lib.data/fftold.c (modified) (6 diffs)

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

@@ -3 +3 @@
 #define FSWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
 
-void fft (float *Data, int N, int isign) {
+void fftold (float *Data, int N, int isign) {
 
   int n,mmax,m,j,istep,i;
@@ -51 +51 @@
 }
 
-void fftold (float *Data, int N, int isign) {
-
-  int n,mmax,m,j,istep,i;
+void fftold_Nindices (float *Data, int N, int isign) {
+
+  int n,ifp1,m,j,ipf2,i;
   double wtemp,wr,wpr,wpi,wi,theta;
   float tempr,tempi, *data;
 
-  data = Data - 1;
+  data = Data;
   n = N << 1;
-  j = 1;
-
-  for (i = 1; i < n; i+=2) {
-    if (j > i) {
-      FSWAP (data[j], data[i]);
-      FSWAP (data[j+1], data[i+1]);
-    }
-    m = n >> 1;
-    while (m >= 2 && j > m) {
-      j -= m;
-      m >>= 1;
-    }
-    j += m;
-  }
-  mmax = 2;
-  while (n > mmax) {
-    istep = 2*mmax;
-    theta = 6.28318530717959 / (isign*mmax);
+
+  ip1 = 2;
+  ip2 = n;
+  i2rev = 0;
+  for (i2 = 0; i2 < ip2; i2 += ip1) {
+    if (i2 < i2ref) {
+      FSWAP (data[i2rev], data[i2]);
+      FSWAP (data[i2rev+1], data[i2+1]);
+    }
+    ibit = ip2 >> 1;
+    while (ibit >= ip1 && i2rev >= ibit) {
+      i2rev -= ibit;
+      ibit >>= 1;
+    }
+    i2rev += ibit;
+  }
+
+  ifp1 = 2;
+  while (ifp1 < ip2) {
+    ifp2 = ifp1 << 1;
+    theta = isign*6.28318530717959 / (2*ifp1/ip1);
     wtemp = sin(0.5*theta);
     wpr = -2.0*wtemp*wtemp;
@@ -82 +85 @@
     wr = 1.0;
     wi = 0.0;
-    for (m = 1; m < mmax; m+=2) {
-      for (i = m; i <= n; i+=istep) {
-        j = i + mmax;
-        tempr = wr*data[j] - wi*data[j+1];
-        tempi = wr*data[j+1] + wi*data[j];
-        data[j] = data[i] - tempr;
-        data[j+1] = data[i+1] - tempi;
-        data[i] += tempr;
-        data[i+1] += tempi;
+    for (i3 = 0; i3 < ifp1; i3+=ip1) {
+      // in the N-d version, this loop is broken into two loops
+      for (i1 = i3; i1 < n; i1+=ipf2) {
+        k1 = i1;
+        k2 = k1 + ifp1;
+        tempr = wr*data[k2] - wi*data[k2+1];
+        tempi = wr*data[k2+1] + wi*data[k2];
+        data[k2] = data[k1] - tempr;
+        data[k2+1] = data[k1+1] - tempi;
+        data[k1] += tempr;
+        data[k1+1] += tempi;
       }
       wr = (wtemp = wr)*wpr - wi*wpi+wr;
       wi = wi*wpr + wtemp*wpi + wi;
     }
-    mmax = istep;
-  }
+    ifp1 = ifp2;
+  }
+
 }
 
 /* convert indices to zero reference */
-void fftN (float *data, int *nn, int ndim, int isign) {
+void fftNold (float *data, int *nn, int ndim, int isign) {
 
   int i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
@@ -138 +144 @@
     while (ifp1 < ip2) {
       ifp2 = ifp1 << 1;
-      theta = isign*6.28318530717959/(ifp2/ip1);
+      theta = isign*6.28318530717959/(2*ifp1/ip1);
       wtemp = sin(0.5*theta);
       wpr = -2.0*wtemp*wtemp;
@@ -182 +188 @@
 */ 
 
-int IsBinary (int N) {
+int IsBinaryOld (int N) {
 
   int i, nbit;
@@ -199 +205 @@
 
 }  
-
-#define FSWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
-
-void fourn (float *data, int *nn, int ndim, int isign) {
-
-  int i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
-  int ibit,idim,k1,k2,n,nprev,nrem,ntot;
-  float tempi,tempr;
-  double theta,wi,wpi,wpr,wr,wtemp;
-
-  ntot = 1;
-  for (idim = 1; idim <= ndim; idim++)
-    ntot *= nn[idim];
-  nprev = 1;
-  for (idim = ndim; idim >= 1; idim--) {
-    n = nn[idim];
-    nrem = ntot/(n*nprev);
-    ip1 = nprev << 1;
-    ip2 = ip1*n;
-    ip3 = ip2*nrem;
-    i2rev = 1;
-    for (i2 = 1; i2 <= ip2; i2+=ip1) {
-      if (i2 < i2rev) {
-        for (i1 = i2; i1 <= i2+ip1-2; i1+=2) {
-          for (i3 = i1; i3 <= ip3;i3+=ip2) {
-            i3rev = i2rev+i3-i2;
-            FSWAP(data[i3],data[i3rev]);
-            FSWAP(data[i3+1],data[i3rev+1]);
-          }
-        }
-      }
-      ibit = ip2 >> 1;
-      while (ibit >= ip1 && i2rev > ibit) {
-        i2rev -=  ibit;
-        ibit >>= 1;
-      }
-      i2rev += ibit;
-    }
-    ifp1 = ip1;
-    while (ifp1 < ip2) {
-      ifp2 = ifp1 << 1;
-      theta = isign*6.28318530717959/(ifp2/ip1);
-      wtemp = sin(0.5*theta);
-      wpr  =  -2.0*wtemp*wtemp;
-      wpi = sin(theta);
-      wr = 1.0;
-      wi = 0.0;
-      for (i3 = 1;i3<=ifp1;i3+=ip1) {
-        for (i1 = i3;i1<=i3+ip1-2;i1+=2) {
-          for (i2 = i1;i2<=ip3;i2+=ifp2) {
-            k1 = i2;
-            k2 = k1+ifp1;
-            tempr = wr*data[k2]-wi*data[k2+1];
-            tempi = wr*data[k2+1]+wi*data[k2];
-            data[k2] = data[k1]-tempr;
-            data[k2+1] = data[k1+1]-tempi;
-            data[k1] += tempr;
-            data[k1+1] += tempi;
-          }
-        }
-        wr = (wtemp = wr)*wpr-wi*wpi+wr;
-        wi = wi*wpr+wtemp*wpi+wi;
-      }
-      ifp1 = ifp2;
-    }
-    nprev *= n;
-  }
-}
-
-#undef FSWAP