在 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)|');
我弄错了频域吗?
我对你的实现分析了很长时间。
您的频域实现看起来不错。为什么?
- 你有高斯 like FFT O/P 用于高斯 like 输入。 (稍后会详细介绍)
- 频谱 bin 递增,即 df = Fs/n。这表示您的频率变换的 X 轴没有任何问题。
- 你的变换的 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))
。代码是正确的。
我在获取高斯 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)|');
我弄错了频域吗?
我对你的实现分析了很长时间。 您的频域实现看起来不错。为什么?
- 你有高斯 like FFT O/P 用于高斯 like 输入。 (稍后会详细介绍)
- 频谱 bin 递增,即 df = Fs/n。这表示您的频率变换的 X 轴没有任何问题。
- 你的变换的 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))
。代码是正确的。