非平稳信号的广义小波分析及其工程应用报告

非平稳信号的广义小波分析及其工程应用报告徐文豪1.窗口傅里叶变换1.1原理分析众所周知，傅里叶变换可以获取信号的全局频谱，但很多时候，我们需要的是信号的瞬时频谱，如机械故障检测、地震信号瞬时属性提取等。

为了达到这个目的，一种很直观的思路就是截取信号的一小段做傅里叶变换，并将得到的频谱作为该小段中心位置的频谱。

这种变换称之为窗口傅里叶变换（WFT），其示意图如下：从数学上描述WFT应为：对平方可积信号()s x和平方可积窗函数()g x，定义窗口傅里叶变换如下：+(,)()()i xS t s x g x t e dxωω∞--∞=-⎰(1) 值得强调的是，式（1）中窗函数()g x的支撑大小对时频谱的有效性有着很大的影响，当窗宽过大或过小时都会使时频谱出现较严重的假象。

作者在实际编程时发现，当信号的长度为L，取一个较小的正值ε作为值零阈值，则取窗宽点数半径为一般能达到较理想的效果。

将()()s x g x t-视为整体，假设2g≠并通过逆傅里叶变换可得逆窗口傅里叶变换如下：221()(,)()2i xs x S t g x t e d dtgωωωπ+∞+∞-∞-∞=-⎰⎰(2) 式（2）对推导WFT值域(,)S tω的性质有着很重要的作用，如重建核方程等。

但由于使用双重积分，其在实现时稍显繁琐，殷勤业的《时频分析及其在工程中的应用》讲义中给出了一个更简单的重构公式如下：*1()(,)2(0)i ts t S t e dgωωωπ+∞-∞=⎰(3) 其推导如下：*()**11(,)()()22()()()()()i i S t e d s g t e d d s g t t d s g t ωμωμτωωτττωππττδμτμμ+∞+∞+∞--∞-∞-∞+∞-∞=-=--=-⎰⎰⎰⎰1.2 程序实现附录1给出了与式（1）对应的WFT 程序，附录2给出了与式（3）对应的IWFT 程序,对如下的四段调频信号：sin(400),[0,256)sin(800),[256,512)()sin(200),[512,768)sin(600),[768,1024)t t ms t t ms s t t t ms t t msππππ∈⎧⎪∈⎪=⎨∈⎪⎪∈⎩ (4)使用WFT 程序和IWFT 程序，得到()s t 的时频谱图和误差图如下：图2 四段跳频信号WFT 时频谱图及重构误差图2. 连续小波变换2.1 原理分析窗口傅里叶变换的缺点在于窗宽是固定的，因而时频谱容易出现假频现象， 如下图所示：图3 WFT 缺点示意图图3中，对快变信号采用大窗或对慢变信号采用小窗都必然会造成假频。

解决这一问题的思路是使窗宽可变，这就是连续小波变换(CWT)的思路。

在CWT 中我们一般把窗函数称为小波函数，其支撑如下所示：图4 小波函数的支撑从数学上描述CWT 应为：对平方可积信号()s x 和平方可积小波()x ψ，定义连续小波变换如下：1(,)()(),0p x bS a b s x dx p a aψ+∞-∞-=>⎰ (5) 其中系数1p a的引入是为了小波在平移伸缩之后1p 范数不变，即111()()pppx bx aaψψ-=(6)对式（5）中的连续小波变换，常见的逆变换公式如下：,23211()(,)()()a b p s t S a b t dadb a d ψψωωω+∞+∞--∞-∞+∞-∞=⎰⎰⎰(7)其中,1()()a b pt bt aaψψ-=称为时间尺度原子。

与WFT 类似，式（7）对推导CWT 值域(,)S a b 的性质有着很重要的作用，如重建核方程等，但其一般只能通过二重数值积分进行计算，速度慢且误差大。

Holschneider 在《Wavelets: An Analysis Tool 》中给出了当范数变量1p =时的一个简单重构公式：10(,)()S a t s t C da a ψ+∞-=⎰(8) 0ˆ()C d ψψωωω+∞=⎰ (9)如果能求出C ψ的值，则可利用数值积分计算式（8）的值。

遗憾的是式（9）是不一定收敛的，例如对常用的Morlet 小波(2()/2ˆωσψω--（)，式（9）在0处发散。

高老师的讲义中给出了当范数1/2p =时的一个简单逆变换公式，但其有效性我还没有搞懂，其计算过程如下：1Re (,)1()Mj S a t s t Cδ==(10)1Re (,0)Mj S a C δ==(11)式（11）中(,0)j S a δ表示单位脉冲信号[]n δ在0时刻尺度为j a 处CWT 的值。

2.2 程序实现附录3给出了与式（8）对应的CWT 程序，附录4给出了与式（10）对应的ICWT 程序，对地球物理中常用的Ricker 信号(取主频为50Hz)，取Gauss 小波(22()/(2)ˆ()2s ωμψω--=，参考刘乃豪师兄提供的资料)使用CWT 程序和ICWT 程序，并利用频率尺度转换关系，得到其时频谱和重构误差如下：图5 50Hz 主频Ricker 信号CWT 时频谱图及重构误差图从图5可以看出，对Ricker 信号，在选取与其相配的小波函数之后，不用考虑窗宽，连续小波变换就能得到很好的时频谱和重构精度。

CWT 具有自适应性的优点可以从其Heisenberg 盒得到进一步验证：图6 连续小波变换的Heisenberg 盒从图6可以看出，CWT 既具有显微镜的功能，又具有望远镜的功能。

即当尺度a较小时，有可能捕捉小范围的信号变化，而当尺度a 较大时，有可能捕捉大范围的信号变化。

然而，福兮祸之所倚，祸兮福之所伏，图6的Heisenberg 盒也反映出了CWT 的一大缺点，即低频可能看不见(频率分辨率过高，而离散尺度有限)，高频可能看不清（频率分辨率过低），这种现象可以从式（4）的四段调频信号反应出来，同样取高斯小波，四段调频信号的CWT 时频谱和重构误差如下：图7四段跳频信号CWT 时频谱图及重构误差图从图7可以很明显的看出，在CWT 时频谱中，信号频率分辨率随着频率的增大而降低。

这里需要指出的是，缓解这一问题的一个有效方法是二进制采样，其能够在低频部分采相对较多的点并在高频部分采相对较少的点。

令一方面，图7也反应出图5的高重构精度并不适合所有信号，其原因可能有两个，一是Ricker 信号的变化较平缓而四段调频信号的变化较剧烈，二是高斯小波可能与四段调频信号不太匹配。

3. 最优对偶标架变换3.1 原理分析标架是指Hilbert 空间H 中的一个完备序列{}ih ，其满足对任意x H ∈，存在,0A B >，使得222,iiA x x hB x ≤≤∑(12)其中x 为x 的由内积诱导的范数。

A 和B 分别称为标架下界和标架上界，特别地，当A B C ==时，称这个标架为紧标架，并称C 为紧标架的冗余度。

已经证明对任意标架{}ih，存在对偶标架{}iγ ，使得对于任意x H ∈，有,i iix x h γ=∑ (13)而当{}ih 为紧标架时，Daubechies 已证明其对偶标架为1i h C ⎧⎫⎨⎬⎩⎭。

对采样间隔为dx ，长度为L 的周期离散信号[]sk (将有限离散信号周期化可以得到更好的边界重构精度)，定义其对偶标架变换及逆变换如下：001,,,0[]((1))L m nm n m n k C s k k dx γ-==-∑ (14)0011,,,00[]((1))M N m n m n m n m n s k C h k dx --===-∑∑(15)其中M 为时间采样点数，N 为波数采样点数，{}00,,()m n m n x γ 为由分析函数()x γ 的平移调制生成的标架，{}00,,()m n m n h x 为由综合函数()h x 的平移调制生成的{}00,,()m n m n x γ的对偶标架，00,,()m n m n h x 的定义如下：000()0,,()(())i n n dN d x m n m n h x h x m m dM dx e ξ+⨯=-+⨯ (16)其中00,m n 分别为时间和波数初始采样点（引入这两个常量是为了在标架相空间中得到更多的空间波数信息），,dM dN 分别为时间和波数采样间隔点数且满足dM M dN N L ⨯=⨯=，2d L dxπξ=⨯为波数离散采样间隔。

通过给定综合函数()hx 的离散采样序列[]h k ，由式（14）和式（15）可推出分析函数()x γ 的离散采样序列[]k γ 所满足的方程组：1*0[][][][]L pk Mk M h k qN W k p q N γδδ--∆=∆+=∑ (17)其中0p M ≤<∆，0q N ≤<∆。

式（17）所对应方程组的系数矩阵H 是dM dN ⨯行L 列的，故若dM dN L⨯=且系数矩阵满秩时可以解出唯一的*γ，然而Daubechies 已经证明这种情况下重构公式是不稳定的，一般称这种情形下的采样为临界采样。

钱世锷在著作《Joint Time –Frequency Analysis 》中提出取dM dN L ⨯<，并在方程组的解空间中寻找与[]hk 最接近的解*optγ ，令方程组的值向量为μ，并假设21h = 则可构造对应的优化问题如下：*21*02[]min []L opt H k k h k γμγγγ-===-∑(18)对式（18）进行简单推导得**1min 1opt H M N γμγγ=⎧⎫∆⎪⎪=-⎨⎬⎪⎪⎩⎭ (19)式（19）意味着要在*H γμ= 的解空间中寻找模最小的解。

由广义逆矩阵理论[21]，可给出所需解的表达式如下：*opt H γμ+=(20)其中H +表示矩阵H 的Moore-Penrose 逆。

对于计算H +，钱世锷书中是通过判断H 是否行满秩并通过SVD 分解将行不满秩的情形转换为行满秩情形，这样不仅计算繁琐还会影响计算精度。

实际上有直接计算矩阵MP 逆的快速算法，如Greville 递推法，而在matlab 中可直接调用库函数pinv 实现。

然而，式（20）只能处理规模相对较小的矩阵，且由于使用复数运算会影响计算精度。

钱世锷书中在以上的基础上通过将大系数矩阵分解为多个小系数矩阵，将复数运算化为实数运算，给出了求解方程组的快速算法：1*01[][][]M l h k l M qN k l M q N γδ-=+∆++∆=∑ (21)其中021q N ≤<∆-，0k M ≤<∆。

这样大复系数矩阵就变成了M ∆个小实系数矩阵，再利用类似推导式（20）的过程即可得到所需解。

对采样间隔为dx 的信号，一般取如下归一化的高斯函数作为综合函数：2()1/422()(),xdx h x e dM N ααπαπ-==⨯ (22)其中α的取值是钱世锷书中给出的[]h k 和[]opt k γ 间取得最小误差的条件。