在Python中证明傅立叶变换操作 [英] Proving Fourier transform operation in Python
问题描述
我在时域有一个表达式
f = -1j*H(t) * exp(-(1j*a+b)*t)
可以使用已知属性进行傅立叶分析(a重组分阶梯功能). FT操作的结果是
which can be Fourier transformed analytically using known properties (H is the Heaviside step function). The result of this FT operation is
F = (w-a-1j*b)/((w-a)**2+b**2)
其中w是频率.
现在，我正在使用本文在Python中对f进行数值傅里叶变换，并确认我确实获得了相同的分析结果F:
Now I'm using the tips in this article to do numerical Fourier transform on f in Python, and confirm that I do get the same analytical result F:
import numpy as np
import matplotlib.pyplot as plt
t = np.linspace(-10,10,1e4) # time
w = np.linspace(-10,10,1e4) # frequency
b = 0.1
a = 1
H = lambda x: 1*(x>0) # heaviside function
# function in time
f = -1j*H(t)*np.exp(-(1j*a+b)*t)
# function in frequency (analytical work)
F = (w-a-1j*b)/((w-a)**2+b**2)
hann = np.hanning(len(t)) # hanning window
# function in frequency (numerical work)
F2 = 2/len(t)*np.fft.fft(hann*f)
plt.figure()
plt.plot(w,F.real,'b',label='analytical')
plt.plot(w,F2.real,'b--',label='fft')
plt.xlabel(r'$\omega$')
plt.ylabel(r'Re($F(\omega)$)')
plt.legend(loc='best')
plt.figure()
plt.plot(w,F.imag,'g',label='analytical')
plt.plot(w,F2.imag,'g--',label='fft')
plt.xlabel(r'$\omega$')
plt.ylabel(r'Im($F(\omega)$)')
plt.legend(loc='best')
plt.show()
但是Python的FFT函数似乎给了我一些完全错误的信息.当绘制F和F2时很明显.
However Python's FFT function seems to give me something completely wrong. This is evident when F and F2 are plotted.
这是情节...
在这些图中并不明显，但是如果放大w=-10和10区域，则可能是由于fft算法引起的小振荡.
It's not obvious in these figures, but if you zoom in around the w=-10 and 10 areas, there are small oscillations, possibly due to the fft algorithm.
推荐答案
FFT算法将计算DFT，该DFT的第一个样本具有原点(在空间和频域上).您需要对信号进行移位(在应用Hanning窗口之后)，以便t = 0是最左边的样本，并且在计算FFT之后，您必须进行反移位.
The FFT algorithm computes the DFT, which has the origin (both spatial and in frequency domain) on the first sample. You need to shift your signal (after applying the Hanning window) so that t=0 is the leftmost sample, and after computing the FFT you have to do the inverse shift.
MATLAB具有ifftshift和fftshift，它们实现了这两个移位. NumPy必须具有类似的功能.
MATLAB has ifftshift and fftshift, which implement those two shifts. NumPy must have similar functions.
代码的另一个问题是计算DFT，并将其绘制在计算出的w给定的位置，但与计算DFT的实际频率无关.
Another issue with your code is that you compute the DFT, and plot it at the locations given by the w that you computed, but is unrelated to the actual frequencies at which the DFT is computed.
这是您的代码，已翻译为MATLAB，并已固定为正确计算F2和w *.我希望这是有用的.需要注意的一件事是您的F与F2不匹配，我相信这不是由于F2中的错误，而是由于您的F计算中的错误.形状相似，但F的缩放比例和镜像不同.
Here is your code, translated to MATLAB, and fixed to properly compute F2 and w *. I hope this is useful. One thing to note is that your F does not match F2, I am confident that this is not due to an error in F2, but an error in your computation of F. The shapes are similar, but F is scaled differently and mirrored.
N = 1e3;
t = linspace(-100,100,N); % time
Fs = 1/(t(2)-t(1));
w = Fs * (-floor(N/2):floor((N-1)/2)) / N; % NOTE proper frequencies
b = 0.1;
a = 1;
H = @(x)1*(x>0); % Heaviside function
% function in time
f = -1j*H(t).*exp(-(1j*a+b)*t);
% function in frequency (analytical work)
F = (w-a-1j*b)./((w-a).^2+b.^2);
% hanning window
hann = 0.5*(1-cos(2*pi*linspace(0,1,N)));
% function in frequency (numerical work)
F2 = fftshift(fft(ifftshift(hann.*f))); % NOTE shifting of origin
figure
subplot(2,1,1), hold on
plot(w,real(F),'b-')
plot(w,real(F2),'r-')
xlabel('\omega')
ylabel('Re(F(\omega))')
legend({'analytical','fft'},'Location','best')
subplot(2,1,2), hold on
plot(w,imag(F),'b-')
plot(w,imag(F2),'r-')
xlabel('\omega')
ylabel('Im(F(\omega))')
legend({'analytical','fft'},'Location','best')
脚注:
*我还更改了颜色，MATLAB的绿色太浅了.
Footnote:
* I also changed the colors, MATLAB's green is too light.
这篇关于在Python中证明傅立叶变换操作的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
