Спектр мощности от FFTW не работает, но в MATLAB он делает

Question

Я пытаюсь сделать анализ спектра мощности с FFTW, используя этот код ниже:

#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; 

} 

Я использую 8000 образцов для FFTW, так что я получаю 4000 значений, которые должны быть спектром мощности. Если я сейчас буду строить данные в MATLAB, график не будет выглядеть как спектр мощности. Вход должен быть правильным, потому что если я раскомментировать эту

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

и комментарии, которые

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

и теперь загрузить вывод программы (которая является входом для быстрого преобразования Фурье) в MATLAB и делать FFT, построенные данные кажутся правильными. Я тестировал его с различными частотами, но при использовании моей собственной программы я всегда получаю странный спектр. Как вы можете видеть, я также попытался добавить окно Hanning перед FFT, но все равно не получилось. Так что я делаю неправильно здесь?

Большое спасибо!

Answer 1

Обычный подозреваемый с подпрограммами FFT является проблемой несоответствия представления. То есть функция FFT заполняет массив с использованием одного типа, и вы интерпретируете этот массив как другой тип.

Вы отлаживаете это, создавая вход синуса. Вы знаете, что это должно дать один ненулевой ввод, и у вас есть разумное ожидание, когда должен быть нуль. Поскольку ваша реализация неверна, фактический FFT вашего синуса будет отличаться, и именно эта разница поможет устранить проблему.

Если вы не можете понять это только из БПФ синуса, попробуйте косинус, синус разных частот и комбинации двух таких простых входов. Косинус - это просто синусоидальный сдвиг фазы, поэтому он должен просто изменить фазу этого единственного ненулевого значения. И FF должен быть линейным, поэтому FF суммы синусов имеет два резких пика.