Naive FFT implementation in C

Intro

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

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

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

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

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:
Input Q:
Output I:
Output Q:

Source

Source located in main branch of git repo

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

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

Thx

Aleksejs - gave tips about js
Dzuris - freedback on code
#developerslv - 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)