Matlab中函数的解析傅立叶变换与FFT

Question

我已经修改了Comparing FFT of Function to Analytical FT Solution in Matlab for this question. I am trying to do FFTs and comparing the result with analytical expressions in the Wikipedia tables中的代码。

我的代码是：

a = 1.223;
fs = 1e5; %sampling frequency
dt = 1/fs;
t = 0:dt:30-dt;     %time vector
L = length(t); % no. sample points
t = t - 0.5*max(t); %center around t=0

y = ; % original function in time
Y = dt*fftshift(abs(fft(y))); %numerical soln

freq = (-L/2:L/2-1)*fs/L; %freq vector
w = 2*pi*freq; % angular freq

F = ; %analytical solution

figure; subplot(1,2,1); hold on
plot(w,real(Y),'.')
plot(w,real(F),'-')
xlabel('Frequency, w')
title('real')
legend('numerical','analytic')
xlim([-5,5])
subplot(1,2,2); hold on;
plot(w,imag(Y),'.')
plot(w,imag(F),'-')
xlabel('Frequency, w')
title('imag')
legend('numerical','analytic')
xlim([-5,5])

如果我研究高斯函数，让

y = exp(-a*t.^2); % original function in time

F = exp(-w.^2/(4*a))*sqrt(pi/a); %analytical solution

在上面的代码中，绘制函数的实部和虚部时看起来很吻合：

但是如果我研究一个衰减指数乘以一个 Heaviside 函数：

H = @(x)1*(x>0); % Heaviside function
y = exp(-a*t).*H(t);

F = 1./(a+1j*w); %analytical solution

然后

为什么会出现差异？我怀疑它与 Y = 行有关，但我不确定原因或方式。

编辑： 我在 Y = dt*fftshift(abs(fft(y))); 中将 ifftshift 更改为 fftshift。然后我也删除了abs。第二张图现在看起来像：

'mirrored' 图表背后的数学原因是什么？我该如何删除它？

Answer 1

问题底部的图没有镜像。如果您使用线而不是点绘制它们，您会看到数字结果具有非常高的频率。绝对分量匹配，但相位不匹配。当发生这种情况时，几乎可以肯定是时域偏移的情况。

确实，您定义了原点在中间的时域函数。 FFT 期望原点位于第一个（最左边）样本。这就是 ifftshift 的用途：

Y = dt*fftshift(fft(ifftshift(y)));

ifftshift将原点移动到第一个样本，为fft调用做准备，fftshift将结果的原点移动到中间，以便显示。

编辑

您的 t 没有 0:

>> t(L/2+(-1:2))
ans =
  -1.5000e-05  -5.0000e-06   5.0000e-06   1.5000e-05

t(floor(L/2)+1)处的样本需要为0，也就是ifftshift移动到最左边样本的样本。（我在那里使用 floor 以防 L 的大小是奇数，这里不是这种情况。）

要生成正确的 t，请执行以下操作：

fs = 1e5; % sampling frequency
L = 30 * fs;
t = -floor(L/2):floor((L-1)/2);
t = t / fs;

我首先生成一个长度正确的整数 t 轴，0 在正确的位置 (t(floor(L/2)+1)==0)。然后我通过除以采样频率将其转换为秒。

有了这个 t、我上面建议的 Y 以及你的代码的其余部分 as-is，我在高斯示例中看到了这个：

>> max(abs(F-Y))
ans =    4.5254e-16

对于另一个函数，我发现差异较大，顺序为 6e-6。这是由于无法对 Heaviside 函数进行采样。您需要在采样函数中使用 t=0，但 H 在 0 处没有值。我注意到实数分量具有相似幅度的偏移量，这是由 t=0 处的样本引起的。

通常，the sampled Heaviside function is set to 0.5 for t=0。如果我这样做，偏移量将被完全移除，实数分量的最大差异将减少 3 个数量级（最大误差发生在非常接近 0 的值上，我看到 zig-zag 模式）。对于虚部，最大误差减少到 3e-6，仍然很大，并且在高频时最大。我将这些错误归因于理想函数和采样函数之间的差异。

您可能应该限制自己使用 band-limited 函数（或接近 band-limited 的函数，例如 Gaussian）。您可能想尝试将 Heaviside 函数替换为具有小 sigma 的误差函数（高斯积分）（sigma = 0.8 * fs 是我考虑进行适当采样的最小 sigma）。 Its Fourier transform is known.

Matlab中函数的解析傅立叶变换与FFT

Analytical Fourier transform vs FFT of functions in Matlab

matlab

fft