FFT：FFTW Matlab FFT2之谜答案

【问题标题】：FFT : FFTW Matlab FFT2 mysteryFFT：FFTW Matlab FFT2之谜
【发布时间】：2015-06-25 13:20:21
【问题描述】：

我继承了带有 fft 子例程的旧 fortran 代码，但我无法追踪该程序的来源。我唯一知道的是调用 ff2prp() 和调用 fft2() 来执行 2D 前向和后向 DFT。为了知道代码在做什么，我对 4x4 2D 数组（矩阵）进行了 DFT，结果与 Matlab 和 FFTW 结果大不相同。

问题：有人可以通过查看输出来判断代码在做什么吗？输入输出都是实数数组

输入数组

 0.20000     0.30000     1.00000     1.20000
 0.00000    12.00000     5.00000     1.30000
 0.30000     0.30000     1.00000     1.40000
 0.00000     0.00000     0.00000     1.50000

使用 fft2() fortran 例程进行前向 FFT 后

 0.16875    -0.01875    -0.05000     0.05625
 0.00000    12.00000     5.00000     1.30000
 0.30000     0.30000     1.00000     1.40000
 0.00000     0.00000     0.00000     1.50000

执行 DCT 的 Matlab 输出：dct2(input)

    6.3750   -0.8429   -3.4250   -2.4922
    2.4620    0.6181   -2.6356   -0.9887
   -4.2750   -0.9798    4.2250    2.2730
   -4.8352   -1.2387    5.0695    3.4819

使用 FFTW 库从 C++ 代码输出。来自 FFTW 的 DCT

(6.3750, 0.00)  (-0.8429, 0.00) (-3.4250, 0.00) (-2.4922, 0.00) 
(2.4620, 0.00)  (0.6181, 0.00)  (-2.6356, 0.00) (-0.9887, 0.00) 
(-4.2750, 0.00) (-0.9798, 0.00) (4.2250, 0.00)  (2.2730, 0.00)  
(-4.8352, 0.00) (-1.2387, 0.00) (5.0695, 0.00)  (3.4819, 0.00)

使用 Matlab 进行正向 FFT - fft2(input)

  25.5000 + 0.0000i  -6.5000 - 7.2000i -10.5000 + 0.0000i  -6.5000 + 7.2000i
  -0.3000 -16.8000i -12.3000 + 4.8000i   0.1000 + 6.8000i  12.1000 + 5.2000i
 -14.1000 + 0.0000i   3.5000 +11.2000i   9.1000 + 0.0000i   3.5000 -11.2000i
  -0.3000 +16.8000i  12.1000 - 5.2000i   0.1000 - 6.8000i -12.3000 - 4.8000i

使用 FFTW 转发 FFT

(25.50, 0.00)   (-6.50, -7.20)  (-10.50, 0.00)  (-6.50, 7.20)   
(-0.30, -16.80) (-12.30, 4.80)  (0.10, 6.80)    (12.10, 5.20)   
(-14.10, 0.00)  (3.50, 11.20)   (9.10, 0.00)    (3.50, -11.20)  
(-0.30, 16.80)  (12.10, -5.20)  (0.10, -6.80)   (-12.30, -4.80)

您可以看到 Matlab 和 FFTW 的输出彼此一致，但与 fortran 代码的输出不一致。我想使用 FFTW，但由于 FFT，结果有所不同。我无法弄清楚 fortran 程序在做什么 FFT。任何人都可以通过查看输出来判断。

【问题讨论】：

如果没有您的代码，我们甚至无法开始诊断您的问题所在。如果没有代码，我们如何重现您的结果并找出问题所在？

标签： c++ matlab fortran fft fftw

【解决方案1】：

据我所知，fft2 似乎已经计算了第一行的 1D FFT（其他 3 个保持不变），结果是按 1/16 缩放并以 r0, r2, r1, i1 格式打包。

换句话说，输出可以在 Matlab 中使用：

input = [0.2 0.3 1 1.2;0 12 5 1.3;0.3 0.3 1 1.4;0 0 0 1.5];
N = size(input,2);
A = fft(input(1,:))/16;
B = reshape([real(A);imag(A)],1,2*N);
B(2) = B(N+1);
output = [B(1:N);A(2:size(input,1),:)];

如果您有理由相信 fft2 应该计算 2D FFT，那么您将数据传递给此例程的方式可能存在一些问题，从而导致结果不正确。此外，针对不同尺寸输入的其他测试用例（或您如何称呼ff2prp）可能会提供有关缩放因子选择的更多见解（例如，是 1/N^2 还是 1/4N，或其他）。

【讨论】：