simple_rfft.c

#include <stdio.h>
#include <math.h>
 
#include <sp/spLib.h>
#include <sp/spMain.h>
 
#define FFT_ORDER 8
#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, fftl2;
    double theta;
    double *real;
    spFFTRec fftrec;
    spPlugin *plugin = NULL;
#ifdef CALCULATE_DISTORTION
    double rx, ix, dist;
    double *real2;
#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);
    fftl2 = fftl / 2;
        
    real = xspAlloc(fftl, double);
 
    for (k = 0; k < fftl; k++) {
        theta = 2.0 * PI * SIN_NCYCLE * (double)k / (double)(fftl - 1);
        real[k] = sin(theta);
    }
 
#ifdef CALCULATE_DISTORTION
    real2 = xspAlloc(fftl, double);
    memcpy(real2, real, fftl * sizeof(double));
#endif
    
    spExecRealFFT(fftrec, real, 0);
    
    printf("%f %f\n", real[0], 0.0);
    for (k = 1; k < fftl2; k++) {
        printf("%f %f\n", real[2 * k], real[2 * k + 1]);
    }
    printf("%f %f\n", real[1], 0.0);
 
#ifdef CALCULATE_DISTORTION
    sprfft(real2, fftl, 0);
    dist = 0.0;
    rx = real[0] - real2[0];
    dist += SQUARE(rx);
    //printf("%ld %f %f %f %f\n", 0L, real[0], 0.0, real2[0], 0.0);
    for (k = 1; k < fftl2; k++) {
        rx = real[2 * k] - real2[2 * k];
        ix = real[2 * k] - real2[2 * k];
        dist += 2.0 * (SQUARE(rx) + SQUARE(ix));
        //printf("%ld %f %f %f %f\n", k, real[2 * k], real[2 * k + 1], real2[2 * k], real2[2 * k + 1]);
    }
    rx = real[1] - real2[1];
    dist += SQUARE(rx);
    //printf("%ld %f %f %f %f\n", k, real[1], 0.0, real2[1], 0.0);
    
    dist = sqrt(dist / (double)fftl);
    printf("dist = %f\n", dist);
    xspFree(real2);
#endif
    
    xspFree(real);
    
    spFreeFFT(fftrec);
    
    return 0;
}