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

1、非平稳信号的广义小波分析及其工程应用报告徐文豪1. 窗口傅里叶变换1.1 原理分析众所周知，傅里叶变换可以获取信号的全局频谱，但很多时候，我们需要的是信号的瞬时频谱，如机械故障检测、地震信号瞬时属性提取等。为了达到这个目的，一种很直观的思路就是截取信号的一小段做傅里叶变换，并将得到的频谱作为该小段中心位置的频谱。这种变换称之为窗口傅里叶变换（WFT），其示意图如下：图1 窗口福利叶变换示意图从数学上描述WFT应为：对平方可积信号和平方可积窗函数，定义 窗口傅里叶变换如下：值得强调的是，式（1）中窗函数的支撑大小对时频谱的有效性有着很大的影响，当窗宽过大或过小时都会使时频谱出现较严重的假象。作者

2、在实际编程时发现，当信号的长度为，取一个较小的正值作为值零阈值，则取窗宽点数半径为一般能达到较理想的效果。将视为整体，假设并通过逆傅里叶变换可得逆窗口傅里叶变换如下：式（2）对推导WFT值域的性质有着很重要的作用，如重建核方程等。但由于使用双重积分，其在实现时稍显繁琐，殷勤业的时频分析及其在工程中的应用讲义中给出了一个更简单的重构公式如下：其推导如下：1.2 程序实现附录1给出了与式（1）对应的WFT程序，附录2给出了与式（3）对应的IWFT程序,对如下的四段调频信号：使用WFT程序和IWFT程序，得到的时频谱图和误差图如下：图2 四段跳频信号WFT时频谱图及重构误差图2. 连续小波变换2.1

3、 原理分析窗口傅里叶变换的缺点在于窗宽是固定的，因而时频谱容易出现假频现象，如下图所示：图3 WFT缺点示意图图3中，对快变信号采用大窗或对慢变信号采用小窗都必然会造成假频。解决这一问题的思路是使窗宽可变，这就是连续小波变换(CWT)的思路。在CWT中我们一般把窗函数称为小波函数，其支撑如下所示：图4 小波函数的支撑从数学上描述CWT应为：对平方可积信号和平方可积小波，定义连续小波变换如下：其中系数的引入是为了小波在平移伸缩之后范数不变，即对式（5）中的连续小波变换，常见的逆变换公式如下：其中称为时间尺度原子。与WFT类似，式（7）对推导CWT值域的性质有着很重要的作用，如重建核方程等，但其一

4、般只能通过二重数值积分进行计算，速度慢且误差大。Holschneider在Wavelets: An Analysis Tool中给出了当范数变量时的一个简单重构公式：如果能求出的值，则可利用数值积分计算式（8）的值。遗憾的是式（9）是不一定收敛的，例如对常用的Morlet小波()，式（9）在0处发散。高老师的讲义中给出了当范数时的一个简单逆变换公式，但其有效性我还没有搞懂，其计算过程如下：式（11）中表示单位脉冲信号在0时刻尺度为处CWT的值。2.2 程序实现附录3给出了与式（8）对应的CWT程序，附录4给出了与式（10）对应的ICWT程序，对地球物理中常用的Ricker信号(取主频为50Hz

5、)，取Gauss小波(，参考刘乃豪师兄提供的资料)使用CWT程序和ICWT程序，并利用频率尺度转换关系，得到其时频谱和重构误差如下：图5 50Hz主频Ricker信号CWT时频谱图及重构误差图从图5可以看出，对Ricker信号，在选取与其相配的小波函数之后，不用考虑窗宽，连续小波变换就能得到很好的时频谱和重构精度。CWT具有自适应性的优点可以从其Heisenberg盒得到进一步验证：图6 连续小波变换的Heisenberg盒从图6可以看出，CWT既具有显微镜的功能，又具有望远镜的功能。即当尺度较小时，有可能捕捉小范围的信号变化，而当尺度较大时，有可能捕捉大范围的信号变化。然而，福兮祸之所倚，祸

6、兮福之所伏，图6的Heisenberg盒也反映出了CWT的一大缺点，即低频可能看不见(频率分辨率过高，而离散尺度有限)，高频可能看不清（频率分辨率过低），这种现象可以从式（4）的四段调频信号反应出来，同样取高斯小波，四段调频信号的CWT时频谱和重构误差如下：图7四段跳频信号CWT时频谱图及重构误差图从图7可以很明显的看出，在CWT时频谱中，信号频率分辨率随着频率的增大而降低。这里需要指出的是，缓解这一问题的一个有效方法是二进制采样，其能够在低频部分采相对较多的点并在高频部分采相对较少的点。令一方面，图7也反应出图5的高重构精度并不适合所有信号，其原因可能有两个，一是Ricker信号的变化较平缓

7、而四段调频信号的变化较剧烈，二是高斯小波可能与四段调频信号不太匹配。3. 最优对偶标架变换3.1 原理分析标架是指Hilbert空间中的一个完备序列，其满足对任意，存在，使得其中为的由内积诱导的范数。和分别称为标架下界和标架上界，特别地，当时，称这个标架为紧标架，并称为紧标架的冗余度。已经证明对任意标架，存在对偶标架，使得对于任意，有而当为紧标架时，Daubechies已证明其对偶标架为。对采样间隔为，长度为的周期离散信号(将有限离散信号周期化可以得到更好的边界重构精度)，定义其对偶标架变换及逆变换如下：其中为时间采样点数，为波数采样点数，为由分析函数的平移调制生成的标架，为由综合函数的平移调

8、制生成的的对偶标架，的定义如下：其中分别为时间和波数初始采样点（引入这两个常量是为了在标架相空间中得到更多的空间波数信息），分别为时间和波数采样间隔点数且满足，为波数离散采样间隔。通过给定综合函数的离散采样序列，由式（14）和式（15）可推出分析函数的离散采样序列所满足的方程组：其中，。式（17）所对应方程组的系数矩阵是行列的，故若且系数矩阵满秩时可以解出唯一的，然而Daubechies已经证明这种情况下重构公式是不稳定的，一般称这种情形下的采样为临界采样。钱世锷在著作Joint Time Frequency Analysis中提出取， 并在方程组的解空间中寻找与最接近的解，令方程组的值向量为

9、，并假设则可构造对应的优化问题如下：对式（18）进行简单推导得式（19）意味着要在的解空间中寻找模最小的解。由广义逆矩阵理论，可给出所需解的表达式如下：其中表示矩阵的Moore-Penrose逆。对于计算，钱世锷书中是通过判断是否行满秩并通过SVD分解将行不满秩的情形转换为行满秩情形，这样不仅计算繁琐还会影响计算精度。实际上有直接计算矩阵MP逆的快速算法，如Greville递推法，而在matlab中可直接调用库函数pinv实现。然而，式（20）只能处理规模相对较小的矩阵，且由于使用复数运算会影响计算精度。钱世锷书中在以上的基础上通过将大系数矩阵分解为多个小系数矩阵，将复数运算化为实数运算，给出

10、了求解方程组的快速算法：其中，。这样大复系数矩阵就变成了个小实系数矩阵，再利用类似推导式（20）的过程即可得到所需解。对采样间隔为的信号，一般取如下归一化的高斯函数作为综合函数：其中的取值是钱世锷书中给出的和间取得最小误差的条件。必须注意到的是，的支撑不应也不必等于，否则当很大时，求解对偶标架将相对比较耗时。由殷勤业讲义3-4节的内容和作者编程时的经验，一般取时可以达到很好的效果，且若让在的值为在的值并在其它点处为，则可以利用式（21）只计算长度为的对偶函数向量。对偶标架变换的重构误差精度与时频采样点数与信号长度之比有关，这可以被称为对偶标架采样冗余度，定义如下：3.2 程序实现附录5给出了与

11、式（21）对应的快速计算最优对偶标架的程序，附录6给出了与式（14）对应的最优对偶标架变换的程序，附录7给出了与式（15）对应的逆最优对偶标架变换程序。同样对式（4）的四段跳频信号，取采样冗余度，使用最优对偶标架变换程序和逆最优对偶标架变换程序得到时频谱和重构误差如下：图8四段跳频信号最优对偶标架时频谱图及重构误差图从图8可以看出，最优对偶标架变换只用了WFT存储量的，就得到了有效的谱图和非常高的重构精度，这暗示了最优对偶标架变换在处理较大数据的潜力。实际上，对长度为131072的音乐片段(许嵩忧伤还是快乐截取，音乐采样频率为44100Hz)，使用冗余度的最优对偶标架变换得到时频图和重构误差如

12、下：图9 音乐片段最优对偶标架变换时频谱图及重构误差图考虑到该段音乐为钢琴片段，图9中的时频谱是比较合理的，且程序运行时间只有26.65s(其中有限支撑最优对偶标架的计算时间仅为0.004s)，再加上非常高的重构精度，完全可以断言最优对偶标架变换具有处理实际数据的能力。最优对偶标架变换另一个比较好的特点是，其自动选择了一个比较好的窗宽(见式（22）)。首先对变化较慢的频率为7.8125Hz的正弦信号，使用冗余度的最优对偶标架变换得到时频谱图和重构误差图如下：图10 慢变信号最优对偶标架变换时频谱图及重构误差图其次，对变化剧烈的单位脉冲信号，仍使用冗余度的最优对偶标架变换，得到其时频谱图和重构误

13、差图如下：图11 快变信号最优对偶标架变换时频谱图及重构误差图从图10和图11可以看出，最优对偶标架自动选择的窗宽可以同时较好地处理慢变信号和快变信号，这给使用带来了较大的便利。最后需要说明的是冗余度和重构精度的关系，事实上求解对偶标架的方程组式（17）的系数矩阵规模为，而。因此越大，对偶标架满足的方程组就越欠定，所得的最优对偶标架就接近紧标架，因而有较高的重构精度(标架正反变换都线性有界(连续)的充要条件是其为紧标架)。而当达到一定程度，使得解空间中已经有很好的对偶标架时，重构精度就应当保持稳定了。上述讨论可以从最优对偶标架变换重构精度随冗余度的变化验证，对长度为1024的100Hz正弦信号

14、，结果如下表所示：表1 最优对偶标架变换重构精度随冗余度变化表13e-321e-1541e-1588e-16168e-16从表1中可以看出，当冗余度时，解空间中应该就已经包含了很接近紧标架的对偶标架了。4. Grossman对偶标架变换4.1 原理分析不同于构造方程组，Grossman给出了一个极其简单的构造性求解对偶标架的方法。为便于该方法的叙述，称和为常数1的函数集为单位1的分割，定义平移算子为，定义调制算子为，则有：定理1：设是对单位1的分割，任意选择，如果令，那么定义了一个Gabor分析标架，相应的综合标架可为，这里定义为。Gross对偶标架变换实现的关键是单位1的分割，但实际上这可以

15、直接通过将一个窗函数离散向量求和归一化实现。4.2 程序实现附录8给出了Grossman对偶标架变换的实现(其需要调用附录6和附录7中的程序)，对主频为50Hz的Ricker信号，取采样冗余度，使用Grossman对偶标架变换程序得到其时频谱图和重构误差如下：图12 50Hz主频Ricker信号Grossman对偶标架变换时频图和重构误差图相对于最优对偶标架变换，Grossman对偶标架变换节省了求解最优对偶标架的时间，然而，当采样冗余度较小时，Grossman对偶标架变换的重构误差精度较低，对长度为1024的100Hz正弦信号，结果如下表所示：表2 Grossman对偶标架变换重构精度随冗余

16、度变化表16e-123e-146e-283e-3165e-6322e-11645e-161285e-162565e-16对比表1和表2并结合计算量可得出结论，当采样冗余度时，应选择最优对偶标架变换，当采样冗余度时，应选择Grossman对偶标架变换。5. 附录附录1：WFT的matlab程序%窗口傅里叶变换程序%s为长度为N的离散信号，w为窗函数指针%S为s的加窗变换,其分量S(i,j)(jlength(hTildeVec) error(采样间隔过大，无法恢复信号！);end%计算所需常量L=length(hTildeVec);N=L/dN; %计算对偶标架频率采样点数M=L/dM; %计算对

17、偶标架时间采样点数%初始化变量gammaTildeVec=zeros(1,L);hKMat=zeros(dN,M);muKBarVec=zeros(dN,1);muKBarVec(1)=1/N;%循环k计算最优对偶标架函数for k=0:dM-1 %构造dN行M列的系数矩阵 for q=0:dN-1 for l=0:M-1 hKMat(q+1,l+1)=hTildeVec(mod(k+l*dM+q*N,L)+1); end end %利用MP逆求解能量最小的解 gammaKBarVec=(pinv(hKMat)*muKBarVec); %注意这里包含共轭了 gammaTildeVec(k+(0

18、:M-1)*dM+1)=gammaKBarVec;endgammaTildeVec=real(gammaTildeVec);end附录6：对偶标架变换的matlab程序%对偶标架变换程序%sTilde为长度为L的离散信号变量命名基于殷勤业讲义%gammaTildeVec为最优对偶标架函数离散向量%dM为对偶标架时间采样间隔点数%dN为对偶标架频率采样间隔点数%m0为空间采样起点，要求其值为0,dM-1间的整数%n0为波数采样起点，要求其值为-L/2,-L/2+dN-1间的整数%cTildeMat为对偶标架系数矩阵function cTildeMat=dualFrameTransform(sTi

19、lde,gammaTildeVec,dM,dN,m0,n0)%判断输入的有效性if length(sTilde)=length(gammaTildeVec) error(要求对偶标架离散向量和信号维数相同！);endif mod(length(sTilde),dM)=0 error(dM必须为整数且能整除信号长度);endif mod(length(sTilde),dN)=0 error(dN必须为整数且能整除信号长度);endif rem(m0,1)=0 | m0dM-1 error(m0必须为0,dM-1-1间的整数);endif rem(n0,1)=0 | n0-length(sTild

20、e)/2+dN-1 error(n0必须为-L/2,-L/2+dN-1间的整数);end%计算所需常量L=length(sTilde);M=L/dM; %计算对偶标架空间采样点数N=L/dN; %计算对偶标架波数采样点数%计算对偶标架系数矩阵cTildeMatcTildeMat=zeros(M,N);for m=0:M-1 cache=fft(sTilde.*conj(gammaTildeVec(mod(0:(L-1)-m0-m*dM,L)+1); cTildeMat(m+1,:)=cache(mod(n0+(0:N-1)*dN,L)+1);endend附录7：对偶标架逆变换的matlab程序

21、%逆对偶标架变换程序%cTildeMat为对偶标架系数矩阵，变量命名基于殷勤业讲义%hTildeVec为综合函数离散向量%dM为对偶标架时间采样间隔点数%m0为空间采样起点，要求其值为0,dM-1间的整数%n0为波数采样起点，要求其值为-L/2,-L/2+dN-1间的整数%sTilde0为对偶标架重构信号function sTilde0=inverseDualFrameTransform(cTildeMat,hTildeVec,dM,dN,m0,n0)%判断输入的有效性if rem(m0,1)=0 | m0dM-1 error(m0必须为0,dM-1-1间的整数);endif rem(n0,1

22、)=0 | n0-length(hTildeVec)/2+dN-1 error(n0必须为-L/2,-L/2+dN-1间的整数);end%计算所需常量M,N=size(cTildeMat);L=length(hTildeVec);unitImaginary=sqrt(-1);%重构信号sTilde0=zeros(1,L);ifftCTildeMat=ifft(cTildeMat,2);for m=0:M-1 sTilde0=sTilde0+hTildeVec(mod(0:L-1)-m0-m*dM,L)+1).*ifftCTildeMat(m+1,mod(0:L-1,N)+1);endsTild

23、e0=N*exp(unitImaginary*n0*(0:L-1)*2*pi/L).*sTilde0;end附录8：Grossman对偶标架变换的matlab程序%主程序clear;clc;%初始化变量dt=0.001;L=210;sTilde=sin(2*pi*100*(0:L-1)*dt).1); %构造正弦信号%sTilde=sin(400*(0:L/4-1)*pi*dt),sin(800*(L/4:2*L/4-1)*pi*dt),sin(200*(2*L/4:3*L/4-1)*pi*dt),sin(600*(3*L/4:L-1)*pi*dt); %构造四段跳频信号%sTilde=(1-

24、2*(pi*(50)*(-L/2:L/2-1)*dt).2).*exp(-(pi*(50)*(-L/2:L/2-1)*dt).2); %ricker信号;%sTilde,fs,=wavread(截取忧伤还是快乐.wav);sTilde=transpose(sTilde(:,1);sTilde=sTilde(1:L);dt=1/fs; %实际音频信号%计算Grossman对偶标架if fix(log2(length(sTilde)=log2(length(sTilde) error(信号的长度必须为2的整数次方);endL=length(sTilde);C1=256; %设置对偶标架变换的冗余度dM=2(ceil(log2(sqrt(L/C1);dN=L/(C1*dM); %模拟紧标架变换deltaX与