使用 fftw 在 C 中生成粉红噪声图像

Question

我正在尝试使用 FFTW 在 C 中生成 2d 粉红噪声 (1/f) 图像

    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);

当我在 matlab 中使用相同的光谱 (real(ifft2(...)) 进行逆傅立叶变换时，我得到了一个典型的粉红噪声图像（左下）。但是 FFTW 返回的图像不是粉红色的噪音（右）：example of pink noise images。 如果我尝试制造棕色噪音 (1/f2)，我会得到更糟糕的结果：example and brown noise images。 有人知道为什么我没有从 FFTW 傅立叶逆变换得到正确的图像吗？我在 matlab 中得到的图像是正确的，所以光谱似乎不是问题。

Answer 1

复数数组 Xf 对于 c2r 转换来说太大了。大小为n0*n1 的实数数组的r2c 变换是大小为n0*(n1/2+1) 的复数数组（参见fftw Real-data DFT Array Format。由于DFT transform 的特定属性，它是有意义的。事实上，如果X 是一个长度为n的实数数组，其DFT变换Xf的辅音Xf[n-k]是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 实数数组。因此，将相应地计算频率。

以下基于您的示例代码生成可以使用 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);
}