功率譜不是工作，而是在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; 

現在程序的輸出（這對FFT輸入）加載到MATLAB和做FFT，繪製的數據似乎是正確的。我用各種頻率對它進行了測試，但在使用我自己的程序時，我總是會得到一個奇怪的頻譜。正如你所看到的，我也嘗試在FFT之前添加一個漢寧窗，但仍然沒有成功。那麼我在這裏做錯了什麼？

非常感謝！

Answer 1

FFT例程的常見嫌疑犯是表示不匹配問題。也就是說，FFT函數使用一種類型填充數組，並將該數組解釋爲另一種類型。

您可以通過創建正弦輸入進行調試。你知道這應該給出一個非零的輸入，並且你有一個合理的期望值應該是零。由於你的實現是錯誤的，你的正弦波的實際FFT會有所不同，正是這種差異將有助於排除故障。

如果你不能從正弦的FFT中找出它，接下來嘗試一個餘弦，不同頻率的正弦和兩個這樣的簡單輸入的組合。餘弦只是一個相移正弦，所以應該只改變那個單一的非零值的相位。 FF應該是線性的，所以正弦和的FF有兩個尖峯。