一、绘制原理:
需要用到的小波工具箱中的三个函数cwt(),centfrq(),scal2frq()。具体参数及用途介绍如下:
(1)COEFS = cwt(S,SCALES,'wname')
该函数实现连续小波变换,其中S为输入信号,SCALES为尺度,wname为小波名称。
(2)FREQ = centfrq('wname')
该函数求以wname命名的母小波的中心频率。
(3)F = scal2frq(A,'wname',DELTA)
该函数能将尺度转换为实际频率,其中A为尺度,wname为小波名称,DELTA为采样周期。
二、尺度与频率之间的关系
设a为尺度,fs为采样频率,Fc为小波中心频率,则a对应的实际频率Fa为Fa=Fc*fs/a,显然,根据采样定理,为使小波尺度图的频率范围为(0,fs/2),尺度范围应为(2*Fc,inf),其中inf表示为无穷大。在实际应用中,只需取尺度足够大即可。
三、尺度序列的确定
由上式可以看出,为使转换后的频率序列是一等差序列,尺度序列必须取为以下形式: c/totalscal, c/(totalscal-1), ...,c/2,c 其中,totalscal是对信号进行小波变换时所用尺度序列的长度(通常需要预先设定好),c为一常数。 而尺度c/totalscal所对应的实际频率应为fs/2,于是可得 c=2*Fc*totalscal 于是可得到所需的尺度序列。
四、时频图的绘制
确定了小波基和尺度后,就可以用cwt求小波系数coefs(系数是复数时要取模),然后用scal2frq将尺度序列转换为实际频率序列f,最后结合时间序列t,用imagesc(t,f,abs(coefs))便能画出小波时频图。
应用例子:下面给出一实际例子来说明小波时频图的绘制。所取仿真信号是由频率分别为50Hz和100Hz的两个正弦分量所合成的信号。
% 小波时频分析
clc
clear all
close all
% 原始信号
fs=1000;
f1=50;
f2=100;
t=0:1/fs:1;
s=sin(2*pi*f1*t)+sin(2*pi*f2*t);
figure
plot(t, s)
% 连续小波变换
wavename='cmor3-3';
totalscal=256;
Fc=centfrq(wavename); % 小波的中心频率
c=2*Fc*totalscal;
scals=c./(1:totalscal);
f=scal2frq(scals,wavename,1/fs); % 将尺度转换为频率
coefs=cwt(s,scals,wavename); % 求连续小波系数
figure
imagesc(t,f,abs(coefs));
set(gca,'YDir','normal')
colorbar;
xlabel('时间 t/s');
ylabel('频率 f/Hz');
title('小波时频图');
运行结果如下:
说明:在这个例子中,最好选用复的morlet小波,其它小波的分析效果不好,而且morlet小波的带宽参数和中心频率取得越大,时频图上反映的时频聚集性越好。