小波变换 mallat教学文案

小波变换m a l l a t实验目的：通过编程实现离散快速小波变换Mallat 算法，从而加深理解二维小波变换的分解与合成，同时，提高编程能力和matlab 的应用，为以后的学习打下基础。

实验原理：1、Mallat 快速算法本实验使用离散快速小波变换快速算法Mallat 算法，算法原理如下1(2)j j k n nc h n k c -=-∑ （1）1(2)j j k n nd g n k c -=-∑ （2）重构算法： 1(2)(2)j j j n k k nnc h n k c g n kd -=-+-∑∑ （3）对于（1）、（2）等效于1j n c -经过冲击响应为[]h n -和[]g n -的数字滤波器，然后再分别进行“二抽取”，Mallat 分解算法的滤波器表示形式如下图C j-1d j (k)C j (k)用滤波器表示如下图d jC j C j-1(k)2、 255*25510lgPSNR MSE='211()*M Nij ij i j ff MSE M N==-=∑∑{}ij f '{}ij f 分别表示原始图像和重建后的图像，1,1i M j N ≤≤≤≤。

3、边界延拓方法有零延拓、周期延拓、对称周期延拓、常数连续延拓等，本实验采用以上四种方法进行原图像的1/8延拓，并进行重构，各种延拓方法所对应的函数为yan0(x)、yancir （x ）、yan(x)、yanc(x),在主程序中，需要某种延拓，便调用某种函数。

实验编程思路：为使程序易于理解，在不考虑算法复杂度的情况下，分解程序采用简洁的循环计算出下一级的分解系数，程序采用的编程思想如下[][][]11100[0][1][2][3][4][5]001[1]00[0][1][2][3]00[1][2][3][4][5]00[0][1]12j j j j j j c c h h h h h h c c h h h h n c n h h h h h h c ---⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎡⎤⎢⎥--⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎣⎦⎣⎦以上矩阵等式左面是进行二抽样的结果，[0][1]2j j nc c -是j 分解的低频部分。

同理，对于j 分解的高频部分有如下矩阵形式：[][][]11100[0][1][2][3][4][5]001[1]00[0][1][2][3]00[1][2][3][4][5]00[0][1]12j j j j j d d g g g g g g d d g g g g n d n g g g g g g d ---⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎡⎤⎢⎥--⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎣⎦⎣⎦ 分解程序： lenx=size(x,2);%x 为一维向量lenh=size(h,2); h=[h,zeros(1,(lenx-lenh))]; g=[g,zeros(1,(lenx-lenh))];r1(1)=sum(h.*x); r2(1)=sum(g.*x);for k=1:1:(lenx/2-1) %循环求出下一级低频和高频分量h=[h(end-1:end),h(1:(end-2))];r1(k+1)=sum(h.*x); g=[g(end-1:end),g(1:1:(end-2))];r2(k+1)=sum(g.*x);endy=[r1,r2];对于重构算法，其等效形式为[][][]1(2)(2)j j j nnc n h n k c k g n kd k -=-+-∑∑上式等号右边部分实质上是对变量k 的数字卷积运算,程序采用频域相乘代替卷积，重建程序为 y=ifft(fft(c3,lenx).*fft(h,lenx))+ifft(fft(d3,lenx).*fft(g,lenx));实验结果及分析：1、多尺度分解与重构图像二维小波变换采用小波采用db3，其峰值信噪比PNSR=230.13db，并对三级分解图像进行归一化，求出0的个数为37626，其所占的百分比为57.41%。

2、延拓重建图像延拓方法周期延拓对称周期延拓零延拓常数连续延拓PNSR 230.18db 230.48db 229.80db 230.13db 从PNSR结果可知，在各种延拓中，对称周期延拓的重建图像结果最好，相比之下零延拓图像效果不如其他方法延拓。

3、不同小波下重构图像的性质用不同小波进行图像重构，所得的重构图像能量分布如下用各种小波进行重构后的图像的均值方差如下表。

小波 db1db2db3db4均值 124.0309 124.0509 124.0509 124.0509 方差 2272227222722272附录：1、主函数程序clc;clear;X=imread('LENA.bmp');%路径X=double(X);% S=yancir(X);A=mallatdec2(X,'db3', 3);image(abs(A));colormap(gray(255));title('3级多尺度分解图像');Y=mallatrec2(A,'db3',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('重建图像');zerosn=numberzeros(A);% Y=Y(33:288,33:288); %当调用延拓图像时，从延拓的重建图像进行截取csize=size(X);sr=csize(1);sc=csize(2);mse=sum(sum( (Y-X).^2,1))/(sr*sc);psnr=10*log(255*255/mse)/log(10)2、分解程序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,ll]=size(X);while t<=level%先进行行小波变换for row=1:hhY(row,1:ll)=mdec1(X(row,1:ll),h,g) ;end%再进行列小波变换for col=1:lltemp=mdec1( Y(1:hh,col)',h,g);Y(1:hh,col)=temp';endt=t+1;hh=hh/2;ll=ll/2;X=Y;endfunction y=mdec1(x,h,g)%输入：x 行数组% h 为低通滤波器% g 为高通滤波器%输出： y 进行一级小波分解后的系数lenx=size(x,2);lenh=size(h,2);h=[h,zeros(1,(lenx-lenh))];g=[g,zeros(1,(lenx-lenh))];r1(1)=sum(h.*x);r2(1)=sum(g.*x);for k=1:1:(lenx/2-1)h=[h(end-1:end),h(1:(end-2))];r1(k+1)=sum(h.*x);g=[g(end-1:end),g(1:1:(end-2))];r2(k+1)=sum(g.*x);endy=[r1,r2];3、重建程序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: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;function y=mrec1(x,h,g)%输入：x 行数组% h 为低通滤波器% g 为高通滤波器%输出： y 进行一级小波重构后值lenx=size(x,2);lenh=size(h,2);h=[h,zeros(1,(lenx-lenh))];g=[g,zeros(1,(lenx-lenh))];;c3=dyadup(x(1,1:lenx*0.5),0,lenx); %内插零use para 0d3=dyadup(x(1,(lenx*0.5+1):lenx),0,lenx); %use para 0y=ifft(fft(c3,lenx).*fft(h,lenx))+ ifft(fft(d3,lenx).*fft(g,lenx));4、延拓程序%零延拓程序function T=yan0(x)[r0,c0]=size(x);for s=1:r0a(s,:)=[zeros(1,c0/8),x(s,:),zeros(1,c0/8)];endc0=c0+c0/4;for t=1:c0T(:,t)=[zeros(r0/8,1)',a(:,t)',zeros(r0/8,1)']';End%连续常数延拓程序function T=yanc(x)[r0,c0]=size(x);for s=1:r0r1(1:c0/8)=x(s,1);r2(1:c0/8)=x(s,c0);a(s,:)=[r1,x(s,:),r2];精品资料endc0=c0+c0/4;for t=1:c0c1(1:r0/8)=a(1,t);c2(1:r0/8)=a(r0,t);T(:,t)=[c1,a(:,t)',c2]';end%对称周期延拓程序function T=yancir(x)[r0,c0]=size(x);for s=1:r0a(s,:)=[x(s,end-(c0/8-1):1:end),x(s,:),x(s,1:1:c0/8)];endc0=c0+c0/4;for t=1:c0T(:,t)=[a(end-(r0/8-1):1:end,t)',a(:,t)',a(1:1:r0/8,t)']'; end%周期延拓程序function T=yan(x)[r0,c0]=size(x);for s=1:r0a(s,:)=[x(s,c0/8:-1:1),x(s,:),x(s,end:-1:end-(c0/8-1))];endc0=c0+c0/4;for t=1:c0T(:,t)=[a(r0/8:-1:1,t)',a(:,t)',a(end:-1:end-(r0/8-1),t)']'; end仅供学习与交流，如有侵权请联系网站删除谢谢11。