在 C++ 中使用 FFTW 改变高斯的中心
Changing the center of a gaussian using FFTW in C++
对于一个项目,我需要使用
使用以原点为中心的高斯的傅里叶变换来传播真实的高斯space
What I want to calculate
这是乳胶代码,因为我还不能包含图片
N(x | \mu, \sigma) = F^{-1}{F{ N(x |0, \sigma)} e^{-i\ 2\pi \mu\omega} \对},
其中 \omega 是傅立叶中的频率 space。
现在我遇到的问题是,在用 fftw 执行 fft 后,我不知道如何计算某些 bin 的频率。这是我正在尝试做的代码。
int main(){
int N = 128; //Number of "pixels" in real space
int N_fft = N/2 + 1;
double *in, *result, *x;
fftw_complex *out;
fftw_plan fw, bw;
in = (double*) fftw_malloc(sizeof(double) * N);
x = (double*) malloc(sizeof(double) * N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N_fft);
result = (double*) fftw_malloc(sizeof(double) * N);
fw = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);
bw = fftw_plan_dft_c2r_1d(N, out, result, FFTW_ESTIMATE);
double min_x = -9.0, max_x = 9.0; //Limits in real space
for (int i=0; i<N; i++){
x[i] = min_x + 2*max_x*i / (N - 1);
in[i] = std::exp(-(x[i]) * (x[i]));
}
for (int i=0; i<N_fft; i++){
out[i][0] = 0.0;
out[i][1] = 0.0;
}
fftw_execute(fw);
double w;
fftw_complex z;
double w_norm;
for (int i=0; i<N_fft; i++){
w = -2*M_PI*i / (max_x - min_x); //WHAT I DON'T KNOW
// Calculating the product with the exponential for translating the gaussian
z[0] = out[i][0]*std::cos(w) - out[i][1]*std::sin(w);
z[1] = out[i][0]*std::sin(w) + out[i][0]*std::cos(w);
out[i][0] = z[0];
out[i][1] = z[1];
}
fftw_execute(bw);
for (int i=0; i<N; i++){
std::cout << x[i] << " " << result[i]/N << " " << std::exp(-x[i] * x[i]) << std::endl;
}
fftw_destroy_plan(fw);
fftw_free(in);
fftw_free(out);
return 0;
}
目前我尝试使用 w-nth = -2*np.pi * 1/(max_x - min_x) * n,它在 python, 但由于某些原因它在 c++
中不起作用
这是我用 c++ 获得的结果
result
这里ref是以0为中心的高斯分布,我得到的应该以1.0为中心,但显然不是这样。
这是乳胶代码,因为我还不能包含图片
(这是乳胶代码,因为我还不能包含图片)
一般来说,越是明显的错误,越需要时间去寻找。
此处已验证。
错误就在这里:
z[1] = out[i][0]*std::sin(w) + out[i][0]*std::cos(w);
应该是:
z[1] = out[i][0]*std::sin(w) + out[i][1]*std::cos(w);
此外,我不知道你为什么不使用 N__ft = N
,但我想这与 fftw 的工作方式有关。
对于一个项目,我需要使用
使用以原点为中心的高斯的傅里叶变换来传播真实的高斯spaceWhat I want to calculate
这是乳胶代码,因为我还不能包含图片
N(x | \mu, \sigma) = F^{-1}{F{ N(x |0, \sigma)} e^{-i\ 2\pi \mu\omega} \对},
其中 \omega 是傅立叶中的频率 space。
现在我遇到的问题是,在用 fftw 执行 fft 后,我不知道如何计算某些 bin 的频率。这是我正在尝试做的代码。
int main(){
int N = 128; //Number of "pixels" in real space
int N_fft = N/2 + 1;
double *in, *result, *x;
fftw_complex *out;
fftw_plan fw, bw;
in = (double*) fftw_malloc(sizeof(double) * N);
x = (double*) malloc(sizeof(double) * N);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N_fft);
result = (double*) fftw_malloc(sizeof(double) * N);
fw = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);
bw = fftw_plan_dft_c2r_1d(N, out, result, FFTW_ESTIMATE);
double min_x = -9.0, max_x = 9.0; //Limits in real space
for (int i=0; i<N; i++){
x[i] = min_x + 2*max_x*i / (N - 1);
in[i] = std::exp(-(x[i]) * (x[i]));
}
for (int i=0; i<N_fft; i++){
out[i][0] = 0.0;
out[i][1] = 0.0;
}
fftw_execute(fw);
double w;
fftw_complex z;
double w_norm;
for (int i=0; i<N_fft; i++){
w = -2*M_PI*i / (max_x - min_x); //WHAT I DON'T KNOW
// Calculating the product with the exponential for translating the gaussian
z[0] = out[i][0]*std::cos(w) - out[i][1]*std::sin(w);
z[1] = out[i][0]*std::sin(w) + out[i][0]*std::cos(w);
out[i][0] = z[0];
out[i][1] = z[1];
}
fftw_execute(bw);
for (int i=0; i<N; i++){
std::cout << x[i] << " " << result[i]/N << " " << std::exp(-x[i] * x[i]) << std::endl;
}
fftw_destroy_plan(fw);
fftw_free(in);
fftw_free(out);
return 0;
}
目前我尝试使用 w-nth = -2*np.pi * 1/(max_x - min_x) * n,它在 python, 但由于某些原因它在 c++
中不起作用这是我用 c++ 获得的结果
result
这里ref是以0为中心的高斯分布,我得到的应该以1.0为中心,但显然不是这样。
这是乳胶代码,因为我还不能包含图片
(这是乳胶代码,因为我还不能包含图片)
一般来说,越是明显的错误,越需要时间去寻找。
此处已验证。
错误就在这里:
z[1] = out[i][0]*std::sin(w) + out[i][0]*std::cos(w);
应该是:
z[1] = out[i][0]*std::sin(w) + out[i][1]*std::cos(w);
此外,我不知道你为什么不使用 N__ft = N
,但我想这与 fftw 的工作方式有关。