计算科学 - Numpy FFT 给了我一个比它应该短的脉冲。不知道我做错了什么 - 吾爱随笔录

Numpy FFT 给了我一个比它应该短的脉冲。不知道我做错了什么

计算科学麻木的傅立叶分析傅里叶变换

2021-12-18 18:19:10

我创建了一个代码（Python，numpy），它在频域中定义了一个超短激光脉冲（脉冲持续时间应为 4 fs），但是当我使用 DFT 执行傅里叶变换时，我在时域中的脉冲实际上比它应该是。

这是我的代码：

import numpy as np
import matplotlib.pyplot as plt

#pulse duration [fs] FWHM -- what I should get after the FT
T0  = 4 

#speed of light [nm/fs]
c = 299792458*10**(-6) 
#central wavelength [nm]                               
wl0 = 800       
#central frequency (angular)[rad/fs]                     
w0 = (2*np.pi*c)/wl0    
#bandwidth [rad/fs] --> from FTL
delta_w = 2*np.pi*0.441/T0         
#bandwidth in [nm]          
delta_wl = delta_w*wl0**2/(2*np.pi*c)        

#angular frequencies [rad/fs]
w = np.linspace(w0-delta_w*8, w0+delta_w*8, 2**9)
#wavelengths [nm]
wl = c/w
#frequencies [PHz]
f = w/(2*np.pi)

#to make the spectrum centered around the carrier frequency
diff_w = w-w0
sigma_w = delta_w/(np.sqrt(8*np.log(2)))
spectrum_w = np.exp(-(diff_w**2)/(2*sigma_w**2))

#phase terms (not relevant here)
phi_w = 0 
def phase(phi_0,phi_1,phi_2,phi_3):
    phi_w = phi_0 + phi_1*(w-w0) + (phi_2*(w-w0)**2)/math.factorial(2) + (phi_3*(w-w0)**3)/math.factorial(3)
    return phi_w

#field in the frequency domain
E_w = np.exp(1j*phase(0,0,0,0))*spectrum_w

#FT:
n = len(E_w)
timestep = 0.01 
fa = 1.0/timestep

t_1 = np.fft.fftfreq(n,d = timestep)
t = np.fft.fftshift(t_1)

field_ft = np.fft.ifft(E_w)
plt.plot(t,field_ft)
plt.show()

new = np.fft.fftshift(field_ft)
plt.xlim(-5,5)
plt.plot(t,new,t,np.abs(new))
plt.show()

我得到的输出是一个比它应该短的脉冲。大约 4 年前，我在这里发现了一个非常相似的问题，但没有收到任何回复。我今天充满希望！

正如之前的海报所说，这应该是一个非常简单明了的代码，但由于我一直遇到的这个问题，它并没有那么简单。

感谢你的帮助！

2个回答

运行你的代码，你的脉搏看起来有点像这样：

（抱歉没有在图中添加单位，我使用的和你一样，即 t 在 fs 中，w 在 rad/fs 中）

因此，FWHM 不正确（应该是 4 fs）并且角载波/中心频率被搞砸了（应该是 ca. 2.3 / ca. 0.37 per fs）。我认为这里有两件事不太对劲：

1）你的w轴的定义

2）频率到时间的转换

根据 FFT 的定义，已经有一个 $2\pi$ 在指数中，我们得到

σ_{t} = \frac{1}{2 σ_{w}}

$\sigma_t = \frac{1}{2\sigma_w}$ 在哪里

σ_{t}

$\sigma_t$ 和

σ_{w}

$\sigma_w$ 分别是 t 空间和 w 空间中强度的标准偏差。

现在，根据维基百科和 $\Delta t$ 和 $\Delta w$ 分别作为 t 和 w 空间中的 FWHM

fac \cdot σ_{t} = Δ t

$\text{fac} \cdot \sigma_t = \Delta t$ 和分别

fac \cdot σ_{w} = Δ w

$\text{fac} \cdot \sigma_w = \Delta w$ 和

fac = \sqrt{8 \cdot log 2}

$\text{fac} = \sqrt{8\cdot \text{log}2}$ . 将所有这些都插入你会得到

Δ w = {fac}^{2} \cdot \frac{1}{2 Δ t} \approx 2 π \cdot 0.441 \cdot \frac{1}{Δ t}

$\Delta w = \text{fac}^2 \cdot \frac{1}{2\Delta t} \approx 2\pi \cdot 0.441 \cdot \frac{1}{\Delta t}$ 正如你所做的那样。然而，在激光物理学中，人们使用场强的半高宽，而不是幅度的半高宽。这就是为什么我们使用

σ_{t} σ_{w} = \frac{1}{2}

$\sigma_t \sigma_w = \frac{1}{2}$ . 这

\frac{1}{2}

$\frac{1}{2}$ 源于强度是幅度的平方这一事实。这是一个小而重要的细节，因为我们必须注意将因子和放在适当的位置。

2

$2$

\sqrt{2}

$\sqrt{2}$

现在对于 1)，FFT 实现将考虑从 w=0 开始的频谱。你定义 w 轴的方式，算法在 E_w 中看到的，是一个以 w 轴为中心的高斯，因为你的轴对称地分布在你的峰值周围。这意味着，在 FFT 之后，您的中心频率明显高于应有的频率。这可以通过让 w 从零开始并达到足够高的值来解决（您选择了 16 个带宽，所以我也坚持使用它）：

#angular frequencies [rad/fs]
w = np.linspace(0, delta_w*16, 2**9)

如上所述，对于高斯，您需要 FWHM 中的。由于整个强度-幅度问题，不要忘记一个因素 $\sigma_w$ $\sqrt{2}$

sigma_w = delta_w/fac # standard deviation of intensity
sigma_w_field = np.sqrt(2)*sigma_w # standard deviation of amplitude

#to make the spectrum centered around the carrier frequency
diff_w = w-w0

我删除了相位函数，因为它在这里似乎没有做任何事情，并且还确保使用 sigma_w_field：

spectrum_w = np.exp(-(diff_w**2)/(2*sigma_w_field**2)) 

#field in the frequency domain
E_w = spectrum_w
plt.figure()
plt.xlabel("w")
plt.ylabel("E_w")
plt.plot(w, E_w)
plt.grid()
plt.show()

从技术上讲，上面的图不是 100% 正确的，因为高斯的左肩可能会重新出现在频谱的高端。这在这里不是问题，因为在 w=0 时，高斯已经衰减得足够多，但是如果你选择更宽的带宽或更低的中心频率，你应该以某种方式处理它（我不知道如何优雅地做到这一点）。这修复了第 1 项）。

对于第 2 项），让我们看一下变量timestep。timestep应该是采样频率的倒数，但是，我认为 0.01 不是正确的值。采样频率的倒数，即采样时间，是 f 空间中的信号长度除以样本数。这是超过的 16 个带宽（这里称为）。然而，这是在 w 空间中并且需要 f 空间，因此需要除以。再次计算强度，所以我们需要另一个将其转换为光谱的带宽： $2^9$ w_sfftfreq $2\pi$ delta_w $\sqrt{2}$

#FT:
n = len(E_w)
w_s =  np.sqrt(2)*16*delta_w/2**9
timestep = w_s/(2*np.pi)
#fa = 1.0/timestep

t_1 = np.fft.fftfreq(n,d = timestep)
t = np.fft.fftshift(t_1)

field_ft = np.fft.ifft(E_w)

new = np.fft.fftshift(field_ft)
plt.figure()
plt.xlabel("t")
plt.ylabel("field")
plt.xlim(-10,10)
plt.plot(t,new, label="new")
plt.plot(t,np.abs(new), label="abs(new)")
plt.plot(t, np.ones(len(new))*0.5*np.max(np.abs(new)), label="half maximum")
plt.grid()
plt.legend()
plt.show()

可以通过扩大范围来改进载波图w。我希望这是您正在寻找的脉搏；）

我不认为 sqrt(2) 在以下情况下是必需的： w_s = np.sqrt(2)*16*delta_w/2**9 当它被移除时，可以获得强度的正确 fwhm 脉冲持续时间。

其它你可能感兴趣的问题

上一篇为什么 OpenFOAM 使用自己的数据结构和线性求解器？下一篇(FEM) 稀疏矩阵存储技术的节点重新排序