matlab时频分析之连续小波变换cwt

（2020年7月更新，第3节绘制了一个实部、虚部的关系图）

1 小波分析简介

和傅里叶变换比，小波变换和短时傅里叶变换都有着相同的优点，就是可以同时在时域和频域观察信号。所以小波变换非定常信号的分析中有很大的作用。

和短时傅里叶变换相比，小波变换有着窗口自适应的特点，即高频信号分辨率高（但是频率分辨率差），低频信号频率分辨率高（但是时间分辨率差），而在工程中常常比较关心低频的频率，关心高频出现的时间，所以近些年用途比较广泛。

2 小波分析基本原理

3 cwt的matlab实现

由于名称一样，在使用中，要想调用不同的版本，只能用输入输出格式来区分。

旧版本cwt用法为：

coefs = cwt(x,scales,'wname') 

新版本cwt用法为：

[wt,f] = cwt(x,wname,fs) 

最明显的区别在于，新版本取消了自定义尺度函数scales的功能。同时新版本还更新替换了一部分小波函数。旧版本支持 ‘haar’，‘db’，‘sym’，‘cmor’，‘mexh’，‘gaus’，‘bior’等小波，新版本支持’morse’， ‘amor’和 ‘bump’小波。

新版本的小波函数用法如下：

% 定义信号信息 fs=2^6; %采样频率 dt=1/fs; %时间精度 timestart=-8; timeend=8; t=(0:(timeend-timestart)/dt-1)*dt+timestart; L=length(t); z=4*sin(2*pi*linspace(6,12,L).*t); %matlab自带的小波变换 %新版本 figure(1) [wt,f,coi] = cwt(z,'amor',fs); pcolor(t,f,abs(wt));shading interp 

% 定义信号信息 fs=2^6; %采样频率 dt=1/fs; %时间精度 timestart=-8; timeend=8; t=(0:(timeend-timestart)/dt-1)*dt+timestart; L=length(t); z=4*sin(2*pi*linspace(6,12,L).*t); %旧版本 wavename='cmor1-3'; %可变参数，分别为cmor的 %举一个频率转尺度的例子 fmin=2; fmax=20; df=0.1; f=fmin:df:fmax-df;%预期的频率 wcf=centfrq(wavename); %小波的中心频率 scal=fs*wcf./f;%利用频率转换尺度 coefs = cwt(z,scal,wavename); figure(2) pcolor(t,f,abs(coefs));shading interp 

对于cmor小波Fc和Fb的选取，可以认为最终结果只与Fc*√Fb的乘积大小有关就行。实际应用中具体值需要根据最终效果选择。

其实cwt的原理很简单，就是用不同尺度的小波逐个窗口去卷积，得到小波系数矩阵。所以根据原理，可以自己编程实现小波变换。

自己编程的cwt函数如下，这里主要算法参考了matlab官方的文档。这里我依然用的是cmor小波作为例子，morlet函数公式可以查询到，也可以用cmorwavf()函数调用：

fs=2^6; %采样频率 dt=1/fs; %时间精度 timestart=-8; timeend=8; t=(0:(timeend-timestart)/dt-1)*dt+timestart; L=length(t); z=4*sin(2*pi*linspace(6,12,L).*t); %定义计算范围和精度 fmin=2; fmax=20; df=0.1; totalscal=(fmax-fmin)/df; f=fmin:df:fmax-df;%预期的频率 wcf=centfrq(wavename); %小波的中心频率 scal=fs*wcf./f; %自己实现的小波函数 coefs2=cwt_cmor(z,1,3,f,fs); figure(3) pcolor(t,f,abs(coefs2));shading interp %后面是函数 function coefs=cwt_cmor(z,Fb,Fc,f,fs) %1 小波的归一信号准备 z=z(:)';%强行变成y向量，避免前面出错 L=length(z); %2 计算尺度 scal=fs*Fc./f; %3计算小波 shuaijian=0.001;%取小波衰减长度为0.1% tlow2low=sqrt(Fb*log(1/shuaijian));%单边cmor衰减至0.1%时的时间长度，参照cmor的表达式 %3小波的积分函数 iter=10;%小波函数的区间划分精度 xWAV=linspace(-tlow2low,tlow2low,2^iter); stepWAV = xWAV(2)-xWAV(1); val_WAV=cumsum(cmorwavf(-tlow2low,tlow2low,2^iter,Fb,Fc))*stepWAV; %卷积前准备 xWAV = xWAV-xWAV(1); xMaxWAV = xWAV(end); coefs = zeros(length(scal),L);%预初设coefs %4小波与信号的卷积 for k = 1:length(scal) %一个scal一行 a_SIG = scal(k); %a是这一行的尺度函数 j = 1+floor((0:a_SIG*xMaxWAV)/(a_SIG*stepWAV)); %j的最大值为是确定的，尺度越大，划分的越密。相当于把一个小波拉伸的越长。 if length(j)==1 , j = [1 1]; end waveinscal = fliplr(val_WAV(j));%把积分值扩展到j区间，然后左右颠倒。f为当下尺度的积分小波函数 %5 最重要的一步 wkeep1取diff(wconv1(ySIG,f))里长度为lenSIG的中间一段 %conv(ySIG,f)卷积。 coefs(k,:) = -sqrt(a_SIG)*wkeep1(diff(conv2(z,waveinscal, 'full')),L); % end end 

4 cwt的边缘效应与影响锥

利用自带的cwt函数，可以很方便的绘制出影响锥。

cwt(z) 

由于小波计算中，小波系数是利用窗口函数和小波卷积而来的，当窗口在信号的边缘时，窗口内会存在一部分没有信号。这时，matlab就把窗口内这部分不完整的信号补零处理，凑够长度。

此时由于信号在边缘被强制补零，导致信号会失真，具体在时频图中表示为频率变宽，信号强度降低。严重的时候，甚至整个低频部分都会出现失真。这就是cwt的边缘效应。

代码如下，依然以cmor小波为例，为了反应相位的影响，这里用的是real()：

%小波变换展示 clear fs=2^6; %采样频率 dt=1/fs; %时间精度 timestart=-8; timeend=8; t=(0:(timeend-timestart)/dt-1)*dt+timestart; L=length(t); z=sin(2*pi*5.*t)+sin(2*pi*9.*t)+sin(2*pi*15.*t); %定义范围 wavename='cmor1-3'; %可变参数 fmin=2; fmax=20; df=0.1; totalscal=(fmax-fmin)/df; f=fmin:df:fmax-df;%预期的频率 wcf=centfrq(wavename); %小波的中心频率 scal=fs*wcf./f; %旧版本 coefs = cwt(z,scal,wavename); figure(2) pcolor(t,f,real(coefs));shading interp tlow2low=sqrt(1*log(1/0.001));%单边cmor衰减至0.1%时的时间长度 tcoi41=tlow2low*scal;%小波一半长度 bur=(tcoi41 
   

5 cwt的重构——icwt

如果用的是matlab新版的默认小波，那么可以直接把cwt之后小波系数icwt即可。

load mtlb; wt = cwt(mtlb); xrec = icwt(wt); 

对于morlet小波，可以直接 sum(real(coefs),1) 来实现重构。但是，这里尺度最好用默认尺度，否则会出现重复叠加导致幅值出现问题。

具体用法可以参见帮助文档。

6 增加cwt的分辨率的wsst

利用wsst可以显著的提高cwt的频率分辨率。具体用法如下：

clear fs=2^6; %采样频率 dt=1/fs; %时间精度 timestart=-8; timeend=8; t=(0:(timeend-timestart)/dt-1)*dt+timestart; L=length(t); z=4*sin(2*pi*linspace(6,12,L).*t); [sst,f] = wsst(z,fs); figure(4) pcolor(t,f,abs(sst));shading interp 

如果想要提高wsst的精度，可以在matlab里选中wsst，然后Ctrl+D进入wsst的代码界面，把里面的na参数，强行改的大一点。

具体方法为在下面代码后，加上一个na=512（这里不一定是512，具体根据na实际值进行适当放大，在我这里的例子中，na=288。如果怕出错，建议备份）。

na = noct*params.nv; na=512 

这样做的优点是，在不影响iwsst的运行结果的基础上，显著的提高分解的精度。