小波分析实验：二维离散小波变换(Mallat快速算法)

小波分析实验：实验2二维离散小波变换（Mallat快速算法）实验目的：在理解离散小波变换原理和Mallat快速算法的基础上，通过编程对图像进行二维离散小波变换，从而加深对二维小波分解和重构的理性和感性认识，并能提高编程能力，为今后的学习和工作奠定基础。

实验工具：计算机，matlab6.5分解算法:重构算法: “"二工必（刃- 2上*［十三g （刃- 2k ）d ［ *分解算法写成矩阵的形式！ （lb g 的长度为4）4[0]如]力⑵ h[3] 0 0 0 '［勺【0】• 记"h[0] h[\]h[2]山⑶ …• ••••・ • •C J=勺【1] • •申[2] h[3] 00 0-.^[0] ^[1]_.勺［乃-1】_>[0] g[l] g ⑵ g[3] 0 • • •e=• 0 •g[0] g[l]g ⑵ • • g[3]■ • •・■ 0• D J =<[i]■•目2] ■g[3]0 0…茎0] 畀］|g[0] g[l] g[2] g[3] 0 0 0 I0 0 g[0] g[l]g[2] S [3] - 0• ••••• • ••・•・・■ • • g[2] g[3] 0 00 ...g[0] g[l]J |_勺4-1[叨］I二・(2»于是Mallat分解公式为矩阵变换？丄Cj- = PC^................. ⑶卩D j = Q D J-L..... .......... ⑷重构算法写成矩阵变换：-C J_I =C$ + Dj------------------------------------ (5) 4M NPPq. 一片『峰值信噪比计算公式：P沁沁逻竺皿E卢H耿V 屈E M ｛皿，00分别表示原始图像和重建图像，且本实验采取的一些小技乐P (I)分SW法…编程时用如下思想：（h, g 的长度为4）“今［1］勺［刀-1］■ V■■丐⑼£［1］ 4刀-1】将数据。

』和低通〔高通）滤波器进行添零到数据长度诳，再2抽祥，相加得到（2） 科瞎小波 系数丄后一半 系数门高通.重构滤波 器卩重构结果加0] 方[1] 农］ 舛3] 00 • •0 o T=纽0] 列1] 風2］疏3］■町] ■列2] ■扯3]• •D 0•…0'smg[l] 列2] 或习0• • • 00 11 =0 烈0] 畀］ 02]烈 3] 0:.°■-g[l]■g[2] ■g[3]• • 0 0■ 0•••0 g[0]J[%[1]勺4-1]怕】编程吋用如下思想：2ii[0]山1]爪)]川］/〔0]0 000 00••••••••0 00•托3] 42]0・・托3]• •• •• •0 ••/(!]0 0 00…灶3] 夙0] 0 00 0•••鼻••♦• • •0 0 00 0■心[0]「心[1]■•十'咕[0「吆[1]•弓4-1]_/j-Jn-l]5・1[0] 5 di]川]■••炷1]址0]勺[1]曲］「77+1实验结果：“多尺度分堺E3像SC1）三级小波分解图Q50100150200250原始图條逅构图像50 15C 200 250 50 100 130 200 250EC2 ）原始凰慷和重构图像峰值信噪比（psm—252.8923db^附录：(1)二维小波分解函数%二维小波分解函数fun cti on Y=mallatdec2(X,w name,level)%输入:X 载入的二维图像像数值；% level 小波分解次(级)数设定值(如果设定值超过最高可分解次数,按最高分解次数分解)% wname 小波名字wavelet name%输岀：丫多极小波分解后的小波系数矩阵[h,g]=wfilters(wname,'d'); %h,g分别为低通和高通滤波器X=double(X);t=1;hh=size(X,2);while t<=level%先进行行小波变换for row=1:hhY(row,1:hh)=mdec1(X(row,1:hh),h,g);end%再进行列小波变换for col=1:hhtemp=mdec1( Y(1:hh,col)',h,g);Y(1:hh,col)=temp';endt=t+1;hh=hh/2;X=Y;end实现%内部子函数，对一行(row)矢量进行一次小波变换，利用fftfun cti on y=mdec1(x,h,g)%输入：x行数组% h为低通滤波器% g为高通滤波器%输岀：y进行一级小波分解后的系数len x=size(x,2);len h=size(h,2);rh=h(e nd:-1:1);rrh=[zeros(1,(le nx-le nh)),rh];rrh=circshift(rrh',1)';rg=g(e nd:-1:1);rrg=[zeros(1,(le nx-le nh)),rg];rrg=circshift(rrg',1)';r1=dyaddow n( ifft(fft(x).*fft(rrh,le nx)),1); %use para 1r2=dyaddow n(ifft(fft(x).*fft(rrg,le nx)),1);y=[r1,r2];(2)二维小波重构函数%二维小波重构函数fun cti on Y=mallatrec2(X,w name,level)%输入:X 载入的小波系数矩阵；,按最高分解次% level 小波分解次(级)数设定值(如果设定值超过最高可分解次数数分解)% wname 小波名字wavelet name%输岀：Y 重构图像矩阵[h,g]=wfilters(wname,'d'); %h,g分别为重构低通滤波器和重构高通滤波器hz=size(X,2);h1=hz/(2A(level-1));while h1<=hz%对列变换for col=1:h1temp=mrec1(X(1:h1,col)',h,g)';X(1:h1,col)=temp;end%再对行变换for row=1:h1temp=mrec1(X(row,1:h1),h,g);X(row,1:h1)=temp;endh1=h1*2;endY=X;%内部子函数，对一行小波系数进行重构fun cti on y=mrec1(x,h,g)%输入：x行数组% h为低通滤波器% g为高通滤波器%输岀：y进行一级小波重构后值len x=size(x,2);r4=dyadup(x(1,(le nx*0.5+1):le nx),O); %use para 0 y=ifft(fft(r3,le nx).*fft(h,le nx))+ ifft(fft(r4,le nx).*fft(g,le nx));(3)测试函数(主函数)%测试函数(主函数)r3=dyadup(x(1,1:le nx*0.5),0); %内插零use para 0clc;clear;X=imread('E:\Libin 的文档实验2 要求\exp2\LENA.bmp');% 路径X=double(X);A = mallatdec2(X,'sym2',3);image(abs(A));colormap(gray(255));ti tle(' 多尺度分解图像’)；Y= mallatrec2(A,'sym2',3);Y=real(Y);figure(2);subplot(1,2,1);image(X);colormap(gray(255));title(' 原始图像');subplot(1,2,2);image(Y);colormap(gray(255));title(' 重构图像');csize=size(X);sr=csize(1);sc=csize(2);mse=sum(sum( (Y-X).A2,1))/(sr*sc);psn r=10*log(255*255/mse)/log(10)小波分析实验：实验1连续小波变换实验目的：在理解连续小波变换原理的基础上，通过编程实现对一维信号进行连续小波变换，(实验中采用的是墨西哥帽小波)，从而对连续小波变换增加了理性和感性的认识，并能提高编程能力，为今后的学习和工作奠定基础。

实验工具：计算机，matlab6.5ZL实验原理:•维连续小波变换公式：°"；（心）二IG2J/（O x y y（:—）处-co1 Qt 7=\^ \2J /⑺*(一^)刃—8 "编程时先离散化,在通过求和代替积分。

卩3+1）心咲Jk 炮A /V /（□）%严'"） 根据衫 〒Cf一丄£ — b"Sr)当小波函第/(Ox 1^1 V(本实验采取的一些小技巧：卩（1）因为实验中的数据通过lo网mta.ma广）（先将数据拷到work文件下）命令后，数据存在了dM矩阵中，所以兀财就是datik）•由于数挺样间隔4为0.03 （常量），所以可以把这个系数忽略，实际上本程序冲-1 ]c-b用的是WAa.b）=\a\2匸产（切;w{^L）J T ay.JrQ）本实验中墨西哥帽函数为：^x/2（x） = （（1~x2）g /|x|-8…0. otherwise〔3）如果按照（◎式计算小波系数矩阵，则算法时间复杂度为。

（沪），3 通的电脑上要运行5分钟左右才能出正确结果。

通过研究（*♦）,可以後对于同一尺度—屮d）与必匕二I）有如下关系：0（上匕）从 a a a二项开始与0（匕）相同，只是第一项不同，（在程序中有说明），所以a把算法时间复杂度降为。

（泌），程序执行时间在5秒钟内• q（4）程序的正确性可由c\vt来验证。

亠（5）本实验还编了个te前程序来说明存在哪些分量，即对正弦波函鐵ZL 样后进行连续小波变换，分析图象结果。

卍实验结果：a(1)实验数据图像和小波系数图像"50 100 150 200 250 300 350400(2) te前数据图像和小波系数图像〜orignal dat5程序附录：(1)墨西哥帽小波函数,按照(***)式编程5fun cti on Y=mexh(x)if abs(x)<=8Y=exp(-x*x/2)*(1-x A 2); elseY=0;End(2) 实验程序，按照(**)式编程，详细过程请参考“本实验采取的一些小技巧” %clc;clear;load('data.mat'); len=len gth(dat);In a=70;a=zeros(1. In a); wfab=zeros(l na,le n); mexhab=zeros(1,le% （尺度a）的长度%小波系数矩阵%离散化小波系数矩阵for s=1: Ina %s表示尺度for k=1:le nmexhab(k)=mexh(k/s);endfor t=1:len % t表示位移wfab(s,t)=(sum(mexhab.*dat))/sqrt(s); 替mexhab=[mexh(-1*t/s),mexhab(1:le n-1)];项并右移%将积分用求和代%mexhab（修改endfigure(1);plot(dat);title(' 原始数据图');figure(2); %小波系数谱image(wfab);colorm ap(pi nk(128));title(' 小波系数图');%surf(wfab);%title('小波系数谱网格图');%pwfab=wfab.*wfab; %%瞬态功率谱%figure(3);%subplot(1,2,1);%surf(pwfab);%title('瞬态功率谱网格图');%subplot(1,2,2);%con tour(pwfab);%title('瞬态功率谱等值线');(3)test 函数。