FFTW 的功率谱不起作用,但在 MATLAB 中确实如此

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

本文关键字:MATLAB 但在 功率 不起作用 FFTW      更新时间:2023-10-16

我正在尝试使用以下代码对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: %sn",
        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: %sn",
        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 occurredn");
  snd_pcm_prepare(handle);
} else if (rc < 0) {
  fprintf(stderr,
          "error from read: %sn",
          snd_strerror(rc));
} else if (rc != (int)frames) {
  fprintf(stderr, "short read, read %d framesn", rc);
}
//    rc = write(1, buffer, size);
//  if (rc != size)
//  fprintf(stderr,
  //        "short write: wrote %d bytesn", 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之前添加一个汉宁窗口,但仍然没有成功。那么我在这里做错了什么?

多谢!

FFT 例程的常见嫌疑人是表示不匹配问题。也就是说,FFT 函数使用一种类型填充数组,您将该数组解释为另一种类型。

您可以通过创建正弦输入来调试此内容。您知道这应该给出单个非零输入,并且您有一个合理的期望零应该在哪里。因为你的实现是错误的,所以你的正弦的实际FFT会有所不同,正是这种差异将有助于解决问题。

如果你不能仅仅从正弦的FFT中弄清楚,接下来尝试一个余弦,不同频率的正弦,以及两个这样的简单输入的组合。余弦只是一个相移正弦,因此应该只改变单个非零值的相位。并且FF应该是线性的,所以正弦和的FF有两个尖锐的峰值。

相关文章: