实验目的：

在理解连续小波变换原理的基础上，通过编程实现对一维信号进行连续小波变换，（实验中采用的是墨西哥帽小波），从而对连续小波变换增加了理性和感性的认识，并能提高编程能力，为今后的学习和工作奠定基础。

实验工具：

计算机，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