workinf dft and fft. Octave test file

2 files changed, 202 insertions, 19 deletions

diff --git a/NaiveFFT/main.c b/NaiveFFT/main.c
index fec893b..f8b56cd 100644
--- a/NaiveFFT/main.c
+++ b/NaiveFFT/main.c
@@ -11,28 +11,31 @@
 #include <string.h>
 #include <math.h>
 #include <complex.h>
+#include <time.h>
 
 //data array
 #define DATA_SIZE (8)
-double data_i[DATA_SIZE] = {1.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f};
+//double data_i[DATA_SIZE] = {0.0f,1.0f,2.0f,3.0f,4.0f,5.0f,6.0f,7.0f};
+//double data_q[DATA_SIZE] = {0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f};
+double data_i[DATA_SIZE] = {1.0f,1.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f};
 double data_q[DATA_SIZE] = {0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f};
+
 //calc the fft
 void fft_if(double *x_i, double *x_q, int n, int inv) {
     
     int i=0,j=0,k=0,m=0, irem=0, sign;
     double tq,ti;
     
-    printf("FFT\n");
     k=n/2;
-    for (j=1,i=1;i<n; i++) {
+    for (j=1,i=0;i<n,j<n; i++) {
         if (i<j) {
             
             //swap values
-            ti = x_i[j];
-            tq = x_q[j];
+            ti = x_i[j-1];
+            tq = x_q[j-1];
             
-            x_i[j] = x_i[i];
-            x_q[j] = x_q[i];
+            x_i[j-1] = x_i[i];
+            x_q[j-1] = x_q[i];
             
             x_i[i] = ti;
             x_q[i] = tq;
@@ -53,7 +56,9 @@ void fft_if(double *x_i, double *x_q, int n, int inv) {
         }
         j = j+k;
     }
-    
+    //print rearranged buffer
+    /*
+    printf("FFT rearranged\n");
     for (i=0;i<DATA_SIZE;i++) {
         printf("%.2f ",x_i[i]);
     }
@@ -62,6 +67,7 @@ void fft_if(double *x_i, double *x_q, int n, int inv) {
         printf("%.2f ",x_q[i]);
     }
     printf("\n");
+    */
     
     //calculate number of stages:
     //m=log2(n)
@@ -81,7 +87,7 @@ void fft_if(double *x_i, double *x_q, int n, int inv) {
     }
     
     // transform
-    for (i=0; i<m; i++) {
+    for (i=1; i<m; i++) {
         int le=pow(2,i);
         int le1=le/2;
         double ui=1.0;
@@ -89,21 +95,24 @@ void fft_if(double *x_i, double *x_q, int n, int inv) {
         double wi = cos(M_PI/le1);
         double wq = sign*sin(M_PI/le1);
         
-        for (j=1; j<le1;j++) {
-            i=j;
-            while (i<n) {
-                int ip = i+le1;
+        printf("le1=%d\n",le1);
+        for (j=0; j<le1;j++) {
+            
+            k=j;
+            while (k<n) {
+                int ip = k+le1;
                 
                 ti = x_i[ip]*ui - x_q[ip]*uq;
                 tq = x_q[ip]*ui + x_i[ip]*uq;
                 
-                x_i[ip] = x_i[i] - ti;
-                x_q[ip] = x_q[i] - tq;
+                x_i[ip] = x_i[k-1] - ti;
+                x_q[ip] = x_q[k-1] - tq;
                 
-                x_i[i] = x_i[i] + ti;
-                x_q[i] = x_q[i] + tq;
+                x_i[k-1] = x_i[k-1] + ti;
+                x_q[k-1] = x_q[k-1] + tq;
                 
-                i = i+le;
+                k = k+le;
+                printf("k=%d\n",k-1);
             }
             double temp = ui*wi - uq*wq;
             uq = uq*wi + ui*wq;
@@ -121,6 +130,124 @@ void fft_if(double *x_i, double *x_q, int n, int inv) {
     
 }
 
+#define complex_mul_re(a_re, a_im, b_re, b_im)  (a_re * b_re - a_im * b_im)
+#define complex_mul_im(a_re, a_im, b_re, b_im)  (a_re * b_im + a_im * b_re)
+//https://github.com/rshuston/FFT-C/blob/master/libfft/fft.c
+void ffti_shuffle_1(double *x_i, double *x_q, uint64_t n) {
+    int Nd2 = n>>1;
+    int Nm1 = n-1;
+    int i,j;
+    
+    
+    for (i = 0, j = 0; i < n; i++) {
+           if (j > i) {
+               double tmp_r = x_i[i];
+               double tmp_i = x_q[i];
+               //data[i] = data[j];
+               x_i[i] = x_q[j];
+               x_q[i] = x_q[j];
+               //data[j] = tmp;
+               x_i[j] = tmp_r;
+               x_q[j] = tmp_i;
+           }
+
+           /*
+            * Find least significant zero bit
+            */
+
+           unsigned lszb = ~i & (i + 1);
+
+           /*
+            * Use division to bit-reverse the single bit so that we now have
+            * the most significant zero bit
+            *
+            * N = 2^r = 2^(m+1)
+            * Nd2 = N/2 = 2^m
+            * if lszb = 2^k, where k is within the range of 0...m, then
+            *     mszb = Nd2 / lszb
+            *          = 2^m / 2^k
+            *          = 2^(m-k)
+            *          = bit-reversed value of lszb
+            */
+
+           unsigned mszb = Nd2 / lszb;
+
+           /*
+            * Toggle bits with bit-reverse mask
+            */
+
+           unsigned bits = Nm1 & ~(mszb - 1);
+           j ^= bits;
+       }
+}
+void fft_1(double *x_i, double *x_q, uint64_t n, uint64_t inv) {
+    uint64_t n_log2;
+    uint64_t r;
+    uint64_t m, md2;
+    uint64_t i,j,k;
+    uint64_t i_e, i_o;
+    double theta_2pi;
+    double theta;
+    
+    double Wm_r, Wm_i, Wmk_r, Wmk_i;
+    double u_r, u_i, t_r, t_i;
+    
+    //find log of n
+    i=n;
+    n_log2 = 0;
+    while (i>1) {
+        i=i/2;
+        n_log2+=1;
+    }
+    
+    if (inv==1) {
+        theta_2pi = -2*M_PI;
+    } else {
+        theta_2pi = 2*M_PI;
+    }
+    
+    
+    for (i=1; i<= n_log2; i++) {
+        m  = 1 << i;
+        md2 = m >> 1;
+        theta = theta_2pi / m;
+        Wm_r = cos(theta);
+        Wm_i = sin(theta);
+        
+        
+        for (j=0; j<n; j+=m) {
+            Wmk_r = 1.0f;
+            Wmk_i = 0.0f;
+            
+            for (k=0; k<md2; k++) {
+                i_e = j+k;
+                i_o = i_e + md2;
+                
+                u_r = x_i[i_e];
+                u_i = x_q[i_e];
+                
+                //t_r = Wmk_r * x_i[i_o] - Wmk_i * x_q[i_o];
+                //t_i = Wmk_r * x_q[i_o] - Wmk_i * x_i[i_o];
+                t_r = complex_mul_re(Wmk_r, Wmk_i, x_i[i_o], x_q[i_o]);
+                t_i = complex_mul_im(Wmk_r, Wmk_i, x_i[i_o], x_q[i_o]);
+                
+                
+                x_i[i_e] = u_r + t_r;
+                x_q[i_e] = u_i + t_i;
+                
+                x_i[i_o] = u_r - t_r;
+                x_q[i_o] = u_i - t_i;
+
+                t_r = complex_mul_re(Wmk_r, Wmk_i, Wm_r, Wm_i);
+                t_i = complex_mul_im(Wmk_r, Wmk_i, Wm_r, Wm_i);
+                
+                Wmk_r = t_r;
+                Wmk_i = t_i;
+            }
+        }
+    }
+}
+
 //dft works fine
 void dft(double *x_i, double *x_q, int n, int inv) {
     double Wn,Wk;
@@ -181,10 +308,13 @@ int main(int argc, const char * argv[]) {
         printf("%.2f ",data_q[i]);
     }
     printf("\n");
-    fft_if(data_i, data_q, DATA_SIZE, 0);
+    ffti_shuffle_1(data_i, data_q, DATA_SIZE);
+    fft_1(data_i, data_q, DATA_SIZE, 0);
     
     //dft(data_i, data_q, DATA_SIZE,0);
     
+    printf("after FFT\n");
+    
     for (i=0;i<DATA_SIZE;i++) {
         printf("%.2f ",data_i[i]);
     }
@@ -194,5 +324,46 @@ int main(int argc, const char * argv[]) {
     }
     printf("\n");
     
+    //big buffer test
+    {
+        time_t timer;
+        char buffer[32];
+        struct tm *tm_info;
+        
+        timer = time(NULL);
+        tm_info = localtime(&timer);
+        
+        strftime(buffer, 32, "%H:%M:%S", tm_info);
+        puts(buffer);
+        
+    }
+    
+    //test big big fft
+    {
+        double *d_i=malloc(sizeof(double)*4*8192+16);
+        double *d_q=malloc(sizeof(double)*4*8192+16);
+        for (i=0;i<4*8192;i++) {
+            d_i[i] = sin(i*1.0/128);
+            d_q[i] = 0.0f;
+        }
+        
+        dft(d_i, d_q, 4*8192,0);
+        free(d_i); d_i = NULL;
+        free(d_q); d_q = NULL;
+    }
+    
+    
+    {
+        time_t timer;
+        char buffer[32];
+        struct tm *tm_info;
+        
+        timer = time(NULL);
+        tm_info = localtime(&timer);
+        
+        strftime(buffer, 32, "%H:%M:%S", tm_info);
+        puts(buffer);
+        
+    }
     return 0;
 }
diff --git a/Test/test_fft.m b/Test/test_fft.m
new file mode 100644
index 0000000..4a11a09
--- /dev/null
+++ b/Test/test_fft.m
@@ -0,0 +1,12 @@
+
+data1 = [1.0 0.0 0.0 0.0  0.0 0.0 0.0 0.0]
+res1 = fft(data1)
+
+data2 = [1.0 1.0 0.0 0.0  0.0 0.0 0.0 0.0]
+res2 = fft(data2)
+
+data3 = [1.0 1.0 i i  0.0 0.0 0.0 0.0]
+res3 = fft(data3)
+
+data4 = [1.0+i 1.0+i 1.0+i 1.0+i  0.0 0.0 0.0 0.0]
+res4 = fft(data4)