小波变换

一、绘制原理:

1.需要用到的小波工具箱中的三个函数cwt(),centfrq(),scal2frq()

COEFS = cwt(S,SCALES,'wname')

该函数实现连续小波变换,其中S为输入信号,SCALES为尺度,wname为小波名称。

 

FREQ = centfrq('wname')

该函数求以wname命名的母小波的中心频率。

 

F = scal2frq(A,'wname',DELTA)

该函数能将尺度转换为实际频率,其中A为尺度,wname为小波名称,DELTA为采样周期。

 

2.尺度与频率之间的关系

设a为尺度,fs为采样频率,Fc为小波中心频率,则a对应的实际频率Fa为

Fa=Fc*fs/a

显然,根据采样定理,为使小波尺度图的频率范围为(0,fs/2),尺度范围应为(2*Fc,inf),其中inf表示为无穷大。在实际应用中,只需取尺度足够大即可。

 

3.尺度序列的确定

由上式可以看出,为使转换后的频率序列是一等差序列,尺度序列必须取为以下形式:

c/totalscal, c/(totalscal-1), ...,c/2,c

其中,totalscal是对信号进行小波变换时所用尺度序列的长度(通常需要预先设定好),c为一常数。

而尺度c/totalscal所对应的实际频率应为fs/2,于是可得

c=2*Fc*totalscal

于是可得到所需的尺度序列。

 

4.时频图的绘制

确定了小波基和尺度后,就可以用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('小波时频图');

image

image

说明:

在这个例子中,最好选用复的morlet小波,其它小波的分析效果不好,而且morlet小波的带宽参数和中心频率取得越大,时频图上反映的时频聚集性越好。

 

与其他时频分析对比,如短时傅里叶变换时频图

image

 

小结:

image

image

image

image

高频时小波效果好,因为小波在高频处分辨率可以自动调整,分辨率高。

代码:

% 时频分析

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);%+randn(1, length(t));

% s = chirp(t,20,0.3,300);

s = chirp(t,20,1,500,'q');

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('小波时频图');

% 短时傅里叶变换时频图

figure

spectrogram(s,256,250,256,fs);

% 时频分析工具箱里的短时傅里叶变换

f = 0:fs/2;

tfr = tfrstft(s');

tfr = tfr(1:floor(length(s)/2), :);

figure

imagesc(t, f, abs(tfr));

set(gca,'YDir','normal')

colorbar;

xlabel('时间 t/s');

ylabel('频率 f/Hz');

title('短时傅里叶变换时频图');

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值