FFTW:前向fft的反向不等于原始函数

我正在尝试使用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网格上的值就是负数。 有谁知道为什么会这样?

冒着涉及死灵法的风险,你应该在fftw docs(这里)中注意到它明确指出fftw没有归一化,因此先前变换的信号的逆变换的结果将是由’n’缩放的原始信号。或信号的长度。

可能是问题。

FFT(正向和反向)具有舍入误差,我认为这就是咬你的原因。 通常,您不应期望零在整个过程中保持为零(尽管对于琐碎的测试用例可能为零)。 在你的测试循环中,是

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

相对于您的数据量大?

作为一个非常简单(病态)的例子,只有1D FFT,大小为2,具有实际值:

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

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

当你在sparseProfile02中除以零样本时,你也可能得到一些NaNs因子。