本文介绍了混淆测试fftw3 - 泊松方程2d测试的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！ 问题描述 29岁程序员，3月因学历无情被辞！ 我无法解释/理解以下现象：要测试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） 我试图用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 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 } } 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 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 transformFrom 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 -lmexecuting : ./a.out > outputand 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 frequency 2PI i/N but also to the frequency 2PI (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's 2PI 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 be Nin double 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 empiricallyThe 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测试的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！ 上岸，阿里云！
