Changing bit_reverse in fft.c

Author: Gathros

chapters/FFT/code/c/fft.c

@@ -4,14 +4,12 @@
 #include <stdlib.h>
 #include <time.h>
 
-#define PI 3.1415926535897932384626
-
 void dft(double complex *X, const size_t N) {
     double complex tmp[N];
     for (size_t i = 0; i < N; ++i) {
         tmp[i] = 0;
         for (size_t j = 0; j < N; ++j) {
-            tmp[i] += X[j] * cexp(-2.0 * PI * I * j * i / N);
+            tmp[i] += X[j] * cexp(-2.0 * M_PI * I * j * i / N);
         }
     }
 
@@ -35,27 +33,30 @@ void cooley_tukey(double complex *X, const size_t N) {
         cooley_tukey(X + N / 2, N / 2);
 
         for (size_t i = 0; i < N / 2; ++i) {
-            X[i + N / 2] = X[i] - cexp(-2.0 * I * PI * i / N) * X[i + N / 2];
+            X[i + N / 2] = X[i] - cexp(-2.0 * I * M_PI * i / N) * X[i + N / 2];
             X[i] -= (X[i + N / 2]-X[i]);
         }
     }
 }
 
 void bit_reverse(double complex *X, size_t N) {
-    double complex temp;
-    unsigned int b;
+    for (int i = 0; i < N; ++i) {
+        int n = i;
+        int a = i;
+        int count = (int)log2((double)N) - 1;
+
+        n >>= 1;
+        while (n > 0) {
+            a = (a << 1) | (n & 1);
+            count--;
+            n >>= 1;
+        }
+        n = (a << count) & ((1 << (int)log2((double)N)) - 1);
 
-    for (unsigned int i = 0; i < N; ++i) {
-        b = i;
-        b = (((b & 0xaaaaaaaa) >> 1) | ((b & 0x55555555) << 1));
-        b = (((b & 0xcccccccc) >> 2) | ((b & 0x33333333) << 2));
-        b = (((b & 0xf0f0f0f0) >> 4) | ((b & 0x0f0f0f0f) << 4));
-        b = (((b & 0xff00ff00) >> 8) | ((b & 0x00ff00ff) << 8));
-        b = ((b >> 16) | (b << 16)) >> (32 - (unsigned int) log2((double)N));
-        if (b > i) {
-            temp = X[b];
-            X[b] = X[i];
-            X[i] = temp;
+        if (n > i) {
+            double complex tmp = X[i];
+            X[i] = X[n];
+            X[n] = tmp;
         }
     }
 }
@@ -65,7 +66,7 @@ void iterative_cooley_tukey(double complex *X, size_t N) {
 
     for (int i = 1; i <= log2((double)N); ++i) {
         int stride = pow(2, i);
-        double complex w = cexp(-2.0 * I * PI / stride);
+        double complex w = cexp(-2.0 * I * M_PI / stride);
         for (size_t j = 0; j < N; j += stride) {
             double complex v = 1.0;
             for (size_t k = 0; k < stride / 2; ++k) {
@@ -93,8 +94,7 @@ int main() {
     }
 
     cooley_tukey(y, 64);
-    //iterative_cooley_tukey(z, 64);
-    dft(z, 64);
+    iterative_cooley_tukey(z, 64);
 
     approx(y, z, 64);