Power spectrum from FFTW not working, but in MATLAB it does

Question

I am trying to do a power spectrum analysis with FFTW using this Code below:

#define ALSA_PCM_NEW_HW_PARAMS_API
#include <iostream>
using namespace std;
#include <alsa/asoundlib.h>
#include <fftw3.h>
#include <math.h> 

float map(long x, long in_min, long in_max, float out_min, float out_max)
{
 return (x - in_min) * (out_max - out_min) / (in_max - in_min) + out_min;
}

float windowFunction(int n, int N)
{
return 0.5f * (1.0f - cosf(2.0f *M_PI * n / (N - 1.0f)));
}

int main() {
//FFTW
int N=8000;
float window[N];
double *in = (double*)fftw_malloc(sizeof(double) * N);
fftw_complex *out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
fftw_plan p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE);

for(int n = 0; n < N; n++)
    window[n] = windowFunction(n, N);

//ALSA
long loops;
int rc;
int size;
snd_pcm_t *handle;
snd_pcm_hw_params_t *params;
unsigned int val;
int dir=0;
snd_pcm_uframes_t frames;
char *buffer;

/* Open PCM device for recording (capture). */
rc = snd_pcm_open(&handle, "default",
                SND_PCM_STREAM_CAPTURE, 0);
if (rc < 0) {
fprintf(stderr,
        "unable to open pcm device: %s\n",
        snd_strerror(rc));
exit(1);
}

/* Allocate a hardware parameters object. */
snd_pcm_hw_params_alloca(&params);

/* Fill it in with default values. */
snd_pcm_hw_params_any(handle, params);

/* Set the desired hardware parameters. */

/* Interleaved mode */
snd_pcm_hw_params_set_access(handle, params,
                  SND_PCM_ACCESS_RW_INTERLEAVED);

/* Signed 16-bit little-endian format */
snd_pcm_hw_params_set_format(handle, params,
                          SND_PCM_FORMAT_S16_LE);

/* One channel (mono) */
snd_pcm_hw_params_set_channels(handle, params, 1);

/* 8000 bits/second sampling rate  */
val = 8000;
snd_pcm_hw_params_set_rate_near(handle, params,
                              &val, &dir);

/* Set period size to 16 frames. */
frames = 16;
snd_pcm_hw_params_set_period_size_near(handle,
                          params, &frames, &dir);

/* Write the parameters to the driver */

rc = snd_pcm_hw_params(handle, params);
if (rc < 0) {
fprintf(stderr,
        "unable to set hw parameters: %s\n",
        snd_strerror(rc));
exit(1);
}

/* Use a buffer large enough to hold one period */
snd_pcm_hw_params_get_period_size(params,
                                  &frames, &dir);
size = frames * 2; /* 2 bytes/sample, 1 channel */
buffer = (char *) malloc(size);
/* We want to loop for 5 seconds */
snd_pcm_hw_params_get_period_time(params,
                                     &val, &dir);
loops = 1000000 / val + 25; //added this, because the first values seem to be useless
int count=0;
while (loops > 0) {
loops--;

rc = snd_pcm_readi(handle, buffer, frames);

int i;
short *samples = (short*)buffer;

for (i=0;i < 16;i++)
{
if(count>24){
//cout << (float)map(*samples, -32768, 32768, -1, 1) << endl;
in[i*count]= /*window[i]*/*(double)map(*samples, -32768, 32768, -1, 1);
}
samples++;
}
count++;
if (rc == -EPIPE) {
  /* EPIPE means overrun */
  fprintf(stderr, "overrun occurred\n");
  snd_pcm_prepare(handle);
} else if (rc < 0) {
  fprintf(stderr,
          "error from read: %s\n",
          snd_strerror(rc));
} else if (rc != (int)frames) {
  fprintf(stderr, "short read, read %d frames\n", rc);
}
//    rc = write(1, buffer, size);
//  if (rc != size)
//  fprintf(stderr,
  //        "short write: wrote %d bytes\n", rc);

}

snd_pcm_drain(handle);
snd_pcm_close(handle);
free(buffer);


//FFTW
fftw_execute(p);

for(int j=0;j<N/2;j++){
//cout << in[j] << endl;
cout << sqrt(out[j][0]*out[j][0]+out[j][1]*out[j][1])/N << endl;
/*if(out[j][1]<0.0){
cout << out[j][0] << out[j][1] << "i" << endl;
}else{
cout << out[j][0] << "+" << out[j][1] << "i" << endl;
}*/
}

fftw_destroy_plan(p);

fftw_free(in); 
fftw_free(out);
fftw_cleanup();
return 0;

}

I am using 8000 Samples for the FFTW, so I get 4000 values back, which should be the power spectrum. If I now plot the data in MATLAB, the plot does not look like a power spectrum. The input must be right, because if I uncomment this

//cout << (float)map(*samples, -32768, 32768, -1, 1) << endl;

and comment that

cout << sqrt(out[j][0]*out[j][0]+out[j][1]*out[j][1])/N << endl;

and now load the output of the program (which is the input for the FFT) into MATLAB and do a FFT, the plotted data seems to be correct. I tested it with various frequencies, but I always get a weird spectrum when using my own program. As you can see, I also tried to add a hanning window before the FFT, but still no success. So what am I doing wrong here?

Thanks a lot!

You should separate your FFT code from your sound acquisition code (put them in separate functions), so that you can test/debug the FFT code in isolation. — Paul R
Strongly second @PaulR’s suggestion. Write a simple stand-alone program that just loads some data from disk, calls FFTW (as in your code above), and saves the power spectrum. Then you can easily compare results with Matlab to see where in the processing chain the C++ code and Matlab start disagreeing. Debugging it like this is too hard. — Ahmed Fasih

MSalters MSalters · Accepted Answer · 2017-01-12T15:00:47

The usual suspect with FFT routines is a representation mismatch issue. That is to say, the FFT function fills in an array using one type, and you interpret that array as another type.

You debug this by creating a sine input. You know that this should give a single non-zero input, and you have a reasonable expectation where the zero should be. Because your implementation is wrong, the actual FFT of your sine will differ, and it's this difference that will help troubleshoot the problem.

If you can't figure it out from just the FFT of a sine, next try a cosine, sine of different frequencies, and combinations of two such simple inputs. A cosine is merely a phase shifted sine, so that should just change the phase of that single non-zero value. And the FF should be linear, so the FF of a sum of sines has two sharp peaks.

Power spectrum from FFTW not working, but in MATLAB it does

1 Answers