2

我对离散 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-')
4

1 回答 1

4

答案基本上是 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，您将看到相位扭曲。off 的效果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宽度不变的同时增加求和中的点数，您将成比例地增加。LxX(k)

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

您可以在此处和此处找到类似问题的其他答案。

于 2013-10-10T14:46:59.077 回答