envelope#
- 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], optional
定義輸入濾波器的頻帶
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, optional
如果不是
None，則輸出將重新取樣為 n_out 個樣本。預設值None將輸出設定為與輸入 z 相同的長度。- squaredbool, optional
如果設定，則傳回包絡線的平方。由於所使用的絕對值函數的非線性性質，平方包絡線的頻寬通常小於非平方包絡線頻寬。也就是說，嵌入式平方根函數通常會產生額外的諧波。預設值為
False。- residualLiteral[‘lowpass’, ‘all’, None], optional
此選項決定了要傳回哪種殘差，即輸入帶通濾波器移除的訊號部分。
'all'傳回除了頻帶bp_in[0]:bp_in[1]的內容之外的所有內容，'lowpass'傳回頻帶< bp_in[0]的內容。如果None，則僅傳回包絡線。預設值：'lowpass'。- axisint, optional
計算包絡線的 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透過 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)|\)，透過 Hilbert 轉換，即 \(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]\)。此外，還描繪了軌跡和管子在 2d 實數和虛數座標平面中的投影。複值訊號的每個點都接觸管子的表面。
左圖顯示了解析訊號，即虛部和實部之間的相位差始終為 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()
