如何在Python中正确提取频谱相位并还原多项式相位

我在时域定义了一个脉冲信号，经傅里叶变换转到频域后，给频域信号添加了指数多项式相位$e^{i\cdot phase}$（其中phase为多项式）。但用numpy的angle函数提取相位时，得到了如图所示的密集尖峰，不确定这个结果是否正确，也不知道怎么还原原本的多项式相位。

代码如下：

import numpy as np
import matplotlib.pyplot as plt
fs = 1e-15
THz = 1e12
nm = 1e-9
c = 3e8
N = 2 ** 13
time_window = 3000 * fs
wavelength = 800 * nm
t = np.linspace(-time_window / 2,time_window / 2, N)
df = np.append(np.linspace(0, N / 2, int(N / 2)),(np.linspace(-N / 2, -1, int(N / 2))))/ time_window
f = c/wavelength + df
dw = 2 * np.pi * df
FWHM = 50 * fs
m = 4 * np.log(2)
A_t = np.exp(-m * t ** 2 * (1 / 2) / FWHM ** 2)
A_w = np.fft.fft(A_t)
GDD = 500 * fs*fs
TOD = 0 * fs*fs*fs
FOD = 0
A_w = np.exp(1j * (GDD / 2.0) * dw**2 +
                                1j * (TOD / 6.0) * dw ** 3+
                                1j * (FOD / 24.0) * dw ** 4) * A_w

fig_1 = plt.figure(1, facecolor='w', edgecolor='k')
ax_1 = fig_1.add_subplot(1, 1, 1)
ax_2 = ax_1.twinx()
ax_1.plot(np.fft.fftshift(f/THz),np.fft.fftshift(np.abs(A_w) ** 2 / max(np.abs(A_w) ** 2)),'b')
ax_2.plot(np.fft.fftshift(f/THz),np.fft.fftshift(np.angle(A_w)),'r')
ax_1.set_ylabel('Intensity / a.u.')
ax_2.set_ylabel('Phase / rad')
ax_1.tick_params(axis='y', colors='b')
ax_2.tick_params(axis='y', colors='r')

plt.xlim(300,450)
plt.show()

相位提取结果图：

问题原因 #

你看到的密集尖峰是**相位卷绕（Phase Wrapping）**导致的，属于正常现象：np.angle函数会将相位值强制限制在$[-\pi, \pi]$区间内，当真实相位的变化幅度超过这个区间时，就会出现跳变尖峰，并非计算错误。

还原多项式相位的方法 #

要得到连续的原始多项式相位，需要对提取的卷绕相位进行解卷绕操作，使用np.unwrap函数即可实现：

修改绘图部分的相位处理代码：

# 提取卷绕相位并解卷绕
wrapped_phase = np.angle(A_w)
unwrapped_phase = np.unwrap(wrapped_phase)
# 绘制解卷绕后的连续相位
ax_2.plot(np.fft.fftshift(f/THz), np.fft.fftshift(unwrapped_phase), 'r')

结果验证 #

解卷绕后，你会得到一条连续的相位曲线，将其与你定义的原始多项式相位（即(GDD/2.0)*dw**2 + (TOD/6.0)*dw**3 + (FOD/24.0)*dw**4）对比，两者会完全一致（仅可能存在傅里叶变换带来的微小线性相位偏移，若需要可进一步去除）。

另外，代码中df的生成可以简化为df = np.fft.fftfreq(N, d=t[1]-t[0])，这是numpy内置的频率轴生成函数，更简洁且不易出错。

内容的提问来源于stack exchange，提问作者Alex