小波图像分解与重构程序存在的问题与解决办法

去年11月发布了一系列有关小波变换和图像处理的文章，把学习小波过程中的心得体会和编写的程序放在网上和大家共享交流。半年来，感谢大家的关注和帮助，在相互的讨论交流中，我不断地从大家提出的问题中拓展自己的知识面，对小波的理论及其应用有了更深入的了解和掌握。的确，分享，才是知识永恒之道！

接下来有新的课题，并且要准备考博，研学征途马上要展开新的探索。在此之前，得先整理下博客文章，根据和大家讨论交流中发现的问题，对博客中的程序进行修正。今天发布的文章，就是对小波图像分解和重构程序作出重要的修正，这个程序是了解小波、使用小波最根本的，希望对大家学习小波有所帮助。

首先指出我原来发布的小波图像分解重构的文章中存在的问题。相关的文章有：

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

其存在的问题有：

（1）程序所用的小波函数只有非标准的Haar小波，其滤波器组为 Lo_D=[1/2 1/2], Hi_D=[-1/2 1/2]，是固化在 mydwt2.m 的程序中的，不能选择其他的小波函数；

（2）非标准的Haar小波，其分解出来的系数矩阵中，高频系数的细节内容（轮廓、边缘等特征）不明显；

（3）函数 mydwt2 中列变换的矩阵对象为输入矩阵，这是错误的，其矩阵对象应该是行变换后的缓存矩阵；

（4）函数 mydwt2 的输出用[LL,HL,LH,HH]表示，不是很规范，应改为[cA,cV,cH,cD]来表示，即一级小波变换输出的系数矩阵有4个部分：平均部分、垂直细节部分、水平细节部分和对角线细节部分。

（5）函数 mydwt2 的输出[LL,HL,LH,HH]中，HL 取的是小波行变换和列变换后所得矩阵的左下角部分，即垂直细节部分，但在最终的显示结果中 HL 是在右上角的位置的；类似地，LH 取行、列变换所得矩阵的右上角部分，即水平细节部分，但最终显示位置是在左下角，当时这样定义，是我从教材中看到小波一级分解后所得矩阵的平均、细节部分分布从上到下、从左到右依次为：平均部分、水平细节部分、垂直细节部分、对角线细节部分。但通过编程实现和仔细分析，我发现其分布次序应该是平均部分、垂直细节部分、水平细节部分、对角线细节部分才对，即水平和细节的位置对调，具体过程下文会给出图示，希望能和大家讨论下哪个次序才是正确的。

（5）函数 mywavedec2 的输出 y 是与输入矩阵 x 相同大小的矩阵，并且已将N级分解后所有的平均、细节系数组合成一体的。实际上，这种定义只对Haar小波有效。因为对于Haar小波来说，其滤波器长度仅为2，输入序列与滤波器系数卷积后，输出序列与输入序列长度是相同的（在Matlab中，若信号长度为n，滤波器长度为2N，则小波分解后cA、cD的长度均为：lnf = floor((n-1)/2)+N，若2N=2，则 lnf*2 = n），从而使小波系数矩阵与图像矩阵的大小相同。若使用其它类型的小波，由于 2N > 2，lnf*2 > n，这样函数 mydwt2 输出的系数矩阵将比输入矩阵大，那么 mywavedec2 的输出 y 将会显示出平均与细节部分边界重叠的问题。这说明了试图把平均和细节系数组合成一体，以显示小波分解的塔式结构，理论上是不可行的，至少对于Matlab来说是不可行的。当然，由公式 lnf = floor((n-1)/2)+N ，我们可以对卷积后的序列进行截取，使卷积输出序列与输入序列长度相同，但通过编程实现发现，截取会找出下一级分解错误，所以函数 mywavedec2 的输出 y 的格式要作适当调整。

（6）原程序中要调用 modmat 函数对图像矩阵进行修剪，使之能被 2 的 N 次方整除，主要是为了生成塔式结构图像而设的，对上述问题修正后，这个 modmat 函数已不需使用了。

针对上述问题，我对程序作了以下几点修正：

（1）函数 mywavedec2 / mywaverec2 的输入参数增加了小波函数'wname'，可以选择任意类型的小波函数来对图像进行小波分解和重构。相应地，调用Matlab函数 wfilters 来求取小波函数的分解滤波器组系数[Lo_D,Hi_D]和重构滤波器组系数[Lo_R,Hi_R]。

（2）函数 mydwt2 的输出应改为[cA,cV,cH,cD]来表示，注意这4个部分的位置分布：cA 在系数矩阵的左上角，cV 在右上角，cH 在左下角，cD 在右下角。

（3）函数 mydwt2 的列变换，其处理对象由输入矩阵 x 改为 行变换后的缓存矩阵 x1，列变换后得到缓存矩阵 x2。输出 cA 取 x2 的左上角1/4部分，cV 取 x2 的右上角1/4部分，cH 取 x2 的左下角1/4部分，cD 取 x2 的右下角1/4部分。

（4）函数 mywavedec2 的输出 y 格式改为与Matlab的 wavedec2 类似，但我的是以“细胞矩阵”的格式，而 wavedec2 是以“向量序列”的格式输出小波系数。小波分解的结果显示方式不再使用塔式结构方式，而是用 subplot 函数把各级平均、细节系数显示出来。

（5）重构函数 mywaverec2 根据分解函数 mywavedec2 的输出格式做了相应的修改，程序更简洁、易懂和合理。并且 mywaverec2 的重构级数可以与分解级数不同，能够实现不同级别的小波重构。

（6）函数 myidwt2 的行、列变换中，卷积后要对卷积输出序列进行裁剪，消除分解滤波器和重构滤波器对序列长度的影响，使重构矩阵的大小与上一分解级的系数矩阵大小一致。

以上是主要的修正点，更详细的修改信息请看下一篇文章给出的程序及其注释。