У меня возникли проблемы с объяснением/пониманием следующего явления: Для тестирования fftw3 я использую двухмерный тестовый пример Пуассона:
laplacian(f(x,y)) = -g(x,y) с периодическими граничными условиями.
После применения преобразования Фурье к уравнению получаем: F(kx,ky) = G(kx,ky) /(kx² + ky²) (1)
если я возьму g(x,y) = sin(x) + sin(y) , (x,y) \in [0,2 \pi], то сразу же f(x,y) = g(x,y)
это то, что я пытаюсь получить с помощью fft:
я вычисляю G из g с помощью прямого преобразования Фурье
Отсюда я могу вычислить преобразование Фурье f с помощью (1).
Наконец, я вычисляю f с обратным преобразованием Фурье (не забывая нормализовать на 1/(nx*ny)).
На практике результаты довольно плохие?
(Например, амплитуда для N = 256 вдвое превышает амплитуду, полученную при N = 512)
Хуже того, если я попробую g(x,y) = sin(x)*sin(y) , кривая не будет иметь даже такой же формы решения.
(обратите внимание, что я должен изменить уравнение; в этом случае я делю лапласиан на два: (1) становится F (kx, ky) = 2 * G (kx, ky) / (kx² + ky²)
Вот код:
/*
* fftw test -- double precision
*/
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fftw3.h>
using namespace std;
int main()
{
int N = 128;
int i, j ;
double pi = 3.14159265359;
double *X, *Y ;
X = (double*) malloc(N*sizeof(double));
Y = (double*) malloc(N*sizeof(double));
fftw_complex *out1, *in2, *out2, *in1;
fftw_plan p1, p2;
double L = 2.*pi;
double dx = L/(N - 1);
in1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
out2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
out1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
in2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(N*N) );
p1 = fftw_plan_dft_2d(N, N, in1, out1, FFTW_FORWARD,FFTW_MEASURE );
p2 = fftw_plan_dft_2d(N, N, in2, out2, FFTW_BACKWARD,FFTW_MEASURE);
for(i = 0; i < N; i++){
X[i] = -pi + i*dx ;
for(j = 0; j < N; j++){
Y[j] = -pi + j*dx ;
in1[i*N + j][0] = sin(X[i]) + sin(Y[j]) ; // row major ordering
//in1[i*N + j][0] = sin(X[i]) * sin(Y[j]) ; // 2nd test case
in1[i*N + j][1] = 0 ;
}
}
fftw_execute(p1); // FFT forward
for ( i = 0; i < N; i++){ // f = g / ( kx² + ky² )
for( j = 0; j < N; j++){
in2[i*N + j][0] = out1[i*N + j][0]/ (i*i+j*j+1e-16);
in2[i*N + j][1] = out1[i*N + j][1]/ (i*i+j*j+1e-16);
//in2[i*N + j][0] = 2*out1[i*N + j][0]/ (i*i+j*j+1e-16); // 2nd test case
//in2[i*N + j][1] = 2*out1[i*N + j][1]/ (i*i+j*j+1e-16);
}
}
fftw_execute(p2); //FFT backward
// checking the results computed
double erl1 = 0.;
for ( i = 0; i < N; i++) {
for( j = 0; j < N; j++){
erl1 += fabs( in1[i*N + j][0] - out2[i*N + j][0]/N/N )*dx*dx;
cout<< i <<" "<< j<<" "<< sin(X[i])+sin(Y[j])<<" "<< out2[i*N+j][0]/N/N <<" "<< endl; // > output
}
}
cout<< erl1 << endl ; // L1 error
fftw_destroy_plan(p1);
fftw_destroy_plan(p2);
fftw_free(out1);
fftw_free(out2);
fftw_free(in1);
fftw_free(in2);
return 0;
}
Я не могу найти (больше) ошибок в своем коде (на прошлой неделе я установил библиотеку fftw3), и я не вижу проблем с математикой, но я не думаю, что это вина fft. Отсюда мое затруднительное положение. У меня закончились идеи и все из Google.
Любая помощь в решении этой головоломки будет принята с благодарностью.
Заметка :
компиляция: g++ test.cpp -lfftw3 -lm
выполнение: ./a.out > вывод
и я использую gnuplot для построения кривых: (в gnuplot) splot «выход» u 1: 2: 4 (для вычисленного решения)
N - 1термины на самом деле должны бытьN? - person Paul R schedule 02.06.2014