實驗?zāi)康?#xff1a;

在理解連續(xù)小波變換原理的基礎(chǔ)上，通過編程實現(xiàn)對一維信號進(jìn)行連續(xù)小波變換，（實驗中采用的是墨西哥帽小波），從而對連續(xù)小波變換增加了理性和感性的認(rèn)識，并能提高編程能力，為今后的學(xué)習(xí)和工作奠定基礎(chǔ)。

實驗工具：

計算機(jī)，matlab6.5

?

程序附錄：

(1) 墨西哥帽小波函數(shù),按照(***)式編程

% mexh.m

function Y=mexh(x)

if abs(x)<=8

Y=exp(-x*x/2)*(1-x^2);

else

??? Y=0;

End

?

(2) 實驗程序，按照(**)式編程，詳細(xì)過程請參考“本實驗采取的一些小技巧”

%

clc;clear;

load('data.mat');

len=length(dat);

lna=70;?????????? % （尺度a）的長度

a=zeros(1,lna);

wfab=zeros(lna,len);?? %小波系數(shù)矩陣

mexhab=zeros(1,len);?? % 離散化小波系數(shù)矩陣

?

for s=1:lna?????????? %s 表示尺度?

??? for k=1:len

??????? mexhab(k)=mexh(k/s);

??? end??

??? for t=1:len??? % t 表示位移

????? wfab(s,t)=(sum(mexhab.*dat))/sqrt(s);?? %將積分用求和代替

????? mexhab=[mexh(-1*t/s),mexhab(1:len-1)];? %mexhab修改第一項并右移

??? end

end

?

figure(1);

plot(dat);

title('原始數(shù)據(jù)圖');

figure(2);? %小波系數(shù)譜

image(wfab);

colormap(pink(128));

title('小波系數(shù)圖');

%surf(wfab);

%title('小波系數(shù)譜網(wǎng)格圖');

%pwfab=wfab.*wfab;? %%瞬態(tài)功率譜

%figure(3);

%subplot(1,2,1);

%surf(pwfab);

%title('瞬態(tài)功率譜網(wǎng)格圖');

%subplot(1,2,2);

%contour(pwfab);

%title('瞬態(tài)功率譜等值線');

?

(3)test函數(shù)。

%test 函數(shù)

clc;clear;

for i=1:200

??? dat(i)=sin(2*pi*i*0.05);???? %正弦波函數(shù)

end

len=length(dat);

lna=40;

wfab=zeros(lna,len);

mexhab=zeros(1,len);

for s=1:lna?????????? %s 表示尺度?

??? for k=1:len

??????? mexhab(k)=mexh(k/s);

??? end??

??? for t=1:len?????????????? % t 表示位移

????? wfab(s,t)=(sum(mexhab.*dat))/sqrt(s);?? %將積分用求和代替

????? mexhab=[mexh(-1*t/s),mexhab(1:len-1)];? %mexhab修改第一項并右移

??? end

end

figure(1);

plot(dat);

title('orignal dat');

figure(2);? %小波系數(shù)譜

image(wfab);

colormap(pink(128));

title('正弦波的小波系數(shù)圖');

（4）用fft實現(xiàn)cwt

%按照圓周卷積定理，原周卷積和線性卷積的關(guān)系L>=M+N-1

%按照圓周卷積的定義，相關(guān)和線性卷積的關(guān)系（原始算法和線性卷積的關(guān)系）

%注意畫圖理解

clc;clear;

t1=cputime;

?

load('data.mat');

len=length(dat);

lna=70;?????????? % a（尺度）的長度

a=zeros(1,lna);??? % a 表示尺度

b=zeros(1,len);???? % b 表示位移

wfab=zeros(lna,len);?? %小波系數(shù)矩陣

mexhab=zeros(1,2*len-1);??

data=[zeros(1,len-1),dat];

Ydata=fft( data ,4*len);

for s=1:lna?????????

??? for k=1:2*len-1

?????? mexhab(k)=mexh((k-len)/s);???

??? end?

?

?? temp=ifft( Ydata.*fft( mexhab,4*len ) ,4*len);

?? wfab(s,:)=real(temp(2*len-1:3*len-2))/sqrt(s); %為什么要取實部而不是取模，我也不是很清楚，可是有種感覺

end

figure(1);

plot(dat);

title('原始數(shù)據(jù)圖');

figure(2);? %小波系數(shù)譜

image(wfab);

colormap(pink(128));

title('小波系數(shù)譜 ');

cputime-t1

?

?

?

4）fft快速計算cwt

%按照圓周卷積的定義，

%注意畫圖理解

clc;clear;

t1=cputime;

?

load('data.mat');

len=length(dat);

lna=70;?????????? % a（尺度）的長度

a=5;

?

data=[dat,zeros(1,len)];

Ydata=fft(dat,2*len);

for s=1:lna?????????

?mexhab=zeros(1,2*len);

?k=[-a*s:1:a*s];

?mexhab(k+len)=mexh2(k./s);

?

?? temp=ifft( Ydata.*fft( mexhab,2*len ) ,2*len);

?? wfab(s,:)=real(temp(len+1:2*len))/sqrt(s); %要取實部而不是取模，呵呵

end

figure(1);

plot(dat);

title('原始數(shù)據(jù)圖');

figure(2);? %小波系數(shù)譜

image(wfab);

colormap(pink(128));

title('小波系數(shù)譜 ');

cputime-t1

?

5）保存為mexh2.m

?

function Y=mexh2(x)

Y=exp(-x.*x/2).*(1-x.^2);

?

Torstan

2005.09.16

轉(zhuǎn)載于:https://www.cnblogs.com/Torstan/archive/2011/08/31/2161427.html

總結(jié)

以上是生活随笔為你收集整理的小波分析实验： 实验1 连续小波变换的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。