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無關！如果w=216和h=216然後`W * H * W * H = 2 176 782 336有符號32位整數的上限爲2 147 483 647你是否面臨溢出

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

：它來自該線...

解決方法是將分母分配到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或者在求和外部循環（sum+=localsum）之前寫入兩個循環並在內部循環上求和（localsum）！

再見，

弗朗西斯