FFTW 的功率谱不起作用,但在 MATLAB 中确实如此
Power spectrum from FFTW not working, but in MATLAB it does
我正在尝试使用以下代码对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(¶ms);
/* 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有两个尖锐的峰值。
相关文章:
- 代码在main()中运行,但在函数中出现错误
- 链接阶段在Ubuntu上失败,但在MacOS上失败
- 对C宏的未定义引用,但在定义它时会出现重新定义错误
- c++17文件系统::recursive_directory迭代器()在mac上没有给出这样的目录,但在windows上
- 断言中的Fold表达式在某些计算机上编译,但在其他计算机上不编译
- 换位表导致测试失败(但在游戏中运行良好)
- 库标题在标题中不可见,但在 cmake build 下.cpp文件中完全可见.为什么?
- 树莓上的 Libtorch 无法加载 pt 文件,但在 ubuntu 上工作
- 在成员dynamic_bitset上使用 boost::from_block_range 时出错,但在本地dynamic
- 编译在我的 Mac 上工作,但在集群 (Linux) 上不起作用
- 我编写了代码将十进制分数转换为其二进制等效数.它编译得很好,但在执行时挂起
- 我的代码运行良好,但在游戏循环中中断
- C++ assigment std::list:<typename>:itrator 在 main 中工作,但在方法中它不起作用
- C++代码在台式机上工作正常,但在笔记本电脑上则不行
- 实现 DFS 在较短的输入下工作正常,但在较大的输入下会抛出分段错误
- 点云库在VS 2019中不起作用,但在VS 2017中确实有效
- C++ Python 模块在 Blender 中崩溃,但在 Python 控制台中不会崩溃
- 分段 Linux Ubuntu 中的 g++ 错误,但在 Windows 中的 g++/MingW 中,在 C++ 中打
- FFTW 的功率谱不起作用,但在 MATLAB 中确实如此
- MATLAB + Mex + OpenCV:正确链接和编译,但在运行时找不到库