这是我的 DFT 实现,使用 {0, 1, 2, 3} 进行测试。
#include <vector>
#include <complex>
#include <numbers>
#include <cmath>
#include <iostream>
std::vector<std::complex<double>> DFT(std::vector<std::complex<double>>&& P)
{
int n = P.size();
if (n == 1) {
return P;
}
std::vector<std::complex<double>> Pe(n / 2);
std::vector<std::complex<double>> Po(n / 2);
for (int i = 0; i < n / 2; ++i) {
Pe[i] = P[2 * i];
Po[i] = P[2 * i + 1];
}
auto ye = DFT(std::move(Pe)), yo = DFT(std::move(Po));
// in place algorithm, use input P to store output
auto wi = std::complex<double>(1.0, 0.0);
auto wn = std::complex<double>(std::cos(2.0 * std::numbers::pi_v<double> / (double)n),
std::sin(2.0 * std::numbers::pi_v<double> / (double)n));
for (int i = 0; i < n / 2; ++i) {
P[i] = ye[i] + wi * yo[i];
P[i + n / 2] = ye[i] - wi * yo[i];
wi = wi * wn;
}
return P;
}
int main()
{
int N = 4;
std::vector<std::complex<double>> p_input(N);
for (int i = 0; i < N; ++i)
{
p_input[i] = {(double)i, 0.0};
}
auto p_output = DFT(std::move(p_input));
for (int i = 0; i < p_output.size(); ++i)
{
std::cout << p_output[i] << std::endl;
}
}
测试结果是
(6,0)
(-2,-2)
(-2,0)
(-2,2)
但是 MATLAB 的答案是
>> fft(0:1:3)
ans =
6.0000 + 0.0000i -2.0000 + 2.0000i -2.0000 + 0.0000i -2.0000 - 2.0000i
我测试了更长的输入长度,并且我的结果总是在奇数索引位置具有相反的顺序。
不知道哪一部分是错误的?我的参考公式是:
P(x): [p_0, p_1, ..., p_{n-1}]
w: [w^0, w^1, ..., w^{n-1}]
Pe(x^2): [p_0, p_2, ..., p_{n-2}]
Po(x^2): [p_1, p_3, ..., p_{n-1}]
ye = [Pe(w^0), Pe(w^2), ..., Pe(w^{n-2})]
yo = [Po(w^0), Po(w^2), ..., Po(w^{n-2})]
P(w^j) = ye[j] + w^j yo[j]
P(w^{j+n/2}) = ye[j] - w^j yo[j]
y = [P(w^0), P(w^1), ..., P(w^{n-1})]
您只是使用了与其实现不同的 n 次单位根,即
cos(2*pi / n) + sin(2*pi / n) * i
而不是
cos(2*pi / n) - sin(2*pi / n) * i
就FFT而言,你使用哪一种并不重要,只要逆FFT与其一致即可。我对 FFT 不太熟悉,不知道按照惯例,其中一种方法是否优于另一种方法。