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

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:

Logarithmic:

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.

Code: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 }

## Bookmarks