Matlab中函数的解析傅里叶变换与FFT [英] Analytical Fourier transform vs FFT of functions in Matlab

问题描述

我已经修改了

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

H = @(x)1*(x>0);% Heaviside 函数y = exp(-a*t).*H(t);F = 1./(a+1j*w);%解析解

然后

为什么会有差异?我怀疑它与 Y = 行有关，但我不确定为什么或如何.

我在 Y = dt*fftshift(abs(fft(y))) 中将 ifftshift 更改为 fftshift;.然后我还删除了 abs.第二个图现在看起来像:

镜像"图形背后的数学原因是什么，我该如何删除它?

解决方案

问题底部的图没有镜像.如果您使用线条而不是点来绘制这些图形，您会看到数字结果的频率非常高.绝对分量匹配，但相位不匹配.发生这种情况时，几乎可以肯定是时域发生了变化.

实际上，您定义了原点在中间的时域函数.FFT 预计原点位于第一个(最左侧)样本处.这就是 ifftshift 的用途:

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

ifftshift 将原点移动到第一个样本，为 fft 调用做准备，fftshift 将结果的原点移动到第一个样本中间，用于展示.

<小时>

编辑

您的 t 没有 0:

<代码>>>t(L/2+(-1:2))答案 =-1.5000e-05 -5.0000e-06 5.0000e-06 1.5000e-05

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

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

fs = 1e5;％ 采样频率L = 30 * fs；t = -地板(L/2):地板((L-1)/2)；t = t/fs;

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

有了这个 t，我上面建议的 Y，以及你的其余代码，我看到了高斯示例:

<代码>>>最大(绝对(F-Y))答案 = 4.5254e-16

对于其他功能，我看到了更大的差异，按 6e-6 的顺序排列.这是由于无法对 Heaviside 函数进行采样.您的采样函数中需要 t=0，但 H 没有 0 值.我注意到实分量具有类似幅度的偏移，这是由 t=0 处的样本引起的.

通常，采样的 Heaviside 函数在 t=0 时设置为 0.5.如果我这样做，偏移将被完全移除，实部的最大差异减少 3 个数量级(最大误差发生在非常接近 0 的值，我看到锯齿形图案).对于虚部，最大误差减少到 3e-6，仍然很大，并且在高频时最大.我将这些错误归因于理想和采样的 Heaviside 函数之间的差异.

您可能应该限制自己使用带限函数(或接近带限的函数，例如高斯函数).您可能想尝试用带有小 sigma 的误差函数(高斯积分)替换 Heaviside 函数(sigma = 0.8 * fs 是我考虑进行适当采样的最小 sigma).其傅立叶变换是已知的.

I have adapted the code in 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.

My code is:

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])

If I study the Gaussian function and let

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

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

in the above code, looks like there is good agreement when the real and imaginary parts of the function are plotted:

But if I study a decaying exponential multiplied with a Heaviside function:

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

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

then

Why is there a discrepancy? I suspect it's related to the line Y = but I'm not sure why or how.

Edit: I changed the ifftshift to fftshift in Y = dt*fftshift(abs(fft(y)));. Then I also removed the abs. The second graph now looks like:

What is the mathematical reason behind the 'mirrored' graph and how can I remove it?

解决方案

The plots at the bottom of the question are not mirrored. If you plot those using lines instead of dots you'll see the numeric results have very high frequencies. The absolute component matches, but the phase doesn't. When this happens, it's almost certainly a case of a shift in the time domain.

And indeed, you define the time domain function with the origin in the middle. The FFT expects the origin to be at the first (leftmost) sample. This is what ifftshift is for:

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

ifftshift moves the origin to the first sample, in preparation for the fft call, and fftshift moves the origin of the result to the middle, for display.

---

Edit

Your t does not have a 0:

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

The sample at t(floor(L/2)+1) needs to be 0. That is the sample that ifftshift moves to the leftmost sample. (I use floor there in case L is odd in size, not the case here.)

To generate a correct t do as follows:

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

I first generate an integer t axis of the right length, with 0 at the correct location (t(floor(L/2)+1)==0). Then I convert that to seconds by dividing by the sampling frequency.

With this t, the Y as I suggest above, and the rest of your code as-is, I see this for the Gaussian example:

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

For the other function I see larger differences, in the order of 6e-6. This is due to the inability to sample the Heaviside function. You need t=0 in your sampled function, but H doesn't have a value at 0. I noticed that the real component has an offset of similar magnitude, which is caused by the sample at t=0.

Typically, the sampled Heaviside function is set to 0.5 for t=0. If I do that, the offset is removed completely, and max difference for the real component is reduced by 3 orders of magnitude (largest errors happen for values very close to 0, where I see a zig-zag pattern). For the imaginary component, the max error is reduced to 3e-6, still quite large, and is maximal at high frequencies. I attribute these errors to the difference between the ideal and sampled Heaviside functions.

You should probably limit yourself to band-limited functions (or nearly-band-limited ones such as the Gaussian). You might want to try to replace the Heaviside function with an error function (integral of Gaussian) with a small sigma (sigma = 0.8 * fs is the smallest sigma I would consider for proper sampling). Its Fourier transform is known.

这篇关于Matlab中函数的解析傅里叶变换与FFT的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！