同步挤压变换(Synchrosqueezing Transform, SST)详解
背景与动机
在信号处理领域,我们经常需要分析包含多个频率成分的复杂信号。传统的傅里叶变换虽然可以将信号从时域转换到频域,但它无法同时提供时间和频率的局部信息。随后发展的短时傅里叶变换(STFT)和连续小波变换(CWT)虽然能够提供时频表示,但存在时频分辨率限制的问题,导致时频图中能量扩散现象,使得信号的瞬时频率难以精确定位。同步挤压变换旨在克服这些局限性,通过"挤压"时频平面上的能量分布,提高时频表示的分辨率和可读性。
基本原理
同步挤压变换的核心思想是利用信号的相位信息来估计瞬时频率,然后将传统时频表示(如CWT或STFT)中的能量重新分配到这些估计的瞬时频率位置上,从而"挤压"能量分布,形成更加集中和清晰的时频表示。
连续小波变换(CWT)回顾
同步挤压变换通常基于连续小波变换实现,所以首先回顾一下CWT的定义。对于信号 f ( t ) f(t) f(t),其连续小波变换为:
W f ( a , b ) = ∫ − ∞ ∞ f ( t ) 1 a ψ ∗ ( t − b a ) d t W_f(a,b) = \int_{-\infty}^{\infty} f(t) \frac{1}{\sqrt{a}} \psi^*\left(\frac{t-b}{a}\right) dt Wf(a,b)=∫−∞∞f(t)a1ψ∗(at−b)dt
其中:
- ψ ( t ) \psi(t) ψ(t) 是母小波函数
- ψ ∗ ( t ) \psi^*(t) ψ∗(t) 表示 ψ ( t ) \psi(t) ψ(t) 的复共轭
- a > 0 a > 0 a>0 是尺度参数,与频率成反比
- b b b 是时间平移参数
- W f ( a , b ) W_f(a,b) Wf(a,b) 是在尺度 a a a 和时间位置 b b b 的小波系数
如果采用傅里叶变换的角度看待CWT,我们可以将其重写为:
W f ( a , b ) = 1 2 π ∫ − ∞ ∞ f ^ ( ξ ) a ψ ^ ∗ ( a ξ ) e i b ξ d ξ W_f(a,b) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \hat{f}(\xi) \sqrt{a}\hat{\psi}^*(a\xi)e^{ib\xi}d\xi Wf(a,b)=2π1∫−∞∞f^(ξ)aψ^∗(aξ)eibξdξ
其中 f ^ ( ξ ) \hat{f}(\xi) f^(ξ) 和 ψ ^ ( ξ ) \hat{\psi}(\xi) ψ^(ξ) 分别是 f ( t ) f(t) f(t) 和 ψ ( t ) \psi(t) ψ(t) 的傅里叶变换。
瞬时频率估计的推导
考虑一个理想的单分量信号 f ( t ) = A ( t ) e i ϕ ( t ) f(t) = A(t)e^{i\phi(t)} f(t)=A(t)eiϕ(t),其中 A ( t ) A(t) A(t) 是幅度调制, ϕ ( t ) \phi(t) ϕ(t) 是相位函数。如果 A ( t ) A(t) A(t) 和 ϕ ′ ( t ) \phi'(t) ϕ′(t) 的变化相对于 ϕ ( t ) \phi(t) ϕ(t) 较慢,则可以证明其小波变换近似为:
W f ( a , b ) ≈ 1 2 A ( b ) a ψ ^ ∗ ( a ϕ ′ ( b ) ) e i ϕ ( b ) W_f(a,b) \approx \frac{1}{2} A(b) \sqrt{a} \hat{\psi}^*(a\phi'(b)) e^{i\phi(b)} Wf(a,b)≈21A(b)aψ^∗(aϕ′(b))eiϕ(b)
对此表达式求关于 b b b 的导数,可以得到:
∂ b W f ( a , b ) ≈ 1 2 ∂ b A ( b ) a ψ ^ ∗ ( a ϕ ′ ( b ) ) e i ϕ ( b ) + 1 2 A ( b ) a ∂ b [ ψ ^ ∗ ( a ϕ ′ ( b ) ) ] e i ϕ ( b ) + 1 2 A ( b ) a ψ ^ ∗ ( a ϕ ′ ( b ) ) i ϕ ′ ( b ) e i ϕ ( b ) \partial_b W_f(a,b) \approx \frac{1}{2} \partial_b A(b) \sqrt{a} \hat{\psi}^*(a\phi'(b)) e^{i\phi(b)} + \frac{1}{2} A(b) \sqrt{a} \partial_b[\hat{\psi}^*(a\phi'(b))] e^{i\phi(b)} + \frac{1}{2} A(b) \sqrt{a} \hat{\psi}^*(a\phi'(b)) i\phi'(b) e^{i\phi(b)} ∂bWf(a,b)≈21∂bA(b)aψ^∗(aϕ′(b))eiϕ(b)+21A(b)a∂b[ψ^∗(aϕ′(b))]eiϕ(b)+21A(b)aψ^∗(aϕ′(b))iϕ′(b)eiϕ(b)
在假设 A ( t ) A(t) A(t) 和 ϕ ′ ( t ) \phi'(t) ϕ′(t) 变化缓慢的条件下,前两项相对于第三项可以忽略,因此:
∂ b W f ( a , b ) ≈ i ϕ ′ ( b ) W f ( a , b ) \partial_b W_f(a,b) \approx i\phi'(b) W_f(a,b) ∂bWf(a,b)≈iϕ′(b)Wf(a,b)
这样就得到了瞬时频率的估计公式:
ω f ( a , b ) = − i ∂ b W f ( a , b ) W f ( a , b ) ≈ ϕ ′ ( b ) \omega_f(a,b) = -i \frac{\partial_b W_f(a,b)}{W_f(a,b)} \approx \phi'(b) ωf(a,b)=−iWf(a,b)∂bWf(a,b)≈ϕ′(b)
这个结果表明,通过小波系数的对数导数,我们可以恢复信号的瞬时频率 ϕ ′ ( b ) \phi'(b) ϕ′(b),这正是同步挤压变换的核心原理。
能量重分配的理论基础
同步挤压变换的能量重分配过程可以从数学上用积分变换的视角理解。考虑从尺度-时间平面 ( a , b ) (a,b) (a,b) 到时间-频率平面 ( b , ω ) (b,\omega) (b,ω) 的变换:
T f ( b , ω ) = ∫ { a : ω f ( a , b ) = ω } W f ( a , b ) ⋅ a − 3 / 2 d a T_f(b,\omega) = \int_{\{a: \omega_f(a,b)=\omega\}} W_f(a,b) \cdot a^{-3/2} da Tf(b,ω)=∫{a:ωf(a,b)=ω}Wf(a,b)⋅a−3/2da
这个表达式可以通过狄拉克δ函数更为优雅地表示:
T f ( b , ω ) = ∫ 0 ∞ W f ( a , b ) ⋅ a − 3 / 2 ⋅ δ ( ω − ω f ( a , b ) ) d a T_f(b,\omega) = \int_{0}^{\infty} W_f(a,b) \cdot a^{-3/2} \cdot \delta(\omega - \omega_f(a,b)) da Tf(b,ω)=∫0∞Wf(a,b)⋅a−3/2⋅δ(ω−ωf(a,b))da
从理论上讲,对于理想的单分量信号 f ( t ) = A ( t ) e i ϕ ( t ) f(t) = A(t)e^{i\phi(t)} f(t)=A(t)eiϕ(t),其同步挤压变换应该在频率 ω = ϕ ′ ( b ) \omega = \phi'(b) ω=ϕ′(b) 处呈现δ函数,即:
T f ( b , ω ) ≈ 1 2 A ( b ) e i ϕ ( b ) δ ( ω − ϕ ′ ( b ) ) T_f(b,\omega) \approx \frac{1}{2} A(b) e^{i\phi(b)} \delta(\omega - \phi'(b)) Tf(b,ω)≈21A(b)eiϕ(b)δ(ω−ϕ′(b))
这表明同步挤压变换能够完美地解析信号的瞬时频率。
高阶同步挤压变换的数学框架
二阶同步挤压变换的精确推导
二阶同步挤压变换的目的是处理频率变化较快的信号。对于这类信号,一阶近似不够精确,需要考虑信号相位的二阶导数。考虑信号 f ( t ) = A ( t ) e i ϕ ( t ) f(t) = A(t)e^{i\phi(t)} f(t)=A(t)eiϕ(t),其中 ϕ ′ ′ ( t ) ≠ 0 \phi''(t) \neq 0 ϕ′′(t)=0。通过Taylor展开,可以得到小波变换在时间 b b b 附近的近似:
W f ( a , t ) ≈ 1 2 A ( b ) a ψ ^ ∗ ( a ϕ ′ ( b ) ) e i ϕ ( b ) e i ϕ ′ ( b ) ( t − b ) e i ϕ ′ ′ ( b ) 2 ( t − b ) 2 W_f(a,t) \approx \frac{1}{2} A(b) \sqrt{a} \hat{\psi}^*(a\phi'(b)) e^{i\phi(b)} e^{i\phi'(b)(t-b)} e^{i\frac{\phi''(b)}{2}(t-b)^2} Wf(a,t)≈21A(b)aψ^∗(aϕ′(b))eiϕ(b)eiϕ′(b)(t−b)ei2ϕ′′(b)(t−b)2
对此表达式进行时间导数,并考虑到二阶项的影响,可以得到更为精确的瞬时频率估计:
ω f ( 2 ) ( a , b ) = ϕ ′ ( b ) + ϕ ′ ′ ( b ) 2 ⋅ ∂ t 2 W f ( a , t ) ∣ t = b ⋅ W f ( a , b ) − ( ∂ t W f ( a , t ) ∣ t = b ) 2 ( ∂ t W f ( a , t ) ∣ t = b ) 2 \omega_f^{(2)}(a,b) = \phi'(b) + \frac{\phi''(b)}{2} \cdot \frac{\partial_t^2 W_f(a,t)|_{t=b} \cdot W_f(a,b) - (\partial_t W_f(a,t)|_{t=b})^2}{(\partial_t W_f(a,t)|_{t=b})^2} ωf(2)(a,b)=ϕ′(b)+2ϕ′′(b)⋅(∂tWf(a,t)∣t=b)2∂t2Wf(a,t)∣t=b⋅Wf(a,b)−(∂tWf(a,t)∣t=b)2
简化后可以表示为:
ω f ( 2 ) ( a , b ) = ω f ( a , b ) + 1 2 ⋅ ∂ b 2 W f ( a , b ) ⋅ W f ( a , b ) − ( ∂ b W f ( a , b ) ) 2 ( ∂ b W f ( a , b ) ) 2 ⋅ ∂ b ω f ( a , b ) \omega_f^{(2)}(a,b) = \omega_f(a,b) + \frac{1}{2} \cdot \frac{\partial_b^2 W_f(a,b) \cdot W_f(a,b) - (\partial_b W_f(a,b))^2}{(\partial_b W_f(a,b))^2} \cdot \partial_b \omega_f(a,b) ωf(2)(a,b)=ωf(a,b)+21⋅(∂bWf(a,b))2∂b2Wf(a,b)⋅Wf(a,b)−(∂bWf(a,b))2⋅∂bωf(a,b)
多分量信号的数学描述
对于包含多个频率成分的信号:
f ( t ) = ∑ k = 1 K f k ( t ) = ∑ k = 1 K A k ( t ) e i ϕ k ( t ) f(t) = \sum_{k=1}^{K} f_k(t) = \sum_{k=1}^{K} A_k(t)e^{i\phi_k(t)} f(t)=k=1∑Kfk(t)=k=1∑KAk(t)eiϕk(t)
其同步挤压变换在理想条件下应该能够分离各个成分:
T f ( b , ω ) ≈ ∑ k = 1 K 1 2 A k ( b ) e i ϕ k ( b ) δ ( ω − ϕ k ′ ( b ) ) T_f(b,\omega) \approx \sum_{k=1}^{K} \frac{1}{2} A_k(b) e^{i\phi_k(b)} \delta(\omega - \phi_k'(b)) Tf(b,ω)≈k=1∑K21Ak(b)eiϕk(b)δ(ω−ϕk′(b))
这表明同步挤压变换可以同时分辨多个瞬时频率轨迹,这是其优于传统时频分析方法的重要特点。
变分模式分解与同步挤压变换的统一框架
变分模式分解(VMD)是另一种处理非平稳信号的方法,通过变分原理将信号分解为一系列带宽有限的内蕴模态函数(IMF)。可以证明,在特定条件下,同步挤压变换和VMD存在数学上的联系。
定义变分问题:
min { u k } , { ω k } ∑ k = 1 K ∥ ∂ t [ ( δ ( t ) + j π t ) ∗ u k ( t ) ] e − j ω k t ∥ 2 2 \min_{\{u_k\}, \{\omega_k\}} \sum_{k=1}^{K} \left\| \partial_t \left[ \left( \delta(t) + \frac{j}{\pi t} \right) * u_k(t) \right] e^{-j\omega_k t} \right\|_2^2 {uk},{ωk}mink=1∑K∥∥∥∥∂t[(δ(t)+πtj)∗uk(t)]e−jωkt∥∥∥∥22
满足约束 ∑ k = 1 K u k = f \sum_{k=1}^{K} u_k = f ∑k=1Kuk=f,其中 u k u_k uk 是分解得到的模态, ω k \omega_k ωk 是中心频率。可以证明,对于线性调频信号,同步挤压变换的峰值轨迹与VMD的最优中心频率 ω k \omega_k ωk 在数学上是等价的:
arg max ω ∣ T f ( b , ω ) ∣ = ω k ( b ) \arg\max_{\omega} |T_f(b,\omega)| = \omega_k(b) argωmax∣Tf(b,ω)∣=ωk(b)
这一等价性为同步挤压变换提供了另一种理论解释。
广义同步挤压变换的数学框架
广义同步挤压变换将经典SST的概念扩展到更广泛的时频表示方法。设 T f \mathcal{T}_f Tf 为任意线性时频变换,定义其相应的瞬时频率估计为:
ω T ( p , q ) = − i ∂ q T f ( p , q ) T f ( p , q ) \omega_{\mathcal{T}}(p,q) = -i \frac{\partial_q \mathcal{T}_f(p,q)}{\mathcal{T}_f(p,q)} ωT(p,q)=−iTf(p,q)∂qTf(p,q)
其中 ( p , q ) (p,q) (p,q) 是时频表示的参数。广义同步挤压变换定义为:
S T ( p , ω ) = ∫ T f ( p , q ) δ ( ω − ω T ( p , q ) ) d q S_{\mathcal{T}}(p,\omega) = \int \mathcal{T}_f(p,q) \delta(\omega - \omega_{\mathcal{T}}(p,q)) dq ST(p,ω)=∫Tf(p,q)δ(ω−ωT(p,q))dq
这一框架可以统一地描述基于不同线性时频表示(如STFT、CWT、S变换等)的同步挤压变换。
同步挤压变换的不确定性原理分析
传统的时频分析方法受不确定性原理的限制,即时间和频率分辨率无法同时达到最优。同步挤压变换虽然在视觉上提高了时频图的锐度,但从理论上讲,它仍然受到基础变换(如CWT或STFT)的不确定性限制。然而,可以证明,对于满足特定条件的信号,同步挤压变换可以在数学上突破这一限制。考虑信号 f ( t ) = A ( t ) e i ϕ ( t ) f(t) = A(t)e^{i\phi(t)} f(t)=A(t)eiϕ(t),如果 A ( t ) A(t) A(t) 和 ϕ ′ ( t ) \phi'(t) ϕ′(t) 满足:
max ( ∣ ∂ t A ( t ) ∣ A ( t ) , ∣ ∂ t ϕ ′ ( t ) ∣ ϕ ′ ( t ) ) ≪ ϕ ′ ( t ) \max\left(\frac{|\partial_t A(t)|}{A(t)}, \frac{|\partial_t \phi'(t)|}{\phi'(t)}\right) \ll \phi'(t) max(A(t)∣∂tA(t)∣,ϕ′(t)∣∂tϕ′(t)∣)≪ϕ′(t)
则其同步挤压变换在频率维度上的分辨率可以达到狄拉克δ函数的级别,从而在数学上"绕过"了不确定性原理的限制。
离散同步挤压变换的数学表述
在实际应用中,我们需要处理离散信号,因此需要离散化同步挤压变换。设离散信号为 f [ n ] , n = 0 , 1 , . . . , N − 1 f[n], n = 0, 1, ..., N-1 f[n],n=0,1,...,N−1,其离散小波变换为:
W f [ j , n ] = ∑ m = 0 N − 1 f [ m ] ψ j , n ∗ [ m ] W_f[j,n] = \sum_{m=0}^{N-1} f[m] \psi_{j,n}^*[m] Wf[j,n]=m=0∑N−1f[m]ψj,n∗[m]
其中 ψ j , n [ m ] = 2 − j / 2 ψ [ 2 − j ( m − n ) ] \psi_{j,n}[m] = 2^{-j/2} \psi[2^{-j}(m-n)] ψj,n[m]=2−j/2ψ[2−j(m−n)] 是离散小波函数。
离散瞬时频率估计可以通过中心差分近似计算:
ω f [ j , n ] = − i W f [ j , n + 1 ] − W f [ j , n − 1 ] 2 ⋅ W f [ j , n ] \omega_f[j,n] = -i \frac{W_f[j,n+1] - W_f[j,n-1]}{2 \cdot W_f[j,n]} ωf[j,n]=−i2⋅Wf[j,n]Wf[j,n+1]−Wf[j,n−1]
离散同步挤压变换定义为:
T f [ n , k ] = ∑ j : ∣ ω f [ j , n ] − ω k ∣ < Δ ω / 2 W f [ j , n ] ⋅ 2 − 3 j / 2 T_f[n,k] = \sum_{j: |\omega_f[j,n] - \omega_k| < \Delta \omega/2} W_f[j,n] \cdot 2^{-3j/2} Tf[n,k]=j:∣ωf[j,n]−ωk∣<Δω/2∑Wf[j,n]⋅2−3j/2
其中 ω k = 2 π k / N , k = 0 , 1 , . . . , N − 1 \omega_k = 2\pi k/N, k = 0, 1, ..., N-1 ωk=2πk/N,k=0,1,...,N−1 是离散频率点。
小波基函数选择的理论分析
小波基函数的选择对同步挤压变换的性能有显著影响。从理论上讲,理想的小波应满足以下条件:
-
频率局部化性好:小波的傅里叶变换 ψ ^ ( ξ ) \hat{\psi}(\xi) ψ^(ξ) 应在某个中心频率 ξ 0 \xi_0 ξ0 附近集中,且衰减迅速;
-
支持紧凑性:在时域或频域具有紧凑支持,以提高计算效率;
-
正交性或双正交性:有利于信号的完美重构;
-
足够的消失矩:能够抑制多项式趋势。
最常用的是Morlet小波和Bump小波,Morlet小波定义为:
ψ σ ( t ) = π − 1 / 4 e − t 2 / ( 2 σ 2 ) ( e i ω 0 t − κ σ ) \psi_{\sigma}(t) = \pi^{-1/4} e^{-t^2/(2\sigma^2)} (e^{i\omega_0 t} - \kappa_{\sigma}) ψσ(t)=π−1/4e−t2/(2σ2)(eiω0t−κσ)
其中 κ σ = e − σ 2 ω 0 2 / 2 \kappa_{\sigma} = e^{-\sigma^2 \omega_0^2/2} κσ=e−σ2ω02/2 是归一化因子, ω 0 \omega_0 ω0 是中心频率, σ \sigma σ 控制时频分辨率。
Bump小波定义为:
ψ ( t ) = F − 1 { e − 1 1 − ( ξ − μ ) 2 / σ 2 χ [ μ − σ , μ + σ ] ( ξ ) } ( t ) \psi(t) = \mathcal{F}^{-1} \left\{ e^{-\frac{1}{1-(\xi-\mu)^2/\sigma^2}} \chi_{[\mu-\sigma, \mu+\sigma]}(\xi) \right\}(t) ψ(t)=F−1{e−1−(ξ−μ)2/σ21χ[μ−σ,μ+σ](ξ)}(t)
其中 χ [ a , b ] \chi_{[a,b]} χ[a,b] 是区间 [ a , b ] [a,b] [a,b] 上的特征函数, μ \mu μ 是中心频率, σ \sigma σ 控制带宽。
从理论上可以证明,对于线性调频信号,Bump小波比Morlet小波具有更好的性能。定义性能度量:
ρ ( ψ , f ) = ∥ T f ( b , ϕ ′ ( b ) ) ∥ 2 2 ∬ ∣ T f ( b , ω ) ∣ 2 d b d ω \rho(\psi, f) = \frac{\left\| T_f(b,\phi'(b)) \right\|_2^2}{\iint |T_f(b,\omega)|^2 db d\omega} ρ(ψ,f)=∬∣Tf(b,ω)∣2dbdω∥Tf(b,ϕ′(b))∥22
可以证明 ρ ( ψ B u m p , f ) > ρ ( ψ M o r l e t , f ) \rho(\psi_{Bump}, f) > \rho(\psi_{Morlet}, f) ρ(ψBump,f)>ρ(ψMorlet,f)。
重构公式的严格推导
同步挤压变换的一个重要特性是其可逆性。对于任意信号 f ( t ) f(t) f(t),可以通过以下步骤严格推导重构公式:
首先考虑小波变换的反变换公式:
f ( t ) = C ψ − 1 ∫ 0 ∞ ∫ − ∞ ∞ W f ( a , b ) 1 a ψ ( t − b a ) d a d b a 2 f(t) = \mathcal{C}_{\psi}^{-1} \int_{0}^{\infty} \int_{-\infty}^{\infty} W_f(a,b) \frac{1}{\sqrt{a}} \psi\left(\frac{t-b}{a}\right) \frac{da db}{a^2} f(t)=Cψ−1∫0∞∫−∞∞Wf(a,b)a1ψ(at−b)a2dadb
其中 C ψ = ∫ 0 ∞ ∣ ψ ^ ( ξ ) ∣ 2 ξ d ξ \mathcal{C}_{\psi} = \int_{0}^{\infty} \frac{|\hat{\psi}(\xi)|^2}{\xi} d\xi Cψ=∫0∞ξ∣ψ^(ξ)∣2dξ。
将同步挤压变换的定义代入,得到:
KaTeX parse error: No such environment: align at position 7: \begin{̲a̲l̲i̲g̲n̲}̲ f(t) &= \mathc…
对于理想的单分量信号,如果 ω f ( a , b ) ≈ ϕ ′ ( b ) \omega_f(a,b) \approx \phi'(b) ωf(a,b)≈ϕ′(b) 对所有 a a a 近似成立,则内层积分可以近似为:
∫ 0 ∞ δ ( ω − ϕ ′ ( b ) ) ⋅ a − 1 ψ ( t − b a ) d a ≈ δ ( ω − ϕ ′ ( b ) ) ⋅ ∫ 0 ∞ a − 1 ψ ( t − b a ) d a \int_{0}^{\infty} \delta(\omega - \phi'(b)) \cdot a^{-1} \psi\left(\frac{t-b}{a}\right) da \approx \delta(\omega - \phi'(b)) \cdot \int_{0}^{\infty} a^{-1} \psi\left(\frac{t-b}{a}\right) da ∫0∞δ(ω−ϕ′(b))⋅a−1ψ(at−b)da≈δ(ω−ϕ′(b))⋅∫0∞a−1ψ(at−b)da
应用傅里叶变换的性质,可以证明:
∫ 0 ∞ a − 1 ψ ( t − b a ) d a = 1 2 π ∫ − ∞ ∞ e i ω ( t − b ) ψ ^ ∗ ( ω ) ω d ω \int_{0}^{\infty} a^{-1} \psi\left(\frac{t-b}{a}\right) da = \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{i\omega(t-b)} \frac{\hat{\psi}^*(\omega)}{\omega} d\omega ∫0∞a−1ψ(at−b)da=2π1∫−∞∞eiω(t−b)ωψ^∗(ω)dω
因此,重构公式简化为:
f ( t ) = C ψ − 1 ⋅ Re { ∫ − ∞ ∞ ∫ 0 ∞ T f ( b , ω ) ⋅ e i ω ( t − b ) d ω d b } f(t) = \mathcal{C}_{\psi}^{-1} \cdot \text{Re} \left\{ \int_{-\infty}^{\infty} \int_{0}^{\infty} T_f(b,\omega) \cdot e^{i\omega(t-b)} d\omega db \right\} f(t)=Cψ−1⋅Re{∫−∞∞∫0∞Tf(b,ω)⋅eiω(t−b)dωdb}
这就是同步挤压变换的标准重构公式。