fftw 输出是否取决于输入的大小？

Question

上周我一直在用 FFTW 编程一些二维卷积，通过将两个信号传递到频域，相乘，然后返回。

令人惊讶的是，只有当输入大小小于固定数字时，我才能得到正确的结果！

我发布了一些工作代码，其中我将简单的初始常数矩阵值 2 作为输入，1 作为空间域上的过滤器。这样，卷积它们的结果应该是第一个矩阵值的平均值的矩阵，即 2，因为它是常数。这是我将宽度和高度的大小分别从 0 更改为 h=215、w=215 时的输出；如果我设置 h=216、w=216 或更大，那么输出就会损坏！！我真的很感激一些关于我可能在哪里犯错误的线索。非常感谢！

#include <fftw3.h>

int main(int argc, char* argv[]) {

int h=215, w=215;

//Input and 1 filter are declared and initialized here
float *in = (float*) fftwf_malloc(sizeof(float)*w*h);
float *identity = (float*) fftwf_malloc(sizeof(float)*w*h);
for(int i=0;i<w*h;i++){
        in[i]=5;
        identity[i]=1;
    }

//Declare two forward plans and one backward    
fftwf_plan plan1, plan2, plan3;

//Allocate for complex output of both transforms
fftwf_complex *inTrans = (fftwf_complex*) fftw_malloc(sizeof(fftwf_complex)*h*(w/2+1));
fftwf_complex *identityTrans = (fftwf_complex*) fftw_malloc(sizeof(fftwf_complex)*h*(w/2+1));

//Initialize forward plans
plan1 = fftwf_plan_dft_r2c_2d(h, w, in, inTrans, FFTW_ESTIMATE);
plan2 = fftwf_plan_dft_r2c_2d(h, w, identity, identityTrans, FFTW_ESTIMATE);

//Execute them
fftwf_execute(plan1);
fftwf_execute(plan2);

//Multiply in frequency domain. Theoretically, no need to multiply imaginary parts; since signals are real and symmetric
//their transform are also real, identityTrans[i][i] = 0, but i leave here this for more generic implementation.

for(int i=0; i<(w/2+1)*h; i++){
    inTrans[i][0] = inTrans[i][0]*identityTrans[i][0] - inTrans[i][1]*identityTrans[i][1];
    inTrans[i][1] = inTrans[i][0]*identityTrans[i][1] + inTrans[i][1]*identityTrans[i][0];
}
//Execute inverse transform, store result in identity, where identity filter lied.
plan3 = fftwf_plan_dft_c2r_2d(h, w, inTrans, identity, FFTW_ESTIMATE);
fftwf_execute(plan3);

//Output first results of convolution(in, identity) to see if they are the average of in.
for(int i=0;i<h/h+4;i++){
    for(int j=0;j<w/w+4;j++){
        std::cout<<"After convolution, component (" <<  i  <<","<< j << ") is " << identity[j+i*w]/(w*h*w*h) << endl;
    }
}std::cout<<endl;

//Compute average of data
float sum=0.0;
for(int i=0; i<w*h;i++)
    sum+=in[i];

std::cout<<"Mean of input was " <<  (float)sum/(w*h)  << endl;
std::cout<< endl;

fftwf_destroy_plan(plan1);
fftwf_destroy_plan(plan2);
fftwf_destroy_plan(plan3);


return 0;
}

Answer 1

您的问题与 fftw 无关！它来自这一行：

std::cout<<"After convolution, component (" <<  i  <<","<< j << ") is " << identity[j+i*w]/(w*h*w*h) << endl;

如果 w=216 和 h=216 则 `w*h*w*h=2 176 782 336。有符号 32 位整数的上限为 2 147 483 647。您正面临溢出...

解决方案是将分母转换为float。

std::cout<<"After convolution, component (" <<  i  <<","<< j << ") is " << identity[j+i*w]/(((float)w)*h*w*h) << endl;

你要面对的下一个麻烦就是这个：

 float sum=0.0;
 for(int i=0; i<w*h;i++)
      sum+=in[i];

请记住，float 有 7 个有用的十进制数字。如果w=h=4000，计算的平均值将低于实际平均值。使用 double 或编写两个循环并在内循环 (localsum) 上求和，然后再对外循环 (sum+=localsum) 求和！

再见，

弗朗西斯