最新多小波变换的矩阵形式教学课件

《最新多小波变换的矩阵形式教学课件》由会员分享，可在线阅读，更多相关《最新多小波变换的矩阵形式教学课件（101页珍藏版）》请在金锄头文库上搜索。

1、多小波变换的矩阵形式多小波变换的矩阵形式多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 预滤波与后滤波的矩阵表示预滤波与后滤波的矩阵表示 从前面的解释中可以看出，假设多小波的重数是从前面的解释中可以看出，假设多小波的重数是r维的，维的，多小波分解算法要求有初始系数（即多小波分解算法要求有初始系数（即r维向量）才能进行计维向量）才能进行计算，而通常的输入数据是一维的信号，或者是函数算，而通常的输入数据是一维的信号，或者是函数f(t)的等间的等间距采样值，为了应用上面的公式，这就需要将一维数据或者距采样值，为了应用上面的公式，这就需要将一维数据或者f

2、(t)的等间距采样值转化为的等间距采样值转化为r维向量，这个转化的过程一般称维向量，这个转化的过程一般称为预处理或者预滤波。对一维数据进行预滤波后，进行小波为预处理或者预滤波。对一维数据进行预滤波后，进行小波变换即小波分解，然后进行逆变换即小波重构，还得将重构变换即小波分解，然后进行逆变换即小波重构，还得将重构后的后的r维向量转化为一维数据，这个过程就称为后滤波或者维向量转化为一维数据，这个过程就称为后滤波或者后处理。一般来说，小波重构是精确重构，则需要后滤波是后处理。一般来说，小波重构是精确重构，则需要后滤波是预滤波的逆过程。一维多小波变换过程如下：预滤波的逆过程。一维多小波变换过程如下：

3、一维信号一维信号预滤波预滤波分解分解重构重构后滤波后滤波一维信号一维信号研究多小波变换的矩阵表示，首先应该知道预滤波与后滤波研究多小波变换的矩阵表示，首先应该知道预滤波与后滤波的矩阵表示形式。的矩阵表示形式。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小

4、小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 %本程序演示对本程序演示对GHM多小波，使用多小波，使用ghmap方法进行预滤波与后滤波的矩阵计算方法方法进行预滤波与后滤波的矩阵计算方法x=8,12,0,5,20,3,9,16,22,6,55,21;xx=8,12;0,5;20,3;9,16;22,6;55,21;disp(原始数据原始数据x=);disp(x);

5、disp(向量化后的数据向量化后的数据xx=);disp(xx);fp=prep1D_appe(x,ghmap);disp(使用使用prep1D_appe函数对函数对x进行预滤波后的数据进行预滤波后的数据fp=);disp(fp);PR,PO=coef_prep(ghmap);a=PR(:,1:2);b=PR(:,3:4);x1=; for i=1:length(xx)-1 x1=x1,a*xx(:,i)+b*xx(:,i+1); end x1=x1,a*xx(:,length(xx)+b*xx(:,1); disp(直接使用矩阵乘积计算，得到的直接使用矩阵乘积计算，得到的xx的预滤波的结果的

6、预滤波的结果x1=);disp(x1); aa=PO(:,1:2);bb=PO(:,3:4); x2=; x2=x2,aa*x1(:,length(x1)+bb*x1(:,1); for i=1:length(x1)-1 x2=x2,aa*x1(:,i)+bb*x1(:,i+1); end disp(对前面预滤波的数据对前面预滤波的数据x1,使用矩阵乘积方法进行后滤波，得到使用矩阵乘积方法进行后滤波，得到x2=);disp(x2); fp1=postp1D_appe(fp,ghmap); disp(使用使用prep1D_appe函数对函数对x进行后滤波，得到数据进行后滤波，得到数据fp1=);

7、disp(round(fp1);本程序名：本程序名：mdwt_test1.m, 使用使用ghmap预滤波器计算预滤波器计算多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 mdwt_test1原始数据原始数据x= 8 12 0 5 20 3 9 16 22 6 55 21向量化后的数据向量化后的数据xx= 8 0 20 9 22 55 12 5 3 16 6 21使用使用prep1D_appe函数对函数对x进行预滤波后的数据进行预滤波后的数据fp= 12.7279 9.7227 10.3414 22.3623 25.7210 35.2670 0 2

8、0.0000 9.0000 22.0000 55.0000 8.0000直接使用矩阵乘积计算，得到的直接使用矩阵乘积计算，得到的xx的预滤波的结果的预滤波的结果x1= 12.7279 9.7227 10.3414 22.3623 25.7210 35.2670 0 20.0000 9.0000 22.0000 55.0000 8.0000对前面预滤波的数据对前面预滤波的数据x1,使用矩阵乘积方法进行后滤波，得到使用矩阵乘积方法进行后滤波，得到x2= 8.0000 0 20.0000 9.0000 22.0000 55.0000 12.0000 5.0000 3.0000 16.0000 6.0

9、000 21.0000使用使用prep1D_appe函数对函数对x进行后滤波，得到数据进行后滤波，得到数据fp1= 8 12 0 5 20 3 9 16 22 6 55 21多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 %本程序演示对本程序演示对GHM多小波，使用多小波，使用ghmorap方法进行预滤波与后滤波的矩阵计算方法方法进行预滤波与后滤波的矩阵计算方法x=8,12,0,5,20,3,9,16,22,6,55,21;xx=8,12;0,5;20,3;9,16;22,6;55,21;disp(原始数据原始数据x=);disp(x); dis

10、p(向量化后的数据向量化后的数据xx=);disp(xx);fp=prep1D_appe(x,ghmorap);disp(使用使用prep1D_appe函数对函数对x进行预滤波后的数据进行预滤波后的数据fp=);disp(fp);PR,PO=coef_prep(ghmorap);a=PR(:,1:2);b=PR(:,3:4);x1=; for i=1:length(xx)-1 x1=x1,a*xx(:,i)+b*xx(:,i+1); end x1=x1,a*xx(:,length(xx)+b*xx(:,1); disp(直接使用矩阵乘积计算，得到的直接使用矩阵乘积计算，得到的xx的预滤波的结果

11、的预滤波的结果x1=);disp(x1); aa=PO(:,1:2);bb=PO(:,3:4); x2=; x2=x2,aa*x1(:,length(x1)+bb*x1(:,1); for i=1:length(x1)-1 x2=x2,aa*x1(:,i)+bb*x1(:,i+1); end disp(对前面预滤波的数据对前面预滤波的数据x1,使用矩阵乘积方法进行后滤波，得到使用矩阵乘积方法进行后滤波，得到x2=);disp(x2); fp1=postp1D_appe(fp,ghmorap); disp(使用使用prep1D_appe函数对函数对x进行后滤波，得到数据进行后滤波，得到数据fp1

12、=);disp(round(fp1);本程序名：本程序名：mdwt_test2.m, 使用使用ghmorap预滤波器计算预滤波器计算多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 mdwt_test2原始数据原始数据x= 8 12 0 5 20 3 9 16 22 6 55 21向量化后的数据向量化后的数据xx= 8 0 20 9 22 55 12 5 3 16 6 21使用使用prep1D_appe函数对函数对x进行预滤波后的数据进行预滤波后的数据fp= 12.8245 5.9335 5.7146 17.9971 11.1835 27.7171

13、 -1.2411 19.2250 6.7448 20.2496 51.5994 5.1272直接使用矩阵乘积计算，得到的直接使用矩阵乘积计算，得到的xx的预滤波的结果的预滤波的结果x1= 12.8245 5.9335 5.7146 17.9971 11.1835 27.7171 -1.2411 19.2250 6.7448 20.2496 51.5994 5.1272对前面预滤波的数据对前面预滤波的数据x1,使用矩阵乘积方法进行后滤波，得到使用矩阵乘积方法进行后滤波，得到x2= 8.0000 0.0000 20.0000 9.0000 22.0000 55.0000 12.0000 5.000

14、0 3.0000 16.0000 6.0000 21.0000使用使用prep1D_appe函数对函数对x进行后滤波，得到数据进行后滤波，得到数据fp1= 8 12 0 5 20 3 9 16 22 6 55 21多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 %本程序演示对本程序演示对Sa4多小波，使用多小波，使用sa4ap方法进行预滤波与后滤波的矩阵计算方法方法进行预滤波与后滤波的矩阵计算方法x=8,12,0,5,20,3,9,16,22,6,55,21;xx=8,12;0,5;20,3;9,16;22,6;55,21;disp(原始数据原始

15、数据x=);disp(x); disp(向量化后的数据向量化后的数据xx=);disp(xx);fp=prep1D_appe(x,sa4ap);disp(使用使用prep1D_appe函数对函数对x进行预滤波后的数据进行预滤波后的数据fp=);disp(fp);PR,PO=coef_prep(sa4ap);x1=; for i=1:length(xx) x1=x1,PR*xx(:,i); enddisp(直接使用矩阵乘积计算，得到的直接使用矩阵乘积计算，得到的xx的预滤波的结果的预滤波的结果x1=);disp(x1); x2=;for i=1:length(x1)-1 x2=x2,PO*x1(

16、:,i);enddisp(对前面预滤波的数据对前面预滤波的数据x1,使用矩阵乘积方法进行后滤波，得到使用矩阵乘积方法进行后滤波，得到x2=);disp(x2);fp1=postp1D_appe(fp,sa4ap);disp(使用使用prep1D_appe函数对函数对x进行后滤波，得到数据进行后滤波，得到数据fp1=);disp(round(fp1);本程序名：本程序名：mdwt_test3.m, 使用使用sa4ap预滤波器计算预滤波器计算多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 mdwt_test3原始数据原始数据x= 8 12 0 5 2

17、0 3 9 16 22 6 55 21向量化后的数据向量化后的数据xx= 8 0 20 9 22 55 12 5 3 16 6 21使用使用prep1D_appe函数对函数对x进行预滤波后的数据进行预滤波后的数据fp= 14.1421 3.5355 16.2635 17.6777 19.7990 53.7401 2.8284 3.5355 -12.0208 4.9497 -11.3137 -24.0416直接使用矩阵乘积计算，得到的直接使用矩阵乘积计算，得到的xx的预滤波的结果的预滤波的结果x1= 14.1421 3.5355 16.2635 17.6777 19.7990 53.7401 2

18、.8284 3.5355 -12.0208 4.9497 -11.3137 -24.0416对前面预滤波的数据对前面预滤波的数据x1,使用矩阵乘积方法进行后滤波，得到使用矩阵乘积方法进行后滤波，得到x2= 8.0000 0 20.0000 9.0000 22.0000 12.0000 5.0000 3.0000 16.0000 6.0000使用使用prep1D_appe函数对函数对x进行后滤波，得到数据进行后滤波，得到数据fp1= 8 12 0 5 20 3 9 16 22 6 55 21多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 通过前面的

19、通过前面的3个程序，可以看出，预滤波和后滤波的计算可以个程序，可以看出，预滤波和后滤波的计算可以用矩阵表示，我们看一个例子。假设我们使用用矩阵表示，我们看一个例子。假设我们使用GHM多小波的多小波的双正交双正交GHMAP预滤波与后滤波方法，对预滤波与后滤波方法，对12个数据的一维信号个数据的一维信号s进行计算，预滤波与后滤波变换矩阵可以表示为进行计算，预滤波与后滤波变换矩阵可以表示为:多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 %本程序演示对本程序演示对GHM多小波，使用多小波，使用ghmap方法进行预滤波与后滤波的矩阵计算方法方法进行预滤波

20、与后滤波的矩阵计算方法s=8,12,0,5,20,3,9,16,22,6,55,21;disp(原始数据原始数据x=);disp(s);PR,PO=coef_prep(ghmap);a=PR(:,1:2);b=PR(:,3:4);c=zeros(size(a);%构造预滤波变换矩阵构造预滤波变换矩阵pr_matpr_mat=a,b,c,c,c,c;c,a,b,c,c,c;c,c,a,b,c,c;c,c,c,a,b,c;c,c,c,c,a,b;b,c,c,c,c,a; s1=pr_mat*s;disp(直接使用矩阵乘积直接使用矩阵乘积pr_mat*s计算，得到计算，得到s的预滤波的结果的预滤波的

21、结果s1=);disp(s1); fp=prep1D_appe(s,ghmap);disp(MWMP中的相同预滤波方法计算得中的相同预滤波方法计算得fp,结果如下，可以看出，将结果如下，可以看出，将s1向向量化后就是量化后就是fp);disp(fp);%构造后滤波变换矩阵构造后滤波变换矩阵po_matx=PO(:,1:2);y=PO(:,3:4); z=zeros(size(x);po_mat=y,z,z,z,z,x;x,y,z,z,z,z;z,x,y,z,z,z;z,z,x,y,z,z;z,z,z,x,y,z;z,z,z,z,x,y;s2=po_mat*s1;disp(直接使用矩阵乘积直接使

22、用矩阵乘积po_mat*s1计算，得到的计算，得到的s1的后滤波的结果的后滤波的结果s2=);disp(round(s2);本程序名：本程序名：mdwt_test4.m, 使用使用ghmap预滤波器计算预滤波器计算多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 mdwt_test4原始数据原始数据x= 8 12 0 5 20 3 9 16 22 6 55 21直接使用矩阵乘积直接使用矩阵乘积pr_mat*s计算，得到计算，得到s的预滤波的结果的预滤波的结果s1=12.7279 0 9.7227 20.0000 10.3414 9.0000 22.

23、3623 22.0000 25.7210 55.0000 35.2670 8.0000MWMP中的相同预滤波方法计算得中的相同预滤波方法计算得fp,结果如下，可以看出，将结果如下，可以看出，将s1向量化后向量化后就是就是fp12.7279 9.7227 10.3414 22.3623 25.7210 35.2670 0 20.0000 9.0000 22.0000 55.0000 8.0000直接使用矩阵乘积直接使用矩阵乘积po_mat*s1计算，得到的计算，得到的s1的后滤波的结果的后滤波的结果s2= 8 12 0 5 20 3 9 16 22 6 55 21 注意，如果将上面程序中的注意，

24、如果将上面程序中的ghm小波换成小波换成sa4小波，使用小波，使用sa4ap方法进行计算，那么预滤波矩阵方法进行计算，那么预滤波矩阵P和后滤波矩阵和后滤波矩阵Q是什么样子是什么样子的形式？的形式？显然，它们都是对角矩阵。显然，它们都是对角矩阵。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 %本程序演示对本程序演示对SA4多小波，使用多小波，使用sa4ap方法进行预滤波与后滤波的矩阵计算方法方法进行预滤波与后滤波的矩阵计算方法s=8,12,0,5,20

25、,3,9,16,22,6,55,21;disp(原始数据原始数据s=);disp(s);PR,PO=coef_prep(sa4ap);x=PR;y=PO;z=zeros(size(x);%构造预滤波变换矩阵构造预滤波变换矩阵pr_matpr_mat=x,z,z,z,z,z;z,x,z,z,z,z;z,z,x,z,z,z;z,z,z,x,z,z;z,z,z,z,x,z;z,z,z,z,z,x; s1=pr_mat*s;disp(直接使用矩阵乘积直接使用矩阵乘积pr_mat*s计算，得到的计算，得到的s的预滤波的结果的预滤波的结果s1=);disp(s1); fp=prep1D_appe(s,sa

26、4ap);disp(MWMP中的相同预滤波方法计算得中的相同预滤波方法计算得fp,结果如下，可以看出，将结果如下，可以看出，将s1向量化后就是向量化后就是fp);disp(fp);%构造后滤波变换矩阵构造后滤波变换矩阵po_matpo_mat=y,z,z,z,z,z;z,y,z,z,z,z;z,z,y,z,z,z;z,z,z,y,z,z;z,z,z,z,y,z;z,z,z,z,z,y;s2=po_mat*s1;disp(直接使用矩阵乘积直接使用矩阵乘积po_mat*s1计算，得到的计算，得到的s1的后滤波的结果的后滤波的结果s2=);disp(round(s2);从前一个幻灯片知道，如果预滤波

27、和后滤波器的长度为从前一个幻灯片知道，如果预滤波和后滤波器的长度为1，则，则预滤波变换矩阵预滤波变换矩阵P和后滤波变换矩阵和后滤波变换矩阵Q都是对角矩阵，而都是对角矩阵，而SA4多小波的预滤波和后滤波器的长度为多小波的预滤波和后滤波器的长度为1。下面的程序名为。下面的程序名为mdwt_test5.m，演示对，演示对12个数据个数据s，使用此方法计算的一段程，使用此方法计算的一段程序。序。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 mdwt_test5原始数据原始数据s= 8 12 0 5 20 3 9 16 22 6 55 21直接使用矩阵乘

28、积直接使用矩阵乘积pr_mat*s计算，得到的计算，得到的s的预滤波的结果的预滤波的结果s1= 14.1421 2.8284 3.5355 3.5355 16.2635 -12.0208 17.6777 4.9497 19.7990 -11.3137 53.7401 -24.0416MWMP中的相同预滤波方法计算得中的相同预滤波方法计算得fp,结果如下，可以看出，将结果如下，可以看出，将s1向量化后向量化后就是就是fp 14.1421 3.5355 16.2635 17.6777 19.7990 53.7401 2.8284 3.5355 -12.0208 4.9497 -11.3137 -2

29、4.0416直接使用矩阵乘积直接使用矩阵乘积po_mat*s1计算，得到的计算，得到的s1的后滤波的结果的后滤波的结果s2= 8 12 0 5 20 3 9 16 22 6 55 21 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 function prmat,pomat=pmatrix(len,pflt)PR,PO=coef_prep(pflt); i,j=size(PR); r=min(i,j);z=zeros(r,r);k=max(i/r,j/r);n=len/r; pr1=PR,zeros(r,len-r*k);po1=zeros(r,

30、len-r*k),PO; prmat=;pomat=;for qq=1:n prmat=prmat;pr1;pr1=circshift(pr1,r); po1=circshift(po1,r);pomat=pomat;po1;end假设进行预滤波或者后滤波变换的数据长度为假设进行预滤波或者后滤波变换的数据长度为len,预滤波方法预滤波方法pflt就是就是MWMP软件包中所定义的方法，则我编写的函数软件包中所定义的方法，则我编写的函数 p, q=pmatrix(len,pflt)返回返回nXn的预滤波矩阵的预滤波矩阵p及后滤波矩阵及后滤波矩阵q, 如果变换数据为如果变换数据为s,则计则计算算s1

31、=p*s，就是对信号，就是对信号s进行了预滤波，结果是进行了预滤波，结果是s1,计算计算s2=q*s1,则则s2就是对经过预滤波的信号就是对经过预滤波的信号s1进行了后滤波，请注意，此程进行了后滤波，请注意，此程序使用了序使用了coef_prep.m中的预滤波与后滤波方法，你要安装中的预滤波与后滤波方法，你要安装MWMP软件包或者将软件包或者将coef_prep.m拷贝到与你运行的程序在同拷贝到与你运行的程序在同一目录下，或者是一目录下，或者是MATLAB的的WORK子目录中。子目录中。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 p,q=pm

32、atrix(10,sa4ap);s=(p*8,12,0,5,20,3,9,16,22,6);(q*s)ans = 8.0 12.0 0 5.0 20.0 3.0 9.0 16.0 22.0 6.0 p,q=pmatrix(10,ghmap);s=(p*8,12,0,5,20,3,9,16,22,6);(q*s) ans = 8.0 12.0 0 5.0 20.0 3.0 9.0 16.0 22.0 6.0 p,q=pmatrix(10,ghmorap);s=(p*8,12,0,5,20,3,9,16,22,6);(q*s) ans = 8.0 12.0 0 5.0 20.0 3.0 9.0 1

33、6.0 22.0 6.0 fp=prep1D_appe(8,12,0,5,20,3,9,16,22,6,ghmorap) fp = 12.8245 5.9335 5.7146 17.9971 8.9024 -1.2411 19.2250 6.7448 20.2496 6.0699 p,q=pmatrix(10,ghmorap);s=(p*8,12,0,5,20,3,9,16,22,6) s = 12.8245 -1.2411 5.9335 19.2250 5.7146 6.7448 17.9971 20.2496 8.9024 6.0699 reshape(s,2,5) ans = 12.82

34、45 5.9335 5.7146 17.9971 8.9024 -1.2411 19.2250 6.7448 20.2496 6.0699下面使用下面使用10个数据，验证程序的运行情况：个数据，验证程序的运行情况：多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 到此为止，我们已经知道，可以使用自编子程序到此为止，我们已经知道，可以使用自编子程序pmatrix，获，获得多小波变换的预滤波变换矩阵和后滤波变换矩阵。我们知得多小波变换的预滤波变换矩阵和后滤波变换矩阵。我们知道，多小波变换的过程一般为道，多小波变换的过程一般为 一维信号一维信号预滤波预滤

35、波分解分解重构重构后滤波后滤波一维信号一维信号如果也能够将小波变换的主要步骤，即如果也能够将小波变换的主要步骤，即 分解分解重构重构写成矩阵变换形式的表达式，那么，结合预滤波变换矩阵和写成矩阵变换形式的表达式，那么，结合预滤波变换矩阵和后滤波变换矩阵的表达式，我们就可以将多小波变换写成一后滤波变换矩阵的表达式，我们就可以将多小波变换写成一个统一的矩阵变换形式。个统一的矩阵变换形式。 基于此，下面我们考虑，进行基于此，下面我们考虑，进行1次小波分解与重构的矩阵表次小波分解与重构的矩阵表示形式。下面的程序，需要使用示形式。下面的程序，需要使用MWMP中的文件中的文件coef.m，请，请安装安装MW

36、MP软件包，或者将文件软件包，或者将文件coef.m拷贝到你目前的程序拷贝到你目前的程序目录中，或者是目录中，或者是MATLAB的的WORK子目录中。子目录中。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 先看看先看看coef.m中多小波滤波器中多小波滤波器H_k,G_k的输入格式：的输入格式：通过比较，如果你通过比较，如果你在论文中看到一个在论文中看到一个新小波，你就知道新小波，你就知道怎么加入到怎么加入到coef.m文件中去了。文件中去了。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 实际计

37、算中，一般多小波的重数实际计算中，一般多小波的重数r=2，s=0, 此时加细方程为此时加细方程为上面的加细方程是小波与多小波一书的写法，此时应为上面的加细方程是小波与多小波一书的写法，此时应为则小波变换的低频变换矩阵则小波变换的低频变换矩阵L与高频变换矩阵与高频变换矩阵H为为(对重数对重数w=2)：多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 function L_matrix,H_matrix=mdwt_matrix(len,flt)%在文件在文件coef.m中，低频滤波器系与高频滤波器系要相同，中，低频滤波器系与高频滤波器系要相同，%如果不

38、相同，要使用如果不相同，要使用0矩阵，补足到长度相同矩阵，补足到长度相同L,H=coef(flt);i,j=size(L); r=min(i,j);z=zeros(r,r);k=max(i/r,j/r);n=len/r;L1=L,zeros(r,len-r*k);H1=H,zeros(r,len-r*k);L_matrix=;H_matrix=;for qq=1:n/2 L_matrix=L_matrix;L1;L1=circshift(L1,2*r); H_matrix=H_matrix;H1;H1=circshift(H1,2*r);endreturn下面的函数调用下面的函数调用L, H=

39、mdwt_matrix(len,flt)，对于，对于coef.m中的中的小波小波flt,预滤波后的数据长度为预滤波后的数据长度为len的数据，就能够返回的数据，就能够返回1次小波次小波变换的低频变换矩阵变换的低频变换矩阵L与高频变换矩阵与高频变换矩阵H。如果经过预滤波后的。如果经过预滤波后的数据即向量为数据即向量为s1，则矩阵乘法运算，则矩阵乘法运算 s2=L*s1, 或者或者 s2=H*s1就是经过就是经过1次小波变换后的低频或者高频系数，自编函数如下：次小波变换后的低频或者高频系数，自编函数如下：多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现

40、 n=16;x=1:n;x1=prep1D_appe(x,ghmorap);x2=dec1D_pe(x1,ghm,1) x2 = 4.7834 11.3154 17.8474 24.0827 0.0000 0.0000 0.0350 4.9339 5.6918 10.3106 14.9644 4.7474 0.0000 0.0000 0.0494 -6.0871 n=16;p,q=pmatrix(n,ghmorap);s=1:n;s1=(p*s);LL,HH=mdwt_matrix(n,ghm);s2=HH*s1;s2 ans = 0.0000 0.0000 0.0000 0.0000 0.0

41、350 0.0494 4.9339 -6.0871 reshape(s2,2,n/4) ans = 0.0000 0.0000 0.0350 4.9339 0.0000 0.0000 0.0494 -6.0871 n=16;p,q=pmatrix(n,ghmorap);s=1:n;s1=(p*s);LL,HH=mdwt_matrix(n,ghm);s2=LL*s1;s2 ans = 4.7834 5.6918 11.3154 10.3106 17.8474 14.9644 24.0827 4.7474 reshape(s2,2,n/4) ans = 4.7834 11.3154 17.8474

42、 24.0827 5.6918 10.3106 14.9644 4.7474 n=16;x=1:n;x1=prep1D_appe(x,sa4ap);x2=dec1D_pe(x1,sa4,1) x2 = 9.0000 17.0000 25.0000 17.0000 -0.0000 -0.0000 -0.0000 0.0000 1.8730 1.8730 1.8730 -13.6190 0 -0.0000 0.0000 -4.0000 n=16;p,q=pmatrix(n,sa4ap);s=1:n;s1=(p*s);LL,HH=mdwt_matrix(n,sa4);s2=LL*s1;s2 ans

43、= 9.0000 1.8730 17.0000 1.8730 25.0000 1.8730 17.0000 -13.6190 reshape(s2,2,n/4) ans = 9.0000 17.0000 25.0000 17.0000 1.8730 1.8730 1.8730 -13.6190 n=16;p,q=pmatrix(n,sa4ap);s=1:n;s1=(p*s);LL,HH=mdwt_matrix(n,sa4);s2=HH*s1;s2 ans = 0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -4.0000 reshape

44、(s2,2,n/4) ans = 0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 -4.0000下面是此函数调用的演示结果：下面是此函数调用的演示结果：多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 假设要变换的初始数据假设要变换的初始数据x=x1,x2,xnT的元素个数是的元素个数是n=2m个个, 在这种情况下，对数据进行第在这种情况下，对数据进行第1次小波变换可表示为：次小波变换可表示为： S1=L1*P*x, D1=H1*P*x其中其中S1、D1的元素个数是的元素个数是n/2个，个，

45、S1是低频系数，是低频系数，D1是高频系是高频系数，数，L1是低频变换矩阵，是低频变换矩阵，H1是高频变换矩阵，是高频变换矩阵，P是预滤波变换是预滤波变换矩阵，具体调用形式为：矩阵，具体调用形式为： p,q=pmatrix(n,pflt); L1,H1=mdwt_matrix(n,flt); S1=(L1*p*x); D1=H1*p*x;看看下面的例子看看下面的例子: n=16;x=rand(n,1); xans = 0.4451 0.9318 0.4660 0.4186 0.8462 0.5252 0.2026 0.6721 0.8381 0.0196 0.6813 0.3795 0.831

46、8 0.5028 0.7095 0.4289 p,q=pmatrix(n,sa4ap); L1,H1=mdwt_matrix(n,sa4);S1=(L1*p*x); S1ans = 1.1462 0.1607 0.9511 0.0112 1.1626 0.1618 1.1897 0.1574 D1=(H1*q*x); D1ans = 0.2400 -0.6664 0.0009 0.0774 -0.3677 0.4197 -0.0936 -0.1867 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 如果变换数据如果变换数据x的长度为的长度为n,对

47、经过对经过 p,q=pmatrix(n,pflt); L1,H1=mdwt_matrix(n,flt); S1=(L1*p*x); D1=H1*p*x;变换的小波系数，小波的逆变换可以写成变换的小波系数，小波的逆变换可以写成 n=16;x=(1:n);p,q=pmatrix(n,sa4ap); L1,H1=mdwt_matrix(n,sa4);S1=L1*p*x;D1=H1*p*x; xx=q*L1;H1*S1;D1; xxans = 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 11.0000 1

48、2.0000 13.0000 14.0000 15.0000 16.0000 下面对下面对x=1,2,16,使用使用sa4小波，先进行小波，先进行sa4ap格式的预滤格式的预滤波，然后正变换波，然后正变换1次，然后进行反变换，最后进行后滤波，请次，然后进行反变换，最后进行后滤波，请注意注意MATLAB中构造矩阵时中构造矩阵时a,b与与a;b的区别：的区别：多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 由上面的推导过程，得出进行由上面的推导过程，得出进行1次小波变换的统一形式：次小波变换的统一形式： n=16;x=(1:n);p,q=pmatri

49、x(n,sa4ap);L1,H1=mdwt_matrix(n,sa4); x1=L1;H1*p*x; x1ans = 9.0000 1.8730 17.0000 1.8730 25.0000 1.8730 17.0000 -13.6190 0 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0 -4.0000 x2=prep1D_appe(x,sa4ap);x3=dec1D_pe(x2,sa4,1)x3 = 9.0000 17.0000 25.0000 17.0000 -0.0000 -0.0000 -0.0000 0.0000 1.8730 1.8730 1.8

50、730 -13.6190 0 -0.0000 0.0000 -4.0000 reshape(x3,1,n)ans = 9.0000 1.8730 17.0000 1.8730 25.0000 1.8730 17.0000 -13.6190 -0.0000 0 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -4.0000下面的程序对下面的程序对x=1,2,16，使用我们的算法与，使用我们的算法与MWMP算法分别计算了算法分别计算了1次，结果是相同的。次，结果是相同的。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 下面

51、考虑进行下面考虑进行2次小波变换的矩阵表示形式，首先要说明的是次小波变换的矩阵表示形式，首先要说明的是函数函数p,q=mdwt_matrix(n,flt)返回的矩阵返回的矩阵p和和q都是都是(n/2)Xn的矩的矩阵，如果原始数据阵，如果原始数据x的长度为的长度为n，则第，则第1次变换时，调用参数是次变换时，调用参数是p1,q1=mdwt_matrix(n,flt)，再进行第再进行第2次变换时，此时只对次变换时，此时只对第第1次变换的低频进行变换，因而数据量减少了一半，从而使次变换的低频进行变换，因而数据量减少了一半，从而使用用p2,q2= mdwt_matrix(n/2,flt), 返回矩阵返

52、回矩阵p2,q2是是(n/4)X(n/2)的大小。进行第的大小。进行第2次小波变换的推导公式如下：次小波变换的推导公式如下：多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 n=16;x=(1:n);p,q=pmatrix(n,sa4ap); L2,H2=mdwt_matrix(n/2,sa4);L1,H1=mdwt_matrix(n,sa4); A=L2*L1;H2*L1;H1;x1=A*p*x;x1ans = 30.8882 6.8465 17.1951 -6.8465 1.5881 -1.0607 9.3663 6.7175 0 -0.000

53、0 -0.0000 0.0000 -0.0000 -0.0000 0 -4.0000 x2=prep1D_appe(x,sa4ap);x3=dec1D_pe(x2,sa4,2);x4=reshape(x3,n,1); x4ans = 30.8882 6.8465 17.1951 -6.8465 1.5881 -1.0607 9.3663 6.7175 -0.0000 0 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -4.0000 y=q*A*x1; yans = 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8

54、.0000 9.0000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000 16.0000 注意，对于正交插值，由于注意，对于正交插值，由于P的逆矩阵等于的逆矩阵等于P的转置，因此对于正交插值有的转置，因此对于正交插值有结果结果Q=P，例如，例如ghmap,clap就不是正交插值就不是正交插值(是双正交插值是双正交插值)，此时，此时Q虽然等虽然等于于P的逆矩阵，但是的逆矩阵，但是Q不等于不等于P的转置，其它方法例如的转置，其它方法例如sa4ap，ghmorap是正是正交插值，此时交插值，此时Q等于等于P的转置。的转置。下面是进行下面是进行2次小波变

55、换正变换后，又进行次小波变换正变换后，又进行2次反变换的演示程序，同时也将变换结果与次反变换的演示程序，同时也将变换结果与MWMP进行了比较，结果相同，进行了比较，结果相同，由此看出，我们编写的程序完全可以代替由此看出，我们编写的程序完全可以代替MWMP软件包。软件包。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 上面是进行上面是进行2次小波变换的整体变换，实际上，你也可以考虑次小波变换的整体变换，实际上，你也可以考虑不进行整体变换，而是直接计算变换的各个子带，例如第不进行整体变换，而是直接计算变换的各个子带，例如第1次次变换的低频部分是变换的

56、低频部分是L1*P*x, 高频是高频是H1*P*x，第，第2次变换后的低次变换后的低频是频是L2*L1*P*x, 高频是高频是H2*L1*P*x，演示程序如下：，演示程序如下： n=16;x=(1:n);p,q=pmatrix(n,sa4ap); L2,H2=mdwt_matrix(n/2,sa4);L1,H1=mdwt_matrix(n,sa4); A=L2*L1;H2*L1;H1;x1=A*p*x;x1 %经过经过2次小波变换的整体结果次小波变换的整体结果ans = 30.8882 6.8465 17.1951 -6.8465 1.5881 -1.0607 9.3663 6.7175 0

57、-0.0000 -0.0000 0.0000 -0.0000 -0.0000 0 -4.0000(L1*p*x) %第第1次变换的低频次变换的低频ans = 9.0000 1.8730 17.0000 1.8730 25.0000 1.8730 17.0000 -13.6190 (H1*p*x) %第第1次变换的高频次变换的高频ans = 0 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0 -4.0000 (L2*L1*p*x) %第第2次变换的低频次变换的低频ans = 30.8882 6.8465 17.1951 -6.8465 (H2*L1*p*x) %

58、第第2次变换的高频次变换的高频ans = 1.5881 -1.0607 9.3663 6.7175多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 下面考虑进行下面考虑进行3次小波变换的矩阵表示形式，推导公式如下：次小波变换的矩阵表示形式，推导公式如下：多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 n=32;x=(1:n);p,q=pmatrix(n,sa4ap); L2,H2=mdwt_matrix(n/2,sa4);L1,H1=mdwt_matrix(n,sa4); L3,H3=mdwt_mat

59、rix(n/4,sa4); B=L3*L2*L1;H3*L2*L1;H2*L1;H1; A=B*p; x1=A*x; x1ans = 92.7399 19.8490 39.2601 -0.4841 8.8490 -7.3750 10.51597.3750 -0.0000 0.1796 -0.0000 0.1796 3.1763 -2.300918.7326 13.2554 0 -0.0000 -0.0000 0.0000 -0.0000-0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 0 -0.0000 0.0000 -0.0000 -8.0000 x

60、2=prep1D_appe(x,sa4ap);x3=dec1D_pe(x2,sa4,3); x4=reshape(x3,n,1); x4ans = 92.7399 19.8490 39.2601 -0.4841 8.8490 -7.3750 10.51597.3750 0.0000 0.1796 0.0000 0.1796 3.1763 -2.300918.7326 13.2554 -0.0000 0 -0.0000 -0.0000 -0.00000.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.00000 0.0000 0.0000 -8.0000

61、 x5=q*B*x1; x5ans =1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.00008.0000 9.0000 10.0000 11.0000 12.0000 13.0000 14.000015.0000 16.0000 17.0000 18.0000 19.0000 20.0000 21.000022.0000 23.0000 24.0000 25.0000 26.0000 27.0000 28.000029.0000 30.0000 31.0000 32.0000下面的程序演示了进行下面的程序演示了进行3次小波分解与重构，并将结果与次小波分

62、解与重构，并将结果与MWMP软件包进行比较，其结果是相同的。软件包进行比较，其结果是相同的。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 关于关于1维多小波变换的子带分割问题维多小波变换的子带分割问题我们知道，我们知道，1维多小波变换每次只对低频子带分解，每次将低维多小波变换每次只对低频子带分解，每次将低频子带进一步分成频子带进一步分成1个低频子带和个低频子带和1个高频子带，这个结论不管个高频子带，这个结论不管对单小波，还是多小波，都是适用的。对单小

63、波，还是多小波，都是适用的。 n=16;x=linspace(1,2,n);x1=x.3;x2=dec1D_pe(x1,d4,2)x2 = 2.8567 5.4068 9.2309 12.9056 -0.1183 -0.1429 0.1488 -3.2119 -0.0183 -0.0205 -0.0226 -0.0248 -0.0270 -0.0292 -0.0313 -2.4779上面是对上面是对16个数据个数据x=1,1+h,1+2h,2，其中，其中h=1/16，使用，使用单小波单小波d4（注：就是小波与多小波一书的（注：就是小波与多小波一书的d2）进行）进行2次次变换，数据被变换为变换，

64、数据被变换为3个部分，红色线是低频，兰色线是次个部分，红色线是低频，兰色线是次高频，绿色线是高频。如果同样的数据，使用高频，绿色线是高频。如果同样的数据，使用SA4进行进行2次多次多小波变换，小波变换，MWMP软件结果为：软件结果为： x3=prep1D_appe(x1,sa4ap);x4=dec1D_pe(x3,sa4,2)x4 = 13.1461 8.3499 0.7737 4.3342 -0.0021 -0.0025 -0.0030 -0.1524 3.6121 -4.1068 -0.5161 3.1652 0.0012 0.0012 0.0012 -1.8786低频部分低频部分次高频部

65、分次高频部分高频部分高频部分多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 p,q=pmatrix(n,sa4ap); L1,H1=mdwt_matrix(n,sa4); L2,H2=mdwt_matrix(n/2,sa4); A=L2*L1;H2*L1;H1; A=L2*L1;H2*L1;H1;x5=A*p*x1; x5ans = 13.1461 3.6121 8.3499 -4.1068 0.7737 -0.5161 4.3342 3.1652 -0.0021 0.0012 -0.0025 0.0012 -0.0030 0.0012 -0.1

66、524 -1.8786 reshape(x5,2,n/2)ans = 13.1461 8.3499 0.7737 4.3342 -0.0021 -0.0025 -0.0030 -0.1524 3.6121 -4.1068 -0.5161 3.1652 0.0012 0.0012 0.0012 -1.8786对于这对于这16个数据，使用我们的变换方法也有同样的结果：个数据，使用我们的变换方法也有同样的结果：低频部分低频部分次高频部分次高频部分高频部分高频部分实际上，对于实际上，对于r=2的多小波变换，上面的每个子带部分都应该的多小波变换，上面的每个子带部分都应该更细致地划分为更细致地划分为2个部

67、分，每个部分的第个部分，每个部分的第1行的频率比第行的频率比第2行的行的频率低。这样，对于频率低。这样，对于1维多小波变换，维多小波变换，1次变换分成了次变换分成了4个子带，个子带，2次变换分成次变换分成6个子带，个子带，3次变换分成次变换分成8个子带。个子带。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 假定有假定有n个数据，这个数据，这n个数据经过个数据经过2次多小波变换，即进行了次多小波变换，即进行了 一维信号一维信号预滤波预滤波2次小波分解次小波分解的过程后，例如信号的过程后，例如信号n=16,MATLAB生成程序为：生成程序为： n=

68、16;x=linspace(1,2,n);x1=x.3;使用使用MWMP进行变换，得到结果：进行变换，得到结果： n=16;x=linspace(1,2,n);x1=x.3;x3=prep1D_appe(x1,sa4ap);x4=dec1D_pe(x3,sa4,2)x4 = 13.1461 8.3499 0.7737 4.3342 -0.0021 -0.0025 -0.0030 -0.1524 3.6121 -4.1068 -0.5161 3.1652 0.0012 0.0012 0.0012 -1.8786低频部分低频部分次高频部分次高频部分高频部分高频部分 p,q=pmatrix(n,sa

69、4ap); L1,H1=mdwt_matrix(n,sa4); L2,H2=mdwt_matrix(n/2,sa4); A=L2*L1;H2*L1;H1; A=L2*L1;H2*L1;H1;x5=A*p*x1; x5ans = 13.1461 3.6121 8.3499 -4.1068 0.7737 -0.5161 4.3342 3.1652 -0.0021 0.0012 -0.0025 0.0012 -0.0030 0.0012 -0.1524 -1.8786或者是我们自编子程序的计算结果：或者是我们自编子程序的计算结果：低频部分低频部分次高频部分次高频部分高频部分高频部分多多 小小 波波

70、变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 x4 = 13.1461 8.3499 0.7737 4.3342 -0.0021 -0.0025 -0.0030 -0.1524 3.6121 -4.1068 -0.5161 3.1652 0.0012 0.0012 0.0012 -1.8786低频部分低频部分次高频部分次高频部分高频部分高频部分ans = 13.1461 3.6121 8.3499 -4.1068 0.7737 -0.5161 4.3342 3.1652 -0.0021 0.0012 -0.0025 0.0012 -0.0030 0.0012 -

71、0.1524 -1.8786低频部分低频部分次高频部分次高频部分高频部分高频部分对上面的对上面的2种变换的不同数据形式，即种变换的不同数据形式，即或者是自编子程序的变换结果：或者是自编子程序的变换结果：我们想将小波变换系数，按照多小波变换更细致的频带划分，我们想将小波变换系数，按照多小波变换更细致的频带划分，将系数排列成由低频到高频的形式，例如数据将系数排列成由低频到高频的形式，例如数据x经过经过2 次变换次变换后，被分成后，被分成6个子带，我们想让这个子带，我们想让这6个子带以从低频到高频的个子带以从低频到高频的顺序，形成一个顺序，形成一个1维数组，对于上面数据，就是维数组，对于上面数据，就

72、是13.1461 8.3499 3.6121 -4.1068 0.7737 4.3342 -0.5161 3.1652 -0.0021 -0.0025 -0.0030 -0.1524 0.0012 0.0012 0.0012 -1.8786多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 调用自编子程序调用自编子程序new_coef=mdwt_coef(old_coef,m)就可以达到要就可以达到要求。假设一维信号求。假设一维信号x经过经过m次多小波变换后的系数是次多小波变换后的系数是old_coef，其中，多小波变换可以使用其中，多小波变换可以使

73、用MWMP的方法，对重数的方法，对重数r=2情况，情况，它生成一个它生成一个2X(n/2)的数组，也可以使用我们编写的方法，此时的数组，也可以使用我们编写的方法，此时形成一个形成一个n个元素的一维数组。此函数对于以上个元素的一维数组。此函数对于以上2种情形下生成种情形下生成的的1维或者维或者2维数组维数组old_coef，根据多小波的变换次数，根据多小波的变换次数m，重新构，重新构造系数排列顺序，将造系数排列顺序，将old_coef的系数如刚才的例子一样，排列成：的系数如刚才的例子一样，排列成：多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多

74、 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 自编子程序自编子程序new_coef=mdwt_coef(old_coef,m)就按照这种小波系就按照这种小波系数排列原则，构造了一个向量，它的大小与小波变换前相同。数排列原则，构造了一个向量，它的大小与小波变换前相同。下面是当下面是当m=2, n=16时，此程序的使用例子。时，此程序的使用例子。 n=16;x=linspace(1,2,n);x1=x.3; x3=prep1D_appe(x1,sa4ap);x4=dec1D_pe(x3,sa4,2)x4 = 13.1461 8.3499 0.7737

75、4.3342 -0.0021 -0.0025 -0.0030 -0.1524 3.6121 -4.1068 -0.5161 3.1652 0.0012 0.0012 0.0012 -1.8786 mdwt_coef(x4,2)ans = 13.1461 8.3499 3.6121 -4.1068 0.7737 4.3342 -0.5161 3.1652 -0.0021 -0.0025 -0.0030 -0.1524 0.0012 0.0012 0.0012 -1.8786 p,q=pmatrix(n,sa4ap); L1,H1=mdwt_matrix(n,sa4); L2,H2=mdwt_ma

76、trix(n/2,sa4); A=L2*L1;H2*L1;H1; A=L2*L1;H2*L1;H1;x5=A*p*x1; x5ans = 13.1461 3.6121 8.3499 -4.1068 0.7737 -0.5161 4.3342 3.1652 -0.0021 0.0012 -0.0025 0.0012 -0.0030 0.0012 -0.1524 -1.8786 mdwt_coef(x5,2)ans = 13.1461 8.3499 3.6121 -4.1068 0.7737 4.3342 -0.5161 3.1652 -0.0021 -0.0025 -0.0030 -0.1524

77、 0.0012 0.0012 0.0012 -1.8786通过例子可以看出，通过例子可以看出，mdwt_coef函数可以将经过函数可以将经过MWMP变换变换或者我们自编子程序变换的小波系数排列成低频到高频的形或者我们自编子程序变换的小波系数排列成低频到高频的形式。另外，此程序只适用于重数式。另外，此程序只适用于重数r=2的多小波。的多小波。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 下面是一维多小波变换的一个自编函数，若下面是一维多小波变换的一个自编函数，若old_wcoef是经过是经过trans_number次多小波变换后的小波系数，变换可

78、以使用次多小波变换后的小波系数，变换可以使用MWMP软件包中的方法，也可以是我们提出的方法，则这些软件包中的方法，也可以是我们提出的方法，则这些系数经过系数经过coef=mdwt_coef(old_coef,trans_number)进行频带划进行频带划分后，分后，coef中系数按照多小波变换更细致的频带划分规则，将中系数按照多小波变换更细致的频带划分规则，将系数排列成由低频到高频的形式，其中按行向量方式是由左到系数排列成由低频到高频的形式，其中按行向量方式是由左到右，或者按列向量方式是由上到下，频率由低到高，最左边或右，或者按列向量方式是由上到下，频率由低到高，最左边或者最上面的者最上面的2

79、个部分是变换后的低频部分，最右边或者最下边个部分是变换后的低频部分，最右边或者最下边的的2个部分是变换后的高频部分，中间是变换后的次高频部分，个部分是变换后的高频部分，中间是变换后的次高频部分，此时，调用函数：此时，调用函数： subband=coef1D(wcoef,trans_number,parti)可以得到多小波可以得到多小波trans_number次变换后的第次变换后的第parti个子带个子带subband，其中，其中parti是整数且在允许范围内，例如是整数且在允许范围内，例如trans_number=1，parti=1,2,3,4。 trans_number=2，parti=1,

80、2,3,4,5,6。trans_number=3，parti=1,2,7,8。 trans_number=4，parti=1,2,9,10。 trans_number=5，parti=1,2,11,12。 trans_number=6，parti=1,2,13,14。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 %程序名：程序名： mdwt_test6.mn=16;x=linspace(1,2,n);x1=x.3;p,q=pmatrix(n,sa4ap); L1,H1=mdwt_matrix(n,sa4);L2,H2=mdwt_matrix(n

81、/2,sa4); A=L2*L1;H2*L1;H1; A=L2*L1;H2*L1;H1;x5=A*p*x1; reshape(x5,2,n/2)x6=mdwt_coef(x5,2)coef1D(x6,2,4)coef1D(x6,2,5)coef1D(x6,2,1) mdwt_test6ans =13.1461 8.3499 0.7737 4.3342 -0.0021 -0.0025 -0.0030 -0.1524 3.6121 -4.1068 -0.5161 3.1652 0.0012 0.0012 0.0012 -1.8786x6 =13.1461 8.3499 3.6121 -4.1068

82、 0.7737 4.3342 -0.51613.1652 -0.0021 -0.0025 -0.0030 -0.1524 0.0012 0.0012 0.0012 -1.8786ans = -0.5161 3.1652ans = -0.0021 -0.0025 -0.0030 -0.1524ans = 13.1461 8.3499下面的程序对下面的程序对16个数据进行个数据进行2次多小波变换，然后使用次多小波变换，然后使用mdwt_coef函数按频率高低排列小波系数，最后使用函数按频率高低排列小波系数，最后使用coef1D函函数提取了数提取了3个不同子带的小波系数：个不同子带的小波系数：多多

83、小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 通过上面的解释，可知，若通过上面的解释，可知，若x是多小波变换前的一维数据，是多小波变换前的一维数据，A是变换矩阵，是变换矩阵，p是预滤波方法，是预滤波方法，m是变换次数，则是变换次数，则y=A*p*x就就是多小波变换系数，但这些系数通过这种计算方法分成了是多小波变换系数，但这些系数通过这种计算方法分成了m+1个子带，其中个子带，其中1个低频子带，个低频子带，m个高频子带。而每个子带个高频子带。而每个子带又可分成又可分成2个子带个子带(针对重数针对重数r=2)，将，将y=A*p*x的结果更细致地的结果更细

84、致地分成分成2(m+1)个子带，前面我们使用了个子带，前面我们使用了dwt_coef自编子程序将自编子程序将每个子带划分为每个子带划分为2个，但如果将此程序用矩阵乘法运算代替，个，但如果将此程序用矩阵乘法运算代替，例如矩阵例如矩阵T，则，则y1=T*A*p*x就可实现就可实现dwt_coef的功能。的功能。 实际上，对于每个经过实际上，对于每个经过y=A*p*x计算后的计算后的(m+1)个子带中的个子带中的一个，如果要进一步划分成为一个，如果要进一步划分成为2个子带，只是将这个子带中的个子带，只是将这个子带中的所有奇数下标的系数放到前面，所有偶数下标的系数放到后所有奇数下标的系数放到前面，所有

85、偶数下标的系数放到后面，就可达到要求。例如：面，就可达到要求。例如：多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 若若n是向量是向量x的个数，的个数，trans_number是小波变换次数，则函数是小波变换次数，则函数mat=smatrix(n,trans_number)返回此小波变换的置换矩阵返回此小波变换的置换矩阵T,从而计算从而计算y=T*A*p*x

86、，就将经过小波变换的，就将经过小波变换的trans_number+1个子带划分为个子带划分为2*(trans_number+1)个子带，若个子带，若n=16, trans_number=1, 则则T的形式为：的形式为： smatrix(16,1)ans = 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0

87、 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

88、 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 若若n=16, trans_number=2, 则则T的形式为：的形式为： smatrix(16,2)ans = 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0

89、 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0

90、 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 %程序名：程序名：mdwt_test7.mn=16;x=linspace(1,2,n);x1=x.3;x3=prep1D_appe(x1,sa4ap);x4=dec1D_pe(x3,sa4,2);x5=reshape(x4,1,n);mdwt_coef(x5,2)p,q=pmatrix(n,sa4ap);L1,H1=mdwt_matrix(n,sa4);L2,H2=mdwt_matrix(

91、n/2,sa4);A=L2*L1;H2*L1;H1;x6=smatrix(n,2)*A*p*x1; x6 mdwt_test7ans = 13.1461 8.3499 3.6121 -4.1068 0.7737 4.3342 -0.5161 3.1652 -0.0021 -0.0025 -0.0030 -0.1524 0.0012 0.0012 0.0012 -1.8786ans = 13.1461 8.3499 3.6121 -4.1068 0.7737 4.3342 -0.5161 3.1652 -0.0021 -0.0025 -0.0030 -0.1524 0.0012 0.0012 0

92、.0012 -1.8786下面的程序磁盘文件名为下面的程序磁盘文件名为mdwt_test7.m，程序分成，程序分成2部分，第部分，第1部分使用部分使用MWMP方法进行方法进行2次小波变换，同时应用次小波变换，同时应用mdwt_coef调整小波系数。第调整小波系数。第2部分使用我们提出的方法进行部分使用我们提出的方法进行2次小波变换，次小波变换，其结果是一致的。其结果是一致的。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 通过上面的解释，可知，若通过上面的解释，可知，若x是多小波变换前的一维数据，是多小波变换前的一维数据，A是变换矩阵，是变换矩阵

93、，p是预滤波方法，是预滤波方法，q是后滤波方法，是后滤波方法，m是变换次数，是变换次数，T是置换矩阵，是置换矩阵，y=T*A*p*x经过经过m次多小波变换的小波系数，次多小波变换的小波系数，若若p为正交插值预滤波矩阵，则为正交插值预滤波矩阵，则x=p*A*T*y, 即即(T*A*p)或者或者q*A*T是逆变换的变换矩阵。若是逆变换的变换矩阵。若p为非正交为非正交(例如双正交插值例如双正交插值预滤波预滤波ghmap,clap方法方法)插值预滤波矩阵插值预滤波矩阵, q*A*T是逆变换的是逆变换的变换矩阵。下面程序是变换矩阵。下面程序是mdwt_test9.m，运行结果如下：，运行结果如下：关于关

94、于1维多小波变换的逆变换的矩阵表示形式：维多小波变换的逆变换的矩阵表示形式：%程序名：程序名：mdwt_test9.mn=16;x1=1:16;p,q=pmatrix(n,sa4ap);L1,H1=mdwt_matrix(n,sa4);L2,H2=mdwt_matrix(n/2,sa4);A=L2*L1;H2*L1;H1;T=smatrix(n,2);x6=T*A*p*x1; x7=p*A*T*x6; x7 mdwt_test9ans = 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 11.0000

95、 12.0000 13.0000 14.0000 15.0000 16.0000多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 一个实际例子：一个实际例子：1维信号的去噪试验维信号的去噪试验 下面信号去噪方面的试验。首先是介绍几个下面信号去噪方面的试验。首先是介绍几个MATLAB下的小下的小波工具箱（波工具箱（Wavelet ToolBox）的）的2个函数，即去噪阈值选择函个函数，即去噪阈值选择函数数thselect.m和小波系数收缩函数和小波系数收缩函数wthresh.m，这些函数已经拷，这些函数已经拷贝到当前目录下。基于贝到当前目录下。基于D

96、onoho所提出的去噪阈值选择函数所提出的去噪阈值选择函数: thr = thselect(x,tptr)对于含噪向量或者含噪矩阵对于含噪向量或者含噪矩阵X=X+，去噪方法，去噪方法tptr, 函数函数thr = thselect(x,tptr)返回其含噪向量或者含噪矩阵返回其含噪向量或者含噪矩阵X的小波收缩算法的小波收缩算法的软阈值或者硬软阈值，其中的软阈值或者硬软阈值，其中是服从是服从N(0,1)分布的高斯白噪声，分布的高斯白噪声，tptr=sqtwolog是根据是根据=sqrt(2*ln(lenth(x)选择阈值选择阈值(VisulShrink方法方法,一般称为一般称为通用阈值通用阈值)

97、, tptr =rigrsure是基于是基于Stein无偏估计方法选择阈值无偏估计方法选择阈值(SureShrink方法方法, Stein无偏估计阈无偏估计阈值值), tptr=heursure是结合以上两种方法的改进是结合以上两种方法的改进(HeurShrink方方法法,启发式阈值启发式阈值), tptr =minimaxi (minimaxi方法方法，最大极小方最大极小方差阈值差阈值)。注意含噪信号是。注意含噪信号是X=X+，如果为，如果为X=X+，即，即X中中噪声符合噪声符合N(0, )分布，则应该将返回值分布，则应该将返回值thr乘上乘上。多多 小小 波波 变变 换换 的的 矩矩 阵阵

98、 形形 式式 及及 其其 软软 件件 实实 现现 小波域内的系数收缩函数：小波域内的系数收缩函数： Y = wthresh(X,SORH,T)对于小波子带对于小波子带X，按照阈值，按照阈值T进行系数收缩，其中进行系数收缩，其中SORH=s是软阈值法，是软阈值法，SORH=t是硬阈值法。解释如下：是硬阈值法。解释如下： 当当SORH=t时，如果时，如果X(i)或者或者X(i,j)小于小于T，则返回的，则返回的Y(i)或或者者Y(i,j)=0，否则返回，否则返回Y(i)或者或者Y(i,j)=X(i,j)。 当当SORH=h时，如果时，如果X(i)或者或者X(i,j)小于小于T，则返回的，则返回的Y

99、(i)或者或者Y(i,j)=0, 否则返回的否则返回的Y(i)或者或者Y(i,j)= signY(i)或者或者Y(i,j)*abs(Y(i)或者或者Y(i,j)-T。 注意，对于上面注意，对于上面T的值，正是函数的值，正是函数thr = thselect(x,tptr)的的计算结果计算结果thr, 再乘上信号含有噪声的标准差的大小。如果噪声再乘上信号含有噪声的标准差的大小。如果噪声的标准差不知道，可以采用统计方法进行估计。的标准差不知道，可以采用统计方法进行估计。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 sigma=3;n=512;t=li

100、nspace(0,6*pi,n);x=(10*sin(t);randn(seed,10000);y=x+sigma*randn(size(x);subplot(1,2,1);plot(t,y);pflt=ghmap;flt=ghm;p,q=pmatrix(n,pflt);method=sqtwolog; L1,H1=mdwt_matrix(n,flt);L2,H2=mdwt_matrix(n/2,flt);L3,H3=mdwt_matrix(n/4,flt);A=L3*L2*L1;H3*L2*L1;H2*L1;H1;T=smatrix(n,3);B=T*A*p; y1=B*y;thr = 0.

101、3*sigma*thselect(y1,method); y2=wthresh(y1,s,thr);BB=T*A*q;y3=BB*y2;subplot(1,2,2);plot(t,y3);fprintf(降噪前信噪比降噪前信噪比SNR=%4.4f, 降噪后信噪比降噪后信噪比SNR=%4.4f. n, snr(x,y), snr(x,y3);下面的程序是下面的程序是mdwt_test21.m,它对一个正弦信号，加入高斯噪它对一个正弦信号，加入高斯噪声后进行声后进行3次小波变换，然后将小波系数看成是一个整体子带，次小波变换，然后将小波系数看成是一个整体子带，在小波域内使用在小波域内使用sqtwol

102、og方法选择一个通用的阈值进行小波的方法选择一个通用的阈值进行小波的系数收缩后，再进行小波逆变换，得到去噪后的信号，程序最系数收缩后，再进行小波逆变换，得到去噪后的信号，程序最后画出了去噪前与去噪后的图像，并计算了去噪前、后的信噪后画出了去噪前与去噪后的图像，并计算了去噪前、后的信噪比。注意，程序中应用了自编函数比。注意，程序中应用了自编函数snr.m计算信噪比，此程序计算信噪比，此程序已经拷贝到当前目录下。已经拷贝到当前目录下。 mdwt_test21降噪前信噪比降噪前信噪比SNR=7.4846, 降噪后信噪比降噪后信噪比SNR=12.9462.多多 小小 波波 变变 换换 的的 矩矩 阵阵

103、 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 我们知道，我们知道，3次小波变换将信号次小波变换将信号y分解为分解为4个大的子带，即有个大的子带，即有y1=L3*L2*L1;H3*L2*L1;H2*L1;H1*p*y, 注意，这个变换注意，这个变换没有乘以置换矩阵没有乘以置换矩阵T，乘以置换矩阵，乘以置换矩阵T的目地是将每个子带的的目地是将每个子带的奇数下标的系数与偶数下标的系数集中到一起，进一步分解奇数下标的系数与偶数下标的系数集中到一起，进一步分解为为2个子带。个子带。 在不乘以置换矩阵在不乘

104、以置换矩阵T的情况下，对于信号的的情况下，对于信号的3次小波变换，次小波变换，可以认为分解为可以认为分解为4个子带，其中，个子带，其中，L3*L2*L1*p*y是低频部分，是低频部分，H3*L2*L1*p*y，H2*L1*p*y，H1*p*y是高频子带，对每个是高频子带，对每个子带单独进行通用阈值的系数收缩算法，然后将收缩后的系子带单独进行通用阈值的系数收缩算法，然后将收缩后的系数重新组合成小波系数矩阵，最后使用这个重新组合的系数数重新组合成小波系数矩阵，最后使用这个重新组合的系数矩阵进行小波逆变换，就得到了去噪后的信号。矩阵进行小波逆变换，就得到了去噪后的信号。 下面的程序下面的程序mdwt

105、_test22.m，实现的就是上面所提及的。，实现的就是上面所提及的。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 sigma=3;n=512;t=linspace(0,6*pi,n);x=(10*sin(t);randn(seed,10000);y=x+sigma*randn(size(x);subplot(1,2,1);plot(t,y);pflt=sa4ap;flt=sa4;p,q=pmatrix(n,pflt);method=minimaxi; L1,H1=mdwt_matrix(n,flt);L2,H2=mdwt_matrix(n/2

106、,flt);L3,H3=mdwt_matrix(n/4,flt);A1=L3*L2*L1*p*y; A2=H3*L2*L1*p*y; A3=H2*L1*p*y; A4=H1*p*y; r=0.75;thr1=r*sigma*thselect(A1,method); thr2=r*sigma*thselect(A2,method);thr3=r*sigma*thselect(A3,method); thr4=r*sigma*thselect(A4,method);A1=wthresh(A1,s,thr1); A2=wthresh(A2,s,thr2);A3=wthresh(A3,s,thr3);

107、 A4=wthresh(A4,s,thr4);y1=A1;A2;A3;A4; A=L3*L2*L1;H3*L2*L1;H2*L1;H1; B=A*q; y2=B*y1; subplot(1,2,2);plot(t,y2);fprintf(降噪前降噪前SNR=%4.4f, 降噪后降噪后SNR=%4.4f.n,snr(x,y),snr(x,y2);下面的程序是下面的程序是mdwt_test22.m，与，与mdwt_test22.m不同之处在不同之处在于：（于：（1）将小波改为）将小波改为sa4。（。（2)将系数收缩方法变为极小极大将系数收缩方法变为极小极大方法方法minimaxi。（。（3）将）将

108、mdwt_test22.m在小波系数域中的一在小波系数域中的一个收缩阈值改为在个收缩阈值改为在4个子带中，每个子带中有一个单独的收缩个子带中，每个子带中有一个单独的收缩阈值。阈值。 mdwt_test22降噪前降噪前SNR=7.4846, 降噪后降噪后SNR=13.8072.多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 对于对于1维信号的维信号的3次小波变换，如果将变换后的小波系数细分为次小波变换，如果将变换后的小波系数细分为8个子带，对每个子带使用

109、单独的阈值进行系数收缩，那么，个子带，对每个子带使用单独的阈值进行系数收缩，那么，怎么做呢？可以使用自编函数怎么做呢？可以使用自编函数coef1D。若。若wcoef是经过是经过m次多次多小波变换后的小波系数，变换使用小波变换后的小波系数，变换使用y=T*A*p*x的形式，则经过的形式，则经过此变换后，系数按照多小波变换更细致的频带划分规则，将系此变换后，系数按照多小波变换更细致的频带划分规则，将系数排列成由低频到高频的形式，此时，调用函数：数排列成由低频到高频的形式，此时，调用函数： subband=coef1D(wcoef, m, parti)可以得到多小波可以得到多小波m次变换后的第次变换

110、后的第parti个子带个子带subband，其中，其中parti是整数且在允许范围内，例如是整数且在允许范围内，例如trans_number=1，parti=1,2,3,4。 trans_number=2，parti=1,2,3,4,5,6。trans_number=3，parti=1,2,7,8。 trans_number=4，parti=1,2,9,10。 trans_number=5。 注意，变换一定是注意，变换一定是y=T*A*p*x形式，形式， y=A*p*x形式不能使形式不能使用此函数，例如用此函数，例如mdwt_test22.m程序就不能使用这个函数。程序就不能使用这个函数。 对

111、取出的子带对取出的子带subband进行系数收缩，然后将收缩后的子进行系数收缩，然后将收缩后的子带按顺序组合并进行逆变换，就得到了去噪后的信号。带按顺序组合并进行逆变换，就得到了去噪后的信号。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 sigma=3;n=512;t=linspace(0,6*pi,n);x=(10*sin(t);randn(seed,10000);y=x+sigma*randn(size(x);subplot(1,2,1);plot(t,y);pflt=ghmap;flt=ghm;p,q=pmatrix(n,pflt);me

112、thod=sqtwolog; L1,H1=mdwt_matrix(n,flt);L2,H2=mdwt_matrix(n/2,flt);L3,H3=mdwt_matrix(n/4,flt);A=L3*L2*L1;H3*L2*L1;H2*L1;H1; T=smatrix(n,3);B=T*A*p; y1=B*y; y2=; r=1;for i=1:8 if i2 sub=coef1D(y1,3,i); thr=r*sigma*thselect(sub,method); sub=wthresh(sub,s,thr); y2=y2; sub; else sub=coef1D(y1,3,i); y2=y

113、2; sub; endendB=T*A*q; y3=B*y2; subplot(1,2,2);plot(t,y3);fprintf(降噪前降噪前SNR=%4.4f, 降噪后降噪后SNR=%4.4f.n,snr(x,y),snr(x,y3);下面的程序名是下面的程序名是mdwt_test23.m，对于，对于3次小波变换的次小波变换的8个子带，个子带，它只对小波变换后的第它只对小波变换后的第3到到8个子带进行了系数收缩。个子带进行了系数收缩。 mdwt_test23降噪前降噪前SNR=7.4846, 降噪后降噪后SNR=17.5336.多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及

114、及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 二维多小波变换的矩阵形式及实现方法二维多小波变换的矩阵形式及实现方法 通过上面的解释，可知，若通过上面的解释，可知，若x是多小波变换前的一维数据，是多小波变换前的一维数据，A是变换矩阵，是变换矩阵，p是预滤波方法，是预滤波方法，m是变换次数，是变换次数，T是置换矩阵，是置换矩阵，则运算则运算y=T*A*p*x经过经过m次多小波变换的结果。次多小波变换的结果。 如果如果x是是nXn的二维数据，可先逐列（逐行）做的二维数据，可先逐列（逐行）做1维变换维变换y=T*A*p

115、*x=B*x(其中其中B=T*A*p) ，然后再逐行（逐列）做，然后再逐行（逐列）做1维维变换变换y*p*A*T=y*B，这样就完成了一次，这样就完成了一次2维小波变换维小波变换,两次两次变换合起来就是变换合起来就是B*x*B。第。第1次次2维小波变换，将维小波变换，将x分成了分成了4X4=16个部分，左上角的个部分，左上角的4个部分是变换后的低频部分，设此个部分是变换后的低频部分，设此部分为部分为x1。第。第2次次2维小波变换只对维小波变换只对x1的这的这4个部分进行变换，个部分进行变换，其它的部分不改变。经过变换后，这其它的部分不改变。经过变换后，这4个部分的每一个部分，个部分的每一个部分

116、，又被分成又被分成4个部分，共个部分，共4X4=16个部分，其中这个部分，其中这16个部分中，左个部分中，左上角的上角的4个子带，是第个子带，是第2次变换后的低频部分，设为次变换后的低频部分，设为x2。第。第3次次变换只对变换只对x2的这的这4个子带进行变换。经过第个子带进行变换。经过第3次小波变换后，次小波变换后，x2的这的这4个子带又分成了个子带又分成了16个部分，其中左上角的个部分，其中左上角的4个部分是低频个部分是低频部分，设为部分，设为x3。第。第4次变换只对次变换只对x3变换，依此类推。变换，依此类推。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件

117、件 实实 现现 因此因此, 若若m是变换次数，是变换次数，A是一维变换矩阵，是一维变换矩阵，p是预滤波方法，是预滤波方法，T是是m次变换时的置换矩阵，则运算次变换时的置换矩阵，则运算T*A*p*x*p*A*T就是就是对二维信号对二维信号x做做m次次2维多小波变换维多小波变换, 若设若设B=T*A*p, 则可知经则可知经过过m次变换后的小波系数矩阵为次变换后的小波系数矩阵为y=B*x*B。下面是对。下面是对16X16的的信号信号x=(xi,j)=(i+j)做做1次小波变换的例子。此程序先使用次小波变换的例子。此程序先使用ghm小波对小波对x做做1次变换，然后对次变换，然后对x使用使用ghm小波，

118、应用我们提出的小波，应用我们提出的方法做方法做1次变换，并对次变换，并对2种变换结果相减后取绝对值种变换结果相减后取绝对值 。从。从x4的的结果几乎为零可看出，结果几乎为零可看出，2种计算方法是等价的。种计算方法是等价的。 %程序名：程序名：mdwt_test8.mn=16;for i=1:n for j=1:n x(i,j)=i+j; end; end; x1=prep2D_appe(x,ghmap);x2=dec2D_pe(x1,ghm,1);p,q=pmatrix(n,ghmap);L1,H1=mdwt_matrix(n,ghm);A=L1;H1;x3=smatrix(n,1)*A*p*

119、x*p*A*smatrix(n,1);x4=x3-x2; sum(sum(abs(x4) mdwt_test8ans = 7.6014e-013多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 上面的程序上面的程序mdwt_test8.m中的中的x有对称性，并且有对称性，并且x是是16X16的尺的尺寸太小，与实际应用的尺寸有差别，现在我们将上面的程序中寸太小，与实际应用的尺寸有差别，现在我们将上面的程序中的的x改动一下，使用改动一下，使用rand函数，生成一个函数，生成一个512X512大小的信号大小的信号x, 其它与上面的程序相同，注意，其它与上

120、面的程序相同，注意，rand(seed,10000)是固定生成是固定生成随机序列的初值，使得程序能够生成固定的随机序列。另外，随机序列的初值，使得程序能够生成固定的随机序列。另外，程序比较了我们的算法与程序比较了我们的算法与MWMP算法的计算速度（单位是秒）算法的计算速度（单位是秒），我们的算法比，我们的算法比MWMP算法快算法快1倍以上。倍以上。%程序名：程序名：mdwt_test11.mn=512;rand(seed,10000);x=rand(n,n);t0=cputime;x1=prep2D_appe(x,ghmap);x2=dec2D_pe(x1,ghm,1);t1=cputime-

121、t0; t2=cputime;p,q=pmatrix(n,ghmap);L1,H1=mdwt_matrix(n,ghm);A=L1;H1;x3=smatrix(n,1)*A*p*x*p*A*smatrix(n,1);t3=cputime-t2; x4=x3-x2; sum(sum(abs(x4)fprintf(我们程序用时我们程序用时:%4.4f, MWMP用时用时:%4.4fn, t3, t1); mdwt_test11ans = 3.0715e-011我们程序用时我们程序用时:3.1560, MWMP用时用时:6.8910 我们提出的算我们提出的算法速度快于法速度快于MWMP方法，方法，但

122、是，由于两但是，由于两种算法计算公种算法计算公式不一样，因式不一样，因此，当进行多此，当进行多次分解时，两次分解时，两种算法的计算种算法的计算结果，在各个结果，在各个子带的边界附子带的边界附近，有可能相近，有可能相同位置的小波同位置的小波系数不相等。系数不相等。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 n=16; for i=1:n for j=1:n x(i,j)=i+j; end; end; x1=prep2D_appe(x,clap);x2=dec2D_pe(x1,cl,2);p,q=pmatrix(n,clap);L1,H1=mdw

123、t_matrix(n,cl);L2,H2=mdwt_matrix(n/2,cl);A=L2*L1;H2*L1;H1;T=smatrix(n,2);x3=T*A*p*x*p*A*T;x4=x3-x2; s=sum(sum(abs(x4);fprintf(MWMP方法与我们方法的差的绝对值之和为方法与我们方法的差的绝对值之和为:%2.4fn,s); mdwt_test10MWMP方法与我们方法的差的绝对值之和为方法与我们方法的差的绝对值之和为:67.8823那么，我们所提出的这种算法，在计算结果上，与那么，我们所提出的这种算法，在计算结果上，与MWMP软软件包的计算结果是完全一致的吗？下面来看一段

124、程序，程序名件包的计算结果是完全一致的吗？下面来看一段程序，程序名是是mdwt_test10.m, 此程序使用此程序使用CL2小波计算，对小波计算，对x=(xi,j)=(i+j),先用先用MWMP软件做了软件做了2次小波变换，然后又用次小波变换，然后又用我们的方法做了我们的方法做了2次小波变换，最后计算了两个结果的差的绝次小波变换，最后计算了两个结果的差的绝对值之和。对值之和。此时，这两种算法并不完全一样，为什么？个人认为，这是此时，这两种算法并不完全一样，为什么？个人认为，这是由于我们的方法不进行边界延拓，而由于我们的方法不进行边界延拓，而MWMP在边界上自动进在边界上自动进行周期延拓，因而

125、在边界附近，结果不相同。行周期延拓，因而在边界附近，结果不相同。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 由于小波变换是一个全局变换，由于边界的延拓值在逐次变换由于小波变换是一个全局变换，由于边界的延拓值在逐次变换过程中的传递性，有时候，可能我们的方法与过程中的传递性，有时候，可能我们的方法与MWMP方法的方法的小波变换系数，相差会更多，或者这在小波变换系数，相差会更多，或者这在2种不同变换结果在相种不同变换结果在相同位置上，可能有更多的系数不相等同位置上，可能有更多的系数不相等, 或者相同位置上的小波或者相同位置上的小波系数，两种算法的计

126、算结果相差较大的值。例如下面的系数，两种算法的计算结果相差较大的值。例如下面的mdwt_test12.m程序，就有这样的特点。程序，就有这样的特点。%程序名：程序名：mdwt_test12.mn=32;rand(seed,10000);x=rand(n,n); x1=prep2D_appe(x,ghmap);x2=dec2D_pe(x1,ghm,2);p,q=pmatrix(n,ghmap);L1,H1=mdwt_matrix(n,ghm);L2,H2=mdwt_matrix(n/2,ghm);A=L2*L1;H2*L1;H1;T=smatrix(n,2);x3=T*A*p*x*p*A*T;x

127、4=x3-x2; s=sum(sum(abs(x4);fprintf(MWMP方法与我们方法的差的绝对值之和为方法与我们方法的差的绝对值之和为:%2.4fn,s); mdwt_test12MWMP方法与我们方法的差的绝对值之和为方法与我们方法的差的绝对值之和为:141.3453多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 而且，滤波器的长度大的小波，这种传递性就越快，传递的强而且，滤波器的长度大的小波，这种传递性就越快，传递的强度可能也越大，例如上面的度可能也越大，例如上面的mdwt_test12.m，只是将刚才使用，只是将刚才使用CL2小波（

128、滤波器的长度为小波（滤波器的长度为3）的）的mdwt_test10.m，改成使用，改成使用ghm小波（滤波器的长度为小波（滤波器的长度为4）。那么，我们的算法的结果与）。那么，我们的算法的结果与MWMP不完全相同，能够使用吗？实际上不影响使用，因为不完全相同，能够使用吗？实际上不影响使用，因为它自成体系，只是与它自成体系，只是与MWMP在计算方式上有些不同，只要能在计算方式上有些不同，只要能精确重构，就能够在实际计算中使用。看看下面的例子。精确重构，就能够在实际计算中使用。看看下面的例子。 下面下面程序程序mdwt_test13.m，是由，是由mdwt_test12.m修改的，此程序使用修改的

129、，此程序使用ghm小波和小波和ghmorap方法，先做方法，先做2次正变换，然后做次正变换，然后做2次反变换，次反变换，最后计算反变换后的数据与原始数据差的绝对值的和。最后计算反变换后的数据与原始数据差的绝对值的和。%程序名：程序名：mdwt_test13.mn=32;rand(seed,10000);x=rand(n,n); p,q=pmatrix(n,ghmorap);L1,H1=mdwt_matrix(n,ghm);L2,H2=mdwt_matrix(n/2,ghm);A=L2*L1;H2*L1;H1;T=smatrix(n,2);B=T*A*p; x1=B*x*B;x2=B*x1*B;

130、 s=sum(sum(abs(x2-x);fprintf(MWMP方法与我们方法的差的绝对值之和为方法与我们方法的差的绝对值之和为:%2.4fn,s); mdwt_test13MWMP方法与我们方法的差的绝对值之和为方法与我们方法的差的绝对值之和为:0.0000多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 刚才的幻灯片，使用了小波的逆变换的具体形式，关于小波刚才的幻灯片，使用了小波的逆变换的具体形式，关于小波逆变换，我们有如下结论：逆变换，我们有如下结论：（1）如果小波变换使用的是正交预滤波方法，例如）如果小波变换使用的是正交预滤波方法，例如g

131、hm小波小波使用的是使用的是ghmorap, sa4小波的小波的sa4ap等等，若设等等，若设x是要变换的是要变换的2维数据，则正变换是维数据，则正变换是y=B*x*B，其中，其中B=T*A*p, T是置换矩是置换矩阵，阵，A是小波变换矩阵，是小波变换矩阵， p是预滤波矩阵。反变换则可简单表是预滤波矩阵。反变换则可简单表示为示为x=B*y*B, 此时后滤波矩阵此时后滤波矩阵q=p。（2）如果小波变换使用的是双正交预滤波方法，例如）如果小波变换使用的是双正交预滤波方法，例如ghm小小波使用的是波使用的是ghmap, cl2小波的小波的clap等等，若设等等，若设x是要变换的是要变换的2维维数据，

132、则正变换是数据，则正变换是y=B*x*B，其中，其中B=T*A*p, T是置换矩阵，是置换矩阵，A是小波变换矩阵，是小波变换矩阵， p是预滤波矩阵。此时，若是预滤波矩阵。此时，若q是后滤波矩是后滤波矩阵，则反变换可表示为阵，则反变换可表示为x=BB*y*BB, 其中其中BB=T*A*q，请注，请注意，意，T*A*q中的中的q是是q的置换矩阵，不是整个矩阵相乘后的的置换矩阵，不是整个矩阵相乘后的转置，即转置，即T*A*q = (T*A)*(q) 。 下面我们先看看下面我们先看看lena图像的小波变换的例子，然后再看看图像的小波变换的例子，然后再看看小波反变换的例子。小波反变换的例子。多多 小小

133、波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 程序程序:mdwt_test14.m使用使用ghm小波，对小波，对lena.bmp(512X512)图图像进行了像进行了2次小波变换，变换结果如下：次小波变换，变换结果如下：请注意，请注意，lena.bmp图像一定图像一定要在当前要在当前目录下。目录下。图中的红图中的红色线是我色线是我后加上去后加上去的，是为的，是为了更方便了更方便地观察经地观察经过过2次变次变换后小波换后小波的各个子的各个子带情况。带情况。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 程序程序

134、mdwt_test15.m与与mdwt_test14.m的不同之处在于，的不同之处在于， mdwt_test15.m对对512X512的的lena.bmp图像，进行了图像，进行了3次多小波变次多小波变换，请注意进行换，请注意进行3次变换的变换矩阵、置换矩阵的表示形式。次变换的变换矩阵、置换矩阵的表示形式。%程序名：程序名：mdwt_test15.mx=imread(lena.bmp);x=double(x);m,n=size(x);p,q=pmatrix(n,ghmap);L1,H1=mdwt_matrix(n,ghm);L2,H2=mdwt_matrix(n/2,ghm);L3,H3=mdw

135、t_matrix(n/4,ghm);A=L3*L2*L1;H3*L2*L1;H2*L1;H1;T=smatrix(n,3);x3=T*A*p*x*p*A*T;colormap(gray(256);imagesc(x3);axis(off);下面是此程序的运行结果，以及下面是此程序的运行结果，以及3次变换后子带划分情况。次变换后子带划分情况。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 如果对二维信号做如果对二维信号做4次小波变换，程序如何写？下面仍然

136、对次小波变换，程序如何写？下面仍然对lena.bmp图像，做图像，做4次小波变换，程序如下：次小波变换，程序如下：%程序名：程序名：mdwt_test16.mx=imread(lena.bmp);x=double(x);m,n=size(x);p,q=pmatrix(n,ghmap);L1,H1=mdwt_matrix(n,ghm);L2,H2=mdwt_matrix(n/2,ghm);L3,H3=mdwt_matrix(n/4,ghm);L4,H4=mdwt_matrix(n/8,ghm);A=L4*L3*L2*L1;H4*L3*L2*L1;H3*L2*L1;H2*L1;H1;T=smatr

137、ix(n,4);x3=T*A*p*x*p*A*T;colormap(gray(256);imagesc(x3);axis(off);程序的运行结果为：程序的运行结果为：多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 下面考虑小波变换的逆变换，先看看双正交插值预滤波情况。下面考虑小波变换的逆变换，先看看双正交插值预滤波情况。首先是对首先是对256X256个数据，使用个数据，使用ghm小波做小波做2次变换，同时使次变换，同时使用双正交插值预滤波方法用双正交

138、插值预滤波方法ghmap, 程序分别使用程序分别使用B=T*A*p的的BB=T*A*(q)进行了逆变换，使用进行了逆变换，使用BB做逆变换是正确的。程做逆变换是正确的。程序的文件名是：序的文件名是：mdwt_test17.m, 程序及结果如下：程序及结果如下：n=256;rand(seed,10000);x=rand(n,n);pflt=ghmap;flt=ghm; p,q=pmatrix(n,pflt);L1,H1=mdwt_matrix(n,flt);L2,H2=mdwt_matrix(n/2,flt);A=L2*L1;H2*L1;H1;T=smatrix(n,2);B=T*A*p;x1=

139、B*x*B;x2=B*x1*B;BB=T*A*(q); x3=BB*x1*BB;s1=sum(sum(abs(x2-x); s2=sum(sum(abs(x3-x);fprintf(s1=%4.4f, s2=%4.4fn,s1, s2); mdwt_test17s1=40894.2233, s2=0.0000多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 n=256;rand(seed,10000);x=rand(n,n);pflt=ghmorap;flt=ghm; p,q=pmatrix(n,pflt);L1,H1=mdwt_matrix(n,

140、flt);L2,H2=mdwt_matrix(n/2,flt);A=L2*L1;H2*L1;H1;T=smatrix(n,2);B=T*A*p;x1=B*x*B;x2=B*x1*B;BB=T*A*(q); x3=BB*x1*BB;s1=sum(sum(abs(x2-x);s2=sum(sum(abs(x3-x);fprintf(s1=%4.4f, s2=%4.4fn,s1,s2);如果将如果将mdwt_test17.m中的语句中的语句pflt=ghmap，改为使用正交，改为使用正交预滤波方法预滤波方法pflt=ghmorap，其它的语句都不变，就变成下，其它的语句都不变，就变成下面的程序面的程

141、序mdwt_test18.m，则对进行小波反变换时，两种方法，则对进行小波反变换时，两种方法即小波系数矩阵乘以矩阵即小波系数矩阵乘以矩阵B或者矩阵或者矩阵BB都是正确的，程序及都是正确的，程序及结果如下：结果如下： mdwt_test18s1=0.0000, s2=0.0000多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 下面考虑使用下面考虑使用CL2小波做变换的情况，还是对小波做变换的情况，还是对256X256个数个数据，使用据，使用CL2小波做小波做3次变换，同时使用次变换，同时使用clap方法进行预滤波方法进行预滤波, 程序分别使用程序分

142、别使用B=T*A*p的的BB=T*A*(q)进行了逆变换，使用进行了逆变换，使用BB做逆变换是正确的。程序的文件名是：做逆变换是正确的。程序的文件名是：mdwt_test19.m, 程程序及结果如下：序及结果如下：n=256;rand(seed,10000);x=rand(n,n);pflt=clap;flt=cl; p,q=pmatrix(n,pflt);L1,H1=mdwt_matrix(n,flt);L2,H2=mdwt_matrix(n/2,flt); L3,H3=mdwt_matrix(n/4,flt);A=L3*L2*L1;H3*L2*L1;H2*L1;H1;T=smatrix(n

143、,3);B=T*A*p;x1=B*x*B;x2=B*x1*B;BB=T*A*(q); x3=BB*x1*BB;s1=sum(sum(abs(x2-x);s2=sum(sum(abs(x3-x);fprintf(s1=%4.4f, s2=%4.4fn,s1,s2); mdwt_test19s1=32181.0682, s2=0.0000多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 对文件对文件mdwt_test19.m中的命令中的命令 pflt=clap;flt=cl;修改为修改为 pflt=sa4ap;flt=sa4;即使用即使用sa4小波做同

144、样的变换，程序是小波做同样的变换，程序是mdwt_test20.m，程序，程序代码及结果如下：代码及结果如下： n=256;rand(seed,10000);x=rand(n,n);pflt=sa4ap;flt=sa4; p,q=pmatrix(n,pflt);L1,H1=mdwt_matrix(n,flt);L2,H2=mdwt_matrix(n/2,flt); L3,H3=mdwt_matrix(n/4,flt);A=L3*L2*L1;H3*L2*L1;H2*L1;H1;T=smatrix(n,3);B=T*A*p;x1=B*x*B;x2=B*x1*B;BB=T*A*(q); x3=BB*

145、x1*BB;s1=sum(sum(abs(x2-x);s2=sum(sum(abs(x3-x);fprintf(s1=%4.4f, s2=%4.4fn,s1,s2); mdwt_test20s1=0.0000, s2=0.0000两种计算方法都是正确的。两种计算方法都是正确的。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 自编子程序自编子程序1：2维多小波变换的低频子带的提取维多小波变换的低频子带的提取多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 下面的程序下面的程序mdwt_test24.m对

146、对lena.bmp(512X512)做做3次变换，次变换，提取其低频部分，程序如下：提取其低频部分，程序如下：pflt=sa4ap;flt=sa4;x=imread(lena.bmp);x=double(x);m,n=size(x);p,q=pmatrix(n,pflt);L1,H1=mdwt_matrix(n,flt);L2,H2=mdwt_matrix(n/2,flt);L3,H3=mdwt_matrix(n/4,flt);A=L3*L2*L1;H3*L2*L1;H2*L1;H1;T=smatrix(n,3);x3=T*A*p*x*p*A*T;w=get_low_subband_coef(

147、x3,3);mm=m/23;subplot(2,2,1);colormap(gray(256);imagesc(w(1:mm/2,1:mm/2);axis(off);subplot(2,2,2);colormap(gray(256);imagesc(w(1:mm/2,mm/2+1:mm);axis(off);subplot(2,2,3);colormap(gray(256);imagesc(w(mm/2+1:mm,1:mm/2);axis(off);subplot(2,2,4);colormap(gray(256);imagesc(w(mm/2+1:mm,mm/2+1:mm);axis(off

148、); mdwt_test24 右图是右图是MATLAB的运行的运行结果，它对应的子带的序结果，它对应的子带的序号分别是号分别是 11 12 21 22多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 自编子程序自编子程序2：2维多小波变换高频子带的提取维多小波变换高频子带的提取多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 %程序名：程序名：mdwt_test25.mpflt=sa4ap;flt=sa4;x=imread(lena.bmp);x=double(x);m,n=size(x);p,q=pm

149、atrix(n,pflt);L1,H1=mdwt_matrix(n,flt);L2,H2=mdwt_matrix(n/2,flt);L3,H3=mdwt_matrix(n/4,flt);A=L3*L2*L1;H3*L2*L1;H2*L1;H1;T=smatrix(n,3);x3=T*A*p*x*p*A*T; k=1;for level=1:3 for i=1:4 for j=1:4 if (i=2)&(j=2) else w=get_subband_coef(x3,level,i,j); subplot(6,6,k);colormap(gray(256);imagesc(w);axis(off

150、); ll=int2str(level);ii=int2str(i);jj=int2str(j);s=strcat(ll,ii,jj,);title(s);k=k+1; endend; end; end;下面的程序，对灰度图像下面的程序，对灰度图像lena.bmp，使用，使用sa4多小波做了多小波做了3次次分解，由于每次分解后，有分解，由于每次分解后，有16个子带，而个子带，而16个子带中的个子带中的4个低个低频子带会进一步分解为频子带会进一步分解为16个子带，因此，每次分解后，会出个子带，因此，每次分解后，会出现现12个高频子带，经过个高频子带，经过3次分解后，总计有次分解后，总计有36个高

151、频子带。程个高频子带。程序序mdwt_test25.m在在MATLAB的同一个图形窗口中，显示了的同一个图形窗口中，显示了这这36个高频子带的具体图形，其中，每个图形的标题个高频子带的具体图形，其中，每个图形的标题ki,j的的意义为：意义为：k表示变换的次数，表示变换的次数，i,j表示第表示第k次变换下，按照函数次变换下，按照函数get_subband_coef将高频子带分块标记，第将高频子带分块标记，第i行第行第j列的那个高列的那个高频子带，程序具体如下，执行结果见下一个幻灯片。频子带，程序具体如下，执行结果见下一个幻灯片。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其

152、其 软软 件件 实实 现现 Lena图图像经过像经过3次小波变次小波变换后的换后的36个高频子个高频子带，其中，带，其中，每个图形每个图形的标题的标题ki,j的意的意义为：义为：k表示变换表示变换的次数，的次数，i,j表示表示第第k次变次变换下第换下第i行第行第j列列的高频子的高频子带带多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 自编子程序自编子程序3：2维多小波变换低频子带的修改维多小波变换低频子带的修改多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 自编子程序自编子程序4：2维多小波变换高频

153、子带的修改维多小波变换高频子带的修改多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 下面是图像增强与边缘检测的一个例子下面是图像增强与边缘检测的一个例子,程序程序mdwt_test26.m, 其想法是将图其想法是将图像变换到小波域像变换到小波域, 让低频小波系数变得比原来大一些让低频小波系数变得比原来大一些,让高频系数减小一些让高频系数减小一些, 这样就增强了图像的前景这样就增强了图像的前景, 弱化的图像的背景弱化的图像的背景, 因而增强了图像因而增强了图像,同时对增强同时对增强前景后图像前景后图像,设置一个阈值设置一个阈值, 让小于此值的像素点

154、为让小于此值的像素点为0, 大于此值的像素点为大于此值的像素点为1, 检测图像的边缘检测图像的边缘. 注意注意, 对于大幅度改变低频小波系数值的小波变换对于大幅度改变低频小波系数值的小波变换,一定要一定要选择正交小波及相应的正交预滤波方法选择正交小波及相应的正交预滤波方法, 否则变换后的图像会失真否则变换后的图像会失真. pflt=sa4ap;flt=sa4;x=imread(cameraman.bmp);x=double(x);m,n=size(x);p,q=pmatrix(n,pflt);L1,H1=mdwt_matrix(n,flt);L2,H2=mdwt_matrix(n/2,flt)

155、;L3,H3=mdwt_matrix(n/4,flt);A=L3*L2*L1;H3*L2*L1;H2*L1;H1;T=smatrix(n,3);x3=T*A*p*x*p*A*T;mm=m/23;w=get_low_subband_coef(x3,3);w11=w(1:mm/2,1:mm/2);w12=w(1:mm/2,mm/2+1:mm);w21=w(mm/2+1:mm,1:mm/2);w22=w(mm/2+1:mm,mm/2+1:mm);r11=1.5;r12=1.2;r21=1.2;r22=1;r=0.7,0.9,1.0;ww=r11*w11,r12*w12; r21*w21,r22*w2

156、2; x3=set_low_subband_coef(x3,3,ww);for level=1:3 for i=1:4 for j=1:4 if (i=2)&(j=2) else w=get_subband_coef(x3,level,i,j);x3=set_subband_coef(x3,level,i,j,r(level)*w); endend; end; end;B=T*A*(q); x4=B*x3*B; subplot(2,2,1);colormap(gray(256);imagesc(x);axis(off);title(原始图像效果原始图像效果);subplot(2,2,2);co

157、lormap(gray(256);imagesc(x4);axis(off);title(图像增强后效果图像增强后效果);x_thresh=130;for i=1:m for j=1:n if x(i,j)x_thresh x(i,j)=0; else x(i,j)=1; end; if x4(i,j)x_thresh x4(i,j)=0; else x4(i,j)=1; end; end; end;subplot(2,2,3);colormap(gray(256);imagesc(x);axis(off);title(原图像的二值化原图像的二值化);subplot(2,2,4);colorm

158、ap(gray(256);imagesc(x4);axis(off);title(增强后的二值化增强后的二值化);例例1：使用：使用小波变换方小波变换方法进行灰度法进行灰度图像的增强图像的增强与边缘检测与边缘检测方法的试验方法的试验多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 左上图是左上图是原始图像。原始图像。 右上图是右上图是增强后的增强后的图像。对图像。对于相同的于相同的阈值，左阈值，左下图是原下图是原始图像二始图像二值化的结值化的结果果, 右下右下图是增强图是增强后图像二后图像二值化的结值化的结果。比较果。比较可知，增可知，增强后图像强

159、后图像有效地消有效地消除了图像除了图像背景中的背景中的无用信息无用信息.多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 pflt=sa4ap;flt=sa4;x=imread(lena.bmp);x=double(x);m,n=size(x);p,q=pmatrix(n,pflt);L1,H1=mdwt_matrix(n,flt);L2,H2=mdwt_matrix(n/2,flt);L3,H3=mdwt_matrix(n/4,flt);A=L3*L2*L1;H3*L2*L1;H2*L1;H1;T=smatrix(n,3);x3=T*A*p*x*

160、p*A*T;mm=m/23;w=get_low_subband_coef(x3,3);w11=w(1:mm/2,1:mm/2);w12=w(1:mm/2,mm/2+1:mm);w21=w(mm/2+1:mm,1:mm/2);w22=w(mm/2+1:mm,mm/2+1:mm);r11=1.1; r12=1.1; r21=1.1; r22=0.9; r=0.1, 0.5, 0.7;ww=r11*w11,r12*w12; r21*w21,r22*w22; x3=set_low_subband_coef(x3,3,ww);for level=1:3 for i=1:4 for j=1:4 if (i

161、=2)&(j=2) else w=get_subband_coef(x3,level,i,j);x3=set_subband_coef(x3,level,i,j,r(level)*w); endend; end; end;B=T*A*(q); x4=B*x3*B; subplot(2,2,1);colormap(gray(256);imagesc(x);axis(off);title(原始图像效果原始图像效果);subplot(2,2,2);colormap(gray(256);imagesc(x4);axis(off);title(图像增强后效果图像增强后效果);x_thresh=135;f

162、or i=1:m for j=1:n if x(i,j)x_thresh x(i,j)=0; else x(i,j)=1; end; if x4(i,j)x_thresh x4(i,j)=0; else x4(i,j)=1; end; end; end;subplot(2,2,3);colormap(gray(256);imagesc(x);axis(off);title(原图像的二值化原图像的二值化);subplot(2,2,4);colormap(gray(256);imagesc(x4);axis(off);title(增强后的二值化增强后的二值化);下面的程序下面的程序mdwt_tes

163、t27.m与与mdwt_test26.m的区别是，只修改的区别是，只修改了程序中用红颜色标注的语句，其它与了程序中用红颜色标注的语句，其它与mdwt_test26.m相同。相同。多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 我们在程序我们在程序mdwt_test26.m中解释过，对于大幅度改变低频小波系数值的小中解释过，对于大幅度改变低频小波系数值的小波变换波变换, 要选择正交小波要选择正交小波, 否则变换后的图像会失真。将否则变换后的图像会失真。将

164、mdwt_test26.m的小的小波更改为波更改为ghm小波，其它的都与小波，其它的都与mdwt_test26.m相同，就是程序相同，就是程序mdwt_test28.m，程序见当前目录的具体文件，程序的执行结果如下：，程序见当前目录的具体文件，程序的执行结果如下：多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 %程序名：程序名：mdwt_test29.mpflt=sa4ap;flt=sa4;x=imread(lena.bmp);x=double(x);m,n=size(x);p,q=pmatrix(n,pflt);L1,H1=mdwt_matri

165、x(n,flt);L2,H2=mdwt_matrix(n/2,flt);L3,H3=mdwt_matrix(n/4,flt);A=L3*L2*L1;H3*L2*L1;H2*L1;H1;T=smatrix(n,3);randn(seed,10000);sigma=100;x1=x+sigma*randn(m,n);x3=T*A*p*x1*p*A*T; mm=m/23; method=sqtwolog; r=0.7;for level=1:3 for i=1:4 for j=1:4 if (i=2)&(j mdwt_test29降噪前信噪比降噪前信噪比SNR=2.4754, 降噪后信噪比降噪后信噪

166、比SNR=16.8510.多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 本讲座结束本讲座结束多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现 多多 小小 波波 变变 换换 的的 矩矩 阵阵 形形 式式 及及 其其 软软 件件 实实 现现