实验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 实验2 二维离散小波变换(Mallat快速算法)实验目的： 在理解离散小波变换原理和Mallat快速算法的基础上，通过编程对图像进行二维离散小波变换，从而加深对二维小波分解和重构的理性和感性认识，并能提高编程能力，为今后的学习和工作奠定基础。 实验工具： 计算机，matlab6.5
附录： （1）二维小波分解函数 %二维小波分解函数
function Y=mallatdec2(X,wname,level) %输入:X 载入的二维图像像数值; % level 小波分解次(级)数设定值(如果设定值超过最高可分解次数,按最高分解次数分解) % wname 小波名字wavelet name %输出：Y 多极小波分解后的小波系数矩阵
[h,g]=wfilters(wname,'d'); %h,g分别为低通和高通滤波器 X=double(X); t=1; hh=size(X,2);
while t<=level %先进行行小波变换 for row=1:hh Y(row,1:hh)=mdec1(X(row,1:hh),h,g) ; end %再进行列小波变换 for col=1:hh temp=mdec1( Y(1:hh,col)',h,g); Y(1:hh,col)=temp'; end t=t+1; hh=hh/2; X=Y; end
%内部子函数，对一行(row)矢量进行一次小波变换，利用fft实现 function y=mdec1(x,h,g) %输入：x 行数组 % h为低通滤波器 % g为高通滤波器 %输出： y 进行一级小波分解后的系数 lenx=size(x,2); lenh=size(h,2);
rh=h(end:-1:1); rrh=[zeros(1,(lenx-lenh)),rh]; rrh=circshift(rrh',1)';
rg=g(end:-1:1); rrg=[zeros(1,(lenx-lenh)),rg]; rrg=circshift(rrg',1)'; r1=dyaddown(ifft(fft(x).*fft(rrh,lenx)),1); %use para 1 r2=dyaddown(ifft(fft(x).*fft(rrg,lenx)),1); y=[r1,r2];
（2）二维小波重构函数 %二维小波重构函数 function Y=mallatrec2(X,wname,level) %输入:X 载入的小波系数矩阵; % level 小波分解次(级)数设定值(如果设定值超过最高可分解次数,按最高分解次数分解) % wname 小波名字wavelet name %输出：Y 重构图像矩阵
[h,g]=wfilters(wname,'d'); %h,g分别为重构低通滤波器和重构高通滤波器
hz=size(X,2); h1=hz/(2^(level-1)); while h1<=hz % 对列变换 for col=1:h1 temp=mrec1(X(1:h1,col)',h,g)'; X(1:h1,col)=temp; end %再对行变换 for row=1:h1 temp=mrec1(X(row,1:h1),h,g); X(row,1:h1)=temp; end h1=h1*2;
end Y=X;
%内部子函数，对一行小波系数进行重构 function y=mrec1(x,h,g) %输入：x 行数组 % h为低通滤波器 % g为高通滤波器 %输出： y 进行一级小波重构后值 lenx=size(x,2);
r3=dyadup(x(1,1:lenx*0.5),0); %内插零use para 0 r4=dyadup(x(1,(lenx*0.5+1):lenx),0); %use para 0 y=ifft(fft(r3,lenx).*fft(h,lenx))+ ifft(fft(r4,lenx).*fft(g,lenx));
|
|
|
