如何使用FFTW库C++计算3D阵列的FFT?

Question

我在MATLAB中有3D数组，如下所示：

val(:,:,1) =
1.0000 + 1.0000i   2.0000 + 2.0000i
3.0000 + 3.0000i   4.0000 + 4.0000i

val(:,:,2) =
5.0000 + 5.0000i   6.0000 + 6.0000i
7.0000 + 7.0000i   8.0000 + 8.0000i

我在 MATLAB 中采用一维 FFT，如下所示：

clear
close all
%3D ARRAY Init
A3D(:,:,1)=[1+1i,2+2i;3+3i,4+4i];
A3D(:,:,2)=[5+5i,6+6i;7+7i,8+8i];
n=size(A3D,1);
res=fft(A3D,n,1)

FFT 的结果是：

val(:,:,1) =
4.0000 + 4.0000i   6.0000 + 6.0000i
-2.0000 - 2.0000i  -2.0000 - 2.0000i

val(:,:,2) =
12.0000 +12.0000i  14.0000 +14.0000i
-2.0000 - 2.0000i  -2.0000 - 2.0000i

C++我编写了如下程序，我将输入转换为行主要格式并发送到fftw_plan_dft_1d以制作 FFT，但我的结果与 MATLAB 结果不同？有人可以帮助我解决这个问题吗？

#include <cstdio>
#include <fftw3.h>
int main()
{
fftw_complex *in, *out;
fftw_plan p;
int N = 8;
in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
for (int i=0;i<8;i++){
in[i][0]=i+1;
in[i][1]=i+1;
}

for(int i = 0; i < N; i++)
{
printf("%ft%fn", in[i][0],in[i][1]);
}
p = fftw_plan_dft_1d(8, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p); /* repeat as needed */
for(int i = 0; i < N; i++)
{
printf("%ft%fn", out[i][0],out[i][1]);
}

fftw_destroy_plan(p);
fftw_free(in); fftw_free(out);
return 0;
}

Answer 1

用户应注意，FFTW 计算非规范化 DFT，其指数的符号由 dir 参数 fftw_create_plan 给出。因此，计算前向后变换(反之亦然(会导致原始数组按 n 缩放。

所以，你应该规范化。下面是一个规范化和执行向后转换以恢复原始值的示例，您会注意到这些值将起作用(但如果不规范化将无法工作(。但是，您期望的值不是产生的值。我不了解 matlab，所以我无法解释那部分。

#include <fftw3.h>
#include <iomanip>
#include <ios>
#include <iostream>
std::ostream& operator<<(std::ostream& os, const fftw_complex& c) {
return os << std::setw(10) << c[0] << " + " << std::setw(10) << c[1] << 'i';
}
void display(const char* title, const fftw_complex* arr, size_t N) {
std::cout << title << ":n";
for(size_t i = 0; i < N; i++) {
std::cout << arr[i] << 'n';
}
}
void normalize(fftw_complex* arr, size_t N) {
for(size_t bin = 0; bin < N; ++bin) {
arr[bin][0] /= static_cast<double>(N);
arr[bin][1] /= static_cast<double>(N);
}
}
int main() {
fftw_complex *in, *out;
fftw_plan p, p2;
std::cout << std::setprecision(6) << std::fixed;
size_t N = 8;
in = fftw_alloc_complex(N);
out = fftw_alloc_complex(N);
for(size_t i = 1; i <= N; i++) {
in[i - 1][0] = static_cast<double>(i);
in[i - 1][1] = static_cast<double>(i);
}
display("in", in, N);
p = fftw_plan_dft_1d(static_cast<int>(N), in, out, FFTW_FORWARD, FFTW_ESTIMATE);
p2 = fftw_plan_dft_1d(static_cast<int>(N), out, in, FFTW_BACKWARD, FFTW_ESTIMATE);
fftw_execute(p);
display("out", out, N);
normalize(out, N);
display("normalized", out, N);
fftw_execute(p2);
display("backward", in, N);
fftw_destroy_plan(p2);
fftw_destroy_plan(p);
fftw_free(in);
fftw_free(out);
}

输出：

in:
1.000000 +   1.000000i
2.000000 +   2.000000i
3.000000 +   3.000000i
4.000000 +   4.000000i
5.000000 +   5.000000i
6.000000 +   6.000000i
7.000000 +   7.000000i
8.000000 +   8.000000i
out:
36.000000 +  36.000000i
-13.656854 +   5.656854i
-8.000000 +   0.000000i
-5.656854 +  -2.343146i
-4.000000 +  -4.000000i
-2.343146 +  -5.656854i
0.000000 +  -8.000000i
5.656854 + -13.656854i
normalized:
4.500000 +   4.500000i
-1.707107 +   0.707107i
-1.000000 +   0.000000i
-0.707107 +  -0.292893i
-0.500000 +  -0.500000i
-0.292893 +  -0.707107i
0.000000 +  -1.000000i
0.707107 +  -1.707107i
backward:
1.000000 +   1.000000i
2.000000 +   2.000000i
3.000000 +   3.000000i
4.000000 +   4.000000i
5.000000 +   5.000000i
6.000000 +   6.000000i
7.000000 +   7.000000i
8.000000 +   8.000000i