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. Но это верно только в некоторые моменты времени; на самом деле значения в сетке dummy_sq отрицательны, когда соответствующие значения в сетке sparseProfile02 равны нулю (но не наоборот). Кто-нибудь знает, почему это происходит?

С риском погружения в некромантию вы должны отметить в документах fftw (здесь), что он явно заявляет, что fftw не нормализуется, поэтому результатом обратного преобразования ранее преобразованного сигнала будет исходный сигнал, масштабируемый ‘n’ или длину сигнала.

Может быть, проблема.

БПФ (прямой и обратный) имеют ошибку округления, и я думаю, что это то, что вас кусает. В общем, вы не должны ожидать, что ноль останется ровно нулевым через ваш процесс (хотя он может быть равен нулю для тривиальных тестовых случаев). В тестовом цикле

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

большой по отношению к величине ваших данных?

Как очень простой (патологический) пример с 1D БПФ размером 2 с реальными значениями:

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

1.0 теряется в 1e20 с двойной точностью FFT.

Также возможно, что вы получаете некоторые NaN для коэффициента, когда делитесь нулевым образцом в sparseProfile02.