多小波变换的矩阵形式.ppt

多小波变换的矩阵形式及其软件实现,我们知道，进行1次多小波变换的分解与重构公式为：,与单小波不同之处在于，公式中的s(n, k)是r维列向量，H(k), G(k)是rXr大小的矩阵。,因此，在使用这个公式前，必须进行预处理，使进行小波变换的一维数据变成rXN维的数据，这称为预滤波。,预滤波后，进行小波分解，然后进行重构，重构后再将rXN维数据变成1维数据，这称为后滤波。,预滤波与后滤波的矩阵表示,从前面的解释中可以看出，假设多小波的重数是r维的，多小波分解算法要求有初始系数（即r维向量）才能进行计算，而通常的输入数据是一维的信号，或者是函数f(t)的等间距采样值，为了应用上面的公式，这就需要将一维数据或者f(t)的等间距采样值转化为r维向量，这个转化的过程一般称为预处理或者预滤波。对一维数据进行预滤波后，进行小波变换即小波分解，然后进行逆变换即小波重构，还得将重构后的r维向量转化为一维数据，这个过程就称为后滤波或者后处理。一般来说，小波重构是精确重构，则需要后滤波是预滤波的逆过程。一维多小波变换过程如下： 一维信号预滤波分解重构后滤波一维信号 研究多小波变换的矩阵表示，首先应该知道预滤波与后滤波的矩阵表示形式。,然后做变换，最后结果中，取每个向量的第1个元素就是不重复预滤波方法的后滤波。这种方法我们不讨论，如果你想讨论的话，我们本次讲座的内容也完全能够处理这种方法。 2、恒等预滤波方法 就是按信号的前后顺序，直接将信号变成向量，不做其它处理，直接进行多小波变换的方法。例如对s=a,b,c,d，有,预滤波方法常见有： 1、重复预滤波方法 比如对重数r=2的GHM小波，如果s=a,b,c,d; 则将s变成向量,恒等预滤波方法适用于平衡多小波。,3、插值预滤波方法 这种方法见此软件目录下的原创文献：Design of prefilters for discrete multiwavelet transforms.pdf，具体想法是针对特定的多小波，将给定的要做小波变换的一维信号f(t)看成是等间隔采样的ft_kj, 然后通过插值方法得到向量c_0k, MWMP软件包中文件coef_prep.m中，有关于GHM的2个预滤波方法，例如对于GHM小波，其双正交插值预滤波方法为：,此软件目录下的子目录“预滤波”和“新加坡国立大学Qingtang Jiang的多小波文章(00年前后)”，都是一些关于预滤波的论文，请参考。,4、Hardin-Roach预滤波方法 插值预滤波方法不能同时保持逼近阶和正交性，而这种方法能同时保持逼近阶和正交性，具体参考此目录下2 篇原创文献“Multiwavelet Prefilters I .pdf”、“Multiwavelet Prefilters II.pdf” 。,预滤波的一般表示形式：,由于在实际应用中，我们使用的多小波都是2重的，即r=2，此时，预滤波器可以写成如下的通用形式：,其中P_0, P_1, , P_L就是长度为L+1的预滤波器。 下面我们给出一些小波的预滤波与后处理的具体滤波器，并且与MWMP软件包中文件coef_prep.m的相同小波进行比较说明，通过这种比较，你就知道，如果你在其它文献或者教材上看到一个新的预滤波器，如何加入文件coef_prep.m 中了。,下面是一些常用的双正交插值预滤波、后滤波的滤波器， 其中P_k是预滤波器，Q_k是后滤波器：,预滤波与后滤波的一般关系,STT多小波，也称SA4多小波，其中P_0是预滤波器，Q_0是后滤波器：,预滤波与后滤波的关系,关于预滤波方法的程序验证：,对GHM小波，其双正交插值预滤波方法如下,请注意：后滤波的公式虽然写成这样形式，但它是预滤波的逆过程，因此，请注意观察下面程序中关于后滤波的实际计算方法，或者看看论文Design of prefilters for discrete multiwavelet transforms,%本程序演示对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); 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的预滤波的结果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=);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 20.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.0000 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); disp(向量化后的数据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的预滤波的结果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=);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 -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.0000 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(原始数据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); end disp(直接使用矩阵乘积计算，得到的xx的预滤波的结果x1=);disp(x1); x2=; for i=1:length(x1)-1 x2=x2,PO*x1(:,i); end disp(对前面预滤波的数据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 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= 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.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,通过前面的3个程序，可以看出，预滤波和后滤波的计算可以用矩阵表示，我们看一个例子。假设我们使用GHM多小波的双正交GHMAP预滤波与后滤波方法，对12个数据的一维信号s进行计算，预滤波与后滤波变换矩阵可以表示为:,%本程序演示对GHM多小波，使用ghmap方法进行预滤波与后滤波的矩阵计算方法 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_mat pr_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的预滤波的结果s1=);disp(s1); fp=prep1D_appe(s,ghmap); disp(MWMP中的相同预滤波方法计算得fp,结果如下，可以看出，将s1向量化后就是fp);disp(fp); %构造后滤波变换矩阵po_mat x=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(直接使用矩阵乘积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.3623 22.0000 25.7210 55.0000 35.2670 8.0000 MWMP中的相同预滤波方法计算得fp,结果如下，可以看出，将s1向量化后就是fp 12.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 ,注意，如果将上面程序中的ghm小波换成sa4小波，使用sa4ap方法进行计算，那么预滤波矩阵P和后滤波矩阵Q是什么样子的形式？显然，它们都是对角矩阵。,%本程序演示对SA4多小波，使用sa4ap方法进行预滤波与后滤波的矩阵计算方法 s=8,12,0,5,20,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_mat pr_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,sa4ap); disp(MWMP中的相同预滤波方法计算得fp,结果如下，可以看出，将s1向量化后就是fp);disp(fp); %构造后滤波变换矩阵po_mat po_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);,从前一个幻灯片知道，如果预滤波和后滤波器的长度为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 直接使用矩阵乘积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.0416 MWMP中的相同预滤波方法计算得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 -24.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,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=p*s，就是对信号s进行了预滤波，结果是s1,计算s2=q*s1,则s2就是对经过预滤波的信号s1进行了后滤波，请注意，此程序使用了coef_prep.m中的预滤波与后滤波方法，你要安装MWMP软件包或者将coef_prep.m拷贝到与你运行的程序在同一目录下，或者是MATLAB的WORK子目录中。, p,q=pmatrix(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 16.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.8245 5.9335 5.7146 17.9971 8.9024 -1.2411 19.2250 6.7448 20.2496 6.0699,下面使用10个数据，验证程序的运行情况：,到此为止，我们已经知道，可以使用自编子程序pmatrix，获得多小波变换的预滤波变换矩阵和后滤波变换矩阵。我们知道，多小波变换的过程一般为 一维信号预滤波分解重构后滤波一维信号 如果也能够将小波变换的主要步骤，即 分解重构 写成矩阵变换形式的表达式，那么，结合预滤波变换矩阵和后滤波变换矩阵的表达式，我们就可以将多小波变换写成一个统一的矩阵变换形式。 基于此，下面我们考虑，进行1次小波分解与重构的矩阵表示形式。下面的程序，需要使用MWMP中的文件coef.m，请安装MWMP软件包，或者将文件coef.m拷贝到你目前的程序目录中，或者是MATLAB的WORK子目录中。,先看看coef.m中多小波滤波器H_k,G_k的输入格式：,通过比较，如果你在论文中看到一个新小波，你就知道怎么加入到coef.m文件中去了。,实际计算中，一般多小波的重数r=2，s=0, 此时加细方程为,上面的加细方程是小波与多小波一书的写法，此时应为 则小波变换的低频变换矩阵L与高频变换矩阵H为(对重数w=2)：,function L_matrix,H_matrix=mdwt_matrix(len,flt) %在文件coef.m中，低频滤波器系与高频滤波器系要相同， %如果不相同，要使用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); end return,下面的函数调用L, H=mdwt_matrix(len,flt)，对于coef.m中的小波flt,预滤波后的数据长度为len的数据，就能够返回1次小波变换的低频变换矩阵L与高频变换矩阵H。如果经过预滤波后的数据即向量为s1，则矩阵乘法运算 s2=L*s1, 或者 s2=H*s1 就是经过1次小波变换后的低频或者高频系数，自编函数如下：, 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.0350 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 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 = 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(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个，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); x ans = 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.8318 0.5028 0.7095 0.4289 p,q=pmatrix(n,sa4ap); L1,H1=mdwt_matrix(n,sa4);S1=(L1*p*x); S1 ans = 1.1462 0.1607 0.9511 0.0112 1.1626 0.1618 1.1897 0.1574 D1=(H1*q*x); D1 ans = 0.2400 -0.6664 0.0009 0.0774 -0.3677 0.4197 -0.0936 -0.1867 ,如果变换数据x的长度为n,对经过 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; xx ans = 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 11.0000 12.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=pmatrix(n,sa4ap);L1,H1=mdwt_matrix(n,sa4); x1=L1;H1*p*x; x1 ans = 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.8730 -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次，结果是相同的。,下面考虑进行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), 返回矩阵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;x1 ans = 30.8882 6.8465 17.1951 -6.8465 1.5881 -1.0607 9.3663 6.7175 0 -0.0000 -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); x4 ans = 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; y ans = 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.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次小波变换正变换后，又进行2次反变换的演示程序，同时也将变换结果与MWMP进行了比较，结果相同，由此看出，我们编写的程序完全可以代替MWMP软件包。,上面是进行2次小波变换的整体变换，实际上，你也可以考虑不进行整体变换，而是直接计算变换的各个子带，例如第1次变换的低频部分是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 -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) %第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_matrix(n/4,sa4); B=L3*L2*L1;H3*L2*L1;H2*L1;H1; A=B*p; x1=A*x; x1 ans = 92.7399 19.8490 39.2601 -0.4841 8.8490 -7.3750 10.5159 7.3750 -0.0000 0.1796 -0.0000 0.1796 3.1763 -2.3009 18.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 x2=prep1D_appe(x,sa4ap);x3=dec1D_pe(x2,sa4,3); x4=reshape(x3,n,1); x4 ans = 92.7399 19.8490 39.2601 -0.4841 8.8490 -7.3750 10.5159 7.3750 0.0000 0.1796 0.0000 0.1796 3.1763 -2.3009 18.7326 13.2554 -0.0000 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 -8.0000 x5=q*B*x1; x5 ans =1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000 16.0000 17.0000 18.0000 19.0000 20.0000 21.0000 22.0000 23.0000 24.0000 25.0000 26.0000 27.0000 28.0000 29.0000 30.0000 31.0000 32.0000,下面的程序演示了进行3次小波分解与重构，并将结果与MWMP软件包进行比较，其结果是相同的。,关于1维多小波变换的子带分割问题,我们知道，1维多小波变换每次只对低频子带分解，每次将低频子带进一步分成1个低频子带和1个高频子带，这个结论不管对单小波，还是多小波，都是适用的。, 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次变换，数据被变换为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,低频部分,次高频部分,高频部分, 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; x5 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 -0.1524 -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个部分，每个部分的第1行的频率比第2行的频率低。这样，对于1维多小波变换，1次变换分成了4个子带，2次变换分成6个子带，3次变换分成8个子带。,假定有n个数据，这n个数据经过2次多小波变换，即进行了 一维信号预滤波2次小波分解 的过程后，例如信号n=16,MATLAB生成程序为： n=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