FFTW 和 fft 与 MatLab

Question

我对离散 fft 有一个奇怪的问题。我知道高斯函数 exp(-x^2/2) 的傅里叶变换又是同一个高斯函数 exp(-k^2/2)。我尝试在 MatLab 和 FFTW 中使用一些简单的代码对其进行测试，但得到了奇怪的结果。

首先，结果的虚部不应该是零（在 MatLab 中）。

其次，实部的绝对值是高斯曲线，但没有绝对值的一半模态具有负系数。更准确地说，每一秒的模式都有一个系数，它应该是负数。

第三，得到的高斯曲线的峰值（在取实部的绝对值之后）不是一而是高得多。它的高度与 x 轴上的点数成正比。但是，比例因子不是 1，而是接近 1/20。

谁能解释我做错了什么？

这是我使用的 MatLab 代码：

    function [nooutput,M] = fourier_test

    Nx = 512;      % number of points in x direction

    Lx = 50;        % width of the window containing the Gauss curve

    x = linspace(-Lx/2,Lx/2,Nx);     % creating an equidistant grid on the x-axis

    input_1d = exp(-x.^2/2);                 % Gauss function as an input
    input_1d_hat = fft(input_1d);            % computing the discrete FFT
    input_1d_hat = fftshift(input_1d_hat);   % ordering the modes such that the peak is centred

    plot(real(input_1d_hat), '-')
    hold on
    plot(imag(input_1d_hat), 'r-')

Answer 1

答案基本上是 Paul R 在他的第二条评论中所建议的，您引入了一个相移（与频率线性相关），因为input_1d_hat 描述的高斯中心实际上位于k>0，其中k+1是input_1d_hat 的索引。相反，如果您将数据居中（例如 input_1d_hat(1) 对应于中心），您将在频域中获得相位校正的高斯：

    Nx = 512;      % number of points in x direction
    Lx = 50;        % width of the window containing the Gauss curve

    x = linspace(-Lx/2,Lx/2,Nx);     % creating an equidistant grid on the x-axis

    %%%%%%%%%%%%%%%%
    x=fftshift(x);   % <-- center
    %%%%%%%%%%%%%%%%

    input_1d = exp(-x.^2/2);                 % Gauss function as an input
    input_1d_hat = fft(input_1d);            % computing the discrete FFT
    input_1d_hat = fftshift(input_1d_hat);   % ordering the modes such that the peak is centered

    plot(real(input_1d_hat), '-')
    hold on
    plot(imag(input_1d_hat), 'r-')

根据 DFT 的定义，如果高斯没有居中，以至于最大值出现在 k=0，您将看到相位扭曲。关闭fftshift的效果是对数据集进行循环移位或左右交换，相当于将峰的中心移到k=0。

至于幅度缩放，这是在 Matlab 中实现的 DFT 定义的问题。来自 FFT 的文档：

For length N input vector x, the DFT is a length N vector X,
with elements
                 N
   X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
                n=1
The inverse DFT (computed by IFFT) is given by
                 N
   x(n) = (1/N) sum  X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.
                k=1

请注意，在前向步骤中，总和未由 N 归一化。因此，如果您增加总和中的点数 Nx，同时保持高斯函数的宽度 Lx 不变你会按比例增加X(k)。

至于信号泄漏到虚频维度，这是由于 DFT 的离散形式，这会导致截断和其他影响，正如 Paul R 再次指出的那样。如果在保持 Nx 的同时减少 Lx恒定，您应该会看到虚维度中的信号量相对于相对于实际维度有所减少（比较光谱，同时保持实际维度中的峰值强度相等）。

您会发现类似问题的更多答案here 和here。