FFTW 根本不会返回除无穷大、接近零的值或负无穷大以外的值答案

【问题标题】：FFTW simply will not return values other than infinity, values that approach zero, or negative infinityFFTW 根本不会返回除无穷大、接近零的值或负无穷大以外的值
【发布时间】：2014-09-11 10:26:51
【问题描述】：

我正在尝试计算 20hz 的正弦波的 DFT。

首先我用 20 赫兹的正弦函数的 10 个周期填充一个 std::vector：

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 的一切。这是输入的所有 5 个周期的完整输出，显示了输入、输出和 hz 值：

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;

}

【问题讨论】：

您的输出函数只输出 500 个结果中的前 50 个。似乎您正试图输出每 5 个结果或其他内容。
请将 N 设为整数类型。
我只输出了输入的第一个完整周期。现在代码输出所有周期，我创建了一个输出要点，显示了周期的输入、输出和相应的 hz 值。

标签： c++ signal-processing fft fftw dft

【解决方案1】：

您分配了一个fftw_complexes 数组，它由两个元素组成——实部和复部——但您只是在初始化每个样本的实部。这可能会留下包含随机数据的复杂组件，导致意想不到的疯狂结果！

如果您不需要处理复杂的样本（很可能），您可能希望使用one of FFTW's real-data DFT functions，它将double 数组作为输入。否则，将所有复杂组件 (in[i][1]) 初始化为零。

此外，您没有查看输出的复杂部分。你可能错过了一些重要的事情。

【讨论】：

好主意，我确信这会解决问题，但无济于事！我确实实施了你的建议。您可以在编辑后的帖子中看到新代码。
检查你得到的转换输出的复杂组件。我认为那里有些东西。
看起来输出的虚部都是乱码，除了第一个周期的 50hz 的随机 -125 值？ gist.github.com/anonymous/ed419f4cfb43887c810e
那是你的人。它显示为一个复杂的分量，因为您的信号在原点周围是反对称的（因为它是正弦波，而不是余弦波）。

【解决方案2】：

您的输入是纯正弦波的 10 个周期。这意味着：

您的输出是纯虚数，由奇正弦函数 (sin(x) = -sin(-x)) 引起在余弦函数的情况下，输出是纯实数因为余弦是偶函数 (cos(x) = cos(-x)
频率线只有两条，一条在正频段一个在负带中。
因为您的输入中有 10 个句点，您的行应该出现在 out[10][1] 和 out[N-10][1]。
因为您没有对输出进行归一化并且样本数为 500，幅度应该是 out[10][1] = -250 和 out[490][1] = +250。
所有其他值都应为零。

【讨论】：