轉(zhuǎn)自：http://blog.21ic.com/user1/5396/archives/2009/57134.html，稍作修改，方便理解。

一、繪制原理：

1.需要用到的小波工具箱中的三個(gè)函數(shù)cwt()，centfrq()，scal2frq()

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

該函數(shù)實(shí)現(xiàn)連續(xù)小波變換，其中S為輸入信號(hào)，SCALES為尺度，wname為小波名稱。

FREQ = centfrq('wname')

該函數(shù)求以wname命名的母小波的中心頻率。

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

該函數(shù)能將尺度轉(zhuǎn)換為實(shí)際頻率，其中A為尺度，wname為小波名稱，DELTA為采樣周期。

2.尺度與頻率之間的關(guān)系

設(shè)a為尺度，fs為采樣頻率，F(xiàn)c為小波中心頻率，則a對(duì)應(yīng)的實(shí)際頻率Fa為

Fa＝Fc*fs/a

顯然，根據(jù)采樣定理，為使小波尺度圖的頻率范圍為(0,fs/2)，尺度范圍應(yīng)為(2*Fc,inf),其中inf表示為無窮大。在實(shí)際應(yīng)用中，只需取尺度足夠大即可。

3.尺度序列的確定

由上式可以看出，為使轉(zhuǎn)換后的頻率序列是一等差序列，尺度序列必須取為以下形式：

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

其中，totalscal是對(duì)信號(hào)進(jìn)行小波變換時(shí)所用尺度序列的長度(通常需要預(yù)先設(shè)定好)，c為一常數(shù)。

而尺度c/totalscal所對(duì)應(yīng)的實(shí)際頻率應(yīng)為fs/2，于是可得

c=2*Fc*totalscal

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

4.時(shí)頻圖的繪制

確定了小波基和尺度后，就可以用cwt求小波系數(shù)coefs（系數(shù)是復(fù)數(shù)時(shí)要取模），然后用scal2frq將尺度序列轉(zhuǎn)換為實(shí)際頻率序列f，最后結(jié)合時(shí)間序列t，用imagesc(t,f,abs(coefs))便能畫出小波時(shí)頻圖。

二、應(yīng)用例子

下面給出一實(shí)際例子來說明小波時(shí)頻圖的繪制。所取仿真信號(hào)是由頻率分別為50Hz和100Hz的兩個(gè)正弦分量所合成的信號(hào)。

% 小波時(shí)頻分析

clc

clear all

close all

% 原始信號(hào)

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)

% 連續(xù)小波變換

wavename='cmor3-3';

totalscal=256;

Fc=centfrq(wavename); % 小波的中心頻率

c=2*Fc*totalscal;

scals=c./(1:totalscal);

f=scal2frq(scals,wavename,1/fs); % 將尺度轉(zhuǎn)換為頻率

coefs=cwt(s,scals,wavename); % 求連續(xù)小波系數(shù)

figure

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

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

colorbar;

xlabel('時(shí)間 t/s');

ylabel('頻率 f/Hz');

title('小波時(shí)頻圖');

說明：

在這個(gè)例子中，最好選用復(fù)的morlet小波，其它小波的分析效果不好，而且morlet小波的帶寬參數(shù)和中心頻率取得越大，時(shí)頻圖上反映的時(shí)頻聚集性越好。

與其他時(shí)頻分析對(duì)比，如短時(shí)傅里葉變換時(shí)頻圖

小結(jié)：

高頻時(shí)小波效果好，因?yàn)樾〔ㄔ诟哳l處分辨率可以自動(dòng)調(diào)整，分辨率高。

代碼：

% 時(shí)頻分析

clc

clear all

close all

% 原始信號(hào)

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)

% 連續(xù)小波變換時(shí)頻圖

wavename='cmor3-3';

totalscal=256;

Fc=centfrq(wavename); % 小波的中心頻率

c=2*Fc*totalscal;

scals=c./(1:totalscal);

f=scal2frq(scals,wavename,1/fs); % 將尺度轉(zhuǎn)換為頻率

coefs=cwt(s,scals,wavename); % 求連續(xù)小波系數(shù)

figure

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

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

colorbar;

xlabel('時(shí)間 t/s');

ylabel('頻率 f/Hz');

title('小波時(shí)頻圖');

% 短時(shí)傅里葉變換時(shí)頻圖

figure

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

% 時(shí)頻分析工具箱里的短時(shí)傅里葉變換

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('時(shí)間 t/s');

ylabel('頻率 f/Hz');

title('短時(shí)傅里葉變換時(shí)頻圖');