我将尝试计算一个20hz正弦波的离散傅里叶变换(DFT)。
首先,我使用正弦函数在20 hz下填充了一个std :: vector中的10个周期:
std::vector<double> sinx;
double samplerate = 1000.0;
double frequency = 20.0;
double num_cycles = 10.0;
for (int i=0; i<samplerate/frequency*num_cycles; i++){
sinx.push_back(sin(2.0*M_PI*((double)i)*frequency/samplerate));
}
double N = sinx.size();
然后我初始化FFTW的输入和输出数组以及FFTW计划:
fftw_complex *in, *out;
fftw_plan p;
in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)* N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
然后,我使用正弦波函数的值来填充输入数组:
for (int i=0; i<N; i++){
in[i][0] = sinx[i];
}
然后我计算FFT:
fftw_execute(p);
计算完成FFT后,我会检查结果并进行清理。
double hz_per_index = samplerate/2.0/N*num_cycles;
for (int i=0; i<N/num_cycles; i++){
std::cout<<"hz: "<<hz_per_index*i<<" out: "<<fabs(out[i][0])<<" "<<out[i][1]<<" in: "<<in[i][0]<<std::endl;
}
我预期会在20hz处看到一个峰值,但是却没有任何峰值。所有返回的值要么趋近于零,要么是无穷大或负无穷大。我的第一个想法是输入正在被覆盖,但在计算后检查输入是正确的。我尝试以不同的模式运行FFTW,包括FFTW_FORWARD、FFTW_BACKWARD、FFTW_ESTIMATE和FFTW_RELATIVE,但没有什么帮助。我尝试在不同的频率、采样率和循环次数下计算不同的正弦函数。
可悲的是,我曾经让它工作过。然后我继续工作,发现它不起作用了,然后重新打开了那个工作的文件,现在它又不起作用了!
还有其他人遇到这个问题吗?
编辑:根据建议的新代码
它仍然不起作用。我尝试将复数的虚部初始化为0,但仍然存在同样的问题,所以这里是利用fftw库的实入/复出函数的代码。
基本上,我得到的一切都是0。 以下是显示输入、输出和hz值的所有5个周期的完整输出:
https://gist.github.com/anonymous/ed419f4cfb43887c810e
double *in;
fftw_complex* out;
fftw_plan p;
std::vector<double> sinx;
double samplerate = 1000.0;
double frequency = 20.0;
double num_cycles = 5.0;
for (int i=0; i<samplerate/frequency*num_cycles; i++){
sinx.push_back(sin(2.0*M_PI*((double)i)*frequency/samplerate));
}
int N = sinx.size();
in = (double*) fftw_malloc(sizeof(double)* N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
p = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);
for (int i=0; i<N; i++){
in[i] = sinx[i];
}
fftw_execute(p);
double hz_per_index = samplerate/2.0/N*num_cycles;
for (int i=0; i<N; i++){
std::cout<<"hz: "<<hz_per_index * (i % (int)(N/num_cycles))<<" out: "<<out[i][0]<<" in: "<<in[i]<<std::endl;
}
fftw_destroy_plan(p);
fftw_free(in); fftw_free(out);
return 0;
}