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)