path: root/md/writeup/naive_fft_implementation_in_c.md

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>
<script src="/demo/naive_fft/index.js"></script>
<script type='text/javascript'>
          var Module = {};
          fetch('/demo/naive_fft/index.wasm')
            .then(response =>
              response.arrayBuffer()
            ).then(buffer => {
              //Module.canvas = document.getElementById("canvas");
              Module.wasmBinary = buffer;
              var script = document.createElement('script');
              script.src = "/demo/naive_fft/index.js";
              script.onload = function() {
                console.log("Emscripten boilerplate loaded.");
                run_dft         = Module.cwrap("dft",            [],[['float'],['float'],'number','number']);
                run_fft_shuffle = Module.cwrap("ffti_shuffle_1", [],[['float'],['float'],'number']);
                run_fft         = Module.cwrap("fft_1",          [],[['float'],['float'],'number','number']);



                var inputDataI    = document.querySelector('.inputDataI');
                var inputDataIerr = document.querySelector('.inputDataIerr');
                var inputDataQ    = document.querySelector('.inputDataQ');
                var inputDataQerr = document.querySelector('.inputDataQerr');
                var outputDataI   = document.querySelector('.outputDataI');
                var outputDataQ   = document.querySelector('.outputDataQ');

                document.querySelector("#calcButton").onclick = function() {
                      console.log("calc FFT");

                      function cArray(size) {
                          var offset = Module._malloc(size * 8);
                          Module.HEAPF64.set(new Float64Array(size), offset / 8);
                          return {
                              "data": Module.HEAPF64.subarray(offset / 8, offset / 8 + size),
                              "offset": offset
                          }
                      }

                      function checkArr(inputData,errorOutput) {
                        var strArr = inputData.value;
                        
                        strArr = strArr.replace(/ +(?= )/g,'');
                        var splArr = strArr.split(" ");
                        var data = Array(splArr.length).fill(0);
                        for (let i=0; i<splArr.length;i++) {
                          try {
                            var j=splArr[i];
                            if ( j != parseFloat(j)) throw new Error("ERROR: not float "+splArr[i]+" at index "+i);
                            data[i] = parseFloat(j)
                          } catch (err) {
                            console.log("Error: ",err.message)
                            errorOutput.innerHTML = err.message;
                            return {"result":false,"data":[]};
                          }

                        }
                        return {"result":true,"data":data};
                      }
                      console.log("checkArr");
                      var arrIok = checkArr(inputDataI,inputDataIerr);
                      console.log("checkArr");
                      var arrQok = checkArr(inputDataQ,inputDataQerr);
                      //max size 8 values, as thats enought for a demo
                      //https://stackoverflow.com/questions/17883799/how-to-handle-passing-returning-array-pointers-to-emscripten-compiled-code
                      if (!arrIok.result) return;
                      if (!arrQok.result) return;

                      var N=8;
                      var myArray_i = cArray(N);
                      var myArray_q = cArray(N);

                      //set test data
                      for(let i=0,j=0;i<N,j<arrIok.data.length;i++,j++) {
                        myArray_i.data[i] = arrIok.data[j];
                      }
                      for(let i=0,j=0;i<N,j<arrQok.data.length;i++,j++) {
                        myArray_q.data[i] = arrQok.data[j];
                      }
                      console.log(myArray_i.data);
                      console.log(myArray_q.data);
                      //DFT
                      run_dft(myArray_i.offset,myArray_q.offset,N,0);
                      //FFT
                      //run_fft_shuffle(myArray_i.offset,myArray_q.offset,N);
                      //run_fft(myArray_i.offset,myArray_q.offset,N,1);
                      
                      console.log(myArray_i.data);
                      console.log(myArray_q.data)
                      function outputResult(array,output) {
                        var arr = Array(array.slice(0,N));
                      
                        arr = arr.map(function(item){
                          return item.map(function(num){
                            return parseFloat(num.toFixed(4));
                          });
                        });
                        output.value = arr;
                      }
                      outputResult(myArray_i.data, outputDataI);
                      outputResult(myArray_q.data, outputDataQ);
                      
                  }

              } //end onload
              document.body.appendChild(script);
            });
</script>

Time to implement DFT/FFT by myself. Just followed Ifeachor book on theory part and 
IOWA FFT source on implementation side.


## Implementation

### DFT 

Formula for most simple Discreate Fourier Transformation. Also most slowest one.

Official DFT formula

$$X(k) = F_{D}[x(nT)]  = $$     

$$ \sum_{n=0}^{N-1} x(nT) e^{-jk\omega nT} , k=0,1,...,N-1 $$

My wording of DFT algorithm:
Go over data array and sumup multiplication of each "frequency" $$e^{-jk\omega nT}$$ over data array $$x(nT)$$




```c
void dft(double *x_i, double *x_q, int n, int inv) {
    double Wn,Wk;
    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

This is more advanced version of Fourier Transformation, with rearrangement of data there is possible to
reduce amount of calculations and that give speed increase. 
DFT speed is 
$$ O(N^2)$$
and FFT becomes as $$ O(N\log N) $$

#### Shuffling data

Rearranging data before running 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;
       }
}
```


#### FFT 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

Quick verification of FFT and DFT with octave 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)

```

### Demo

Here hookedup demo of compiled wasm fft code location of demo source http://git.main.lv/cgit.cgi/NaiveFFT.git/tree/Build/index.html?h=main

Two boxes for input and two boxes for output. Values set as vectors of complex numbers. First box for real
part, second box for imaginary part or in case of DSP IQ sample values.


Input I:  
<textarea class="inputDataI">1 1 0</textarea><sup class="inputDataIerr"></sup></br> 
Input Q:  
<textarea class="inputDataQ">0 0 0</textarea><sup class="inputDataQerr"></sup></br> 
Output I:
<textarea class="outputDataI"></textarea></br> 
Output Q:
<textarea class="outputDataQ"></textarea></br><button id="calcButton" type="button">Cals</button>


### Source

Source located in main branch of git repo

Browse source http://git.main.lv/cgit.cgi/NaiveFFT.git

```bash
git clone http://git.main.lv/cgit.cgi/NaiveFFT.git
git checkout main
```

### Build

#### Linux

```bash
cd NaiveFFT/Build
make
```

#### Macos

Open with Xcode and build

## Thx

[Aleksejs](https://github.com/ivanovsaleksejs) - gave tips about js  
[Dzuris]() - freedback on code  
[#developerslv](discord.gg/e9t3a9NPBH) - having patiens for listening to js nonsence from js-newbie  

## 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)