FFT en C C++

Après avoir comparé les FFT en python et en C C++, il n’y a pas photo… La library C la plus puissante est FFTW. Le rapport de puissance entre 2 programmes équivalents python/C est largement en faveur de C que ce soit en charge processeur ou utilisation de la RAM. Ces critères sont très importants car sur ma raspberry le programme python remplissait la RAM jusqu’à la bloquer.

Voici le site de FFTW.

Et voici comment j’ai réalisé mes FFT en C en sortit de mon RTLSDR:

#include <fftw3.h>

#define CAPTURE 1000 //taille de la FFT

fftw_complex * in;
fftw_complex * out;
fftw_plan p;


void setFFT() {
    in = (fftw_complex * ) fftw_malloc(sizeof(fftw_complex) * CAPTURE + 1);
    out = (fftw_complex * ) fftw_malloc(sizeof(fftw_complex) * CAPTURE + 1);
    p = fftw_plan_dft_1d(CAPTURE, in , out, FFTW_FORWARD, FFTW_ESTIMATE);
}



double * fft(uint8_t * data, int * count, double * power_spectrum) {

    double * power_spectrum_tmp = (double * ) malloc(sizeof(double) * CAPTURE + 1);
    fftw_complex *cx_power_spectrum = (fftw_complex * ) malloc(sizeof(fftw_complex) * CAPTURE + 1);
    int w = 0;
    int x = 0;

    //Transform to complex data, rtlsdr receive 2 numbers for each iteration.
    //One and another with phase a 90°
    //So real data lenght is 2x complex data lenght
    for (int i = 0; i < CAPTURE * 2; i+=2) {
        in[x][0] = (data[i]/(255.0/2.0))-1.0; 
        in[x][1] = (data[i+1]/(255.0/2.0))-1.0;
        x++;
    }

    fftw_execute(p);

    // c = conj(c)*c
    for (int i = 0; i< CAPTURE; i++){
        double complex complex_tmp = out[i][0] + out[i][1]*I;
        complex_tmp = conj(complex_tmp) * complex_tmp;
        cx_power_spectrum[w][0] = creal(complex_tmp) ;
        cx_power_spectrum[w][1] =  cimag(complex_tmp);
        w++;
    }
    
    //complex to double
    for (int i = 0; i < CAPTURE; i++)
        power_spectrum_tmp[i] = cx_power_spectrum[i][0];

    //reorganize data
    for (int i = 0; i < CAPTURE/2; i++)
        power_spectrum[i] = power_spectrum_tmp[i+(CAPTURE/2)];
    for (int i = 0; i < CAPTURE/2; i++)
        power_spectrum[i+(CAPTURE/2)] = power_spectrum_tmp[i];

    free(power_spectrum_tmp);
    free(cx_power_spectrum);


    return power_spectrum;
}