FFT-C

fft.c

#include <stdio.h>
#include <math.h>
#include <time.h>

//algorithm taken from
//https://rosettacode.org/wiki/Fast_Fourier_transform#C
//simplified version to remove double complex numbers dependency to see exactly what the code does
//(N log(N)) - complexity

typedef float fftn;
fftn PI;
fftn sampleRate = 44100;
int fftSize = 4096;
int repeats = 100;


void printData(const char *s, fftn buf[], int len) {
    printf("%s", s);
    for (int i = 0; i < len; i += 2)
        printf("%i => (%g, %g) \n", i / 2, buf[i], buf[i + 1]);
}

void _fft(fftn buf[], fftn out[], int n, int step) {
    if (step < n) {
        _fft(out, buf, n, step * 2);
        _fft(out + (2 * step), buf + (2 * step), n, step * 2);

        for (int _i = 0; _i < n; _i += 2 * step) {
            float v1Re = 0.0f;
            float v1Im = -PI * _i / n;

            float r = expf(0.0);
            v1Re = r * cosf(v1Im);
            v1Im = r * sinf(v1Im);

            float v2Re = out[(_i + step) * 2];
            float v2Im = out[((_i + step) * 2) + 1];

            //Complex(re * x.re - im * x.im, re * x.im + im * x.re)
            float rRe = (v1Re * v2Re) - (v1Im * v2Im);
            float rIm = (v1Re * v2Im) + (v1Im * v2Re);

            int i1 = _i;
            int i2 = (_i + n);
            //v1
            int _i1 = _i * 2;
            int _i2 = _i1 + 1;

            buf[i1] = out[_i1] + rRe;
            buf[i1 + 1] = out[_i2] + rIm;
            //v2
            buf[i2] = out[_i1] - rRe;
            buf[i2 + 1] = out[_i2] - rIm;
        }
    }
}

void fft(fftn buf[], int n) {
    fftn out[n * 2];
    for (int i = 0; i < n * 2; i++) out[i] = buf[i];
    _fft(buf, out, n, 1);
}

int main() {
    PI = atan2f(1, 1) * 4;
    double PId = atan2(1, 1) * 4;
    int dataLen = fftSize * 2;
    float buf[fftSize * 2];// = {1, 1, 1, 1, 0, 0, 0, 0};


    for (int i = 0; i < fftSize; i++) {
        float step = (float) i / (float) sampleRate;
        double c = cos(3.0 * step * 2.0 * PId);
        int iv = i * 2;
        //double to float to get same value as java, working with floats, is quite different
        buf[iv] = (fftn) c;
        buf[iv + 1] = 0;
    }

    clock_t start, end;
    double cpu_time_used;
    //printData("Data: ", buf, dataLen);
    start = clock();

    for (int i = 0; i < repeats; i++) {
        fft(buf, fftSize);
    }

    end = clock();

    cpu_time_used = ((double) (end - start));
    printf("FFT takes %f ms in avarage (%i iterations)", cpu_time_used / repeats, repeats);
    if (0) {
        printData("\nFFT Result:\n", buf, dataLen);
    }

    return 0;
}