CSDN 一维小波分解与重构python 一维连续小波变换

实验目的：

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

实验工具：

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