[Help-gsl] mixed-radix fft

[Sean Parsons]

Hi, 

I was using the mixed-radix routines to compute the FFT of some data.

However I noticed that the "halfcomplex" output of gsl_fft_real_transform is 
different according to whether the no. of samples (nsample) is a power of two 
or not. With nsample a power of 2, I get a sensible result, but with nsample 
not a power of two I get rubbish.

I have illustrated this with the program below (also attached). I input a sum 
of three sine waves. With nsample = 4096 I get three nice peaks at the correct 
frequencies, but with nsample = 5000, I get just a sequence of small values.

Could you help me? I realise that upon the reverse transform, I get my original 
data, whether nsample = 4096 or 5000. But I need to know how the forward 
transformed data is represented so I can apply filters.

Thank you very much,

Sean Parsons

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
#include <cstdlib>
#include <iostream>
#include <cmath>
#include "gsl/gsl_math.h"
#include "gsl/gsl_fft_real.h"
#include "gsl/gsl_fft_halfcomplex.h"

using namespace std;

int main(int argc, char *argv[])
{
    int nsample = 4096;      //no. of samples
    double rsample = 1;      //sampling rate

    double f1, f2, f3, m1, m2, m3;      //sine components: frequencies and 
magnitudes
    f1 = 0.01;
    f2 = 0.05;
    f3 = 0.1;
    m1 = 10;
    m2 = 15;
    m3 = 40;

    //create input data (sum of sine waves)
    double *data = (double*) new double[nsample];
    double time = 0;
    for(int i = 0;i<nsample;i++) {
        time = (double)i / rsample;
        data[i] = (m1 * sin(time*f1*2*M_PI))
                        + (m2 * sin(time*f2*2*M_PI))
                        + (m3 * sin(time*f3*2*M_PI));
    }

    //allocate wavetable + workspace
    gsl_fft_real_wavetable *realWT = gsl_fft_real_wavetable_alloc (nsample);
    gsl_fft_real_workspace *realWS = gsl_fft_real_workspace_alloc (nsample);

    //forward transform
    gsl_fft_real_transform(data,1,nsample,realWT,realWS);

    //free wavetable + workspace
    gsl_fft_real_wavetable_free (realWT);
    gsl_fft_real_workspace_free (realWS);

    //output real component of transform (magnitude) according to FFTPACK 
storage.
    double *copy = (double*) new double[nsample/2 +1];
    copy[0] = pow(data[0],2);
    for (int i = 1; i <= floor(double(nsample)/2.0); i++){
            copy[i] = pow(data[(i*2)-1],2);
            cout<<"freq = "<<(double)i*rsample/nsample<<"\t\t mag = 
"<<copy[i]<<endl;
    }

    system("PAUSE");
    return EXIT_SUCCESS;
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
 Sean Parsons
McMaster University
Intestinal Disease Research Program
1200 Main St. West, Room 3N6-9,
Hamilton,
Ontario, L8N 3Z5.
Tel: 905-966-2540
address@hidden
address@hidden

main_fft3.cpp
Description: Binary data