c++ - 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;

}

我为 FFTW 使用了 8000 个样本，所以我得到了 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 有两个尖峰。