Matlab中函数的解析傅里叶变换与FFT的比较分析

Question

Matlab中函数的解析傅里叶变换与FFT的比较分析

6

我已经根据Matlab中比较函数FFT和解析FT解法的代码进行了修改以回答这个问题。我正在尝试使用FFT并将结果与维基百科表格中的解析表达式进行比较。

我的代码如下：

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。第二个图现在看起来像：

什么是“镜像”图形背后的数学原因，我该如何去除它？

- Medulla Oblongata

傅里叶变换虽然与离散傅里叶变换（通过FFT算法实现）密切相关，但它们并不相同。因此，在某些特定条件下，您可能会得到非常接近的结果，但往往会出现明显的差异（就像您观察到的那样）。 - SleuthEye

通过使用 abs 函数，您可以消除虚部，因此它将始终为 0。我认为您需要使用 fftshift 而不是 ìfftshift。 - ViG

@Zep 这只是一个 MWE - 最终我需要对更复杂的函数进行 FFT 和 IFFT，这些函数没有解析形式。因此，我正在做类似于 FFT{ IFFT[F(w)]*IFFT[G(w)] / IFFT[J(w)] } 的事情，其中我只知道频率下的 F、G、J 向量。 - Medulla Oblongata

@Medulla Oblongata 我删除了我的评论，因为我认为 ifftshift 会移动并进行逆变换，但实际上它只是移动。 - Zep

我猜第一个实验中的高斯函数之所以有效，是因为使用了 abs 函数？它去除了相位成分，这个计算是错误的。在第二个图中，你可以看到 abs 的效果：没有虚部，只有正实部。当你在第三个图中去掉 abs 时，你会发现相位完全错误，但幅值是正确的。 - Cris Luengo

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Cris Luengo · Accepted Answer

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

事实上，您定义了以中间为原点的时间域函数。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和你现有的代码，对于高斯示例，我看到的结果是：

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

对于另一个函数，我看到了更大的差异，大约为6e-6的数量级。这是由于无法采样海维赛德函数所致。在采样函数中，您需要t=0，但H在0时没有值。我注意到实部有一个类似大小的偏移量，这是由t=0处的采样引起的。

通常，采样后的海维赛德函数在t=0处设置为0.5。如果我这样做，偏移完全被去除，并且实部最大差异降低了3个数量级（最大误差发生在非常接近0的值，其中我看到了一种锯齿状模式）。对于虚部，最大误差降至3e-6，仍然相当大，并且在高频率时最大。我将这些误差归因于理想和采样后的海维赛德函数之间的差异。

你应该限制自己使用带限函数（或几乎带限的函数，如高斯函数）。你可以尝试用误差函数（高斯函数的积分）替换Heaviside函数，并使用小的sigma值（sigma = 0.8 * fs是适当采样的最小sigma值）。它的傅里叶变换已知。