用Python证明Fourier变换运算

2024-04-25 07:45:37 发布

您现在位置:Python中文网/ 问答频道 /正文

我有一个时域表达式

f = -1j*H(t) * exp(-(1j*a+b)*t)

它可以用known propertiesH是Heaviside阶跃函数)进行解析变换。这种FT操作的结果是

^{pr2}$

其中w是频率。

现在我使用this article中的技巧在Python中对f进行数值傅里叶变换,并确认我确实得到了相同的分析结果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和{}时,这一点很明显。

编辑:以下是绘图。。。

在这些图中并不明显,但是如果放大w=-1010区域,会有小的振荡,可能是由于fft算法造成的。


Tags: 函数inimportfftplotasnpfunction
1条回答
网友
1楼 · 发布于 2024-04-25 07:45:37

FFT算法计算DFT,它在第一个样本上有原点(包括空间域和频域)。你需要移动你的信号(在应用汉宁窗口之后),使t=0是最左边的样本,在计算FFT之后,你必须进行逆移位。在

MATLAB有ifftshift和{},它们实现了这两个移位。NumPy必须有类似的功能。在

代码的另一个问题是计算DFT,并将其绘制在您计算的w给定的位置,但与计算DFT的实际频率无关。在

这是您的代码,翻译成MATLAB,并修复了正确计算F2和{}*的代码。我希望这是有用的。需要注意的一点是你的FF2不匹配,我相信这不是由于F2中的错误,而是由于你计算{}的错误。形状相似,但F的缩放和镜像方式不同。在

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

output of code

脚注:
*我也改变了颜色,MATLAB的绿色太浅了。

相关问题 更多 >