FFTW：正向 fft 的逆函数不等于原始函数答案

【问题标题】：FFTW: Inverse of forward fft not equal to original functionFFTW：正向 fft 的逆函数不等于原始函数
【发布时间】：2011-04-12 21:02:25
【问题描述】：

我正在尝试使用 FFTW 来计算快速求和，但遇到了一个问题：

int numFreq3 = numFreq*numFreq*numFreq;

FFTW_complex* dummy_sq_fft = (FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*numFreq3 );

FFTW_complex* dummy_sq = (FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*numFreq3 );

FFTW_complex* orig = (FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*numFreq3 );

FFTW_plan dummyPlan = FFTW_plan_dft_3d( numFreq, numFreq, numFreq,
                  orig, dummy_sq_fft,
                  FFTW_FORWARD, FFTW_MEASURE );

FFTW_plan dummyInvPlan = FFTW_plan_dft_3d( numFreq, numFreq, numFreq,
                      dummy_sq_fft, dummy_sq,
                      FFTW_BACKWARD, FFTW_MEASURE );

for(int i= 0; i < numFreq3; i++) {
  orig[ i ][ 0 ] = sparseProfile02[ 0 ][ i ][ 0 ];
  //img. part == 0
  orig[ i ][ 1 ]  = sparseProfile02[ 0 ][ i ] [ 1 ];
}

FFTW_execute(dummyPlan);
FFTW_execute(dummyInvPlan);

int count = 0; 
for(int i=0; i<numFreq3; i++) {
  double factor = dummy_sq[ i ][ 0 ]/sparseProfile02[ 0 ][ i ][ 0 ];

  if(factor < 0) {
    count++;
  }
}

std::cout<<"Count "<<count<<"\n";

FFTW_free(dummy_sq_fft);
FFTW_free(dummy_sq);
FFTW_destroy_plan(dummyPlan);
FFTW_destroy_plan(dummyInvPlan);

（这里 sparseProfile02[0] 是 FFTW_complex* 类型，只包含正实数据。）

既然我们有 dummy_sq = IFFT(FFT(sparseProfile02[ 0 ]))，我们必须有 dummy_sq = n^3*sparseProfile02。但这仅在某些时候是正确的；事实上，只要 sparseProfile02 网格上的对应值为零（反之亦然），dummy_sq 网格上的值就是负数。有谁知道为什么会这样？

【问题讨论】：

请注意，引号格式（使用“>”）和代码格式（使用四个空格缩进）之间存在差异。另见stackoverflow.com/editing-help
fftw3 inverse transform not work的可能重复

标签： c fftw

【解决方案1】：

冒着涉足死灵术的风险，您应该在 fftw 文档（此处）中指出，它明确指出 fftw 不会归一化，因此先前转换的信号的逆变换结果将是原始信号按比例缩放'n' 或信号的长度。

可能是问题。

【讨论】：

【解决方案2】：

FFT（正向和反向）存在舍入误差，我认为这就是困扰您的原因。一般来说，您不应该期望零在您的过程中保持完全为零（尽管对于琐碎的测试用例它可能为零）。在您的测试循环中，是

fabs(dummy_sq[i][0] - numFreq*numFreq*numFreq*sparseProfile02[0][i][0])

相对于您的数据量大？

作为一个非常简单（病态）的示例，只有大小为 2 的 1D FFT 和实际值：

ifft(fft([1e20, 1.0])) != [2e20, 2.0]

1.0 在双精度 FFT 的 1e20 中丢失。

当您在 sparseProfile02 中除以零样本时，您也可能会得到一些 NaN。

【讨论】：