使用 numpy.fft 的时间带宽积（高斯宽度）
Time-bandwidth product (width of Gaussian) using numpy.fft
我想使用 Python 将模拟激光脉冲从时域傅里叶变换到频域。我从高斯函数开始，因为已知 "time-bandwidth product"（时域宽度乘以频域宽度）是 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()
编辑 2：似乎魔鬼在物理中，而不是数学中，请参阅 Dan 的 self-answer。高斯平方确实将 half-maxima 的位置移动了 1 / sqrt(2) 每个，所以一切都很好。我仍然要谦虚而毫无保留地向 RP Photonics 道歉。编辑 2 结束。
我很确定您正在寻找的 "bug" 是 0.44，您 link 的参考看起来不是 100% 可靠。
所以让我们自己计算一下会发生什么。傅里叶变换有不同的定义； this one 似乎是 numpy 所坚持的。在这个约定中，高斯和它的傅立叶变换的标准差的乘积是 1 / (2pi)。具有 SD sigma 的 zero-mean Gaussian 的半最大值位于 +/- sigma sqrt(2 log 2)。因此 FWHM 的乘积是 1 / (2pi) 8 log 2 = 4 / pi log 2 = 0.8825...
换句话说：您观察到的是正确的。
编辑：为了公平对待 RP Photonics，他们不一定是错误的，他们可能只是使用了傅立叶变换的另一种定义。
@PaulPanzer 走在正确的轨道上！在比较两个高斯函数的 FWHM 时，我们确实希望找到 0.88 作为 time-bandwidth 乘积。
但为什么大多数参考文献 [1,2,3] 都说 0.44 是 time-bandwidth 激光脉冲的乘积？关键是我们实际观察到的是电场（E）的强度（I），其中I=E^2。因此，实际上，比较 强度 剖面的宽度最有意义，而不是 电场 剖面的宽度。当我们比较强度分布时，我们发现 time-bandwidth 乘积确实是 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 是一个很棒的资源。是激光和光子学领域的主要教材之一。
我想使用 Python 将模拟激光脉冲从时域傅里叶变换到频域。我从高斯函数开始，因为已知 "time-bandwidth product"（时域宽度乘以频域宽度）是 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()
编辑 2：似乎魔鬼在物理中，而不是数学中，请参阅 Dan 的 self-answer。高斯平方确实将 half-maxima 的位置移动了 1 / sqrt(2) 每个，所以一切都很好。我仍然要谦虚而毫无保留地向 RP Photonics 道歉。编辑 2 结束。
我很确定您正在寻找的 "bug" 是 0.44，您 link 的参考看起来不是 100% 可靠。
所以让我们自己计算一下会发生什么。傅里叶变换有不同的定义； this one 似乎是 numpy 所坚持的。在这个约定中，高斯和它的傅立叶变换的标准差的乘积是 1 / (2pi)。具有 SD sigma 的 zero-mean Gaussian 的半最大值位于 +/- sigma sqrt(2 log 2)。因此 FWHM 的乘积是 1 / (2pi) 8 log 2 = 4 / pi log 2 = 0.8825...
换句话说：您观察到的是正确的。
编辑：为了公平对待 RP Photonics，他们不一定是错误的，他们可能只是使用了傅立叶变换的另一种定义。
@PaulPanzer 走在正确的轨道上！在比较两个高斯函数的 FWHM 时，我们确实希望找到 0.88 作为 time-bandwidth 乘积。
但为什么大多数参考文献 [1,2,3] 都说 0.44 是 time-bandwidth 激光脉冲的乘积？关键是我们实际观察到的是电场（E）的强度（I），其中I=E^2。因此，实际上，比较 强度 剖面的宽度最有意义，而不是 电场 剖面的宽度。当我们比较强度分布时，我们发现 time-bandwidth 乘积确实是 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 是一个很棒的资源。是激光和光子学领域的主要教材之一。