我想使用Python将时域到频域的模拟激光脉冲进行傅立叶变换。我从高斯函数开始,因为当用术语定义宽度时,已知“时间带宽积”(时域中的宽度乘以频域中的宽度)为0.44。高斯的半峰全宽。
但是,当使用numpy.fft.fft
时,我发现时间带宽乘积为0.88,是应有的两倍。
这是我的代码(前几行的最小示例,其余的只是作图):
import numpy as np
import matplotlib.pyplot as plt
fwhm = 40e-15 # using a 40 femtosecond pulse
t = np.linspace(-500e-15, 500e-15, 2000)
Et = np.exp( -t**2 / (2*(fwhm / 2.35482)**2) ) # gaussian function
Ef = np.abs(np.fft.fftshift( np.fft.fft(Et) )) # take the fourier transform
f = np.fft.fftshift( np.fft.fftfreq(Ef.shape[0],t[1]-t[0]) ) # generate the frequencies
fwhm_fft = 2 * np.abs( f[ np.argmin(np.abs(0.5*np.max(Ef)-Ef)) ] ) # find the fwhm of the frequnecy-domain signal
print 'Observed time-bandwidth product: %.3f'%(fwhm*fwhm_fft)
# just making plots from here onwards:
fig, axs = plt.subplots(2,1, figsize=(6,8))
axs[0].set_title('Time domain')
axs[0].plot(t,Et)
axs[0].axvline(-fwhm*0.5, color='r', alpha=0.5, label='Full-width at half-maximum (FWHM) = %.1f fs'%(fwhm*1e15))
axs[0].axvline( fwhm*0.5, color='r', alpha=0.5)
axs[0].set_ylim(0,1.3)
axs[0].set_xlabel('Time (sec)')
axs[1].set_title('Frequency domain')
axs[1].plot(f,Ef)
axs[1].axvline(-0.44/fwhm*0.5, color='r', alpha=0.5, label='FWHM should be %.1f THz'%(0.44/fwhm*1e-12) )
axs[1].axvline( 0.44/fwhm*0.5, color='r', alpha=0.5)
axs[1].axvline(-fwhm_fft*0.5, color='g', alpha=0.5, label='FWHM is actually %.1f THz'%(fwhm_fft*1e-12) )
axs[1].axvline( fwhm_fft*0.5, color='g', alpha=0.5)
axs[1].set_xlim(-5e13,5e13)
axs[1].set_ylim(0,120)
axs[1].set_xlabel('Frequency (Hz)')
for ax in axs:
ax.legend(fontsize=10)
ax.set_ylabel('Electric field intensity (arbitrary units)')
plt.tight_layout()
plt.savefig('time-bandwidth-product.png', dpi=200)
plt.show()
最佳答案
@PaulPanzer在正确的轨道上!当比较两个高斯函数的FWHM时,我们确实希望找到0.88作为时间带宽乘积。
但是,为什么大多数引用 [1,2,3]都说0.44是激光脉冲的时间带宽积?关键是我们实际观察到的是电场强度(E)的强度(I),其中I = E ^ 2。因此,实际上,比较强度分布的宽度而不是电场分布最有意义。当比较强度分布时,我们发现时间带宽积确实为0.44。
修改后的代码:
import numpy as np
import matplotlib.pyplot as plt
fwhm = 40e-15 # using a 40 femtosecond pulse
t = np.linspace(-1000e-15, 1000e-15, 4000)
It = np.exp( -t**2 / (2*(fwhm / 2.35482)**2) ) # Intensity in the time domain
Et = np.sqrt(It) # E-field in the time domain
Ef = np.abs(np.fft.fftshift( np.fft.fft(Et) )) # FT to get E-field in frequency domain
If = Ef**2 # Intensity in the frequnecy domain
f = np.fft.fftshift( np.fft.fftfreq(Ef.shape[0],t[1]-t[0]) ) # generate the frequencies
fwhm_fft = 2 * np.abs( f[ np.argmin(np.abs(0.5*np.max(If)-If)) ] ) # find the fwhm of the frequency-domain signal
print 'Observed time-bandwidth product: %.3f'%(fwhm*fwhm_fft)
# just making plots from here onwards:
fig, axs = plt.subplots(2,1, figsize=(6,8))
axs[0].set_title('Time domain')
axs[0].plot(t,It)
axs[0].axvline(-fwhm*0.5, color='r', alpha=0.5, label='Full-width at half-maximum (FWHM) = %.1f fs'%(fwhm*1e15))
axs[0].axvline( fwhm*0.5, color='r', alpha=0.5)
axs[0].set_xlim(-150e-15, 150e-15)
axs[0].set_ylim(0,1.3)
axs[0].set_xlabel('Time (sec)')
axs[1].set_title('Frequency domain')
axs[1].plot(f,If)
axs[1].axvline(-0.44/fwhm*0.5, color='r', alpha=0.5, label='FWHM should be %.1f THz'%(0.44/fwhm*1e-12) )
axs[1].axvline( 0.44/fwhm*0.5, color='r', alpha=0.5)
axs[1].axvline(-fwhm_fft*0.5, color='g', alpha=0.5, ls='dashed', label='FWHM is actually %.1f THz'%(fwhm_fft*1e-12) )
axs[1].axvline( fwhm_fft*0.5, color='g', alpha=0.5, ls='dashed')
axs[1].set_xlim(-2.0e13,2.0e13)
axs[1].set_ylim(0,30000)
axs[1].set_xlabel('Frequency (Hz)')
for ax in axs:
ax.legend(fontsize=10)
ax.set_ylabel('Electric field intensity (arbitrary units)')
plt.tight_layout()
plt.savefig('time-bandwidth-product.png', dpi=200)
plt.show()
PS:RP-Photonics是很棒的资源。它是激光和光子学领域的主要教科书之一。