【发布时间】:2021-01-20 15:58:27
【问题描述】:
在我的程序中,我有一个执行快速傅里叶变换的函数。我知道有很好的免费实现,但这是一个学习的东西,所以我不想使用这些。我最终找到了 this comment 的以下实现(它起源于 FFT 的意大利条目):
void transform(complex<double>* f, int N) //
{
ordina(f, N); //first: reverse order
complex<double> *W;
W = (complex<double> *)malloc(N / 2 * sizeof(complex<double>));
W[1] = polar(1., -2. * M_PI / N);
W[0] = 1;
for(int i = 2; i < N / 2; i++)
W[i] = pow(W[1], i);
int n = 1;
int a = N / 2;
for(int j = 0; j < log2(N); j++) {
for(int k = 0; k < N; k++) {
if(!(k & n)) {
complex<double> temp = f[k];
complex<double> Temp = W[(k * a) % (n * a)] * f[k + n];
f[k] = temp + Temp;
f[k + n] = temp - Temp;
}
}
n *= 2;
a = a / 2;
}
free(W);
}
到目前为止,我已经做出了很多改变,但这是我的出发点。我所做的更改之一是不缓存旋转因子,因为我决定先看看是否需要它。现在我决定要缓存它们。这个实现的方式似乎是它有这个长度为N/2 的数组W,其中每个索引k 的值都是。我不明白的是这个表达:
W[(k * a) % (n * a)]
请注意,n * a 始终等于 N/2。我知道这应该等于,我可以看到,这是依赖的。我也知道可以在这里使用模数,因为旋转因子是循环的。但有一件事我不明白:这是一个长度为 N 的 DFT,但只有 N/2 旋转因子被计算过。数组的长度不应该是N,模数应该是N吗?
【问题讨论】:
-
基本上它在 2 行之后的作用
标签: c++ algorithm math fft dft