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

1非平稳信号的广义小波分析及其工程应用报告徐文豪1.窗口傅里叶变换1.1原理分析众所周知，傅里叶变换可以获取信号的全局频谱，但很多时候，我们需要的是信号的瞬时频谱，如机械故障检测、地震信号瞬时属性提取等。为了达到这个目的，一种很直观的思路就是截取信号的一小段做傅里叶变换，并将得到的频谱作为该小段中心位置的频谱。这种变换称之为窗口傅里叶变换（WFT），其示意图如下：图1窗口福利叶变换示意图从数学上描述WFT应为：对平方可积信号()sx和平方可积窗函数()gx，定义窗口傅里叶变换如下：+(,)()()ixStsxgxtedx(1)值得强调的是，式（1）中窗函数()gx的支撑大小对时频谱的有效性有着很大的影响，当窗宽过大或过小时都会使时频谱出现较严重的假象。作者在实际编程时发现，当信号的长度为L，取一个较小的正值作为值零阈值，则取窗宽点数半径为L一般能达到较理想的效果。将()()sxgxt视为整体，假设20g并通过逆傅里叶变换可得逆窗口傅里叶变换如下：221()(,)()2ixsxStgxteddtg(2)式（2）对推导WFT值域(,)St的性质有着很重要的作用，如重建核方程等。但由于使用双重积分，其在实现时稍显繁琐，殷勤业的《时频分析及其在工程中的应用》讲义中给出了一个更简单的重构公式如下：*1()(,)2(0)itstStedg(3)其推导如下：00.511.500.511.522.5xA窗口傅里叶变换示意图2*()**11(,)()()22()()()()()iiStedsgteddsgttdsgt1.2程序实现附录1给出了与式（1）对应的WFT程序，附录2给出了与式（3）对应的IWFT程序,对如下的四段调频信号：sin(400),[0,256)sin(800),[256,512)()sin(200),[512,768)sin(600),[768,1024)ttmsttmsstttmsttms(4)使用WFT程序和IWFT程序，得到()st的时频谱图和误差图如下：图2四段跳频信号WFT时频谱图及重构误差图2.连续小波变换2.1原理分析窗口傅里叶变换的缺点在于窗宽是固定的，因而时频谱容易出现假频现象，如下图所示：图3WFT缺点示意图3图3中，对快变信号采用大窗或对慢变信号采用小窗都必然会造成假频。解决这一问题的思路是使窗宽可变，这就是连续小波变换(CWT)的思路。在CWT中我们一般把窗函数称为小波函数，其支撑如下所示：图4小波函数的支撑从数学上描述CWT应为：对平方可积信号()sx和平方可积小波()x，定义连续小波变换如下：1(,)()(),0pxbSabsxdxpaa(5)其中系数1pa的引入是为了小波在平移伸缩之后1p范数不变，即111()()pppxbxaa(6)对式（5）中的连续小波变换，常见的逆变换公式如下：,23211()(,)()()abpstSabtdadbad(7)其中,1()()abptbtaa称为时间尺度原子。与WFT类似，式（7）对推导CWT值域(,)Sab的性质有着很重要的作用，如重建核方程等，但其一般只能通过二重数值积分进行计算，速度慢且误差大。Holschneider在《Wavelets:AnAnalysisTool》中给出了当范数变量1p时的一个简单重构公式：10(,)()SatstCdaa(8)0ˆ()Cd(9)如果能求出C的值，则可利用数值积分计算式（8）的值。遗憾的是式（9）是不一定收敛的，例如对常用的Morlet小波(2()/2ˆ=2e（）)，式（9）在0处发散。高老师的讲义中给出了当范数1/2p时的一个简单逆变换公式，但其有效性我还没有搞懂，其计算过程如下：1Re(,)1()MjjjSatstCa(10)41Re(,0)MjjjSaCa(11)式（11）中(,0)jSa表示单位脉冲信号[]n在0时刻尺度为ja处CWT的值。2.2程序实现附录3给出了与式（8）对应的CWT程序，附录4给出了与式（10）对应的ICWT程序，对地球物理中常用的Ricker信号(取主频为50Hz)，取Gauss小波(22()/(2)ˆ()2s，参考刘乃豪师兄提供的资料)使用CWT程序和ICWT程序，并利用频率尺度转换关系，得到其时频谱和重构误差如下：图550Hz主频Ricker信号CWT时频谱图及重构误差图从图5可以看出，对Ricker信号，在选取与其相配的小波函数之后，不用考虑窗宽，连续小波变换就能得到很好的时频谱和重构精度。CWT具有自适应性的优点可以从其Heisenberg盒得到进一步验证：图6连续小波变换的Heisenberg盒从图6可以看出，CWT既具有显微镜的功能，又具有望远镜的功能。即当尺度a较小时，有可能捕捉小范围的信号变化，而当尺度a较大时，有可能捕捉大范围的信号变化。5然而，福兮祸之所倚，祸兮福之所伏，图6的Heisenberg盒也反映出了CWT的一大缺点，即低频可能看不见(频率分辨率过高，而离散尺度有限)，高频可能看不清（频率分辨率过低），这种现象可以从式（4）的四段调频信号反应出来，同样取高斯小波，四段调频信号的CWT时频谱和重构误差如下：图7四段跳频信号CWT时频谱图及重构误差图从图7可以很明显的看出，在CWT时频谱中，信号频率分辨率随着频率的增大而降低。这里需要指出的是，缓解这一问题的一个有效方法是二进制采样，其能够在低频部分采相对较多的点并在高频部分采相对较少的点。令一方面，图7也反应出图5的高重构精度并不适合所有信号，其原因可能有两个，一是Ricker信号的变化较平缓而四段调频信号的变化较剧烈，二是高斯小波可能与四段调频信号不太匹配。3.最优对偶标架变换3.1原理分析标架是指Hilbert空间H中的一个完备序列ih，其满足对任意xH，存在,0AB，使得222,iiAxxhBx(12)其中x为x的由内积诱导的范数。A和B分别称为标架下界和标架上界，特别地，当ABC时，称这个标架为紧标架，并称C为紧标架的冗余度。已经证明对任意标架ih，存在对偶标架i，使得对于任意xH，有,iiixxh(13)而当ih为紧标架时，Daubechies已证明其对偶标架为1ihC。对采样间隔为dx，长度为L的周期离散信号[]sk(将有限离散信号周期化可6以得到更好的边界重构精度)，定义其对偶标架变换及逆变换如下：001,,,0[]((1))LmnmnmnkCskkdx(14)0011,,,00[]((1))MNmnmnmnmnskChkdx(15)其中M为时间采样点数，N为波数采样点数，00,,()mnmnx为由分析函数()x的平移调制生成的标架，00,,()mnmnhx为由综合函数()hx的平移调制生成的00,,()mnmnx的对偶标架，00,,()mnmnhx的定义如下：000()0,,()(())inndNdxmnmnhxhxmmdMdxe(16)其中00,mn分别为时间和波数初始采样点（引入这两个常量是为了在标架相空间中得到更多的空间波数信息），,dMdN分别为时间和波数采样间隔点数且满足dMMdNNL，2dLdx为波数离散采样间隔。通过给定综合函数()hx的离散采样序列[]hk，由式（14）和式（15）可推出分析函数()x的离散采样序列[]k所满足的方程组：1*0[][][][]LpkMkMhkqNWkpqN(17)其中0pM，0qN。式（17）所对应方程组的系数矩阵H是dMdN行L列的，故若dMdNL且系数矩阵满秩时可以解出唯一的*，然而Daubechies已经证明这种情况下重构公式是不稳定的，一般称这种情形下的采样为临界采样。钱世锷在著作《JointTime–FrequencyAnalysis》中提出取dMdNL，并在方程组的解空间中寻找与[]hk最接近的解*opt，令方程组的值向量为，并假设21h则可构造对应的优化问题如下：*21*02[]min[]LoptHkkhk(18)对式（18）进行简单推导得**1min1optHMN(19)式（19）意味着要在*H的解空间中寻找模最小的解。由广义逆矩阵理论[21]，可给出所需解的表达式如下：7*optH(20)其中H表示矩阵H的Moore-Penrose逆。对于计算H，钱世锷书中是通过判断H是否行满秩并通过SVD分解将行不满秩的情形转换为行满秩情形，这样不仅计算繁琐还会影响计算精度。实际上有直接计算矩阵MP逆的快速算法，如Greville递推法，而在matlab中可直接调用库函数pinv实现。然而，式（20）只能处理规模相对较小的矩阵，且由于使用复数运算会影响计算精度。钱世锷书中在以上的基础上通过将大系数矩阵分解为多个小系数矩阵，将复数运算化为实数运算，给出了求解方程组的快速算法：1*01[][][]MlhklMqNklMqN(21)其中021qN，0kM。这样大复系数矩阵就变成了M个小实系数矩阵，再利用类似推导式（20）的过程即可得到所需解。对采样间隔为dx的信号，一般取如下归一化的高斯函数作为综合函数：2()1/4212()(),xdxhxedMNdx(22)其中的取值是钱世锷书中给出的[]hk和[]optk间取得最小误差的条件。必须注意到的是，[]hx的支撑hL不应也不必等于L，否则当L很大时，求解对偶标架将相对比较耗时。由殷勤业讲义3-4节的内容和作者编程时的经验，一般hL取N时可以达到很好的效果，且若让[]hk在0:1hL的值为()hx在(/2:/21)hhLLdx的值并在其它点处为0，则可以利用式（21）只计算长度为hL的对偶函数向量。对偶标架变换的重构误差精度与时频采样点数MN与信号长度L之比有关，这可以被称为对偶标架采样冗余度C，定义如下：MNLCLdMdN(23)3.2程序实现附录5给出了与式（21）对应的快速计算最优对偶标架的程序，附录6给出了与式（14）对应的最优对偶标架变换的程序，附录7给出了与式（15）对应的逆最优对偶标架变换程序。同样对式（4）的四段跳频信号，取采样冗余度2C，使用最优对偶标架变换程序和逆最优对偶标架变换程序得到时频谱和重构误差如下：8图8四段跳频信号最优对偶标架时频谱图及重构误差图从图8可以看出，最优对偶标架变换只用了WFT存储量的1/512，就得到了有效的谱图和非常高的重构精度，这暗示了最优对偶标架变换在处理较大数据的潜力。实际上，对长度为131072的音乐片段(许嵩《忧伤还是快乐》截取，音乐采样频率为44100Hz)，使用冗余度64C的最优对偶标架变换得到时频图和重构误差如下：图9音乐片段最优对偶标架变换时频谱图及重构误差图考虑到该段音乐为钢琴片段，图9中的时频谱是比较合理的，且程序运行时间只有26.65s(其中有限支撑最优对偶标架的计算时间仅为0.004s)，再加上非常高的重构精度，完全可以断言最优对偶标架变换具有处理实际数据的能力。最优对偶标架变换另一个比较好的特点是，其自动选择了一个比较好的窗宽(见式（22）)。首先对变化较慢的频率为7.8125Hz的正弦信号，使用冗余度64C的最优对偶标架变换得到时频谱图和重构误差图如下：9图10慢变信号最优对偶标架变换时频谱图及重构误差图其次，对变化剧烈的单位脉冲信号，仍使用冗余度64C的最优对偶标架变换，得到其时频谱