小波变换课件 第4章 小波变换的实现技术.docx

小波变换课件 第4章 小波变换的实现技术第4章 小波变换的实现技术 4.1 Mallat算法 %=h%、g%=g%n、h=hn、g=gn为实系数双双正交小波变换的Mallat算法：设hn%，g正交小波滤波器。h%是小波分析滤波器，h，g是小波综合滤波器。h表示h的逆序，即hn=h-n。若输入信号为an，它的低频部分和高频部分以此为an-1和dn-1，小波分解与重构的卷积算法： n-1n%)ìïa=D(a*h ín-1n%)=D(a*gïîdn-1 an=(Uan-1)*h+Ud(*)g 先进行输入信号和分析滤波器的巻积，再隔点采样，以形成低频和高频信号。对于有限的数据量，经过多次小波变化后数据量大减，因此需对输入数据进行处理。 4.1.1 边界延拓方法 下面给出几种经验方法。 1. 补零延拓 是假定边界以外的信号全部为零，这种延拓方式的缺点是，如果输入信号在边界点的值与零相差很大，则零延拓意味着在边界处加入了高频成分，造成很大误差。实际应用中很少采用。 .,0,s0,s1,s2,.,sn-1,0,0,. 14424432.简单周期延拓 将信号看作一个周期信号,即sk+n=sk。简单周期延拓后的信号变为 .,sn-1,s0,s1,s2,.,sn-1,s0,s1,s2,.,sn-1,s0,s1,s2,.,sn-1,s0,.42443144244314424431424314这种延拓方式的不足之处在于，当信号两端边界值相差很大时，延拓后的信号将存在周期性的突变，也就是说简单周期延拓可在边界引入大量高频成分，从而产生较大误差。 3. 周期对称延拓 这种方法是将原信号在边界上作对称折叠，一般分二 1）当与之做卷积的滤波器为奇数时，周期延拓信号为 2）当与之做卷积的滤波器为偶数时，周期延拓信号为 sn-1,sn-1,sn-2,.,s1,s0,s0,s1,s2,.,sn-1,sn-1,sn-2,.,s1,s0,s0,.424431442443144244314.sn-2,sn-1,sn-2,.,s3,s2,s1,s0,s1,s2,.,sn-2,sn-1,144424443144424443sn-2,.,s2,s1,s0,s1,.14424434. 光滑常数延拓 在原信号两端添加与端点数据相同的常数。 .,s0,s0,s0,s1,s2,.,sn-1,sn-1,sn-1,.142431442443142435. 平滑延拓 在原信号两端用线性外插法补充采样值，即沿着信号两端包络线的一阶导数方向增加采样值。 l 实际应用时，在变换前对输入信号进行边界进行延拓，使之变成无限长的信号，变换后，、 在尽可能不丢失信息的情况下，适当截取部分变换系数作为低频信号和高频信号，以保证小波分解后的数据总量保持不变。见下图。 为实现完全重构，先对有限长序列进行延拓，然后再插值和滤波，对滤波后的信号相加，再 适当截取，以恢复原信号。见下图。 4.1.3 用小波处理函数信号的基本步骤 1. 初始化 l 对于时间t的连续信号f(t)，选择适当的j=J，使得2J大于信号的抽样频率(不同 的应用决定了不同的抽样率)。 l 设信号在最高初始分辨率级J下的光滑逼近为fjPjf(t)ÎVj，则有 fj(t)=åckJkfj,k(t) -1由式, 既N1/2an»f(NLn)可得 ck»2J-J/2f(k/2) J在实际应用中，由原信号确定的k的范围是有限的，譬如信号的持续时间为0£t£1，则k的范围为0£k£2J-1。 2. 小波分解 应用Mallat 算法，得到离散信号的小波变换c0,d0,d1,.,dJ-1，相应地，得到f(t)的分 辨率表示： fJ(t)=f0(t)+d0(t)+d1(t)+.+dJ-1(t) 其中f0(t)ÎV0，d0ÎW0，d1ÎW1，.， dj(t)ÎWj。具体地 f0(t)=c0f0,0(t)=c0f(t),dj(t)=00kÎkådkyj,k(t),j=0,1,.,J-1 j实际应用中，可以根据需要控制分解的级数，不一定达到j=0级。 3. 小波系数处理 针对不同的应用目标，对小波系数c0,d0,d1,.dJ-1进行处理获得新的小波系数%0,d%1,.d%J-1。譬如，在进行信号的数据压缩时，将绝对值小于某一阈值的系数置为%0,dc零，保留剩余的系数，用于重构信号；而在去噪时，将绝对值小于某一阈值的系数置为零，用于重构去噪信号。 4. 小波重构 %0,d%1,.d%J-1，%J。%0,d对处理后的小波系数c重构出分辨率j=J时的离散信号c一般%J是cJ的逼近信号。进而可以得到fJ(t)或f(t)的重构信号fJ(t)。 地，c对于离散信号的小波处理过程为 ，设bn是一个离散输入信号，采样间隔为N-1，其中N=2。可将bn与f(t)=LåanfL,n(t)ÎVL联系起来，使bn为f(t)的均匀采样，即f(t)f(Nn)。根据式可得NL1/2Lan»bn。由此可获得f(t)在最高分辨率L下的初始系数序列an。然后，利用Mallat算法对该序列进行小波分解、对小波系数处理以及处理后的系数进行小波重构等。 4.1.4 应用举例 例4.1对单位区间上一个连续信号，将信号离散化为28个采样值，相应的逼近信号记为 f8(t)。用Haar小波对信号进行3级小波分解，写出信号f8(t)的多分辨表示，并画出该信 号在不同分辨率下的逼近信号f7(t)、f6(t)和f5(t)的图形。 假设信号为f(t)=sin(4pt)+12cos(10pt)，它在V8中的投影记为f8(t)，则f8(t)的图 形见图42a。用Haar小波对信号进行3级小波分解，其多分辨表示为 f8(t)=f7(t)+d7(t) =f6(t)+d6(t)+d7(t) =f5(t)+d5(t)+d6(t)+d7(t) 其中f7(t)、f6(t)、f5(t)波形如图42b、c、d所示。 图4-2 一个函数的多分辨逼近函数 例4.2 对于例4-2中的信号及逼近信号。若用正交尺度函数和正交小波函数进行小波分 析解，则可得到： f(t)=f0+d0+d1+d2+,.,+d7 0其中，f0(t)=c0f(t)ÎV0，dj(t)=ådkÎkjkyj,k(t)ÎWj，j=0,1,.,7 1）用Haar尺度函数和Haar小波函数分解信号，令绝对值最小的80%和90%的系数为0，对信号进行小波压缩。画出相应的重构信号波形，并求出相应的相对误差。 2）用Daubechice尺度函数和Daubechice小波 (如db2)分解信号，令绝对值最小的80%和90%的系数为0，对信号进行小波压缩。画出相应的重构信号波形，并求出相应的相对差。 3）比较1）和2）的压缩效果。 4）用FFT在相同条件下压缩信号，所得的相对误差如何。 求解过程如下： 1） 用Haar小波函数分解信号，令绝对值最小的80%和90%的系数为0，得到重构信号图形如图4-3a所示。所得的均方差为0.7991；相对误差为0.0050。 如果令绝对值最小的90%的系数为0，得到重构信号图形如图4-3b所示。 在这种情况下，得到均方误差为2.9559，相对误差为0.0185。 图4-3 用Haar小波压缩后的重构信号 2） 用Daubechice小波函数(如Db2)分解信号，令绝对值最小的80%的系数为0，得到的重构信号波形如图4-4a所示，得到均方误差为0.0277，相对误差为0.00017。如果令绝对值最小的90%的系数为0，得到重构信号图形如图4-4b所示，这时得到均方误差为0.2159，相对误差为0.0014。 图4-4 用Daubechice小波压缩后的重构信号 3） 比较1）和2）的压缩效果可以看出，对于同一种小波，保留更多的变换后的系数可以得到更好的重构信号。对于相同比例的保留系数，用Daubechice小波分解信号，再重构，得到的效果好于Haar小波。这是因为Daubechice小波的连续性较好，更适合处理连续性较好的信号。 4) 用FFT对信号进行处理。令绝对值最小的80%的系数为0，则重构信号的图形如图4-5a所示。得到均方误差为0.0025，相对误差为1.59´10-5。 图4-5 用FFT压缩后的重构信号 注释：上两例是说明用小波处理信号的基本过程和在压缩中的应用，并不是想与FFT比较谁的效果好。实际上，这里采用的信号周期性很好，用傅立叶变换处理更有优势。一般地，小波更适合处理突变信号，而傅立叶变换更适合处理周期信号。 4.2 多孔算法 l Mallat算法存在的问题是 数据逐级减少问题。 原因是逐级二抽样，每经过一 级小波分解，数据减少一半，因 此，随着分辨率减少，低频分量 的数据越来越少。 l 多孔算法 两通道Mallat算法等价的z变换如图47表示。 记H%(z)=h%(z-1)，G%(z)=g%(z-1) l 二分树算法Mallat算法的小波分解迭代过程如图4-8所示。其中， x%j-1(z)=xj(z)H%(z)，d%j-1(z)=xj(z)G%(z) l z变换的等效易位性质： x(z)h(z)y(z)x(z)¯2Ûh(z2)y(z)¯2因为 左边 y(z)=1/22x(1z/2+)x-(z1)h (z)右边 y(z)=12x(1z/2)h(z+)-x(1/z2)h (z)图4-9 Mallat 算法的一种等效形式 %如果不考虑每个分支的最后的抽样环节，则dj-1%(z)，dj-2%(z)，dj-2%(z)，dj-3(z)，.相当于xj(z)中各点的小波变换全部计算出来，这叫非抽样小波变换。如图4-9所示。 h(z2)表示在滤波器hn的任意两点间插入2j-1个零所得到的滤波器Z变换，所以，非%，g%个相邻点之间插入2j-1个零再与低频信号做卷积，抽样小波变换就是把滤波器h故称j%每两个%(z4)是将g%(z2)是将h%每两个点之间插入三个零得到的新滤波器，H多孔算法。G点之间插入一个零得到的新滤波器，这样就把每一级的抽取移到了最后，保证了总数据不会%，g%(z4)和H%(z)是使h%中补零 ，增加了逐级减少，有效地实现了Mallat算法。由于G2空隙，故称多孔算法。 设原始信号x=x0,x1,.,x2J-1%J=x，根据下面的分解算法： 的长度为2J，记xJ%j-1%(z2-1)ì%j(z)×Gïd(z)=x íj-1Jj2-1%(z%(z)=x%(z)×H)ïîx%k的伪码程序： %0，d计算各抽样点处的小波变换xj=JWhilej>0%j-1=x%J-j%j-1*hd%xj-1%*h=xjJ-jj=j-1endofwhile从多孔算法的分解过程）可知， J-jJ-jJ-jj22ìd%j-1(z)g(z2%(z)×G(z)=x)g(z)ï íJ-jJ-jJ-jj-12j22%(zï%(z)h(z%(z)×H)=x)h(z)îx于是，由完全重构条件）可得 %j(z)=x12%xj-1(z)h(z2J-j%)+dj-1(z)g(z2J-j) 4.3 小波变换的提升实现 优点 可以实现更快速的小波变换算法，一般比Mallat算法快2倍。 可以实现完全的同址运算。 能很好地克服小波变换的边界问题。 提升算法小波变换的描述简单，可以避免使用傅里叶变换。 在时域或空域直接实现小波构造，既工程师可以按着自己的要求来构造不同的小波，不再紧紧依赖数学家。 4.3.0小波变换的提升实现 l 回顾 Haar小波变换 按着Haar小波滤波器组，两个数a,b的平均与细节分别为 s=a+b2d=b-a 于是，a,b的小波变换为s,d。如果a，b高度相关，则d很小。由s,d恢复a,b的计算公式如下： a=s-d/2 b=s+d/2 n对于长度为2的信号s=sn,l|0£l<2，将求平均与细节运算应用到每对数据nna=sn,2l，b=sn,2l+1(l=0,1,2,.,2n-1-1)上，记 ìssïn,2l+sn,2l+1ín-1,l=2 ïîdn-1,l=sn,2l+1-sn,2l这里是将信号序列sn的偶序号点和奇序号点相加取平均得平均值sn-1,l；将奇序号点减去偶序号点得到差值dn-1,l。sn-1,l和dn-1,l各有2n-1个样本，看作是信号sn的概貌和细节。当l遍历0到2n-1-1之间的所有整数时，sn-1,l组成sn-1序列，可以对sn-1进行类似的分解，得到s2-2n-2,l和dn-2,l，它们各有2n-个样本，sn-2,l组成了sn序列。这样的分解可以进行n次，最后的概貌信号只有一个点s0=sn0,0，它是信号s的均值或直流分量。最后的概貌信号加上各级细节信号，正好是2n个。上述变换正是Haar小波变换。多级Haar小波变换的分解过程和重构过程如下图。 图1-21 多级小波变换过程示意图 l Haar小波变换的提升实现 小波变换的原位实现-不增加额外空间 d=b-a s =a+d/2 (s=a+d/2=a+b/2- a/2= (a+b)/2) l 提升算法中的信号分解 提升过程分为三步：分裂、预测和更新。 考虑信号s=sj,l|0l<2提升算法分三步完成 jj一级变换后®j-1j-1，ds 概貌信号细节信号Split(分裂) 将信号简单地分为两部分，有多种方法，这里将信号序列按奇、偶分为两个子集： j-1ìïoddj-1=sj,2l+1|0£l<2 íj-1even=s|0£l<2ïj-1j,2lî只将完整信号序列分成二部分，不做其他处理，故称懒小波变换。这里Split是分裂算子。 这种分解方法可表示为 (evenj-1,oddj-1):=Split(sj) 注释：将信号分解为成两个序列的方法有很多，将信号序列按奇、偶分为两个子集是一种方法；或将前一半和后一半分成两个子集也是一种方法；或将相邻两个数之和分给一个子集，而相邻两个数之差分给另一个子集又是一种方法。总之，不同的分裂方法，相当于采用不同的小波基。实际应用中，最常用的是懒小波最为第一级分解。 Predict(预测) 若原信号sj具有局部相关性，则子集evenj-1和oddj-1也具有相关性。因此，只要知道其中任一个，就可以以合理地精度预测另一个，通常用偶子集预测奇子集。在Haar小波变换的情况下，预测误差为dj-1,l=sj,2l+1-sj,2l，特别地，如果原信号是一个常量，则所有预测误差均为零。一般情况，定义预测算子P，且 dj-1=oddj-1-P(evenj-1) j-1P(evenj-1)表示用evenj-1值的某种运算或某种组合来预测odd的值。在Haar小波变换的情况下，简单地选择偶次项sj,2l，去减奇次项，得差值dj-1,l。预测误差表示信号的细节信息。这一步骤在提升算法中被成为“对偶提升”。当信号的相关性较大时，预测是非常有效的，若信号为常数时，dj-1,l恒为零。 Update 低频概貌信号的一个关键性质是，它与原信号应具有相同的平均值，即 s=2-jj2-1åsj,l，并且不随j变化。这能确保最后的变换系数是原信号的总平均。Updatel=0操作可保证该性质成立，既用细节信息子集dj-1来更新偶序号子集evenj-1，既 sj-1=evenj-1+U(dj-1j-1) 式中U为更新算子，表示对d测误差信号dj-1的某种运算或某种组合。对于Haar小波，可以简单地用预更新偶子集信号evenj-1，既sj-1=sj,2l+dj-1,l/2。 在提升算法中，更新称为“原始提升”，故才有“提升算法”一词。 容易验证， 2-(j-1)2j-1-12j-1-1åsj-1,l=l=0å-j(sj,2l+dj-1,l/2) j-1jl=02-1=2å(sj,2l+sj,2l+1)=2-j2-1l=0åsj,l l=0以上三步操作相当于对信号sj进行一级小波变换，将它分解为低频信号sj-1和细节信号dj-1。一般地，对于提升算法存在分裂算子Split、预测算子P和更新算子U，使 (evenj-1,oddj-1)=Split(s)dj-1j=oddj-1-P(evenj-1) =evenj-1+U(dj-1sj-1)上述所有操作可以实现原位实现，即偶数位置用平均值重写，奇数位置用细节重写。 伪码实现： (oddj-1,evenj-1):=Split(s)oddj-1-=P(evenj-1)evenj-1+=U(oddj-1)j提升方案实现小波分解的最大优点在于将小波分解成了几个简单的基本步骤，且每个步骤都非常容易找到它的逆变换。 l 提升算法中信号的重构过程 相反地，从低频信号s(1)反更新 j-1和细节信号dj-1恢复原信号s的提升算法为 j给定sj-1和dj-1，可由下式恢复出偶序号序列 j-1 evenj-1=s反预测 -U(dj-1) 用反更新计算出的偶序号序列evenj-1和给定的细节d oddj-1=d合并 j-1j-1，可通过下式预测出奇序号序列 +P(evenj-1) 通过反更新和反预测步骤，分别获得偶序号序列和奇序号序列，将它们合并即可恢复原信号。这一步骤称为invert Lazy wavelet transform，记作 s=Merge(evenj-1,oddj-1)逆懒小波变换 对于Haar小波，恢复 sj,2l=sj-1,l-dj-1,l/2 sj,2l+1=dj-1,l+sj,2l j4.3.1小波分解与重构的多相位表示 %、g%、h、g为双正交小波滤波器组，对应讨论有限冲激响应的双正交滤波器的情况。设h二通道Mallat算法的等价Z变换如图4-7所示。 从图3-2到图3-3，是将时域表示成z域，图中z-1与时域中时序反转相对应。 l %、h、g%、g的约束条件 两通道分析滤波器和综合滤波器在理想重构条件下，h设x=xn是一个信号序列， z变换为x(z)=åxzn-n。当xn是一个有限序列时，称x(z)是一个Laurent多项式。 设小波滤波器h=hn的Z变换为h(z)=åhnzn-n。可推出h(z)=h(z)。同理， -1%的巻积x*h%的Z变换等于它们Z变换的乘积，既 %(z)=h%(z-1)。另外，由于x=x与hhn(x*h%)(z)=x(z)h%(z)=x(z)h%(z%(z-1)=令x(z)h-1) 于是，图3-2的滤波器组的等价的Z变换形式如图3-3。 åckzk-k，则根据二元下采样的定义，有 a(z)=%å(D(x*h)kzk-k=åc2kzk-k从而 a(z)=1%(z-1/2)x(z1/2)+h%(-z-1/2)x(-z1/2)ù 的结果 éhû2ë又y0(z)=a(z2)h(z)，故 y0(z)=1%-11%-1h(z)h(z)x(z)+h(-z)h(z)x(-z) 221212类似地，有 y1(z)=%(z-1)g(z)x(z)+g%(-z-1)g(z)x(-z) g因此， y(z)=y0(z)+y1(z) 1%-11%-1%(z-1)g(z)x(z)+h%(-z-1)g(z)x(-z) h(z)h(z)+g(-z)h(z)+g224444424444413由x(-z)引起的混叠=x(z)之前第1个求和项应等于1，是确保y(z)=x(z)的一个条件；必须消除由x(-z)引起的混叠，既第2个求和项应消除。于是，滤波器组对任何输入信号实现精确重构，下式是二通道滤波器组完全重构条件，既PR (Perfect Reconstruction)条件。 %-1ì%(z-1)g(z)=2ïh(z)h(z)+g í-1-1%(-z)g(z)=0ïîh(-z)h(z)+gl PR条件的矩阵形式： %(z-1)éhê-1ëh(-z)%(z-1)ùéh(z)ùé2ùgú=êú -1úêg(z)%g(-z)ûëûë0û%(z-1)éh当矩阵ê-1h(-z)ë%(z-1)ùg%(z-1)g%(-z-1)g%(-z-1)-h%(z-1)¹0时，综合的行列式D(z)=h-1ú%(-z)ûg%、g%确定。于是， 滤波器h、g完全由分析滤波器héh(z)ù2=êúëg(z)ûD(z)ég%(-z-1)ùêú -1%êë-h(-z)úû所以， h(z)=2V(z)%(-z-1)，g(z)=g%(-z-1)ù é-hûV(z)ë2将式代入式，并注意到D(z)=-D(-z)，可得 %(z-1)+h(-z)h%(-z-1)=2 h(z)h表明，当D(z)¹0时，完全重构条件等价于重构条件式和式同时成立。 l %,h,g%,g的完全重构条件： 有限长度滤波器h-1对于有限长滤波器，根据定义，D(z)是Z的Laurent多项式，而由式知，D(z)也是Z的Laurent多项式，因此，D(z)必是一个单项式。又因为D(z)=-D(-z)，故D(z) 是一个奇数次的单项式，既 D(z)=-2az2l+1,aÎR,lÎZ 代入式，整理得 %(-z-1) %(z)=az-(2l+1)h(-z-1)，g(z)=a-1z-(2l+1)hg ¨ 一种取法是a=-1，l=0，于是式变为 %(-z-1) %(z)=-z-1h(-z-1)，g(z)=-z-1h g¨ 另一种取法是a=1，l=0，于是式变为 %(-z-1) %(z)=z-1h(-z-1)，g(z)=z-1h g%,h,g%,g的完全重构条件为 按式，有限长度滤波器h%(z-1)+h(-z)h%(-z-1)=2ìh(z)hï%(z)=-z-1h(-z-1) ígï-1%-1îg(z)=-zh(-z)l 多相表示的基本思想： ¨ 一个多相表示的例子 对于多项式 X(z)=1+1.5z取M=2，可写成 X0(z)=1+2.2zX1(z)=1.5+4z-1-1-1+2.2z-2+4z-3+2.2z-4+1.5z-5+z-6+2.2z+1.5z-2-2+z-3X(z)=(1+2.2z-2+2.2z-4+z-6)+z(1.5+4z-1-2+1.5z-4) 用z2替换z，得 =X0(z)+z2-1X1(z) 2这就是多项式的两相表示。 ¥¨ 设滤波器h(z)=k=-¥¥åhkz-k，将其分裂成z的偶次幂和奇次幂二部分： h(z)=k=-¥å¥h2kz-2k¥+k=-¥åh2k+1z-(2k+1)¥=k=-¥åh2kz-2k¥+k=-¥åh2k+1z-2kz-1 =k=-¥å¥h2kz-2k¥+z-1k=-¥åh2k+1z¥-2k =k=-¥¥åh2k(z)2-k+z-1k=-¥¥åh2k+1z(2-k) 定义，he(z)=k=-¥åh2k(z)-k和h0(z)=¥k=-¥¥åh2k+1(z)-k则 he(z)=从而有 2k=-¥åh2k(z)2-k，h0(z)=2k=-¥åh2k+1(z)2-k2) h(z)=he(z+-1z0h(z ) 2其中，he(z)包含了h(z)的所有偶系数，而ho(z)包含了h(z)的所有奇系数。更进一步，任一给定整数M，可将h(z)分解成模M不同余数次幂的M部分。即： ¥¥¥ h(z)k=-¥M-1åhkMz-kM+k=-¥åhKM+1z-KM+z-M-k=-¥åhKM+M-1z-KM可简洁地写为：h(z)å-kzhl(z-lM)-第一类多相表示； l=0¥式中，hl(z)开始。 n=-¥åhlkz是第一类多相表示h(z)的元素，对于因果序列，求和下限可从0定义h和g的多相位矩阵为 éhe(z)p(z)=êëho(z)ge(z)ùú go(z)û于是， 1éh(z)p(z)=ê2ëg(z)2Th(-z)ùúg(-z)ûé1êë1zùú -zû%和g类似地，h%的多相位矩阵p%(z)为 %(z)éhe%(z)=êp%êëho(z)2T%e(z)ùgú %o(z)úgû%(z)1éh%(z)=ê同理， p2ëg%(z)%h-(zù)é1úê%g-(zû)ë1zùú -zû因此，小波重构的完全条件式可以写成 %(z)=I2´2 p(z)p-1T%(z-1)T表示矩阵p%(z-1)的转置矩阵，I2´2为单位阵。 其中，pl 小波分解与重构的多相位表示 解释图4-11和图3-3的等价性。设输入序列xn的Z变换为x(z),则流程图(图4-12)的作用相当于对xn进行懒小波变换 , 即抽取xn的偶序列和奇序列,Z变换分别为 e(z)=o(z)=1212x(zz1/21/2)+x(-z1/2-1/2)x(z)-z1/2x(-z1/2)%(z-1)T作用后的低频和高频部分分别为a(z)和d(z)，则 设这两个Z变换经p%(z)p-1Tée(z)ùéa(z)ùú êú=êd(z)o(z)ûëûë可以验证a(z)=此验证。 1%-1/21/2%(-z-1/2)x(-z1/2)与图3-3的结果一样。d(z)也可作h(z)x(z)+h24.3.2 Laurent多项式的Euclidean算法 l %、g%都是有限长的小波滤波器，所以p(z)和p%(z)的行列式detp(z)和 由于h、g、h%(z-1)T=I2´2知，detp(z)及其倒数都是Laurent%(z)都是Laurent多项式。由式p(z)pdetp多项式，故detp(z)为Z的单项式。设detp(z)=1，我们根据小波分解和重构的多相位表示，通过对多相位矩阵p(z)进行因子分解，给出小波变换的实现算法。 l 有限冲激响应滤波器FIR可以描述为一系数集h=hk:kÎz，范围为kb£k£ke。 FIR滤波器h的Z变换是一个Laurent多项式h(z)，既 h(z)=åhkzkbke-k于是，h(z)的次数为 h(z)=ke-kb 因此，滤波器的长度等于其相应Laurent多项式的次数加1。 l 两个Laurent多项式的带余数除法 对于任何两个Laurent多项式a(z)和b(z)，其中，b(z)¹0，a(z)³b(z)，一定存在Laurent多项式的q(z)和r(z)，使a(z)=b(z)q(z)+r(z)成立。其中， q(z)=a(z)-b(z),r(z)<b(z)或r(z)=0。两个Laurent多项式的商和余数不是唯一的。 例如，对于a(z)=z-1+6+z，b(z)=4+4z，则对于以下几种情况： 1）q(z)=1/4(z-1+5)，r(z)=-4z。 2）q(z)=1/4(z-1+1)，r(z)=4 3）q(z)=1/4(5z-1+1)，r(z)=-4z-1 都满足a(z)=b(z)q(z)+r(z)，且q(z)=a(z)-b(z),r(z)<b(z)。 l Euclidean算法如下： 利用带余数除法，可以给出Laurent多项式的Euclidean算法。 对于任何两个Laurent多项式a(z)和b(z)，其中，b(z)¹0，且a(z)³b(z)。设a0(z)=a(z)，b0(z)=b(z)，从i=0开始进行如下的递归运算： ai+1(z)=bi(z) bi+1(z)=ai(z)%bi(z) 其中，%是表示取“余数”运算。则an(z)=gcd(a(z),b(z)，且an(z)是一个Laurent多项式，下标n是使bn(z)=0的最小数，“gcd”表示取最大公因子。 假设bi+1(z)<bi(z)，则存在m使得bm(z)=0，因此，算法在n=m+1步骤结束，其中，n£b(z)+1。若记qi+1(z)=ai(z)/bi(z)，其中，“/“表示取商运算，则有 1éan(z)ùêú=Õë0ûi=né0êë1ùéa(z)ùúêú -qi(z)ûëb(z)û1这等价于 néa(z)ùé-qi(z)=êúÕêb(z)ëûi=1ë11ùéan(z)ùúêú 0ûë0û显然，an(z)同时整除a(z)和b(z)。如果an(z)是一个单项式，则a(z)和b(z)是互为素数的。 例4.3 a(z)=z-1+6+z，b(z)=4+4z，则由第一步带余除法，可得 a(z)=4+4zb(z)=4q1(z)=1/4z-1 +1/4下一步带余数除法，给出 a(z)=4b2(z)=0q2(z)=1+z所以，a(z)和b(z)是互为素数的，且 éz-1+6+zùé1/4z-1+1/4êú=ê1ë4+4zûë1ùé1+zúê0ûë11ùé4ùúêú 0ûë0û辗转除法的步数为n=2=b(z)+1。 4.3.3 多相位矩阵的因子分解 下面的定理奠定了小波提升实现的基础。 l 定理 4.1 若p的行列式等于1，既detp(z)=1，则总存在Laurent多项式ui(z)和 pi(z)(1£i£m)及非零常数k，使得 mp(z)=Õi=1é1êë0ui(z)ùé1úê1ûëpi(z)0ùékúê1ûë00ùú 1/kû其中，pm(z)=0. 定理证明略。主要介绍提升因子ui(z)和pi(z)(1£i£m)的计算方法。 首先，对he(z)和ho(z)应用Euclidean算法，可得到 éhe(z)ùêú=ëho(z)ûm-1Õi=1é1êë0ui(z)ùé1úê1ûëpi(z)0ùékùúêú 1ûë0û记 éhe(z)p(z)=êëho(z)0ge(z)ùú=0go(z)û0mÕi=1é1êë0ui(z)ùé1úê1ûëpi(z)0ùékúê1ûë00ùú 1/kû注意到 é1êë02ks(z)ùékúê1ûë00ùékú=ê1/kûë0é1ë0ks(z)ùékú=ê1/kûë00ùú1/kûé1êë0s(z)ùú 1û则由式可得，p(z)=p(z)ê0s(z)ù2s(z)=u(z)/k，其中。 mú1ûm-1若记Q(z)=Õi=1é1êë0ui(z)ùé1úê1ûëpi(z)0ùéa(z)ú=ê1ûëkb(z)0c(z)/kùú，则 d(z)/kûékêë00ùéka(z)ú=ê1/kûëkb(z)c(z)/kùú d(z)/kûéhe(z)p(z)=êëho(z)0ge(z)ùú=Q(z)0go(z)û于是，有 é1êë0s(z)ùú=1û(p0(z)-1éka(z)p(z)=êëkb(z)c(z)/kùéhe(z)úêd(z)/kûëho(z)ge(z)ùú go(z)ûc(z)édzh(z)-h0(z)e=êkkêêë-kb(z)he(z)+kh0(z)a(z)ùg0(z)ú kkú-kb(z)