Reaarange source

3 files changed, 294 insertions, 273 deletions

diff --git a/NaiveFFT/fft.c b/NaiveFFT/fft.c
new file mode 100644
index 0000000..cfcf3fa
--- /dev/null
+++ b/NaiveFFT/fft.c
@@ -0,0 +1,271 @@
+//
+//  fft.c
+//  NaiveFFT
+//
+//  Created by Jacky Jack on 26/08/2021.
+//
+
+#include "fft.h"
+
+
+
+
+//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;
+    
+    k=n/2;
+    for (j=1,i=0;i<n,j<n; i++) {
+        if (i<j) {
+            
+            //swap values
+            ti = x_i[j-1];
+            tq = x_q[j-1];
+            
+            x_i[j-1] = x_i[i];
+            x_q[j-1] = x_q[i];
+            
+            x_i[i] = ti;
+            x_q[i] = tq;
+            printf("i=%d j=%d\n",i,j);
+            //find middle point
+            k = n/2;
+            while (k<j) {
+                j = j-k;
+                k=k/2;
+            }
+        } else {
+            printf("i=%d j=%d\n",i,j);
+            k = n/2;
+            while (k<j) {
+                j = j - k;
+                k = k/2;
+            }
+        }
+        j = j+k;
+    }
+    //calculate number of stages:
+    //m=log2(n)
+    m=0;
+    irem = n;
+    while(irem>1) {
+        irem = irem/2;
+        m = m+1;
+    }
+    
+    //FFT or IFFT
+    if (inv==1) {
+        sign = 1;
+    } else {
+        sign = -1;
+    }
+    
+    // transform
+    for (i=1; i<m; i++) {
+        int le=pow(2,i);
+        int le1=le/2;
+        double ui=1.0;
+        double uq=0.0;
+        double wi = cos(M_PI/le1);
+        double wq = sign*sin(M_PI/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[k-1] - ti;
+                x_q[ip] = x_q[k-1] - tq;
+                
+                x_i[k-1] = x_i[k-1] + ti;
+                x_q[k-1] = x_q[k-1] + tq;
+                
+                k = k+le;
+                printf("k=%d\n",k-1);
+            }
+            double temp = ui*wi - uq*wq;
+            uq = uq*wi + ui*wq;
+            ui = temp;
+        }
+    }
+    
+    //if inverse
+    if (inv==1) {
+        for (i=0;i<n;i++) {
+            x_i[i] = x_i[i]/n;
+            x_q[i] = x_q[i]/n;
+        }
+    }
+    
+}
+
+#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;
+    //static array
+    //double Xi[DATA_SIZE],Xq[DATA_SIZE];
+    //dynamic array
+    double *Xi, *Xq;
+    double c,s;
+    int i,j;
+    
+    Xi = malloc(n*sizeof(double));
+    Xq = malloc(n*sizeof(double));
+    
+    Wn = 2*M_PI/n;
+    
+    if (inv==1) {
+        Wn=-Wn;
+    }
+    
+    for (i=0;i<n;i++) {
+        Xi[i] = 0.0f;
+        Xq[i] = 0.0f;
+        
+        Wk = i*Wn;
+        for (j=0;j<n;j++) {
+            c = cos(j*Wk);
+            s = sin(j*Wk);
+            
+            //i - real, q - imaginary
+            Xi[i] = Xi[i] + x_i[j]*c + x_q[j]*s;
+            Xq[i] = Xq[i] - x_i[j]*s + x_q[j]*c;
+        }
+        
+        if (inv==1) {
+            Xi[i] = Xi[i]/n;
+            Xq[i] = Xq[i]/n;
+        }
+    }
+    
+    for (i=0;i<n;i++) {
+        x_i[i] = Xi[i];
+        x_q[i] = Xq[i];
+    }
+    
+}
diff --git a/NaiveFFT/fft.h b/NaiveFFT/fft.h
new file mode 100644
index 0000000..f1ee783
--- /dev/null
+++ b/NaiveFFT/fft.h
@@ -0,0 +1,20 @@
+//
+//  fft.h
+//  NaiveFFT
+//
+//  Created by Jacky Jack on 26/08/2021.
+//
+
+#ifndef fft_h
+#define fft_h
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+void fft_if(double *x_i, double *x_q, int n, int inv);
+void ffti_shuffle_1(double *x_i, double *x_q, uint64_t n);
+void fft_1(double *x_i, double *x_q, uint64_t n, uint64_t inv);
+void dft(double *x_i, double *x_q, int n, int inv);
+
+#endif /* fft_h */
diff --git a/NaiveFFT/main.c b/NaiveFFT/main.c
index 1d71ac5..718e5d0 100644
--- a/NaiveFFT/main.c
+++ b/NaiveFFT/main.c
@@ -14,286 +14,16 @@
 #include <time.h>
 #include <stdint.h>
 
-//data array
+#include "fft.h"
+
+//test data array
 #define DATA_SIZE (8)
 //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;
-    
-    k=n/2;
-    for (j=1,i=0;i<n,j<n; i++) {
-        if (i<j) {
-            
-            //swap values
-            ti = x_i[j-1];
-            tq = x_q[j-1];
-            
-            x_i[j-1] = x_i[i];
-            x_q[j-1] = x_q[i];
-            
-            x_i[i] = ti;
-            x_q[i] = tq;
-            printf("i=%d j=%d\n",i,j);
-            //find middle point
-            k = n/2;
-            while (k<j) {
-                j = j-k;
-                k=k/2;
-            }
-        } else {
-            printf("i=%d j=%d\n",i,j);
-            k = n/2;
-            while (k<j) {
-                j = j - k;
-                k = k/2;
-            }
-        }
-        j = j+k;
-    }
-    //print rearranged buffer
-    /*
-    printf("FFT rearranged\n");
-    for (i=0;i<DATA_SIZE;i++) {
-        printf("%.2f ",x_i[i]);
-    }
-    printf("\n");
-    for (i=0;i<DATA_SIZE;i++) {
-        printf("%.2f ",x_q[i]);
-    }
-    printf("\n");
-    */
-    
-    //calculate number of stages:
-    //m=log2(n)
-    m=0;
-    irem = n;
-    while(irem>1) {
-        irem = irem/2;
-        m = m+1;
-    }
-    
-    
-    //FFT or IFFT
-    if (inv==1) {
-        sign = 1;
-    } else {
-        sign = -1;
-    }
-    
-    // transform
-    for (i=1; i<m; i++) {
-        int le=pow(2,i);
-        int le1=le/2;
-        double ui=1.0;
-        double uq=0.0;
-        double wi = cos(M_PI/le1);
-        double wq = sign*sin(M_PI/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[k-1] - ti;
-                x_q[ip] = x_q[k-1] - tq;
-                
-                x_i[k-1] = x_i[k-1] + ti;
-                x_q[k-1] = x_q[k-1] + tq;
-                
-                k = k+le;
-                printf("k=%d\n",k-1);
-            }
-            double temp = ui*wi - uq*wq;
-            uq = uq*wi + ui*wq;
-            ui = temp;
-        }
-    }
-    
-    //if inverse
-    if (inv==1) {
-        for (i=0;i<n;i++) {
-            x_i[i] = x_i[i]/n;
-            x_q[i] = x_q[i]/n;
-        }
-    }
-    
-}
-
-#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;
-    //static array
-    //double Xi[DATA_SIZE],Xq[DATA_SIZE];
-    //dynamic array
-    double *Xi, *Xq;
-    double c,s;
-    int i,j;
-    
-    Xi = malloc(n*sizeof(double));
-    Xq = malloc(n*sizeof(double));
-    
-    Wn = 2*M_PI/n;
-    
-    if (inv==1) {
-        Wn=-Wn;
-    }
-    
-    for (i=0;i<n;i++) {
-        Xi[i] = 0.0f;
-        Xq[i] = 0.0f;
-        
-        Wk = i*Wn;
-        for (j=0;j<n;j++) {
-            c = cos(j*Wk);
-            s = sin(j*Wk);
-            
-            //i - real, q - imaginary
-            Xi[i] = Xi[i] + x_i[j]*c + x_q[j]*s;
-            Xq[i] = Xq[i] - x_i[j]*s + x_q[j]*c;
-        }
-        
-        if (inv==1) {
-            Xi[i] = Xi[i]/n;
-            Xq[i] = Xq[i]/n;
-        }
-    }
-    
-    for (i=0;i<n;i++) {
-        x_i[i] = Xi[i];
-        x_q[i] = Xq[i];
-    }
-    
-}
 
 int main(int argc, const char * argv[]) {
     int i;