1、1使用 MATLAB 设计小波变换程序中的若干问题在使用 MATLAB 完成小波变换程序和通过阈值来压缩图像的过程中，我和许多同学都是边学边用，是从一个接一个的问题中逐步理解小波和 MATLAB 编写程序的。因此我愿意就个人遇到和解决问题的经验与大家讨论，希望能够对遇到同样问题的人有所帮助。在清华大学林福宗老师倡导的网上互动的学习方式中，老师同学的开诚布公的讨论，尤其是林老师启发大家对出现 问题 采取的态度和做法，对我今后成长为一名合格的清华的研究生意义重大，谨以此文表示对他指导关心的敬意！ 1. 内容简介 本文分为三部分: 如何使用 MATLAB 设计小波标准与标准分解；如何完成使用小波变换
2、压缩图像；仍需探讨的问题。每一部分主要以问题和例子的形式讲述，为了便于参照，附录部分给出部分源代码供大家参考指正。我个人在完成作业的时候，走了许多弯路，最后，反复读老师的第 3 章讲义的第 3-5 节，明白了小波变换的一些非常基础性的知识，正是如此，我主要是按照第三章的例子，做出使用 MATLAB 进行 Haar 变换的实例，这至少对完成非标准 Haar 小波的 3 级分解与合成没有任何问题，而且计算速度远远比使用一维变换快，甚至超过 dwt2 和 wavedec2。再者，是通过这些简短的例子的学习，也能从实践的角度理解小波变换的基础知识和道理，我觉得掌握小波的知识要比简单使用 dwt2 直接
3、分析出结果更重要。为了让使用一维变换或打算采用卷积完成分解重构程序的同学有所参考，也给出了我以前写的代码和设想供有兴趣者参考。在图像压缩的任务中，我个人建议采用 wdencmp，原因是简单可靠，能够生成老师要求小波压缩与重构的演示及PNG 文件。在最后部分着重对使用小波压缩之后重构图像的 PNG 文件会变大问题进行了简单的猜测性的解释和讨论，希望感兴趣的同学能深入研究。2. 如何使用 MATLAB 设计小波标准与标准分解 1. 使用 Haar 和 Db9 小波编写图像文件的标准和非标准分解重构程序 实质上，任务书中规定首先要编写种程序之一：哈尔（haar）小波的标准分解重构程序；Daubech
4、ies 9 小波的标准分解重构程序；哈尔（haar）小波的非标准分解重构程序；Daubechies 9 小波的非标准分解重构程序。如果采用标准方法，需要生成图的系列图像。如果采用非标准方法则需要生成图的系列图像。我个人编写了非标准的分解重构程序，并总结网上同学们公开的方法，认为可以有种方法实现。但前提是必须仔细阅读林老师推荐的补充教材第三章“小波与小波变换 ”的 3.3、3.4 和 3.5 节（页）。如果能使用 MATLAB 简单地实践一下教材的例子，至少完成哈尔小波的标准与非标准程序相当简单。为此我主要介绍如何实现老师教材给我们的例子是如何在MATLAB 中实现的。然后，再引申到使用 MAT
5、LAB 的函数实现方法。为了描述问题简便，使用黑体表示需在 MATLAB 的命令窗口（command window）的输入部分。如果不熟悉下面出现的函数功能和使用方法，请再命令窗口中使用 help 函数名，如help dwt2，能够得到英文的功能和用法说明。2. 哈尔小波变换的 MATLAB 实例 先介绍 MATLAB 的矩阵(Matrix) 和向量(Vector)的赋值方法：在 command window 的符号下输入：一维向量： = 9 7 3 5 = 2, 5, 8, 9, 7, 4, -1, 1二维矩阵： A = 64 2 3 61 60 6 7 57 9 55 54 12 13 5
6、1 50 1617 47 46 20 21 43 42 2440 26 27 37 36 30 31 3332 34 35 29 28 38 39 2541 23 22 44 45 19 18 4849 15 14 52 53 11 10 568 58 59 5 4 62 63 1和现在就是例 3.1例 3.2一维向量, A 就是3.5.1 中的图像矩阵。在进行 Haar 小波的一维与二维变换（分解）之前，首先介绍一下 MATLAB 的矩阵运算。MATLAB 是 Matrix Laboratory 的合写，它对矩阵运算之功能堪称一流。由于使用矩阵描述问题更象数学表达式，所以编写的程序不仅高效，
7、更易读。所以 MATLAB 程序应该尽量使用矩阵直接描述。 如 M 是一 4X4 的矩阵，则B=I*M 就完成了使用两个 for 循环编写乘法的程序。分析老师的例 3.1，那么可以得到如下结果：（9+7）*1/2 （3+5）* 1/2 （9-7）*1/2 （3-5）* 1/2如果把其看作矩阵方式的乘法，那么令 M 为其表述求取平均和差值的系数矩阵，则输入：M = 1/2 0 1/2 021/2 0 -1/2 0 0 1/2 0 1/20 1/2 0 -1/2把上述的 M 输入到 MATLAB 中,然后使用：C=I*M得到的结果就是 8 4 1 -1。这就是非规范化 Haar 小波的第一级分解系
8、数。注意 C 的前两项 8 4 就是近似系数（Approximation Coefficients ）, 后两项1 _1就是细节系数（Detail Coefficients ）. M 就是 Haar 小波的非规范化（non-normalization）系数矩阵。如果我们执行 ：M*M_(M_代表的是 M 的转置矩阵)，你会发现得到对角线为 0.5 其它为 0 的 4X4 矩阵，所以如果让它变成单位矩阵（对角线为1），必须把原来的 1/2 增大，变为 1/sqrt(2), 所以执行：N=M*sqrt(2)后得到新的系数矩阵 0.7071 0 0.7071 00.7071 0 -0.7071 00
9、 0.7071 0 0.70710 0.7071 0 -0.7071再执行N * N_就会得到单位矩阵。N 就是 Haar 小波的规范化（normalization）系数矩阵。现在执行:C*M_将得到 4.5 3.5 1.5 2.5，显然是因为 M*M_的对角矩阵的值为 0.5。现在使用 Cn=I*N得到 11.3137 5.6569 1.4142 -1.4142 规范化的一维 Haar 小波一级分解系数。I1=Cn*N_得到 9 3 7 7 规范化的一维 Haar 小波逆变换（重构）结果。可见得到这样规范化矩阵的好处是为了逆变换时不需要额外的运算，只是需要分解系数矩阵乘以规范化系数矩阵的转置
10、。如果需要对分解系数矩阵进行第 2 级分解，那么此时只对 C或 Cn 的前 2 项的低频系数进行，此时的 Haar 的非规范化与规范化系数矩阵 M1 和 N1 分别为：M1=1/2 1/21/2 1/2N1=1/2 1/21/2 1/2 那么第 2 级分解可以用以下方式进行：C1=C(1:2)*M1Cn1=Cn(1:2)*N1结果是 C1=6 2 , Cn1= 12.0 4.0。注意, C1(1:2) 表示取出C 中第到第个向量，符号“：”的含义是从到的意思，这是 MATLAB 截取数据的方法。逆变换：C1*M1*2 C(3:4)*M*2Cn1*N1 Cn(3:4)*N这里采用了中括号的赋值方
11、式直接将第级的重构合并到第级的系数中再重构原始向量 I。从上面的实验可以看出，除了输入哈尔小波系数矩阵 M 和 N 比较麻烦，整个变换和重构过程异常简单 ,这就启示我们只要编写一个 Haar 系数生成矩阵的函数，整个的分解问题即将化简。1. 问题 1：如何生成 Haar 小波系数矩阵？ MATLAB 的函数可以如以下形势构造：在 MATLAB 中使用 File? New? M-file,输入function M= haarmatrix (rows,cols, flag)%矩阵的行数 rows, 列数 cols, 非规范化 flag=否则就是规范化矩阵 % M 为函数返回的 Haar 系数矩阵M
12、=zeros(rows,cols); %首先制作 rows X cols 的矩阵 if flag= 0 sv=0.5; %non-normalizedelsesv=1.0/sqrt(2.0); % normalizedendhalf=rows/(2); 矩阵从中间分开，左半部为近似系数部分，右半部为细节系数部分for i= 1:2:half*2M(i,ceil(i/2) =sv;M(i+1,ceil(i/2)=sv;M(i,ceil(i/2)+half) =sv;M(i+1,ceil(i/2)+half)=-sv;3end M=sparse(M) %生成稀疏矩阵，提高运算速度注意分号的作用是不
13、在 command 窗口上显示 M 的值。在输入后，请保存与函数名同名的文件名 haarmatrix.m。现在来做例 3.2： = 2, 5, 8, 9, 7, 4, -1, 1N= haarmatrix(8,8,1)C1=F*NN1= haarmatrix(4,4,1)C2= C1(1:4)*N1N2= haarmatrix(2,2,1)C3= C2(1:2)*N2按图 3-21 小波分解的层次结构合成最终系数C=C3, C2(3:4), C1(5:8)得到的 C 就是老师让我们使用第三章 17 页伪代码编写的一维小波变换分解的结果。其逆变换（重构）代码如下：Fr= C3*N2 ,C2(3:
14、4)*N1 ,C1(5:8) * N 注意理解最内部的中扩号实质上是使用 C3 逆变换的结果（低频） 与 C2 的后位细节系数一起逆变换得倒 C1 的前位的低频近似系数，然后与 C1 的后位细节系数构造出原始向量 Fr。为了验证我们的一维 Haar 变换程序，我们使用 MATLAB 提供的一维离散小波变换函数 dwt 来进行同样的计算：第一级分解ca1, cd1=dwt(F, haar )第级分解ca2,cd2=dwt(ca1, haar )第级分解ca3,cd3=dwt(ca2, haar )按图 3-21 小波分解的层次结构合成最终系数c=ca3 cd3 cd2 cd1或直接使用 MATL
15、AB 提供的维多级分解程序c=wavedec(F,3,haar)比较 c 与 C ( 注意 MATLAB 中大小写变量的区别)，可知我们成功构造了 Haar 小波分解合成程序。也许这是你可能觉得刚才的系数矩阵比起 MATLAB 提供的函数来说还是复杂了，但如果把上述过程函数化，做一个 matrixDWT 和 matrixIDWT 函数，使用一样简单，而且关键在于，我们使用矩阵方式解决维图像的小波非标准变换的演示过程来说会带来更大便利。2. 问题 2 如何进行图像的 1 级非标准分解和重构？ 如果我们利用 2.2 节开始时已经输入的 18 页的图像矩阵 A，那么它的第一级二维变换的非标准分解过程
16、只需步： (为了方便核对 MATLAB 的 dwt2 分解结果，使用规范化系数矩阵)M1=haarmatrix(8,8,1) ; AR1=A*M1; % 这就是进逐行分解的图像AC1=M1*AR1; %这就是进逐列分解的图像，也是一次维变换的系数矩阵使用 MATLAB 的 dwt2 来验证ca ch cv cd=dwt2(A,haar)c=ca cv ; ch cd %将个系数合并在一起，注意 cv 与 ch 的位置与第四章图 4-04 区别。对比 ca3, ch3,cv3, cd3 与 AC3，可知矩阵算法正确。细心的同学也许能发现我们的矩阵位置和第章图- (a) 本文图 1解释的不同，这是
17、因为 dwt2 采用的算法与我们不同，但这并不影响小波分解系数的使用。详细解释见附录 1。图 1、多级分解的各个系数矩阵4其重构过程只需要两步：（参照第三章 24 页的公式，但是我们只采用单级重构）AR1=M1*AC1;A=AR1*M1 ;3. 问题 3 如何进行图像的 1 级非标准分解和重构？如果是进行级非标准分解，那么只需要对前一次系数矩阵的的左上角的近似系数图-04（a）中的 a 进行，如对是级分解：M1=haarmatrix(8,8,1) ; AR1=A*M1; % 这就是第 1 次逐行分解的图像AC1=M1*AR1; %这就是第次逐列分解的图像，M2=haarmatrix(4,4,)
18、; AR2=AC1(1:4,1:4)*M2; %第级行分解，只对左上角低频近似系数部分AC2=M2*AR2; %第级列分解M3=haarmatrix(2,2,); AR3=AC2(1:2,1:2)*M3; %第 3 级行分解，只对左上角低频近似系数部分AC3=M3*AR3; %第级列分解执行 MATLAB 提供的维多级小波分解函数 wavedec2C,L=wavedec2(A,3,haar); ca3=appcoef2(C,L,haar,);ch3,cv3,cd3 = detcoef2(all,C,L,3)AC3 相当于上图中(b)左上角的最小的a3 v3 h3 d3，AC2 是c2 v2 h
19、2 d2，AC1 是c1 c2 c3 c4。那么如果合成图 4-04(c)这样的图像只需要把 AC3 替代 AC2 的 a2 ,在把 AC2 替代AC1 的 a1 就可以，相应的代码：AC2(1:2,1:2)=AC3;AC1(1:4,1:4)=AC2;如何由 3 级分解系数矩阵重构出原始矩阵？参照 1 级非标准重构,注意此时从第 3 级向第 1 级重构，即首先由 A3 中的 a h v d 四个系数重构出 A2 的低频部分 a,然后重复这种方法：AR3=M3*AC3; A3=AR3*M3； AC2(1:2,1:2)=A3； %重构的 A3 就是 A2 的低频近似系数，如图 4-04(b)的a
20、h v dAR2=M2*AC2;A2=AR2*M2；AC1(1:4,1:4)=A2； %重构的 A2 就是 A1 的低频近似系数，如图 4-04(a)的 aAR1=M1*AC1;A1=AR1*M1； %A1 就是 A注意：上述分解与重构过程与老师第 3 章中的分解重构过程在一级分解的情况下相同，在第 2 级分解时就发生了变化。上述分解过程使用了与 MATLAB 的 wavedec2 近似的多分辨率分析(multi-resolution analysis), 产生的矩阵系数除了 h 与 v 的位置对调外，二者结果完全相同。然而采用第 3 章 21-24 页介绍的方法（即使是使用 23 页的规范化
21、系数矩阵）不能生成与wavedec2 相同的矩阵，正是因为如此，采用设置阈值来压缩图像时会与使用 wavedec2 的结果不同。原因分析及对应的 24 页的算法的演示程序请参见附录 2。（请林老师务必核对）。3. 三种实现非标准分解的方法： 第一种：2.2.3 介绍的矩阵 A 变成 256X256 的照片矩阵， 1:2 的矩阵下标变成 1:64, 1:4 变成 1:128, 其它不变，保存系数 AR1,AC1,AR2,AC2,AR3,AC3 为 PNG 格式，就是图 3-26 所需要的照片，只不过下 1 级没有显示上一级的细节系数（h v d）部分。使用合成技术类似这样的语句 c=ca cv
22、; ch cd应该很方便得出。当然采用附录 2 中的程序可以更直接，更快速地实现老师的范例。第二种：按照第 3 章伪代码形式使用 MATLAB 一维离散小波变换函数 dwt 编写。dwt 的用法: ca cd = dwt(I, wname); I 是你要分解的一维向量， wname 是字符串，使用单引号扩起来，如db9或haar。返回系数 ca 是低频近似系数， cd 是高频细节系数。 由于 db9 的滤波系数长度为 18，dwt 需要扩展以减少边界误差，重构函数 idwt 能够保证恢复到原来的尺寸。所以可以采用不同的扩展边界模式来使用，缺省值是sym,它在分解后会使 ca 和 cd 的长度变
23、为 (LI + LF-1)/2, LI 是向量 I 的长度，LF 是小波滤波系数的长度。也可以改变边界扩展模式，如使用per,具体方法是ca cd=dwt(I,wname,mode,per).它的 ca 和 cd 的长度为 LI/2,但是它可能带来边界误差。我编写了 2 维非标准分解 nsdwt2 和重构 nsidwt2 函数，可参见附录 3。第三种：直接使用 MATLAB 卷积函数 conv2，即仿照 dwt2 的函数编写。此处不再详细讨论，有兴趣者参见附录 15这三种方式在 MATLAB 下计算速度最快的是第一种，最慢的是第二种，第三种速度与第一种速度比相差不大。在笔者PIII733 上运
24、行分解 wbarb 图像过程所用时间（使用 TIC,TOC）分别为 0.10, 0.28, 7.23 秒。4. 常见问题：如何读取、转换和保存图像文件？X, map = imread(文件名) ; %文件名必须带后缀。使用 size(X)来判断图像矩阵的维数，如果是 256 256 就是 2 维的，那么 map 就是它的色板文件，可以使用 imshow(X,map)来显示它。如果 256 256 3 就是 3 维真彩色图像（True color Image），使用 imshow(X)来显示。一般我们把 2 维的彩色图片称为伪彩色图像（Indexed color image）。一般来说，不管图像
25、文件是不是 PNG 格式， 最好首先使用 MATLAB 的 imwrite 方式重新保存：对于真彩色文件，imwrite（X, 新文件名） ; 对于伪彩色文件，imwrite（X, map, 新文件名 ）; 由 imread 读出的图像文件采用 uint8 的类型，必须转换成 double 类型之后方可进行运算。转换方法： X=doubl(X);经过小波变换和重构的 X 是 double 类型，在保存图像前需要转换成 uint8 类型，如果处理不好，会产生截断误差，许多同学反映阈值为 0 时的文件大小和原始文件大小不一致可能就是出现在这里，转换方法：Y=uint8(round(X); 这里的
26、round 是就近取整，这样 7.99999999 就会变成 8 而不会被截断变为 7。如何处理 Indexed color 文件？原则说，只要可以转换成 double 类型我们就可以开始小波变换。但是由于伪彩色图片的色板 map 保存了每个象素值对应的r g b 的颜色值，因此，如果 map 中的颜色值不连续，就不能保证小波转换的效果。可以使用 imshow(X,map);colorbar;来显示图片的色板是否连续。如果不连续，可以使用 MATLAB 提供的方法转换成灰度图像处理。这是林老师在 “教师答疑”中贴出的MATLAB 帮助文件 Image Conversion，我把它转换成一个小程
27、序，可以修改使用。Anna 指出 ind2gray 和 ind2rgb 也能完成同样的任务。X,map=imread(你的伪彩色图片文件);X=double(X);if min(min(X)=0X=X+1;endimage(X)title(Original Color Indexed Image)colormap(map); colorbarR = map(X,1); R = reshape(R,size(X);G = map(X,2); G = reshape(G,size(X);B = map(X,3); B = reshape(B,size(X);Xrgb = 0.2990*R + 0.
28、5870*G + 0.1140*B;n = 255; % Number of shades in new indexed imageX = round(Xrgb*(n-1) + 1;map2 = gray(n);image(X), title(Processed Gray Scale Indexed Image)colormap(map2), colorbarbaboon= X;map = map2;imwrite(X,map,gray.png); %保存为灰度图像，文件名为gray.png如何处理真彩色文件？有几种方法，在 MATLAB 的 Image Conversion 中有叙述，但许多
29、同学采用了单独处理它的 3 个颜色的 2 维矩阵。设 X为三维矩阵(256,256,3) , X(:,:,1)代表红颜色的 2 维矩阵 X(:,:,2)代表绿 2 维矩阵 , X(:,:,3)代表绿 2 维矩阵。X, map=imread(真彩色文件);r=double(X(:,:,1); %r 是 256 x 256 的红色信息矩阵g=double(X(:,:,2); %g 是 256 x 256 的绿色信息矩阵b=double(X(:,:,3); %b 是 256 x 256 的兰色信息矩阵%开始对 r g b 分别作小波变换，并分别形成各自的小波系数矩阵%把 3 个变换的系数矩阵合并成图
30、像文件Y(:,:,1)=uint8(round( r );Y(:,:,2)=uint8(round( g );Y(:,:,1)=uint8(round( b );可以参照 3.1 的 simplecmp.m 程序来使用。3. 如何完成使用小波变换压缩图像的任务 6由于压缩的目的在于比较使用不同小波压缩后重构图像的失真度视觉效果和使用 PNG 文件保存时的文件大小，如果我们自己编写的小波变换程序存在小的遗漏，可能对压缩结果判断错误，所以我认为至少应当使用 MATLAB 为我们提供的压缩函数记录结果，以便与自己设计的阈值处理程序进行比较。一般来说，两种处理结果应该差别很小，甚至无差别。1. 直接使
31、用 wdencmp: 这是我使用 wdencmp 编写的简单压缩处理程序 simplecmp，可以直接再运行时输出图像、0 的个数、PNG 文件的大小以及计算时间。如果你的图片为 a.png,不管是真彩色、伪彩色，使用：simplecmp(a.png,haar,3) %进行 haar 小波的三级分解压缩合成simplecmp(a.png,db9,3) %进行 db9 小波的三级分解压缩合成下面是具体程序清单。我只对高频细节系数部分进行硬阈值设置。function simplecmp(fname,wname,level)rgb,map=imread(fname);fig=figure;color
32、bar;axis on;axis equal;set(fig,position,10 20 790 580,name,压缩程序演示)if length(size(rgb)=3 rgbcmp(rgb,wname,level);elseindcmp(rgb,map,wname,level);end%-function rgbcmp(rgb,wname,level) %这是压缩真彩色图像rgb=double(rgb);THR=0 5 10 20;otxt=sprintf(%s_zeros.txt,wname);fid=fopen(otxt,w);fprintf(fid,%20s%15s%15s%15
33、sn,文件名 ,大小,阈值,0 数);for T=THRTIC; rgb(:,:,1),cxc,lxc,perf0,perfl2=wdencmp(gbl,rgb(:,:,1),wname,level,T,h,1);num0=length(find(abs(cxc) 1 dofor row 1 to h doDecompositionStep(C row, 1 . . . h)end for% C is changed to Average and Difference% A-Left half, D-right halffor col 1 to h doDecompositionStep(C
34、1 . . . h, col)10end for% C is changed to ca ch,cv,cd, ca-Top Left,% ch-bottom left,cv-top right, cd -bottom righth=h/2end whileend procedure1）两个 DecompositionStep 都处理(C) ，那么 C在 row 处理后，原来的 C 已经改变为行处理后的。应当声明。2）我仔细阅读了教材，并使用 MATLAB 编写上述算法，发现老师没有说明产生了几个系数。这是我分析 MATLAB 的dwt,dwt2,idwt2 之后并亲自编写自己的 nsdwt2
35、和 nsidwt2 时发现的。您的原文“首先对图像的每一行计算像素对的均值和差值，然后对每一列计算像素对的均值和差值。”设平均产生结果为 A,差为 D,试想每一行的平均值是 h/2 个 A,差值是 h/2 个 D,那么仍然存储在 C 中（两者都假设下采样），此时仍然使用 DecompositionStep(C 1 . . . h, col),就会出现 h/4个 AA,h/4 个 AD(A 的均值与差值）,h/4 个 DA,h/4 个DD（A 的均值与差值），这实际上就是说，2 维与 1 维只有A 与 D 两个系数相比，产生了 4 个系数，这也就是使用ca ch cv cd=dwt2(C,wna
36、me)的得到的低频系数，水平、垂直、对角系数。对于最后生成的矩阵 C,左上 ca,左下 ch,右上 cv,右下 cd。虽然您在前面章节介绍了左上角矩阵 AA 为低频，其余均为高频细节系数，但如果使用 dwt 编写上述算法时，很容易只考虑低频部分而丢掉 cv,cd。-dwt2 源程序:前面的源代码略去首先扩展边界，symsizeEXT=(18-1), per18/2y = wextend(2D,dwtEXTM,x,sizeEXT,sizeEXT);在进行低频 Lo_D 系数的每行的过滤得到 z，这相当于行变换，只不过没有下采样，就是上面的 A 部分z = wconv(row,y,Lo_D);然后
37、再用行变换得低频矩阵 z 再一次进行列方向低频滤波卷积下采样，得到的低频部分就是 ca=a,对 z 得高频滤波就是 ch=h,再对比上面对 AA,ADa = convdown(z,Lo_D,sizeKEPT,shift);h = convdown(z,Hi_D,sizeKEPT,shift);在进行高频 Hi_D 系数的每行的过滤得到 z，这相当于行变换，只不过没有下采样，就是上面的 D 部分z = wconv(row,y,Hi_D);然后再用行变换得低频矩阵 z 再一次进行列方向 Lo_D 滤波卷积下采样，得到的高频中的低频部分就是 cv=v,对 z 得高频滤波就是 cd=d 对角部分,再对
38、比上面对 DA,DDv = convdown(z,Lo_D,sizeKEPT,shift);d = convdown(z,Hi_D,sizeKEPT,shift);如果你想使用 dwt2 确得到图 3-26 那样的图像，比使用 dwt要快的 dwt2 中的关键代码的解释： type dwt2.m 与下面的源代码解释相对照，稍加修改，然后换一个函数名称，作业就完成了。可行的方法：对第一个 z 进行下采样 a1=dyaddown(wkeep(z,sizeKEPT)就是图 3-26 的第一次行分解的低频部分，对第 2 个的 z 进行上述过程得到 d1 就是图3-26 的高频部分。然后合成图像：img
39、3_26=uint8(round(a1 d1);这个速度绝对与 dwt2 一致，因为它就是 dwt2.2、采用第 3 章 21-24 页介绍的方法不能生成与 wavedec2 相同的矩阵。原来认为老师第三章的 haar 的 2 级，3 级小波矩阵就只对低频（ca) 部分分解,实际上不是，ch 部分也参加了计算。但是因为矩阵是正交的，重构没有问题，但压缩可能就会出现问题。举例：4X4 A 与 2 级系数 MA= a11 a12 a13 a14 M= 1/2 1/2 0 0a21 a22 a23 a24 1/2 -1/2 0 0a31 a32 a33 a34 0 0 1 0a41 a42 a43
40、a44 0 0 0 1 此时 a11 a12 a21 a22 为 ca; a31 a32 a41 a42 是 ch;A*M= a11*1/2+a12*1/2 a11*1/2-a12*1/2 a13 a14 a21*1/2+a22*1/2 a21*1/2-a22*1/2 a23 a24 -a31*1/2+a32*1/2 a31*1/2-a31*1/2 a33 a34 a41*1/2+a42*1/2 a41*1/2-a41*1/2 a43 a44 由此，ch 参加了运算。下面是根据老师第 3 章 21-24 页编写的标准和非标准分解重构程序。可以看出它分解的高频部分越来越亮 0 的数目不是越来越多
41、。但如果将 X 的赋值 语句的注释符号%去掉，可以看到与老师 3.5 节一样的结果。function showhaardecrec% This program is completely according Professor Linfuzongs additional materials% chapter 3,section 3.5 2-D haar wavelet transform and synthesis.% Here we use the .mat file wbarb which should be in the directory wavedemo% It use function makehaarmatrix make the haar wavelet transform matrix at 3 different level.% and then 1)showing standard decomposition/reconstruction process and display the result
