包络线#
- scipy.signal.envelope(z, bp_in=(1, None), *, n_out=None, squared=False, residual='lowpass', axis=-1)[源代码]#
计算实值或复值信号的包络线。
- 参数:
- zndarray
实值或复值输入信号,假设由
n
个样本组成,采样间隔为T
。 z 也可以是多维数组,时间轴由 axis 定义。- bp_intuple[int | None, int | None], 可选
定义输入滤波器频带
bp_in[0]:bp_in[1]
的 2 元组。截止频率指定为1/(n*T)
的整数倍,允许的频率范围为-n//2 <= bp_in[0] < bp_in[1] <= (n+1)//2
。None
条目分别替换为-n//2
或(n+1)//2
。默认值(1, None)
会移除平均值以及负频率分量。- n_outint | None, 可选
如果不是
None
,则输出将重采样为 n_out 个样本。默认值None
将输出设置为与输入 z 相同的长度。- squaredbool, 可选
如果设置,则返回包络线的平方。由于所使用的绝对值函数的非线性性质,平方包络线的带宽通常小于非平方包络线的带宽。也就是说,嵌入的平方根函数通常会产生额外的谐波。默认值为
False
。- residualLiteral[‘lowpass’, ‘all’, None], 可选
此选项确定返回哪种残差,即输入带通滤波器移除的信号部分。
'all'
返回除频带bp_in[0]:bp_in[1]
内容之外的所有内容,'lowpass'
返回频带< bp_in[0]
的内容。如果None
,则仅返回包络线。默认值:'lowpass'
。- axisint, 可选
计算包络线的 z 的轴。默认值为最后一个轴。
- 返回:
- ndarray
如果参数 residual 为
None
,则返回一个与输入 z 具有相同形状的数组z_env
,其中包含其包络线。否则,返回一个形状为(2, *z.shape)
的数组,其中包含沿第一个轴堆叠的数组z_env
和z_res
。它允许解包,即z_env, z_res = envelope(z, residual='all')
。残差z_res
包含输入带通滤波器移除的信号部分,具体取决于参数 residual。请注意,对于实值信号,返回实值残差。因此,将忽略 bp_in 的负频率分量。
另请参阅
hilbert
通过希尔伯特变换计算解析信号。
注释
任何复值信号 \(z(t)\) 都可以用实值的瞬时幅度 \(a(t)\) 和实值的瞬时相位 \(\phi(t)\) 来描述,即 \(z(t) = a(t) \exp\!\big(j \phi(t)\big)\)。包络线定义为幅度 \(|a(t)| = |z(t)|\) 的绝对值,同时也是信号的绝对值。因此,\(|a(t)|\) “包络”了所有具有幅度 \(a(t)\) 和任意相位 \(\phi(t)\) 的信号类别。对于实值信号,\(x(t) = a(t) \cos\!\big(\phi(t)\big)\) 是类似的公式。因此,可以通过使用希尔伯特变换将 \(x(t)\) 转换为解析信号 \(z_a(t)\) 来确定 \(|a(t)|\),即 \(z_a(t) = a(t) \cos\!\big(\phi(t)\big) + j a(t) \sin\!\big(\phi(t) \big)\),它会产生一个具有相同包络线 \(|a(t)|\) 的复值信号。
该实现基于计算输入信号的 FFT,然后在傅里叶空间中执行必要的操作。因此,需要考虑典型的 FFT 警告。
该信号假定是周期性的。由于吉布斯现象,信号开始和结束之间的不连续性会导致不希望的结果。
如果信号长度是质数或很长,则 FFT 会很慢。此外,内存需求通常高于基于 FIR/IIR 滤波器的类似实现。
带通滤波器的截止频率的频率间隔
1 / (n*T)
对应于scipy.fft.fftfreq(len(z), T)
生成的频率。
如果需要不进行带通滤波的复值信号 z 的包络线,即
bp_in=(None, None)
,则包络线对应于绝对值。因此,使用np.abs(z)
而不是此函数会更有效。虽然基于解析信号计算包络线 [1] 是实值信号的自然方法,但也经常使用其他方法。最流行的替代方法可能是所谓的“平方律”包络检波器及其相关方法 [2]。它们并非总是为所有类型的信号计算出正确的结果,但对于大多数类型的窄带信号,通常是正确的,并且计算效率更高。此处给出的包络线定义在瞬时幅度和相位受关注的情况下很常见(例如,如 [3] 中所述)。还存在其他概念,这些概念依赖于包络线的通用数学思想 [4]:一种实用的方法是确定所有上限和下限信号峰值,并使用样条插值来确定曲线 [5]。
参考
[1]“解析信号”,维基百科,https://en.wikipedia.org/wiki/Analytic_signal
[4]“包络(数学)”,维基百科,https://en.wikipedia.org/wiki/Envelope_(mathematics)
示例
下图展示了一个具有可变频率和低频漂移的信号的包络。为了将漂移与包络分离,使用了 4 Hz 的高通滤波器。输入带通滤波器的低通残余部分用于确定一个非对称的上限和下限,以包围该信号。由于所得包络的平滑性,它从 500 个样本被下采样到 40 个样本。请注意,瞬时幅度
a_x
和计算出的包络x_env
并非完全相同。这是因为信号不是完全周期性的,以及x_carrier
和x_drift
之间存在一些频谱重叠。因此,它们不能被带通滤波器完全分离。>>> import matplotlib.pyplot as plt >>> import numpy as np >>> from scipy.signal.windows import gaussian >>> from scipy.signal import envelope ... >>> n, n_out = 500, 40 # number of signal samples and envelope samples >>> T = 2 / n # sampling interval for 2 s duration >>> t = np.arange(n) * T # time stamps >>> a_x = gaussian(len(t), 0.4/T) # instantaneous amplitude >>> phi_x = 30*np.pi*t + 35*np.cos(2*np.pi*0.25*t) # instantaneous phase >>> x_carrier = a_x * np.cos(phi_x) >>> x_drift = 0.3 * gaussian(len(t), 0.4/T) # drift >>> x = x_carrier + x_drift ... >>> bp_in = (int(4 * (n*T)), None) # 4 Hz highpass input filter >>> x_env, x_res = envelope(x, bp_in, n_out=n_out) >>> t_out = np.arange(n_out) * (n / n_out) * T ... >>> fg0, ax0 = plt.subplots(1, 1, tight_layout=True) >>> ax0.set_title(r"$4\,$Hz Highpass Envelope of Drifting Signal") >>> ax0.set(xlabel="Time in seconds", xlim=(0, n*T), ylabel="Amplitude") >>> ax0.plot(t, x, 'C0-', alpha=0.5, label="Signal") >>> ax0.plot(t, x_drift, 'C2--', alpha=0.25, label="Drift") >>> ax0.plot(t_out, x_res+x_env, 'C1.-', alpha=0.5, label="Envelope") >>> ax0.plot(t_out, x_res-x_env, 'C1.-', alpha=0.5, label=None) >>> ax0.grid(True) >>> ax0.legend() >>> plt.show()
第二个例子提供了复值信号的几何包络解释:以下两个图显示了复值信号为蓝色 3D 轨迹,包络为橙色圆管,其直径可变,即为 \(|a(t)| \exp(j\rho(t))\),其中 \(\rho(t)\in[-\pi,\pi]\)。此外,还描绘了轨迹和圆管在二维实部和虚部坐标平面上的投影。复值信号的每个点都接触到圆管的表面。
左图显示了一个解析信号,即虚部和实部之间的相位差始终为 90 度,从而形成螺旋轨迹。可以看出,在这种情况下,实部也具有预期的包络,即表示瞬时幅度的绝对值。
右图显示了解析信号的实部被解释为复值信号,即虚部为零。那里的结果包络不像解析情况那样平滑,并且在实平面中没有恢复瞬时幅度。如果
z_re
作为实值信号传递,即z_re = z.real
而不是z_re = z.real + 0j
,结果将与左图相同。原因是实值信号被解释为复值解析信号的实部。>>> import matplotlib.pyplot as plt >>> import numpy as np >>> from scipy.signal.windows import gaussian >>> from scipy.signal import envelope ... >>> n, T = 1000, 1/1000 # number of samples and sampling interval >>> t = np.arange(n) * T # time stamps for 1 s duration >>> f_c = 3 # Carrier frequency for signal >>> z = gaussian(len(t), 0.3/T) * np.exp(2j*np.pi*f_c*t) # analytic signal >>> z_re = z.real + 0j # complex signal with zero imaginary part ... >>> e_a, e_r = (envelope(z_, (None, None), residual=None) for z_ in (z, z_re)) ... >>> # Generate grids to visualize envelopes as 2d and 3d surfaces: >>> E2d_t, E2_amp = np.meshgrid(t, [-1, 1]) >>> E2d_1 = np.ones_like(E2_amp) >>> E3d_t, E3d_phi = np.meshgrid(t, np.linspace(-np.pi, np.pi, 300)) >>> ma = 1.8 # maximum axis values in real and imaginary direction ... >>> fg0 = plt.figure(figsize=(6.2, 4.)) >>> ax00 = fg0.add_subplot(1, 2, 1, projection='3d') >>> ax01 = fg0.add_subplot(1, 2, 2, projection='3d', sharex=ax00, ... sharey=ax00, sharez=ax00) >>> ax00.set_title("Analytic Signal") >>> ax00.set(xlim=(0, 1), ylim=(-ma, ma), zlim=(-ma, ma)) >>> ax01.set_title("Real-valued Signal") >>> for z_, e_, ax_ in zip((z, z.real), (e_a, e_r), (ax00, ax01)): ... ax_.set(xlabel="Time $t$", ylabel="Real Amp. $x(t)$", ... zlabel="Imag. Amp. $y(t)$") ... ax_.plot(t, z_.real, 'C0-', zs=-ma, zdir='z', alpha=0.5, label="Real") ... ax_.plot_surface(E2d_t, e_*E2_amp, -ma*E2d_1, color='C1', alpha=0.25) ... ax_.plot(t, z_.imag, 'C0-', zs=+ma, zdir='y', alpha=0.5, label="Imag.") ... ax_.plot_surface(E2d_t, ma*E2d_1, e_*E2_amp, color='C1', alpha=0.25) ... ax_.plot(t, z_.real, z_.imag, 'C0-', label="Signal") ... ax_.plot_surface(E3d_t, e_*np.cos(E3d_phi), e_*np.sin(E3d_phi), ... color='C1', alpha=0.5, shade=True, label="Envelope") ... ax_.view_init(elev=22.7, azim=-114.3) >>> fg0.subplots_adjust(left=0.08, right=0.97, wspace=0.15) >>> plt.show()