混淆测试fftw3 - 泊松方程2d测试 [英] Confusion testing fftw3 - poisson equation 2d test
问题描述
我无法解释/理解以下现象:
要测试fftw3我使用2d poisson测试用例:
laplacian(f ,y))= -g(x,y)与周期边界条件。在应用傅立叶变换到方程式后,我们得到:F(kx,ky)=如果我取g(x,y)= sin(x)+ sin(y),则G(kx,ky)/(kx 2 + ky 2)(1)
< (x,y)\in [0,2 \pi]我立即有f(x,y)= g(x,y)
我试图用fft获取:
我用g进行正向傅里叶变换计算g
最后,我用后向傅里叶变换计算f(不用忘记归一化1 /(nx *)),我们可以计算f的傅立叶变换。(b)
(例如,幅度为更糟糕的是,如果我尝试g(x,y)= sin(x)* sin(y),则N = 256是N = 512时获得的振幅的两倍)
< ,该曲线甚至没有相同的解的形式。
(注意,我必须改变方程; i在这种情况下除以2的拉普拉斯算子:(1)变为F(kx,ky)= 2 * G(kx,ky)/(kx 2 + ky 2)
这是代码:
/ *
* 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 X [i] = -pi + i * dx;
for(j = 0; 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]); //第二个测试用例
in1 [i * N + j] [1] = 0;
}
}
fftw_execute(p1); // FFT forward
for(i = 0; i for(j = 0; j < N * j + 1e-16);(3)N + j = 1。
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) //第二测试用例
// in2 [i * N + j] [1] = 2 * out1 [i * N + j] [1] /(i * i + j * j + 1e- ;
}
}
fftw_execute(p2); // FFT backward
//检查计算结果
double erl1 = 0 .;
for(i = 0; i for(j = 0; 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; //>输出
}
}
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的错。因此我的困境。
任何帮助解决这个难题都会非常感激。
注意:
执行:./ a.out> output
并且我使用gnuplot绘制曲线:
(在gnuplot中)splotoutputu 1:2:4对于计算解)
这里有一些小修改:
-
您需要考虑所有小频率,包括负频率!索引
i 对应于频率
2PI i / N
,但也对应于频率2PI )/ N
。在傅里叶空间中,数组的结束与开始一样重要!在我们的例子中,我们保持最小的频率:对于数组的前半部分,它是2PI i / N
,在后半部分是2PI(iN)/ N。 p> -
当然,如Paul所说,
N-1
应该是N
在
double dx = L /(N-1);
=>c $ c>
N-1
不对应于连续的周期性信号。 -
缩放...我是凭经验做的
我获得的结果更接近预期结果,对于这两种情况。这是代码:
/ *
* 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 * output1,* in2,* out2,* in1;
fftw_plan p1,p2;
double L = 2. * pi;
double dx = L /(N);
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 X [i] = --pi + i * dx;
for(j = 0; 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]); //第二个测试用例
in1 [i * N + j] [1] = 0;
}
}
fftw_execute(p1); // FFT forward
for(i = 0; i for(j = 0; j < N; j ++){
double fact = 0;
in2 [i * N + j] [0] = 0;
in2 [i * N + j] [1] = 0;
if(2 * i fact =((double)i * i);
} else {
fact =((double)(N-i)*(N-i));
}
if(2 * j fact + =((double)j * j);
} else {
fact + =((double)(N-j)*(N-j));
}
if(fact!= 0){
in2 [i * N + j] [0] = out1 [i * N + j] [0]
in2 [i * N + j] [1] = out1 [i * N + j] [1]
} else {
in2 [i * N + j] [0] = 0;
in2 [i * N + j] [1] = 0;
}
// in2 [i * N + j] [0] = out1 [i * N + j] [0]
// in2 [i * N + j] [1] = out1 [i * N + j] [1]
// in2 [i * N + j] [0] = out1 [i * N + j] [0] *(1.0 /(i * i + 1e-16)+ 1.0 / 1e-16)+1.0 /((Ni)*(Ni)+ 1e-16)+1.0 /((Nj)*(Nj)+ 1e-16)
// in2 [i * N + j] [1] = out1 [i * N + j] [1] *(1.0 /(i * i + 1e-16)+ 1.0 / 1e-16)+1.0 /((Ni)*(Ni)+ 1e-16)+1.0 /((Nj)*(Nj)+ 1e-16)
// in2 [i * N + j] [0] = 2 * out1 [i * N + j] [0] /(i * i + j * j + 1e-16) //第二测试用例
// in2 [i * N + j] [1] = 2 * out1 [i * N + j] [1] /(i * i + j * j + 1e- ;
}
}
fftw_execute(p2); // FFT backward
//检查计算结果
double erl1 = 0 .;
for(i = 0; i for(j = 0; 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; //>输出
// cout<< i<<<< j<< sin(X [i])* sin(Y [j])<< out2 [i * N + j] [0] /(N * N)< endl; //>输出
}
}
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;
}
这段代码远非完美,它既不优化也不美观。但它几乎给出了预期的结果。
再见,
I am having trouble explaining/understanding the following phenomenon: To test fftw3 i am using the 2d poisson test case:
laplacian(f(x,y)) = - g(x,y) with periodic boundary conditions.
After applying the fourier transform to the equation we obtain : F(kx,ky) = G(kx,ky) /(kx² + ky²) (1)
if i take g(x,y) = sin (x) + sin(y) , (x,y) \in [0,2 \pi] i have immediately f(x,y) = g(x,y)
which is what i am trying to obtain with the fft :
i compute G from g with a forward Fourier transform
From this i can compute the Fourier transform of f with (1).
Finally, i compute f with the backward Fourier transform (without forgetting to normalize by 1/(nx*ny)).
In practice, the results are pretty bad?
(For instance, the amplitude for N = 256 is twice the amplitude obtained with N = 512)
Even worse, if i try g(x,y) = sin(x)*sin(y) , the curve has not even the same form of the solution.
(note that i must change the equation; i divide by two the laplacian in this case : (1) becomes F(kx,ky) = 2*G(kx,ky)/(kx²+ky²)
Here is the code:
/*
* 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;
}
I can't find any (more) mistakes in my code (i installed the fftw3 library last week) and i don't see a problem with the maths either but i don't think it's the fft's fault. Hence my predicament. I am all out of ideas and all out of google as well.
Any help solving this puzzle would be greatly appreciated.
note :
compiling : g++ test.cpp -lfftw3 -lm
executing : ./a.out > output
and i use gnuplot in order to plot the curves : (in gnuplot ) splot "output" u 1:2:4 ( for the computed solution )
Here are some little points to be modified :
You need to account for all small frequencies, including the negative ones ! Index
i
corresponds to the frequency2PI i/N
but also to the frequency2PI (i-N)/N
. In the Fourier space, the end of the array matters as much as the beginning ! In our case, we keep the smallest frequency : it's2PI i/N
for the first half of the array, and 2PI(i-N)/N on the second half.Of course, as Paul said,
N-1
should beN
indouble dx = L/(N - 1);
=>double dx = L/(N );
N-1
does not correspond to a continious periodic signal. It woud be hard to use it as a test case...Scaling...I did it empirically
The result i obtain is closer to the expected one, for both cases. Here is the code :
/*
* 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 );
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++){
double fact=0;
in2[i*N + j][0]=0;
in2[i*N + j][1]=0;
if(2*i<N){
fact=((double)i*i);
}else{
fact=((double)(N-i)*(N-i));
}
if(2*j<N){
fact+=((double)j*j);
}else{
fact+=((double)(N-j)*(N-j));
}
if(fact!=0){
in2[i*N + j][0] = out1[i*N + j][0]/fact;
in2[i*N + j][1] = out1[i*N + j][1]/fact;
}else{
in2[i*N + j][0] = 0;
in2[i*N + j][1] = 0;
}
//in2[i*N + j][0] = out1[i*N + j][0];
//in2[i*N + j][1] = out1[i*N + j][1];
// in2[i*N + j][0] = out1[i*N + j][0]*(1.0/(i*i+1e-16)+1.0/(j*j+1e-16)+1.0/((N-i)*(N-i)+1e-16)+1.0/((N-j)*(N-j)+1e-16))*N*N;
// in2[i*N + j][1] = out1[i*N + j][1]*(1.0/(i*i+1e-16)+1.0/(j*j+1e-16)+1.0/((N-i)*(N-i)+1e-16)+1.0/((N-j)*(N-j)+1e-16))*N*N;
//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<< 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;
}
This code is far from being perfect, it is neither optimized nor beautiful. But it gives almost what is expected.
Bye,
这篇关于混淆测试fftw3 - 泊松方程2d测试的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!