小波分析算法的公式与C语言实现

　　一、小波分析算法的计算

　　1、Mallat算法［经典算法］

　　在小波理论中，多分辨率分析是一个重要的组成部分。多分辨率分析是一种对信号的空间分解方法，分解的最终目的是力求构造一个在频率上高度逼近L2（R）空间的正交小波基，这些频率分辨率不同的正交小波基相当于带宽各异的带通滤波器。因此，对于一个能量有限信号，可以通过多分辨率分析的方法把其中的逼近信号和细节信号分离开，然后再根据需要逐一研究。多分辨率分析的概念是S.Mallat在构造正交小波基的时候提出的，并同时给出了著名的Mallat算法。Mallat算法在小波分析中的地位相当于快速傅立叶变换在经典傅立叶变换中的地位，为小波分析的应用和发展起到了极大的推动作用。

　　MALLAT算法的原理

　　在对信号进行分解时，该算法采用二分树结构对原始输入信号x（n）进行滤波和二抽取，得到第一级的离散平滑逼近和离散细节逼近，再采用同样的结构对进行滤波和二抽取得到第二级的离散平滑逼近和离散细节逼近，再依次进行下去从而得到各级的离散细节逼近对，即各级的小波系数。重构信号时，只要将分解算法中的步骤反过来进行即可，但要注意，此时的滤波器与分解算法中的滤波器不一定是同一滤波器，并且要将二抽取装置换成二插入装置才行。

　　2、多孔算法

　　多孔算法是由M.shen于1992年提出的一种利用Mallat算法结构计算小波变换的快速算法，因在低通滤波器和高通滤波器中插入适当数目的零点而得名。它适用于的二分树结构，与Mallat算法的电路实现结构相似。先将Mallat算法的电路实现的基本支路作一下变形。令的z变换为，下两条支路完全等价，只不过是将插值和二抽取的顺序调换一下罢了。图中其它的上下两条支路也为等效支路，可仿照上面的方法证明。这样，我们便可由Mallat算法的二分树电路结构得出多孔算法的电路级联图，原Mallat算法中的电路支路由相应的等效支路所取代，所以整个电路形式与Mallat算法非常相似。如果舍去最后的抽取环节们实际上相当于把所有点的小波变换全部计算出来。

　　3、基干FFT的小波快速算法

　　Mallat算法是由法国科学家StephaneG.Mallat提出的计算小波分解与重构的快速算法，能大大降低小波分解与重构的计算量，因此在数字信号处理和数字通信领域中得到了广泛的应用。但是如果直接采用该算法计算信号的分解和重构，其运算量还是比较大。主要体现在信号长度较大时，与小波滤波器组作卷积和相关的乘加法的计算量很大，不利于信号的实时处理。故有必要对该算法作进一步的改进。众所周知，FFT是计算离散傅里叶变换（DFT）的一种快速算法，如能将它和Mallat算法结合在一起，势必会进一步降低小波分解和重构的计算量，事实证明这一想法是可行的。

　　基于FFT的小波变换快速算法是通过离散傅里叶变换建立起FFT和mallat算法之何的桥梁，从而将、FFT引入到小波变换中来，达到改小波变换快速算法及硬件实现的研究进Mallat算法的目的。

　　当信号长度较小时，FFT算法效率不及直接算法;随着长度的增加，特别是对于长度是2的幕次方的信号，FFT算法比直接算法更适用，能大大降低计算t。当信号是长序列信号时，小波分解与重构中，滤波器要补很多的零，这对信号的实时计算很不利，我们可以采用长序列快速相关卷积算法对信号进行分段后再运用FFT算法，提高运算速度。

　　4、基于算术傅里叶变换的小波变换快速算法

　　算术傅里叶变换（AFT）是1988年由Tufts和Sadasiv提出的一种用Mobius反演公式计算连续函数傅里叶系数的方法。它具有乘法运算t仅为O（N）算法简单、并行性好的优点。根据DPT和连续函数傅里叶系数的关系，可以用AFT计算DFT。同直接算法相比，APT方法可以将DFT的计算时间减少90%，尤其是对于含有较大素因子，特别是其长度本身为素数的DFT，它的速度比传统的FFT更快。另一方面，Mallat算法的分解和重构算法也可由DFT来计算，从而将AFT与Mallat算法联系了起来，从而为小波变换快速算法开辟了新的途径。

　　对于尺度

　　1）为j的快速分解算法步骤如下：

　　1）选定滤波器系数h（n）和g（n），再根据FFT的性质2，用N点的AFT分别计算出H（k）和G（k），分别取共扼，进而得到H*（k），G*（k）。

　　2）在已知cj（n）的情况下，用N点的AFT求出其DFTCj（k）

　　3）分别计算出H*（k）Cj（k），G*（k）Cj（k），即C’j（k）和D’j（k）

　　4）用N点的AFT求出C’j+1（k）和D’j+1（k）IDFT，得到C’j+1（n）和D’j+1（n）IDFT，再分别对它

　　们作二抽取，就可求出Cj+1（n）和Dj+1（n）。

　　在进行分解计算时，H（k）G（k）只要计算一次即可。重复步骤（2）一（4）可实现下一尺度小波分解，直到达到规定的尺度为止。不过要注意：尺度增加一个级别，信号长度减半。

　　2）对于尺度为j+1的快速重构算法为：

　　1）对Cj+1（n）和Dj+1（n）进行二插值，得到C’j+1（n）和D’j+1（n）;

　　2）用N点的AFT分别求出h（n）、g（n）的DFTH（k）和G（k）

　　3）用N点的AFT分别求出C’j+1（n）和D’j+1（n）的DFTC’j+1（k）和D’j+1（k）;

　　4）根据（17）式求出Cj（k），再用N点的AFT进行IDFT，可求出cj（n）。

　　5、基于Hermite插值的小波变换模极大值重构信号快速算法

　　信号在不同尺度上的小波变换模极大值包含了信号中的重要信息，因此研究如何由小波变

　　换模极大值重构信号是很有意义的。论文提出了一种基于Hermite插值多项式由二进小波变换模极大值重构信号的快速算法。数值试验表明，与S.Mallat提出的经典交替投影算法相比，该算法可以在保证重构质量的前提下简化计算过程，提高计算效率，计算所需时间与交替投影算法相比大大减少，是一种实用性较强的信号重构算法。

　　Hermite插值［11］方法是一种具有重节点的多项式插值方法，由于它要求在节点处满足相应的导数条件，因此也称为切触差值。由于小波系数模极大值点的导数为零，这与Hermite插值对节点的导数要求不谋而合，因此我们选用Hermite插值多项式作为改进的插值方法。

　　6、强奇异积分方程小波Petrov-Galerkin快速算法

　　通过构造具有高阶消失矩、小支集和半双正交性质的分片多尺度小波基底，给出第2类强奇异积分方程的小波Petrov-Galerkin快速算法，并证明该算法收敛阶达到最佳，条件数有界，计算复杂性几乎最佳。构造配置泛函的思想，构造分片多项式空间Xn上2列具有半双正交性的小波基，其中一列具有高阶消失矩性质。

　　二、小波分析算法的C语言实现

　　小波变换程序：

　　voidDB4DWT（doubleData［］，intn）

　　{

　　if（n》=4）

　　}

　　inti，j;

　　intbalf=n》》1;

　　double*tmp=newdouble［n］;

　　i=0;

　　for（j=0;j《half;j++）

　　{

　　tmp［j］=Data［（2*j）%n］*h0+

　　Data（2*j+l）%n］*h1+

　　Data［（2*j+2）%n］*h2+

　　Data［（2*j+3）%n］*h3;

　　tmp［j+half］=Datal（2*j）%n］*g0+

　　Data［（2*j+I）%n］*g1+

　　Datal（2*j+2）%n］*g2+

　　Data［（2*j+3）%n］*g3;

　　}

　　for（i=0;i《n;i++）

　　{

　　Data［i］=tmp［i］;

　　}}}

　　提升算法的程序

　　小波的分解算法的提升过程如下：

　　

　　使用提升算法的程序:

　　voidDB4LiftDWT（doubleData［］，intn）

　　应用研究

　　inti，half=n》》1;

　　double*pS=newdouble［half］;

　　//临时变量存放平滑系数

　　double*pD=new

　　double［half］;

　　//临时变量存放细节系数

　　for（i=0;i《half;it+）/赋值

　　{

　　pS［j］=Data［2*i］：//even，

　　pD［i］=Data［2*i+1］：//odd

　　}

　　for（i=0;i《half;i++）11DB4变换Update1

　　{

　　pS［i］=pS［i］+pD［i］*sqrt_3;

　　}

　　for（i=0;i《half;i++）//DB4的predict

　　{

　　pD［i］=pD［i］-sqrt_3*p$［i14

　　（sqt_3-2）*pS［（-1）》=0？（i-1）：（halfti-1）14;

　　//边界是采用周期延拓

　　}

　　for（i=0;i《half;i++）//DB4的update2

　　{

　　pS［i］=pS［i］-pD［（i+1）%half］;

　　/1边界采用周期延拓

　　}

　　for（i=0;i《half;i++）

　　{

　　p$［i］=（sqt_3-1）*pS［iV（sqnt_2）;/比例系数

　　pD［i］=（sqrt_3+1）*pD［i］（sqrt_2）;

　　}

　　for（i=0;i《half;i++）

　　{

　　Data［i］=pS［i］;//将平滑系数放回原数组

　　}

　　for（i=half;i《length;i++）

　　{

　　Data［i］=pD［i-half］;

　　/将细节系数放回原数组

　　}

　　deleteOpD;

　　delete［DpS;

　　在目述程序中用到了两个临时变量数组，虽然程

　　序结构清楚，容易阅读但消耗较大的内存空间，这里给

　　出另一种方法不用辅助数组，带来节约内存，但耗时。

　　voidsplit（doubleData［］，intn）

　　{intstart=1;

　　intend=n-1;

　　while（start《end）

　　{for（inti=start;i《end;i=i+2）

　　{doubletmp=Data［i］;

　　Data［i］=Data［i+1］;

　　Data［i+l］=tmp;}

　　start=start+1;

　　end=end-1;}}

　　对原数组调用split函数进行处理后，则Data数组的前半部分是偶数索引值，后半部分是奇数索引值。再对提升算法的程序进行一些修正就可以起到节约内存的目的。具体的实现限于篇幅不多介绍。利用提升算法的小波程序结构清楚，与Mallat算法相比运算量也小［3］。同时它的反变换也很易实现，这里由于篇幅限制，不对反变换作多的介绍。

　　本文给出了一维小波变换的C的实现，可在基础上实现：二维的小波变换，希望对要自己实现小波变换的读者有一定的帮助。