我在解释/理解以下现象时遇到困难：
要测试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
由此我可以用（1）计算f的傅立叶变换。
最后，我使用后向傅立叶变换来计算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中）散点图“输出” u 1：2：4（对于计算出的解）
最佳答案
以下是一些需要修改的地方：
您需要考虑所有小的频率，包括负频率！索引i对应于频率2PI i/N，但也对应于频率2PI (i-N)/N。在傅立叶空间中，数组的结尾和开头一样重要！在我们的例子中，我们保持最小的频率：数组的前半部分是2PI i/N，下半部分是2PI（i-N）/ N。
当然，正如保罗所说，N-1应该在N => double dx = L/(N - 1); double dx = L/(N );中与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 *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;
}
这段代码远非完美，它既没有优化也不美观。但是它几乎提供了预期的结果。
再见