【问题标题】:Changing the center of a gaussian using FFTW in C++在 C++ 中使用 FFTW 更改高斯中心
【发布时间】:2021-12-17 06:29:30
【问题描述】:

对于一个项目,我需要使用以原点为中心的高斯傅里叶变换在真实空间中传播高斯,使用

What I want to calculate

这是乳胶代码,因为我还不能包含图片

N(x | \mu, \sigma) = F^{-1}{F{ N(x |0, \sigma)} e^{-i\ 2\pi \mu\omega} \right} ,

其中 \omega 是傅里叶空间中的频率。

现在我遇到的问题是,在使用 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 为中心,但这显然没有发生。

这是乳胶代码,因为我还不能包含图片

(这里是乳胶代码,因为我还不能包含图片)

【问题讨论】:

    标签: c++ fftw


    【解决方案1】:

    一般来说,错误越明显,找到它所需的时间就越长。

    此处已验证。

    错误就在这里:

    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 的工作方式有关。

    【讨论】:

    • 谢谢!我不应该在凌晨 2 点发布这个。当我发现错误时,我来到这里删除问题,但出于对您的尊重,我将把它留在这里。
    • @TerryHintz 确实勾选了这个答案,所以这个问题可以标记为已解决:)
    猜你喜欢
    • 2014-04-02
    • 1970-01-01
    • 2012-06-06
    • 2021-05-22
    • 2011-12-06
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多