Added FFT implementation post

1 files changed, 241 insertions, 0 deletions

diff --git a/md/writeup/naive_fft_implementation_in_c.md b/md/writeup/naive_fft_implementation_in_c.md
new file mode 100644
index 0000000..a6ca9dc
--- /dev/null
+++ b/md/writeup/naive_fft_implementation_in_c.md
@@ -0,0 +1,241 @@
+title:Naive FFT implementation in C
+keywords:c,linux,fir,dsp,octave,fft
+
+# Naive FFT implementation in C
+
+## Intro
+
+<script>
+MathJax = {
+  tex: {
+    inlineMath: [['$', '$'], ['\\(', '\\)']]
+  }
+};
+</script>
+
+<script type="text/javascript" id="MathJax-script" async
+  src="/js/tex-chtml.js">
+</script>
+
+Time to implement DFT/FFT by myself. Just followed Ifeachor book on theory part and 
+IOWA FFT source on implementation side.
+
+
+## Implementation
+
+### DFT 
+
+
+$$X(k) = F_{D}[x(nT)]  = $$     
+
+$$ \sum_{n=0}^{N-1} x(nT) e^{-jk\omega nT} , k=0,1,...,N-1 $$
+
+
+```c
+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];
+    }
+    
+}
+```
+### FFT
+#### Shuffling arrays 
+
+```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;
+       }
+}
+```
+
+
+#### Implementation
+
+```c
+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 = 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;
+            }
+        }
+    }
+```
+
+
+
+## Octave verification code
+
+```matlab
+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)
+
+```
+
+### Source
+
+Source located in main branch of git repo
+
+```
+git clone http://git.main.lv/cgit.cgi/NaiveFFT.git
+git checkout main
+```
+
+### Build
+
+#### Linux
+
+```
+cd NaiveFFT/Build
+make
+```
+
+#### Macos
+
+Open with Xcode and build
+
+
+## Links
+
+
+[01] https://rosettacode.org/wiki/Fast_Fourier_transform#C  
+[02] https://www.math.wustl.edu/~victor/mfmm/fourier/fft.c  
+[03] https://www.strauss-engineering.ch/libdsp.html  
+[04] http://www.iowahills.com/FFTCode.html  
+[05] https://github.com/rshuston/FFT-C/blob/master/libfft/fft.c  
+[06] http://www.guitarscience.net/papers/fftalg.pdf  
+[07] https://root.cern/doc/master/FFT_8C_source.html  
+[08] https://lloydrochester.com/post/c/example-fft/  
+[09] https://github.com/mborgerding/kissfft/blob/master/kiss_fft.c  
+[10] https://community.vcvrack.com/t/complete-list-of-c-c-fft-libraries/9153  
+[11] https://octave.org/doc/v4.2.1/Signal-Processing.html  
+[12] https://en.wikipedia.org/wiki/Fast_Fourier_transform  
+[13] http://www.iowahills.com/FFTCode.html  
+[14] Digital Signal Processing: A Practical Approach by (Emmanuel Ifeachor, Barrie Jervis)  
+
+