實驗目的：

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

實驗工具：

計算機，matlab6.5

?

程序附錄：

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

% mexh.m

function Y=mexh(x)

if abs(x)<=8

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

else

??? Y=0;

End

?

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

%

clc;clear;

load('data.mat');

len=length(dat);

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

a=zeros(1,lna);

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('原始數據圖');

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

image(wfab);

colormap(pink(128));

title('小波系數圖');

%surf(wfab);

%title('小波系數譜網格圖');

%pwfab=wfab.*wfab;? %%瞬態功率譜

%figure(3);

%subplot(1,2,1);

%surf(pwfab);

%title('瞬態功率譜網格圖');

%subplot(1,2,2);

%contour(pwfab);

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

?

(3)test函數。

%test 函數

clc;clear;

for i=1:200

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

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);? %小波系數譜

image(wfab);

colormap(pink(128));

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

（4）用fft實現cwt

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

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

%注意畫圖理解

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);?? %小波系數矩陣

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('原始數據圖');

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

image(wfab);

colormap(pink(128));

title('小波系數譜 ');

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('原始數據圖');

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

image(wfab);

colormap(pink(128));

title('小波系數譜 ');

cputime-t1

?

5）保存為mexh2.m

?

function Y=mexh2(x)

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

?

Torstan

2005.09.16

轉載于:https://www.cnblogs.com/Torstan/archive/2011/08/31/2161427.html

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。