生成用C粉红噪声图像与FFTW

问题描述：

fftw_complex * Xf = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*nrows*ncolumns); 
    fftw_plan ift = fftw_plan_dft_c2r_2d(nrows,ncolumns,Xf,image,FFTW_BACKWARD|FFTW_ESTIMATE); 
    for (int rr=0; rr<nrows; rr++) { 
     for (int cc=0; cc<ncolumns; cc++) { 
      if (rr<=nrows/2) { 
       u = 1.0*rr/nrows; 
      } 
      else { 
       u = 1.0*(rr-nrows)/nrows; 
      } 
      if (cc<=ncolumns/2) { 
       v = 1.0*cc/ncolumns; 
      } 
      else { 
       v = 1.0*(cc-ncolumns)/ncolumns; 
      } 
      // 1/f power spectrum 
      w = pow(u,2)+pow(v,2); 
      if (w!=0) { 
       Sf = pow(w,-1/2); 
      } 
      else { 
       Sf = 0; 
      } 
      // random phase 
      phi = 1.0*rand()/RAND_MAX; 
      // complex spectrum 
      Xf[rr+nrows*cc][0] = sqrt(Sf) * cos(2*pi*phi); 
      Xf[rr+nrows*cc][1] = sqrt(Sf) * sin(2*pi*phi); 
     } 
    } 
    fftw_execute(ift);

以产生2D粉红噪声（1/f）在C图像当我做逆傅立叶在Matlab使用相同的变换频谱（真实（ifft2（...）），我得到了一个典型的粉红噪声图像（左下方），但FFTW返回的图像不是粉红噪声（右）：example of pink noise images。如果我试图制造棕色噪音（1/f2），我得到的东西更糟：example and brown noise images。有没有人有任何想法，为什么我不能从FFTW逆傅里叶变换得到正确的图像？我在matlab中得到的图像是正确的，所以频谱似乎没有成为问题。

答

复杂阵列Xf对于c2r转换而言太大。大小n0*n1的实阵列的R2C变换是大小n0*(n1/2+1)的一系列复杂的（见FFTW Real-data DFT Array Format。这是有道理由于DFT transform的特定属性。事实上，如果X是长度n的实阵列的conponant Xf[n-k]其DFT变换Xf中是Xf[k]复共轭。其结果是，时间和存储器可以通过删除上述复数排列的一半被保存。

通过调用fftw_plan ift = fftw_plan_dft_c2r_2d(nrows,ncolumns,Xf,image,FFTW_BACKWARD|FFTW_ESTIMATE); FFTW创建一个计划到nrows*(ncolumns/2+1)复数排列Xf变换回为nrow*ncolumn reals数组，因此，频率将被相应地计算。

以下基于您的示例代码可生成可在paraview中计算的VTK图像。它是由编译gcc main.c -o main -lfftw3 -lm -Wall

#include<stdlib.h> 
#include<complex.h> 
#include<math.h> 
#include<fftw3.h> 

#define PI 3.14159265358979323846 

int main(void){ 

    int nrows=400; 
    int ncolumns=1000; 

    double* image=malloc(nrows*ncolumns*sizeof(double)); 

    fftw_complex * Xf = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*nrows*(ncolumns/2+1)); 
    fftw_plan ift = fftw_plan_dft_c2r_2d(nrows,ncolumns,Xf,image,FFTW_BACKWARD|FFTW_ESTIMATE); 

    int rr; 
    int cc; 
    double u,v,w,phi,Sf; 

    for (rr=0; rr<nrows; rr++) { 
     for (cc=0; cc<ncolumns/2+1; cc++) { 
      if (rr<=nrows/2) { 
       u = 1.0*rr/nrows; 
      } 

      else { 
       u = 1.0*(rr-nrows)/nrows; 
      } 

      v = 1.0*cc/ncolumns; 


      // 1/f power spectrum 
      w = pow(u,2)+pow(v,2); 
      if (w!=0) { 
       // Sf = pow(w,-1./2); 
       Sf = pow(w,-1.); 
      } 
      else { 
       Sf = 0; 
      } 
      // random phase 
      phi = 1.0*rand()/RAND_MAX; 
      // complex spectrum 
      //Xf[rr*(ncolumns/2+1)+cc][0] = sqrt(Sf) * cos(2*pi*phi); 
      //Xf[rr*(ncolumns/2+1)+cc][1] = sqrt(Sf) * sin(2*pi*phi); 
      Xf[rr*(ncolumns/2+1)+cc]=sqrt(Sf)*(cos(2*PI*phi)+I*sin(2*PI*phi)); 
     } 
    } 
    fftw_execute(ift); 

    // writing to VTK file 

    FILE* file=fopen("image.vtk","w"); 
    fprintf(file,"# vtk DataFile Version 2.0\n"); 
    fprintf(file,"pinknoise\n"); 
    fprintf(file,"ASCII\n"); 
    fprintf(file,"DATASET STRUCTURED_POINTS\n"); 
    fprintf(file,"DIMENSIONS %d %d 1\n",nrows,ncolumns); 
    fprintf(file,"ASPECT_RATIO 1 1 1\n"); 
    fprintf(file,"ORIGIN 0 0 0\n"); 
    fprintf(file,"POINT_DATA %d\n",nrows*ncolumns); 
    fprintf(file,"SCALARS image double 1\n"); 
    fprintf(file,"LOOKUP_TABLE default\n"); 
    for (cc=0; cc<ncolumns; cc++) { 
     for (rr=0; rr<nrows; rr++) { 

      fprintf(file,"%lf ",image[rr*(ncolumns)+cc]); 
     } 
    } 
    fclose(file); 

    fftw_destroy_plan(ift); 
    fftw_free(Xf); 
    free(image); 

    return(0); 
}

是的，它确实有效。我误解了FFTW文档。谢谢， – baca

生成用C粉红噪声图像与FFTW

相关推荐