在 MATLAB 中对高斯函数使用 FFT 和正确重新缩放

Using FFT and rescaling properly on a Gaussian in MATLAB

我在获取高斯 exp(-2(x^2 + y^2) 的傅里叶变换的正确输出时遇到问题。我从分析上知道傅里叶变换应该是 0.25 * exp(-0.125(k ^2 + j^2)), 其中k和j是傅立叶变量,但是FFT输出的绝对值的输出与这个不匹配。这是解析解:

与我的输出相比:

这是相关代码:

Fs = 10;                                        %space frequency
Dx = 1 / Fs;                                    %sampling period
x = -10: Dx : 10;                               %space vector on domain [-10, 10], one dimensional
L = length(x);                                  %length of signal
[X, Y] = meshgrid(x, x);                        %2D domain defined by the space vector x

U =  exp(-2 *( X.^2 + Y.^2));                                      
n = 2^nextpow2(L);                              %padding for FFT to optimize performance

U_hat = fft2(U, n, n);                          %FFT of the Gaussian

Dk = Fs*((-n/2):((n/2) - 1)) / n;               %space frequency domain (ie, the fourier domain) in one dimension, shifted, rescaled by n
P = abs(fftshift(U_hat) / n);                   %power spectrum (ie | X_hat | = sqrt(X_hat * complex_conjugate(X_hat)), shifted
                                                                                                                 
[K1, K2] = meshgrid(Dk, Dk);                    %Fourier domain

mesh(K1, K2, P);
title('Gaussian Pulse in Frequency Domain');
xlabel('Frequency (k_1)');
ylabel('Frequency (k_2)');
zlabel('|P(f)|');                  

我弄错了频域吗?

我对你的实现分析了很长时间。 您的频域实现看起来不错。为什么?

  1. 你有高斯 like FFT O/P 用于高斯 like 输入。 (稍后会详细介绍)
  2. 频谱 bin 递增,即 df = Fs/n。这表示您的频率变换的 X 轴没有任何问题。
  3. 你的变换的 Y 轴是 fft 的 fftshift 的 abs 除以 n,以标准化变换。你没有乘以大多数人会采用的 2,因为你也对负频率感兴趣。

所以实现是正确的,但是最后的结果是错误的... 这可能只意味着一件事。您的高斯信号看起来只是高斯信号,但实际上不是!

高斯波应具有以下方程。

Click for Gaussian signal equation

如果你的均值为0,方差和sd为1。它就变成了

U = (1/(2pi))exp(-( X.^2 + Y.^2)/ 2) --- (1)

即使您要忽略 (1) 的第一部分,即 1/2*pi,因为它是一个简单的幅度缩放,您的等式 *exp(-2 ( X.^2 + Y.^2))还是不遵循高斯波的基本方程

请重新检查您的输入。

感谢 Cris Luengo 在评论中的回答 - 问题中的解析解不正确。 Mathematica 报告傅立叶变换为 0.25 * exp(-0.125(k^2 + j^2)),但实际变换为 (pi / 2) * exp(-( (pi^2) / 2) * (k^2 + j^2))。代码是正确的。