uOpt

May 21st, 2010, 07:46 PM

I am trying to get spectrum data (the numbers, not the graphics) out of FFTW, like this (from Audacity):

http://wavehh.dyndns.org:3734/tmp/fft_auda.png

After studying fftw's documentation and all Google links I could find I am still unable to interpret the output of FFTW correctly. I just don't know where to find which frequency in there.

Here is what I got. It is obviously center-plotted but let's leave that aside, it doesn't look any good if you shift it.

Linear X scale:

http://wavehh.dyndns.org:3734/tmp/fft1.png

Logarithmic:

http://wavehh.dyndns.org:3734/tmp/fft2.png

Any idea what I am doing wrong here? The code is a routine from a larger program, I can make a runnable thing available if you like.

struct frames {

int n_frames;

int n_files;

char **files;

frames(): n_frames(-1), n_files(-1){}

};

typedef double fftw_input_type;

typedef double fftw_output_type;

typedef double fftw_spectrum_type;

static int fft(const struct frames &frames)

{

int n_points = frames.n_frames;

n_points = 50000; // short sample

fprintf(stderr, "In fft on %d samples\n", n_points);

fftw_input_type *in;

fftw_output_type *out;

fftw_complex *comout;

fftw_plan plan;

fftw_plan cplan;

fftw_spectrum_type *spectrum;

in = (fftw_input_type *) fftw_malloc(sizeof(fftw_input_type) * n_points);

out = (fftw_output_type *) fftw_malloc(sizeof(fftw_output_type) * n_points);

comout = (fftw_complex *) fftw_malloc(sizeof(*comout) * n_points);

spectrum = (fftw_spectrum_type *)malloc(sizeof(*spectrum) * n_points / 2 + 1);

plan = fftw_plan_r2r_1d(n_points, in, out, FFTW_R2HC, FFTW_FORWARD);

cplan = fftw_plan_dft_r2c_1d(n_points, in, comout, FFTW_FORWARD);

channeldata *data1 = (channeldata*)frames.files[0];

for (int i = 0; i < n_points; i++) {

// fixme, channel 0 only

in[i] = (double)data1[i].channels[0];

}

fftw_execute(plan); /* repeat as needed */

fftw_execute(cplan); /* repeat as needed */

spectrum[0] = out[0] * out[0];

int n_spectrum_points = (n_points + 1) / 2;

for (int k = 1; k < n_spectrum_points; k++) /* (k < N/2 rounded up) */

spectrum[k] = out[k] * out[k] + out[n_points - k] * out[n_points - k];

if (n_points % 2 == 0) /* N is even */

spectrum[n_points / 2] = out[n_points / 2] * out[n_points / 2]; /* Nyquist freq. */

#if 0

for (int i = 0; i < n_spectrum_points; i++) {

printf("Spectrum: %f\n", spectrum[i]);

}

#endif

double correction = (double)sample_frequency / (double)n_points / 2;

double accu = 0;

double freq1 = -1;

double freq2 = -1;

const int every_nth = n_points / 1000; // want this number of lines

// rotate output array by half

int max = n_points - 2;

for (int k = 0; k < max; k++) {

int index_in_output = k + max / 2;

if (index_in_output > max)

index_in_output -= max;

double absval = sqrt(out[index_in_output] * out[index_in_output]);

double cc = (double)k * correction;

double maybe_freq = (double)k / (double)n_points;

if (k % every_nth == 0) {

if (k != 0) {

freq2 = cc;

// FIXME - last value documented to be lost

printf("%f %f %f %f ( %f - %f ) %f\n"

, (freq2 + freq1) / 2.0

, absval / every_nth

, sqrt(comout[index_in_output][0] * comout[index_in_output][0])

, sqrt(comout[index_in_output][1] * comout[index_in_output][1])

, freq1, freq2

, maybe_freq

);

}

freq1 = cc;

accu = 0;

}

accu += absval;

}

fftw_destroy_plan(plan);

fftw_free(in);

fftw_free(out);

return 0; // success, this goes to exit

}

http://wavehh.dyndns.org:3734/tmp/fft_auda.png

After studying fftw's documentation and all Google links I could find I am still unable to interpret the output of FFTW correctly. I just don't know where to find which frequency in there.

Here is what I got. It is obviously center-plotted but let's leave that aside, it doesn't look any good if you shift it.

Linear X scale:

http://wavehh.dyndns.org:3734/tmp/fft1.png

Logarithmic:

http://wavehh.dyndns.org:3734/tmp/fft2.png

Any idea what I am doing wrong here? The code is a routine from a larger program, I can make a runnable thing available if you like.

struct frames {

int n_frames;

int n_files;

char **files;

frames(): n_frames(-1), n_files(-1){}

};

typedef double fftw_input_type;

typedef double fftw_output_type;

typedef double fftw_spectrum_type;

static int fft(const struct frames &frames)

{

int n_points = frames.n_frames;

n_points = 50000; // short sample

fprintf(stderr, "In fft on %d samples\n", n_points);

fftw_input_type *in;

fftw_output_type *out;

fftw_complex *comout;

fftw_plan plan;

fftw_plan cplan;

fftw_spectrum_type *spectrum;

in = (fftw_input_type *) fftw_malloc(sizeof(fftw_input_type) * n_points);

out = (fftw_output_type *) fftw_malloc(sizeof(fftw_output_type) * n_points);

comout = (fftw_complex *) fftw_malloc(sizeof(*comout) * n_points);

spectrum = (fftw_spectrum_type *)malloc(sizeof(*spectrum) * n_points / 2 + 1);

plan = fftw_plan_r2r_1d(n_points, in, out, FFTW_R2HC, FFTW_FORWARD);

cplan = fftw_plan_dft_r2c_1d(n_points, in, comout, FFTW_FORWARD);

channeldata *data1 = (channeldata*)frames.files[0];

for (int i = 0; i < n_points; i++) {

// fixme, channel 0 only

in[i] = (double)data1[i].channels[0];

}

fftw_execute(plan); /* repeat as needed */

fftw_execute(cplan); /* repeat as needed */

spectrum[0] = out[0] * out[0];

int n_spectrum_points = (n_points + 1) / 2;

for (int k = 1; k < n_spectrum_points; k++) /* (k < N/2 rounded up) */

spectrum[k] = out[k] * out[k] + out[n_points - k] * out[n_points - k];

if (n_points % 2 == 0) /* N is even */

spectrum[n_points / 2] = out[n_points / 2] * out[n_points / 2]; /* Nyquist freq. */

#if 0

for (int i = 0; i < n_spectrum_points; i++) {

printf("Spectrum: %f\n", spectrum[i]);

}

#endif

double correction = (double)sample_frequency / (double)n_points / 2;

double accu = 0;

double freq1 = -1;

double freq2 = -1;

const int every_nth = n_points / 1000; // want this number of lines

// rotate output array by half

int max = n_points - 2;

for (int k = 0; k < max; k++) {

int index_in_output = k + max / 2;

if (index_in_output > max)

index_in_output -= max;

double absval = sqrt(out[index_in_output] * out[index_in_output]);

double cc = (double)k * correction;

double maybe_freq = (double)k / (double)n_points;

if (k % every_nth == 0) {

if (k != 0) {

freq2 = cc;

// FIXME - last value documented to be lost

printf("%f %f %f %f ( %f - %f ) %f\n"

, (freq2 + freq1) / 2.0

, absval / every_nth

, sqrt(comout[index_in_output][0] * comout[index_in_output][0])

, sqrt(comout[index_in_output][1] * comout[index_in_output][1])

, freq1, freq2

, maybe_freq

);

}

freq1 = cc;

accu = 0;

}

accu += absval;

}

fftw_destroy_plan(plan);

fftw_free(in);

fftw_free(out);

return 0; // success, this goes to exit

}