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

问题描述

我已经修改了>比较FFT的代码这个问题对Matlab的FT分析解析函数有作用.我正在尝试进行FFT，并将结果与​​ Wikipedia表中的分析表达式进行比较.

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.

我的代码是:

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:

但是，如果我研究衰减指数乘以Heaviside函数的话:

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

然后

为什么会有差异?我怀疑它与Y =行有关，但是我不确定原因或方式.

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

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

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.

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

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将原点移动到第一个样本，为fft调用做准备，而fftshift将结果的原点移动到中间，以进行显示.

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.

修改

您的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的大小为奇数，而不是这里的情况.)

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

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

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;

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

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.

有了这个t，我上面建议的Y以及其他代码，我在高斯示例中看到了这一点:

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

对于其他功能，我发现差异更大，顺序为6e-6.这是由于无法采样Heaviside功能.您的采样函数中需要t = 0，但H的值不为0.我注意到实分量的偏移量相似，这是由t = 0处的样本引起的.

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.

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

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.

您可能应该将自己限制在频带受限的函数(或高频带等近频带受限的函数)中.您可能想用一个具有较小sigma(sigma = 0.8 * fs是我认为适当采样所考虑的最小sigma)的误差函数(高斯积分)替换Heaviside函数. 其傅立叶变换是已知的.

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屋！