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

Question

我正在尝试使用 FFTW 在 C 中生成二维粉红噪声 (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. It makes sense due to a particular property of the 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);
}

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

generating pink noise images in C with fftw

c

fftw