小波图像分解程序：function coef=mywavedec2(x,N,wname) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 函数MYWA VEDEC2() 对输入矩阵x 进行dim 层分解，得到相应的分解系数矩阵y% 输入参数：x ——输入矩阵% N ——分解级数% wname ——分解所用的小波函数% 输出参数：coef ——分解系数矩阵，其结构如下：% coef = {cA_N;cV_N;cH_N;cD_N;cV_N-1;cH_N-1;cD_N-1;……;cV_1;cH_1;cD_1}% Copyright by Zou Yuhua ( chenyusiyuan ), original : 2007-11-10, modified: 2008-06-04 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 求出小波函数的滤波器组系数向量[Lo_D,Hi_D] = wfilters(wname,'d');% 画出原始图像imshow(x);title('Original Image');% 标明图像大小[r,c]=size(x);xlabel(['Size : ',num2str(r),'*',num2str(c)]);% 将矩阵x的数据格式转换为适合数值处理的double格式xd=double(x);coef=[];for i=1:N[cA,cV,cH,cD]=mydwt2(xd,Lo_D,Hi_D);% 第i 级小波分解xd=cA;% 将第i 级分解得到的低频系数矩阵作为第i+1 级分解的源矩阵outmp={cV;cH;cD};% 将第i 级分解得到的高频系数矩阵cV,cH,cD存入细胞矩阵outmp % 注意细胞矩阵的赋值是用大括号“{}”的，而普通矩阵赋值是用方括号“[]”% 细胞矩阵不要求其中的子矩阵的行列数都相同coef=[outmp;coef];% 将细胞矩阵outmp 存入输出矩阵coef，coef将由空矩阵变为细胞矩阵% 注意这里的方括号不能用大括号取代% 否则，使用大括号会将初始的coef空矩阵也作为细胞矩阵的子矩阵% 而且，在迭代中coef 将是一个不断嵌套的细胞矩阵，不便于后续处理和读取% 上面这个语句是一种有效的在迭代过程中保存数据的方法% 设待存数据为data，可以是单个数、向量或矩阵% 保存数据的矩阵为mat，初始为空矩阵：mat=[]% 则可按以下格式保存迭代过程产生的数据% mat=[mat;data];% 方括号内的分号“;”表示数据data 是按“列”排序的方式存入矩阵mat% mat=[mat,data];% 方括号内的逗号“,”表示数据data 是按“行”排序的方式存入矩阵mat% data 也可以在mat 前嵌入，即mat=[data;mat] 或mat=[data,mat]end% 迭代结束后，矩阵coef 中保存的是各级分解中的高频系数矩阵% 故需将迭代后得到的矩阵cA，即第dim 级低频矩阵存入矩阵coefcoef=[cA;coef];% 最后，小波系数矩阵coef 的结构如下% coef = {cA_N;cV_N;cH_N;cD_N;cV_N-1;cH_N-1;cD_N-1;……;cV_1;cH_1;cD_1}% 画出各级低频、高频系数矩阵% 首先建立一个名为“Wavelet Decomposition -- Wavelet Type: , Levels: ”的图像窗口figure('Name',['Wavelet Decomposition -- Wavelet Type: ',wname,' , Levels: ',num2str(N)]);% 图像的第1行显示低频系数，置中，左右两个subplot为空subplot(N+1,3,2);yt=uint8(coef{1});[yrow,ycol]=size(yt);imshow(yt);title( ['Approximation A',num2str(N)]);xlabel(['Size : ',num2str(yrow),'*',num2str(ycol)]);% 第2-(N+1)行显示各级高频系数titllist={['Vertical Detail V'];['Horizontal Detail H'];['Diagonal Detail D']};pn=2;% pn 是子图的显示序号for pr=1:Nfor pc=1:3subplot(N+1,3,pn+2);yt=[];% 为了使高频细节内容（轮廓、边缘）更清晰，将高频系数增加100灰度值yt=uint8(coef{pn})+100;[yrow,ycol]=size(yt);imshow(yt);title([ titllist{pc},num2str(N-pr+1)]);xlabel(['Size : ',num2str(yrow),'*',num2str(ycol)]);% 每行的第一个图像的Y轴，显示该行高频系数对应的分解级别if mod(pn+2,3)==1ylabel(['Level ',num2str(N-pr+1)]);endpn=pn+1;endendfunction [cA,cV,cH,cD]=mydwt2(x,Lo_D,Hi_D) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 函数MYDWT2() 对输入的r*c维矩阵x 进行二维小波分解，输出四个分解系数子矩阵[LL,HL,LH,HH]% 输入参数：x ——输入矩阵，为r*c维矩阵。
% Lo_D,Hi_D ——小波分解的滤波器组系数向量% 输出参数：cA,cV,cH,cD ——是小波分解系数矩阵的四个相等大小的子矩阵% cA：低频部分分解系数；cV：垂直方向分解系数；% cH：水平方向分解系数；cD：对角线方向分解系数。
% Copyright by Zou Yuhua ( chenyusiyuan ), original : 2007-11-10, modified: 2008-06-04%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [row,col]=size(x);% 读取输入矩阵的大小for j=1:row% 首先对输入矩阵的每一行序列进行一维离散小波分解tmp1=x(j,:);[ca1,cd1]=mydwt(tmp1,Lo_D,Hi_D,1);% tmp1 长度为row ，滤波器长度为lnf ，则[ca1,cd1] 的总长为( row + lnf -1 )x1(j,:)=[ca1,cd1];% 将分解系数序列存入缓存矩阵x1 中end[row1,col1]=size(x1);% row1=row + lnf -1, col1=col+lnf-1for k=1:col1% 再对缓存矩阵x1 的每一列序列进行一维离散小波分解tmp2=x1(:,k);[ca2,cd2]=mydwt(tmp2,Lo_D,Hi_D,1);x2(:,k)=[ca2,cd2]' ;% 将分解所得系数存入缓存矩阵x2 中% 注意不要遗漏了上一行代码中的转置符号“ ’”。
Matlab 6.5 及以下较低的版本不支% 持行、列向量的相互赋值，故要把行向量[ca2,cd2]转置为列向量，再存入x2 的相应列end[row2,col2]=size(x2);cA=x2(1:row2/2,1:col2/2);% cA是矩阵x2的左上角部分cV=x2(1:row2/2,col2/2+1:col2);% cV是矩阵x2的右上角部分cH=x2(row2/2+1:row2,1:col2/2);% cH是矩阵x2的左下角部分cD=x2(row2/2+1:row2,col2/2+1:col2);% cD是矩阵x2的右下角部分function [cA,cD] = mydwt(x,lpd,hpd,dim) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 函数[cA,cD]=MYDWT(X,LPD,HPD,DIM) 对输入序列x进行一维离散小波分解，输出分解序列[cA,cD]% 输入参数：x——输入序列；% lpd——低通滤波器；% hpd——高通滤波器；% dim——小波分解层数。
% 输出参数：cA——平均部分的小波分解系数；% cD——细节部分的小波分解系数。
% Copyright by Zou Yuhua ( chenyusiyuan ), original : 2007-11-10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% cA=x;% 初始化cA，cDcD=[];for i=1:dimcvl=conv(cA,lpd);% 低通滤波，为了提高运行速度，调用MATLAB提供的卷积函数conv()dnl=downspl(cvl);% 通过下抽样求出平均部分的分解系数cvh=conv(cA,hpd);% 高通滤波dnh=downspl(cvh);% 通过下抽样求出本层分解后的细节部分系数cA=dnl;% 下抽样后的平均部分系数进入下一层分解cD=[cD,dnh];% 将本层分解所得的细节部分系数存入序列cDendfunction y=downspl(x) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 函数Y=DOWMSPL(X) 对输入序列进行下抽样，输出序列Y。
% 下抽样是对输入序列取其偶数位，舍弃奇数位。
例如x=[x1,x2,x3,x4,x5]，则y=[x2,x4].% Copyright by Zou Yuhua ( chenyusiyuan ), original : 2007-11-10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% N=length(x);% 读取输入序列长度M=floor(N/2);% 输出序列的长度是输入序列长度的一半（带小数时取整数部分）i=1:M;y(i)=x(2*i);。