我正在尝试使用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因子。
需要了解更多c/c++开发分享FFTW:前向fft的反向不等于原始函数,也可以关注C/ C++技术分享栏目—计算机技术网(www.ctvol.com)!
以上就是c/c++开发分享FFTW:前向fft的反向不等于原始函数相关内容,想了解更多C/C++开发(异常处理)及C/C++游戏开发关注计算机技术网(www.ctvol.com)!)。
本文来自网络收集,不代表计算机技术网立场,如涉及侵权请联系管理员删除。
ctvol管理联系方式QQ:251552304
本文章地址:https://www.ctvol.com/c-cdevelopment/979499.html