simple_fft.c

#include <stdio.h>
#include <math.h>
 
#include <sp/spLib.h>
#include <sp/spMain.h>
 
#define FFT_ORDER 10
#define SIN_NCYCLE 5
 
#define CALCULATE_DISTORTION
 
#if 1
#define FFT_PLUGIN_NAME "oourafft"
#elif 0
#define FFT_PLUGIN_NAME "fftw"
#elif 0
#define FFT_PLUGIN_NAME "cufft"
#else
#define FFT_PLUGIN_NAME ""
#endif
 
static char plugin_name[SP_MAX_PATHNAME] = FFT_PLUGIN_NAME;
 
int spMain(int argc, char **argv)
{
    long k;
    long fftl;
    double theta;
    double *real, *imag;
    spFFTRec fftrec;
    spPlugin *plugin = NULL;
#ifdef CALCULATE_DISTORTION
    double rx, ix, dist;
    double *real2, *imag2;
#endif
 
    if (!spStrNone(plugin_name)) {
        if ((plugin = spLoadFFTPlugin(plugin_name)) == NULL) {
            fprintf(stderr, "Can't open FFT plugin: %s\n", plugin_name);
        } else {
            fprintf(stderr, "FFT plugin: %s\n", plugin_name);
        }
    }
    
    fftrec = spInitFFTByPlugin(plugin, FFT_ORDER, SP_FFT_DEFAULT_PRECISION);
    fftl = spGetFFTLength(fftrec);
        
    real = xspAlloc(fftl, double);
    imag = xspAlloc(fftl, double);
 
    for (k = 0; k < fftl; k++) {
        theta = 2.0 * PI * SIN_NCYCLE * (double)k / (double)(fftl - 1);
        real[k] = sin(theta);
        //printf("%ld %f\n", k, real[k]);
        imag[k] = 0.0;
    }
 
#ifdef CALCULATE_DISTORTION
    real2 = xspAlloc(fftl, double);
    imag2 = xspAlloc(fftl, double);
    memcpy(real2, real, fftl * sizeof(double));
    memcpy(imag2, imag, fftl * sizeof(double));
#endif
    
    spExecFFT(fftrec, real, imag, 0);
    
    for (k = 0; k < fftl; k++) {
        printf("%f %f\n", real[k], imag[k]);
    }
 
#ifdef CALCULATE_DISTORTION
    spfft(real2, imag2, fftl, 0);
    dist = 0.0;
    for (k = 0; k < fftl; k++) {
        rx = real[k] - real2[k];
        ix = imag[k] - imag2[k];
        dist += (SQUARE(rx) + SQUARE(ix));
        //printf("%ld %f %f %f %f\n", k, real[k], imag[k], real2[k], imag2[k]);
    }
    dist = sqrt(dist / (double)fftl);
    printf("dist = %f\n", dist);
    xspFree(real2);
    xspFree(imag2);
#endif
    
    xspFree(real);
    xspFree(imag);
    
    spFreeFFT(fftrec);
    
    return 0;
}