小波图像分解 Matlab 程序 - V2.0版

注意！程序有新的修正了，详细请见如下文章：

%----------------------------------------------------------%

小波分解重构 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级分解