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