小波分析实验： 实验1 连续小波变换

实验目的：

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

实验工具：

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