包络提取方法解读

包络提取方法解读

包络提取是信号处理中的一项重要技术，用于从复杂信号中提取其幅度变化的轮廓。简单来说，包络就是连接所有波峰的平滑曲线。

基本概念与数学基础

对于一个复数信号 $z(t)=x(t)+jy(t)$，其包络定义为模值：

$$|z(t)|=\sqrt{x^2(t)+y^2(t)}$$

对于实信号 $x(t)$，我们通常希望构造一个解析信号 $z(t)$，使得：

$$z(t)=A(t)e^{j\varphi (t)}=A(t)[\cos (\varphi (t))+j\sin (\varphi (t))]$$

其中 $A(t)=|z(t)|$ 是瞬时幅度（包络），$\varphi (t)=\arg (z(t))$ 是瞬时相位。

瞬时频率定义为相位的时间导数：

$${\omega }_i(t)=\frac{d\varphi (t)}{dt}=\frac{d}{dt}\arctan \left(\frac{y(t)}{x(t)}\right)$$

通过链式法则展开：

$${\omega }_i(t)=\frac{x(t)\frac{dy(t)}{dt}-y(t)\frac{dx(t)}{dt}}{x^2(t)+y^2(t)}$$

希尔伯特变换法

希尔伯特变换是包络提取中最经典和理论最完善的方法。对于实信号 $x(t)$，其希尔伯特变换定义为：

x^(t)=H{x(t)}=1πP.V.∫−∞∞x(τ)t−τdτ\hat{x}(t) = \mathcal{H}\{x(t)\} = \frac{1}{\pi} \text{P.V.} \int_{-\infty}^{\infty} \frac{x(\tau)}{t-\tau} d\taux^(t)=H{x(t)}=π1​P.V.∫−∞∞​t−τx(τ)​dτ

其中 P.V. 表示柯西主值积分。这个积分的数学意义可以通过极限来理解：

$$\hat{x}(t)=\underset{ϵ\to 0^{+}}{\lim }\frac{1}{\pi }\left[{\int }_{-\infty }^{t-ϵ}\frac{x(\tau )}{t-\tau }d\tau +{\int }_{t+ϵ}^{\infty }\frac{x(\tau )}{t-\tau }d\tau \right]$$

在频域中，希尔伯特变换对应于：

X^(ω)=F{x^(t)}=−j⋅sgn(ω)⋅X(ω)\hat{X}(\omega) = \mathcal{F}\{\hat{x}(t)\} = -j \cdot \text{sgn}(\omega) \cdot X(\omega)X^(ω)=F{x^(t)}=−j⋅sgn(ω)⋅X(ω)

其中符号函数为：

$$\text{sgn}(\omega )=\{\begin{array}{ll}1,& \omega >0\\ 0,& \omega =0\\ -1,& \omega <0\end{array}$$

解析信号的构造为：

$$z(t)=x(t)+j\hat{x}(t)$$

其频谱具有单侧性质：

$$Z(\omega )=\{\begin{array}{ll}2X(\omega ),& \omega >0\\ X(0),& \omega =0\\ 0,& \omega <0\end{array}$$

瞬时幅度（包络）为：

$$A(t)=|z(t)|=\sqrt{x^2(t)+{\hat{x}}^2(t)}$$

瞬时相位为：

$$\varphi (t)=\arctan \left(\frac{\hat{x}(t)}{x(t)}\right)$$

希尔伯特变换具有重要的数学性质。对于正弦函数 $x(t)=A\cos ({\omega }_{0}t+{\varphi }_0)$，其希尔伯特变换为：

$$\hat{x}(t)=A\sin ({\omega }_{0}t+{\varphi }_0)$$

这表明希尔伯特变换将余弦函数转换为同频率的正弦函数。

对于复杂信号，希尔伯特变换的卷积形式为：

$$\hat{x}(t)=x(t)\ast h(t)=x(t)\ast \frac{1}{\pi t}$$

其中 $h(t)=\frac{1}{\pi t}$ 是希尔伯特变换的冲激响应。

包络检波法

包络检波的数学模型基于整流和滤波操作。对于AM调制信号：

$$x(t)=A(t)[1+m(t)]\cos ({\omega }_{c}t+{\varphi }_c)$$

其中 $A(t)$ 是缓变包络，$m(t)$ 是调制信号，${\omega }_c$ 是载波频率。

全波整流后的信号为：

$$x_r(t)=|x(t)|=|A(t)[1+m(t)]||\cos ({\omega }_{c}t+{\varphi }_c)|$$

利用三角函数的傅里叶级数展开：

$$|\cos ({\omega }_{c}t+{\varphi }_c)|=\frac{2}{\pi }-\frac{4}{\pi }\sum _{n=1}^{\infty }\frac{\cos (2n{\omega }_{c}t+2n{\varphi }_c)}{4n^2-1}$$

整流信号包含直流分量和高频分量：

$$x_r(t)=\frac{2A(t)[1+m(t)]}{\pi }-\frac{4A(t)[1+m(t)]}{\pi }\sum _{n=1}^{\infty }\frac{\cos (2n{\omega }_{c}t+2n{\varphi }_c)}{4n^2-1}$$

通过截止频率为 ${\omega }_L$ 的低通滤波器，滤波器的传递函数为 $H(j\omega )$，则输出为：

y(t)=F−1{Xr(ω)⋅H(jω)}y(t) = \mathcal{F}^{-1}\{X_r(\omega) \cdot H(j\omega)\}y(t)=F−1{Xr​(ω)⋅H(jω)}

对于理想低通滤波器：

$$H(j\omega )=\{\begin{array}{ll}1,& |\omega |\le {\omega }_L\\ 0,& |\omega |>{\omega }_L\end{array}$$

如果选择 ${\omega }_{max}<{\omega }_L<2{\omega }_c$，其中 ${\omega }_{max}$ 是 $m(t)$ 的最高频率分量，则可以有效提取包络。

对于实际的一阶RC低通滤波器：

$$H(j\omega )=\frac{1}{1+j\omega RC}$$

其时域冲激响应为：

$$h(t)=\frac{1}{RC}{e}^{-\frac{t}{RC}}u(t)$$

包络估计为：

$$\hat{A}(t)={\int }_0^{\infty }{x}_r(t-\tau )\frac{1}{RC}{e}^{-\frac{\tau }{RC}}d\tau$$

这可以写成微分方程的形式：

$$RC\frac{d\hat{A}(t)}{dt}+\hat{A}(t)=x_r(t)$$

峰值检测法

峰值检测法的数学基础是寻找信号的局部极值。对于连续信号 $x(t)$，局部最大值点 $t_0$ 满足：

$$\frac{dx(t)}{dt}{|}_{t=t_0}=0\text{且}\frac{d^{2}x(t)}{d{t}^2}{|}_{t=t_0}<0$$

对于离散信号 $x[n]$，峰值检测条件为：

$$x[n_0]>x[n_0-1]\text{且}x[n_0]>x[n_0+1]$$

更严格的峰值检测可以使用二阶差分：

$${\Delta }^{2}x[n]=x[n+1]-2x[n]+x[n-1]$$

峰值点处有 ${\Delta }^{2}x[n_0]<0$。

检测到峰值点 $(t_i,x_i)$ 后，使用三次样条插值重构包络。三次样条函数在区间 $[t_i,t_{i+1}]$ 上定义为：

$$S_i(t)=a_i(t-t_i{)}^3+b_i(t-t_i{)}^2+c_i(t-t_i)+d_i$$

边界条件和连续性条件给出参数方程组：

$$\left[\begin{array}{cccc}1& 0& 0& 0\\ 1& h_i& h_i^2& h_i^3\\ 0& 1& 0& 0\\ 0& 1& 2h_i& 3h_i^2\end{array}\right]\left[\begin{array}{c}{d}_i\\ c_i\\ b_i\\ a_i\end{array}\right]=\left[\begin{array}{c}{x}_i\\ x_{i+1}\\ S_i^{\prime }(t_i)\\ S_i^{\prime }(t_{i+1})\end{array}\right]$$

其中 $h_i=t_{i+1}-t_i$。

自然样条的边界条件为 $S_0^{\prime \prime }(t_0)=S_{n-1}^{\prime \prime }(t_n)=0$，导致三对角线性方程组：

$$\left[\begin{array}{ccccc}2(h_0+h_1)& h_1& 0& \cdots & 0\\ h_1& 2(h_1+h_2)& h_2& \cdots & 0\\ ⋮& ⋮& \ddots & \ddots & ⋮\\ 0& 0& \cdots & h_{n-2}& 2(h_{n-2}+h_{n-1})\end{array}\right]\left[\begin{array}{c}{M}_1\\ M_2\\ ⋮\\ M_{n-1}\end{array}\right]=\left[\begin{array}{c}6f[t_0,t_1,t_2]\\ 6f[t_1,t_2,t_3]\\ ⋮\\ 6f[t_{n-2},t_{n-1},t_n]\end{array}\right]$$

其中 $M_i=S_i^{\prime \prime }(t_i)$，$f[t_i,t_{i+1},t_{i+2}]$ 是二阶差商。

滑动窗口法

滑动窗口法通过局部统计量估计包络。对于长度为 $L$ 的窗口，RMS包络定义为：

$$A_{RMS}[n]=\sqrt{\frac{1}{L}\sum _{k=n-L+1}^n{x}^2[k]}$$

这等价于信号平方通过移动平均滤波器：

$$A_{RMS}^2[n]=x^2[n]\ast w[n]$$

其中 $w[n]=\frac{1}{L}$ 对 $0\le n\le L-1$，否则为0。

在频域中，移动平均滤波器的频率响应为：

$$W(e^{j\omega })=\frac{1}{L}\cdot \frac{\sin (\omega L/2)}{\sin (\omega /2)}\cdot e^{-j\omega (L-1)/2}$$

其幅度响应为：

$$|W(e^{j\omega })|=\frac{1}{L}\left|\frac{\sin (\omega L/2)}{\sin (\omega /2)}\right|$$

对于指数加权移动平均，权重函数为：

$$w[k]=\alpha (1-\alpha {)}^k,k=0,1,2,\dots$$

其中 $0<\alpha <1$ 是遗忘因子。指数加权RMS包络为：

$$A_{EWMA}^2[n]=\alpha x^2[n]+(1-\alpha )A_{EWMA}^2[n-1]$$

这对应于一阶IIR滤波器：

$$H(z)=\frac{\alpha }{1-(1-\alpha )z^{-1}}$$

自适应窗口长度可以基于信号的局部统计特性。定义局部方差为：

$${\sigma }^2[n]=\frac{1}{L_0}\sum _{k=n-L_0+1}^n(x[k]-\bar{x}[n]{)}^2$$

其中 $\bar{x}[n]=\frac{1}{L_0}{\sum }_{k=n-L_0+1}^{n}x[k]$ 是局部均值。

自适应窗口长度为：

$$L[n]=L_{min}+(L_{max}-L_{min})\cdot \exp \left(-\frac{{\sigma }^2[n]}{{\sigma }_0^2}\right)$$

其中 ${\sigma }_0^2$ 是归一化参数。

经验模态分解法

EMD的数学基础是将信号自适应地分解为固有模态函数（IMF）。对于信号 $x(t)$，EMD分解为：

$$x(t)=\sum _{i=1}^n{c}_i(t)+r_n(t)$$

其中 $c_i(t)$ 是第 $i$ 个IMF，$r_n(t)$ 是残余分量。

筛选过程的数学描述如下。设上包络为 $e_{max}(t)$，下包络为 $e_{min}(t)$，均值包络为：

$$m(t)=\frac{e_{max}(t)+e_{min}(t)}{2}$$

第 $k$ 次筛选的结果为：

$$h_k(t)=h_{k-1}(t)-m_{k-1}(t)$$

其中 $h_0(t)=x(t)$。

IMF的判定准则包括：

1. 极值点数量条件：$|N_{ext}-N_{zero}|\le 1$，其中 $N_{ext}$ 是极值点数量，$N_{zero}$ 是零交叉点数量。

2. 均值条件：$\frac{1}{T}{\int }_0^T|m(t){|}^{2}dt<ϵ$，其中 $T$ 是信号长度，$ϵ$ 是阈值。

3. 标准差条件：$SD={\sum }_{t=0}^T\frac{|h_{k-1}(t)-h_k(t){|}^2}{h_{k-1}^2(t)}<\delta$

包络插值通常使用三次样条。对于 $n$ 个极值点 $(t_i,x_i)$，三次样条插值的矩阵形式为：

$$\mathbf{AM}=\mathbf{B}$$

其中 $\mathbf{A}$ 是三对角矩阵：

$$A_{i,j}=\{\begin{array}{ll}{h}_{i-1},& j=i-1\\ 2(h_{i-1}+h_i),& j=i\\ h_i,& j=i+1\\ 0,& \text{其他}\end{array}$$

$\mathbf{M}=[M_1,M_2,\dots ,M_{n-1}{]}^T$ 是二阶导数向量，$\mathbf{B}$ 是右端向量。

集合经验模态分解（EEMD）通过添加白噪声改善模态混叠问题：

$$x_i(t)=x(t)+n_i(t)$$

其中 $n_i(t)$ 是第 $i$ 次实现的白噪声，方差为 ${\sigma }^2$。最终的IMF为：

$$c_j(t)=\frac{1}{N}\sum _{i=1}^N{c}_{j,i}(t)$$

其中 $c_{j,i}(t)$ 是第 $i$ 次噪声实现下的第 $j$ 个IMF。

小波变换法

连续小波变换（CWT）提供了时频分析的强大工具。对于信号 $x(t)$ 和母小波 $\psi (t)$，CWT定义为：

$$W_x(a,b)=\frac{1}{\sqrt{a}}{\int }_{-\infty }^{\infty }x(t){\psi }^{\ast }\left(\frac{t-b}{a}\right)dt$$

其中 $a>0$ 是尺度参数，$b\in \mathbb{R}$ 是平移参数。

母小波必须满足可容许性条件：

$$C_{\psi }={\int }_{-\infty }^{\infty }\frac{|\Psi (\omega ){|}^2}{|\omega |}d\omega <\infty$$

其中 $\Psi (\omega )$ 是 $\psi (t)$ 的傅里叶变换。这个条件隐含了 $\Psi (0)=0$，即小波函数的零均值特性。

常用的Morlet小波定义为：

$$\psi (t)={\pi }^{-1/4}{e}^{-t^2/2}{e}^{j{\omega }_{0}t}$$

其傅里叶变换为：

$$\Psi (\omega )={\pi }^{1/4}\sqrt{2}{e}^{-(\omega -{\omega }_0{)}^2/2}$$

Mexican hat小波（二阶高斯导数）为：

$$\psi (t)=\frac{2}{\sqrt{3}{\pi }^{1/4}}(1-t^2)e^{-t^2/2}$$

小波变换的模值 $|W_x(a,b)|$ 提供了信号在时间-尺度平面上的能量分布。对于包络提取，脊线算法寻找使小波系数模值最大的尺度曲线：

$$a_{ridge}(b)=\arg \underset{a}{\max }|W_x(a,b)|$$

脊线检测可以通过梯度方法实现：

$$\frac{\partial }{\partial a}|W_x(a,b)|=0$$

这导致：

$$\text{Re}\left\{W_x(a,b)\frac{\partial W_x^{\ast }(a,b)}{\partial a}\right\}=0$$

其中：

$$\frac{\partial W_x(a,b)}{\partial a}=-\frac{1}{2a^{3/2}}{\int }_{-\infty }^{\infty }x(t)\left[{\psi }^{\ast }\left(\frac{t-b}{a}\right)+\frac{t-b}{a}{\psi }^{\prime \ast }\left(\frac{t-b}{a}\right)\right]dt$$

基于脊线的瞬时频率估计为：

$${\omega }_i(b)=\frac{{\omega }_{\psi }}{2\pi a_{ridge}(b)}$$

其中 ${\omega }_{\psi }$ 是母小波的中心频率。

小波包络估计为：

$$A(t)=|W_x(a_{ridge}(t),t)|$$

数学形态学方法

数学形态学基于集合论和拓扑学，为信号处理提供了非线性滤波工具。对于一维信号 $f:\mathbb{Z}\to \mathbb{R}$ 和结构元素 $B\subset \mathbb{Z}$，基本形态学操作定义为：

腐蚀：$(f\ominus B)(n)={\inf }_{b\in B}f(n+b)$

膨胀：$(f\oplus B)(n)={\sup }_{b\in B}f(n-b)$

开运算：$f\circ B=(f\ominus B)\oplus B$

闭运算：$f\bullet B=(f\oplus B)\ominus B$

这些操作具有重要的代数性质。腐蚀和膨胀是对偶操作：

$$(f\ominus B{)}^c=f^c\oplus \hat{B}$$

其中 $f^c$ 是 $f$ 的补集，B^={−b:b∈B}\hat{B} = \{-b : b \in B\}B^={−b:b∈B} 是 $B$ 的反射。

开运算和闭运算是幂等的：

$$(f\circ B)\circ B=f\circ B$$
$$(f\bullet B)\bullet B=f\bullet B$$

并且满足对偶性：

$$(f\circ B{)}^c=f^c\bullet \hat{B}$$

形态学梯度定义为膨胀和腐蚀的差：

$${\nabla }_{B}f=(f\oplus B)-(f\ominus B)$$

这提供了边缘检测的工具。

对于包络提取，交替序列滤波器（ASF）特别有用：

$${\text{ASF}}_n^{\circ \bullet }(f)=((\cdots ((f\circ B_1)\bullet B_1)\circ B_2)\bullet B_2\cdots )\circ B_n)\bullet B_n$$

其中 $B_1\subset B_2\subset \cdots \subset B_n$ 是递增的结构元素序列。

形态学包络可以定义为：

上包络：$e_{upper}(n)=(f\bullet B)(n)$

下包络：$e_{lower}(n)=(f\circ B)(n)$

信号的形态学包络为：

$$A(n)=\frac{e_{upper}(n)+|e_{lower}(n)|}{2}$$

结构元素的优化可以通过最小化重构误差来实现：

$$B_{opt}=\arg \underset{B}{\min }\sum _n|f(n)-{\text{Recon}}_B(f)(n){|}^2$$

其中 ${\text{Recon}}_B(f)$ 是基于结构元素 $B$ 的形态学重构。

滤波器组方法

滤波器组方法将信号分解到多个频率子带，每个子带独立进行包络提取。完美重构滤波器组满足：

$$\sum _{k=0}^{M-1}{H}_k(z)G_k(z)=c{z}^{-n_0}$$

其中 $H_k(z)$ 是分析滤波器，$G_k(z)$ 是综合滤波器，$c$ 是常数，$n_0$ 是延迟。

对于均匀DFT滤波器组，分析滤波器为：

$$H_k(z)=H(z{W}_M^{-k})$$

其中 $W_M=e^{j2\pi /M}$，$H(z)$ 是原型滤波器。

子带信号为：

$$x_k[n]=\sum _{m}x[m]h_k[n-m]$$

每个子带的希尔伯特变换为：

$${\hat{x}}_k[n]=\sum _m{x}_k[m]\frac{\sin (\pi (n-m)/2)}{\pi (n-m)/2}$$

子带包络为：

$$A_k[n]=\sqrt{x_k^2[n]+{\hat{x}}_k^2[n]}$$

加权合成的总包络为：

$$A[n]=\sum _{k=0}^{M-1}{w}_k[n]A_k[n]$$

自适应权重基于子带能量：

$$w_k[n]=\frac{E_k[n]}{{\sum }_{i=0}^{M-1}{E}_i[n]}$$

其中瞬时能量为：

$$E_k[n]=\alpha E_k[n-1]+(1-\alpha )(x_k^2[n]+{\hat{x}}_k^2[n])$$

多相实现提高了计算效率。多相分解为：

$$H_k(z)=\sum _{l=0}^{M-1}{z}^{-l}{E}_{k,l}(z^M)$$

其中 $E_{k,l}(z)$ 是多相分量。