matlab 如何对图像进行9/7小波分解

小波分解重构 V2.0 版程序存在的问题分析
http://blog.csdn.net/chenyusiyuan/archive/2008/07/09/2628911.aspx

小波图像分解 Matlab 程序 - V3.0版
http://blog.csdn.net/chenyusiyuan/archive/2008/07/09/2630153.aspx

小波图像重构 Matlab 程序 - V3.0版
http://blog.csdn.net/chenyusiyuan/archive/2008/07/09/2630365.aspx
%----------------------------------------------------------%

本文给出了小波图像分解程序的修正代码，并对一些细节问题进行了图示讨论。
修正前的小波图像分解与重构程序，请看如下文章：
相关的文章有：
1、自己动手编写小波信号分解与重构的Matlab程序
http://blog.csdn.net/chenyusiyuan/archive/2007/11/13/1881781.aspx
2、用自编的程序实现小波图像分解与重构
http://blog.csdn.net/chenyusiyuan/archive/2007/11/13/1881940.aspx
下面是针对上述文章中存在的问题而修改的小波图像分解程序。

function coef=mywavedec2(x,N,wname)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 函数 MYWAVEDEC2() 对输入矩阵 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 级低频矩阵存入矩阵 coef
coef=[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:N
for pc=1:3
subplot(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)==1
ylabel(["Level ",num2str(N-pr+1)]);
end
pn=pn+1;
end
end

function [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-1
for 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，cD
cD=[];
for i=1:dim
cvl=conv(cA,lpd);
% 低通滤波，为了提高运行速度，调用MATLAB提供的卷积函数conv()
dnl=downspl(cvl); % 通过下抽样求出平均部分的分解系数
cvh=conv(cA,hpd); % 高通滤波
dnh=downspl(cvh); % 通过下抽样求出本层分解后的细节部分系数
cA=dnl; % 下抽样后的平均部分系数进入下一层分解
cD=[cD,dnh]; % 将本层分解所得的细节部分系数存入序列cD
end

function 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);

—— 图示讨论
1、小波分解的行、列变换过程（使用Haar小波）
% 行变换代码
[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

行变换的结果图示：
可见，行变换将图像矩阵分为左右两部分，左边是平均系数，右边是细节系数，并且由图可见细节系数是垂直性的，属于 vertical detail。
% 列变换代码
[row1,col1]=size(x1); % row1=row + lnf -1, col1=col+lnf-1
for k=1:col1 % 再对缓存矩阵 x1 的每一列序列进行一维离散小波分解
tmp2=x1(:,k);
[ca2,cd2]=mydwt(tmp2,Lo_D,Hi_D,1);
x2(:,k)=[ca2,cd2]"; % 将分解所得系数存入缓存矩阵 x2 中
end

列变换的结果图示：
列变换后,所得矩阵就是一级小波变换的结果，可分为4部分：左上角的平均系数 cA、右上角的垂直细节系数 cV、左下角的水平细节系数 cH、右下角的对角线细节系数 cD。则 mydwt2 的输出序列是 [cA,cV,cH,cD]。不过，我不大理解的是，一般教材和Matlab的说明文档都是把系数序列按这样的次序列出的：[cA,cH,cV,cD] ，即先水平后垂直，在显示时，水平细节在右上角，垂直细节在左下角。
2、小波分解的结果
（1）Haar 小波，2级分解
（2）Bior3.7 小波，2级分解