Daubechies小波分析的部分Matlab实现

注：本实验在原本操作过程中在重构部分存在着问题，经思考检验调试之后，排除了隐藏的问题，实验得到成功进行同时达到了预期的结果。

第一部分

Daubechies 尺度序列，二进点上的尺度函数图像及小波函数图像

第一部分实验的目标有两个：

1，给出各阶Daubechies小波的低通滤波器，即尺度函数两尺度序列。具体的程序算法由《小波导论》中7.3节中的引理7.16，定理7.17给现，其中定理7.17的证明过程实际上给出了一个显式的计算过程，这是程序实现的重要依据。

2，根据两尺度关系，由第一步骤中得到的两尺度序列逐步计算得到Daubechies 尺度函数及小波函数的二进点值，进而给出近似图像。

第一步骤

相关函数程序解释：

1、sinj(n)，实现将引理7.16中的(sin)2n转化为形如cosk，n即表示(sin)的次数。 22k2、sinj(n)中涉及函数cosn(n), cosn(n)实现将cos转化为

cos(kw)，n表示三角函数次数。

kk程序编写的思路就是先将(sin2)2n转化为cosk再转化为cos(kw)，然后根据定理7.17的证明过程计算出两尺度序列，其中用于提取多

项式系数的函数利用了三角函数的正交性，使用了积分函数int()。daucoef(n)给出n阶daubechies尺度函数的两尺度序列。具体的原程序代码见附

件，部分程序运行结果如下：

各阶（2--12）Daubechies 尺度函数的两尺度序列如下：（规范化因子为2）

表1 2~12阶Daubechies 尺度函数的两尺度序列 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 2 3 4 0.325803 1.010946 0.892200 5 0.226419 0.853944 1.024327 6 0.157742 0.699504 1.062264 0.445831 7 0.110099 0.560791 1.031148 0.664372 8 0.076956 0.442467 0.955486 0.827817 9 0.053850 0.344834 0.855349 0.929546 10 0.037717 0.266122 0.745575 0.973628 0.397638 11 0.026438 0.203742 0.636254 0.969708 0.582606 12 0.018544 0.154950 0.533661 0.929419 0.729574 0.683013 0.470467 1.183013 1.141117 0.316987 0.650365 -0.183013 -0.190934 -0.039575 0.195767 0.000000 0.049817 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.043616 0.046504 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -0.120832 -0.264507 -0.342657 -0.319987 -0.203514 -0.022386 0.188370 0.109703 0.137888 0.100846 0.114003 0.000668 0.182076 -0.045601 -0.183518 -0.316835 -0.401659 -0.414752 -0.353336 -0.229492 -0.063306 -0.136954 -0.277110 -0.387821 -0.447144 0.210068 0.180127 0.131603 0.093400 0.211866 -0.033629 0.258064 -0.014987 -0.008827 0.038923 0.004717 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000783 0.006756 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -0.017792 -0.044664 -0.053782 -0.024564 0.043453 0.017750 0.019772 0.012369 0.000355 0.031624 -0.023440 -0.062350 -0.095647 -0.100967 -0.065733 0.007580 -0.041659 -0.093959 -0.136376 0.046970 0.044315 0.029473 0.015343 0.058755 -0.001524 0.000608 0.000500 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -0.002548 -0.006888 -0.006680 0.005100 0.000955 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.002613 0.001973 0.002818 -0.000554 -0.006055 -0.015179 -0.021729 -0.017280 -0.004725 -0.018160 0.006970 0.009491 -0.000166 0.000326 0.000056 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -0.000356 -0.000970 -0.000436 0.003180 -0.000165 -0.001263 -0.003082 0.000132 0.000000 0.000000 0.000000 0.000000 0.000352 0.000009 0.000550 -0.000034 0.000018 -0.000002 -0.000019 0.000077 0.000006 0.000000 0.000000 -0.000049 -0.000125

第二步骤：

给出n阶daubechies尺度函数图像，由函数dauinterp(n,max1)实现，即给出n阶daubechies尺度函数的max1次迭代的函数值，所得到的点都是函数的二进点。Dauinterpwa(n,max1)则给出小波函数的二进点上的函数值。 这一部分内容中，比较难以实现的地方在于下标的变换，这常常引起一些数据的混乱或错误，相关的数学推导总是难以避免的，由于过程过于繁复，这里不给出详细的推导过程及解释。 相关程序运行结果如下，分别给出了大图与小图，方便查阅比较。 以下罗列的二进点迭代图形都是经过8次迭代得到的图形，即每个单位区间上有256个点。如果增加迭代次数的话当然可以使图形更加光滑，但这同时也会增加程序运行的时间。

各阶（2--10）Daubechies 尺度函数图形如下：

图1 2~10阶Daubechies 尺度函数图像

分别用二进点迭代8次得到的消失矩为3、5、7的Daubechies尺度函数大图：

Daubechies 3 尺度函数图像

Daubechies 5 尺度函数图像

Daubechies 7 尺度函数图像

从图像直接观察可知，Daubechies 尺度函数随着其阶数的增大而变得越来越光滑，而且在其支撑区间上，总会在某一个区间段上，其值会变得很小，甚至于充分接近0。因此，实际中起主要作用的也就是其函数值较大的区间段而已。 相应地，我们画出对应的Daubechies小波函数的图像，综合图如下：

图2 2~10阶Daubechies 小波函数图像

下面单独给出消失矩为3、5、7、9的Daubechies小波函数图像 消失矩为3的Daubechies小波函数图像

消失矩为5的Daubechies小波函数图像

消失矩为7的Daubechies小波函数图像

消失矩为9的Daubechies小波函数图像

以上图像都是根据两尺度方程进行迭代，得到小波函数及尺度函数在二进点上的值，这里的迭代只进行了8次，由于迭代的次数增多，相应地就要增加指数倍的运算量，然而在图形上的改变却是微小的，所以，在只想要看到图形的情况下，也就没有太多次的迭代了。

而从图形来看，无论是小波函数，还是尺度函数，都有一种收敛的趋势，也就是说，它们在图形上有着一定的一致性，并且向着某种形状慢慢接近，而这种“终极图形”也许就像是一个振荡微小的波浪。

第二部分

用Cascade迭代方法得到的Daubechies尺度函数图像及小波函数图像

一般地，我们会使用Cascade 算法迭代过程来产生尺度函数及小波函数的一个逼近函数，而初始选择的函数一般不会作特别的要求，所以我们总可以选择合适的函数来作为初始函数。实验结果表明，若用Haar函数作为初始函数，由于Haar函数不是连续的，因此其收敛速度很慢。因此，我们使用1阶B样条函数作为初始函数，迭代产生的尺度函数图像具有较好的光滑性（至少在视觉上是这样的），而且收敛速度明显比用Haar函数要快了许多。实现这一过程的是函数：daucascade(n,max1,t,p),及daucascadewa(n,max1,t,p),n表示n阶小波，max1表示迭代次数，t为变量，p为两尺度序列，因此，使用以上两个函数之前，要先用p=daucoef(n)进行赋值。例如：test501.m给出了daubechies6的前四次的迭代图形；test502.m给出了daubechies2的第四次的迭代图形。

下面是部分实验结果：

例1：用Cascade迭代法进行4次迭代得到的Daubechies 2、3、6尺度函数图像

图2-1 Cascade 4次迭代得到的Daubechies 2尺度函数图像

图2-2 Cascade 4次迭代得到的Daubechies 3尺度函数图像

图2-3 Cascade 4次迭代得到的Daubechies 6尺度函数图像test501.m

另外，Cascade算法并生的函数值在尺度函数的支撑区间之外自动地等于0，以db2 尺度函数为例：test502.m

图2-4 Cascade算法产生的db2尺度函数在支撑区间[0，3]之外全为0

下面再给出2-7阶的Daubechies 尺度函数的图像（Cascade迭代4次）：test503.m

然而，用Cascade 算法得到尺度函数的图像虽然要比用二进点迭代法来得直接，但是，当计算的阶数过高，迭代次数较大，并且所取点较密时，需要的时间也是很长的，甚至需要比二进点迭代长数倍的时间，原因就是计算过程中不断重复的迭代。因此，如果想要得到高阶的Daubechies 函数的图像，可以使用二进点迭代法，而且，相对而言，二进迭代得到的函数值都是精确的，而Cascade迭代则仅仅是一个逼近。因此，实际应用中，或许应用二进迭代会是一个更好的选择。

下面给出Daubechies小波函数的几个迭代图像

例2：用Cascade迭代法进行3次迭代得到的Daubechies 2、3、6小波函数图像test504.m

图2-5 Cascade 3次迭代得到的Daubechies 2小波函数图像

图2-6 Cascade 3次迭代得到的Daubechies 3小波函数图像

图2-7 Cascade 3次迭代得到的Daubechies 6小波函数图像

下面再给出用Cascade迭代得到的2-7阶Daubechies 小波函数图像test505.m

图2-8 Cascade 3次迭代得到的Daubechies 2-7小波函数图像

用Cascade算法的好处在于得到的是一个函数，我们可以在任意点上得到Daubechies小波函数的一个近似值，而且这种给出可以是独立的，而用二进点迭代法却必须得到所有低层次的二进点以后才能用之来出所求的值，而且当所求点不是二进点时，还必须进行转化，取近似，之后才能进行运算，这是Cascade算法所不要求的。所以单纯从一个点来说，我们可以迭代更多的次数以得到一个更加精确的值。而如果是要得到小波函数或是尺度函数的一个大致图形，那么，使用二进迭代在程序运行时间上会更加经济，因为Cascade算法的迭代次数太多，因而使运算次数翻倍地指数地增长。至少，给出图2-8的程序运行时间就要达到20分钟左右。而二进迭代2—10阶小波函数也只用了几分钟而已。然而，实际上看，图形的整体效果却是没有太大差别的。

第三部分

用Cascade 迭代的分解重构

前面提到，用Cascade算法可以得到一个函数，这个函数使得我们可以在任意点上取值，而如果用二进迭代，那么我们就不得不面对卷积了，而且由于二进点取到了值都是确定了的二进点，所以，在使用上还是有一定的复杂性的，在这次作业中没有对这种做法进行试验，相应地只使用Cascade 迭代算法进行试验，而且，我们知道，Cascade算法在迭代3-4次之后，对于本身的小波函数或是尺度函数已经达到了很好的逼近效果，所以我们可以采用Cascade算法进行试验，而且，我们可以预测到，这样处理所得到的结果在准确性上的误差也是不会太大的。

下面是我们对练习4-9的一个顶级为8的分解图：

这个分解图是使用Daubechies2 的Cascade算法4次迭代得到的尺度函数进行分解的。这里必须提到的一点是，我们可以清楚地看到，在图形的边缘，也就是接近横坐标为0的地方，总是有一个跳动，这是延拓所产生的副作用。而我们所用的延拓是对称延拓。或许应用其它的延拓方法，可以得到更好的图形。

图2-9 使用Daubechies2 的Cascade算法4次迭代的练习4-9分解

这里我们与Haar小波分解作一下对比：

图2-10 使用Haar小波 的练习4-9分解

接着是其各阶小波分量的图形：

图2-11 使用Daubechies2 的Cascade算法4次迭代的练习4-9小波分量图形

比较Haar小波分解

图2-12使用Haar小波 的练习4-9小波分解

从图形上看，在总体的逼近上，Daubechies小波逼近的效果要比Haar小波好，因为Haar小波的间断性，使得它在于Daubechies小波分解作对比时有些相形见绌，但是由于Haar小波不需要延拓，所以即使在边缘处，其误差也不会太大。或者说，由于其插值性，限制了其误差的波动。而Daubechies小波不但要做迭代逼近，还要做边缘延拓，这一步步的处理极大地损伤了Daubechies小波的优越性，实在是杯具了。

但即便如此，我们依然看到，Daubechies小波的分解效果依然比Haar小波好很多，其小波系数要比Haar小波系数小得多。相信在技术上加以完善，Daubechies小波的优越性将会得到极大的释放。

根据重构算法，给出重构算法原程序如下：

s=10; p=daucoef(2) syms soft

[a1,b1]=daurecon(2,8,soft,0) t=0:0.001:1;

f=exp(-t.^2/10).*(sin(2*t)+2*cos(4*t)+0.4*sin(t).*sin(50*t)); while s>1 g=0;

for k=1:(2^(s-2))

g=g+a1(s-1,k)*daucascade(2,4,2^(s-2)*t-k+1,p); end

subplot(3,3,s-1),plot(t,g,t,f),title(['the recon in V',int2str(s-2)]) grid on s=s-1 end

原来的程序在运行之后得到的不良结果如下，经检查调试之后，发现了问题出现在重构函数daurecon()的一个语句中，而并不是之前预测的边界对于整体的影响。在检查过程中，我认为边界点的值只能在边缘上形成一定影响，并不能对总体构成威胁，因此集中精神对函数daurecon()作了检查，最后排除了这个问题。

原程序零延拓结果如下：

原程序对称延拓结果如下：

修改后程序运行结果，应用daubechies2的重构图：test508.m 零延拓结果：

由原来的试验可知，对称延拓的效果并没有零延拓好，因此，修改后程序不再对对称延拓作试验，我们再用daubechies5作一次试验，以作比较。

Daubechies3重构图：

从图形上比较容易看出，daubechies2的重构效果要比daubechies3好，这种优势体现在重构图形的边缘两侧，显然daubechies2在边缘的图开效果要比daubechies3好很多。我猜想，由于分解算法在一端需要延拓，而重构在另外一端需要延拓，这样分解重构分别产生了一端的误差，这样，才会使得重构图形在两端有较大的误差。

因此，一方面，如果小波阶数越大，是有误差的区间将会相应加大；另一方面，如果想减少或消除这种误差，必要的边缘或延拓技术是必不可少的，因此，这依然存在一个很费脑筋的问题，这里就不再进一步思考了。

而由以上数据与图形结果终于可以证明这个耗时良久的实验可是算是取得成功了。高阶的daubechies小波分解重构需要在高性能的电脑上进行试验，否则将会因为程序运行时间过长而浪费了时间和精力。

实验完毕。