计算科学 - 如何使用 FFTW 计算 2D r2c 变换后的幅度？ - 吾爱随笔录

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

int main(){
    int N=2000;
    int i,j;
    FILE *filepointer;
    filepointer=fopen("2DDFT_spacetime.plt","w");
    //double in[N][N];
    double *in;
    fftw_complex *out;
    fftw_plan p;
    double fx=13.0;
    double fz=9.0;
    double x[N];
    double xstart=0.0;
    double xend=5.0/fx;
    double z[N];
    double zstart=0.0;
    double zend=5.0/fz;
    double dx=(xend-xstart)/(N-1);
    double dz=(zend-zstart)/(N-1);
    x[0]=xstart;
    z[0]=zstart;
    in = (double*) malloc(sizeof(double) * N * N); //allocates input array
    out = fftw_alloc_complex(N*((int)floor(N/2)+1)); //wrapper function ;allocates output array
    p = fftw_plan_dft_r2c_2d(N,N, in, out, FFTW_MEASURE);
    for(i=1;i<N;i++) {
        x[i]=x[i-1]+dx;
    }
    for(i=1;i<N;i++) {
        z[i]=z[i-1]+dz;
    }
    for(i=0;i<N;i++) {
        for(j=0;j<N;j++) {
            in[i*N+j]=cos(2*M_PI*fx*x[i]+2*M_PI*fz*z[j]);
        }
    }

    fftw_execute(p);

    fprintf(filepointer,"TITLE =\"FFTW\"\nVARIABLES=\"Wavenumber-x\"\n\"Wavenumber-z\"\n\"Amplitude\"\nZONE T=\"Amplitude\"\n I=%d, J=%d, K=1, ZONETYPE=Ordered\n DATAPACKING=POINT\n DT=(SINGLE SINGLE SINGLE)\n",N,(int)floor(N/2)+1);
    for(j=0;j<(int)floor(N/2)+1;j++) { 
        for(i=0;i<N;i++) { 
            fprintf(filepointer," %.9E %.9E %.9E\n",i/(xend-xstart),j/(zend-zstart),sqrt(pow(out[i*((int)floor(N/2)+1)+j][0],2)+pow(out[i*((int)floor(N/2)+1)+j][1],2)));
        }
    }
    fftw_destroy_plan(p);
    free(in); 
    fftw_free(out);
    fclose(filepointer);
    return(1);
}

我正在使用上述程序初始化 2000x2000 数据元素的 2D 字段（例如在 xz 平面中），并使用 FFTW 将数据从真实平面转换为波数平面，据我所知，波数平面具有 2000x1001 FFTW 文档中的元素。数据由函数给出：

i n = \cos (2 π f_{x} + 2 π f_{z})

$in=\cos(2\pi f_x+2\pi f_z)$

在哪里 $f_x$ 是 x 方向的频率， $f_z$ 是 z 方向的频率。我希望我选择的波数的最大幅度为 1.0。但是，我得到数字 1999959 作为最大幅度。我知道我必须对其进行归一化以获得最大幅度，但是如何计算我必须将最大幅度除以达到 Amplitude=1.0 的数字？如果数字是 2000x1001，这是否意味着由于某些数值错误，我无法获得与 1.0 完全相同的幅度？