Mallat算法及C语言实现（一维DB小波分解与重构）

附：理论部分来源于书籍与网络资源，如有侵权，请于文章作者联系，立即处理！！！

一、Mallat原理分析

1.1矩阵变换角度看小波分解与重构：

DB4有4个系数：c0、c1、c2和c3。可以有以下变换矩阵，将它作用于一列数据矢量的左边。即（注意行列对齐，手打不易）

其中，空白处为0,。注意，此矩阵结构，第一行产生一个数据与滤波系数c0、c1、c2、c3卷积的分量，依次类推，第3行、第5行和其余奇数行的结果一样。如果偶数行以这种形式出现，正负交替，那么矩阵将是循环的，也就是一般的卷积，可以用FFT方法计算（注意，最后两行像具有周期边界条件下的卷积那样环绕起来）。但是，偶数行并不是以系数c0、c1、c2、c3进行卷积，而是以系数c3，-c2，c1，-c0进行卷积。整个矩阵的作用就是进行两个相关的卷积，然后各去掉一半数值，将剩下的各一半融合在一起。

将滤波器c0、c1、c2、c3看成是一个光滑滤波器，称它为H,H有点像4个点的移动平均。而将滤波器c3，-c2，c1，-c0称为G，由于G中滤波器带有负号，因此，G不能看成一个光滑滤波器。在信号处理的文献中，H和G被称为求积镜像滤波器。实际上，c值的选取是使G对光滑的数据矢量尽可能产生零响应。通过以上运算，导致H的输出在去掉以后，精确地代表了数据的“近似”信息。G的输出同样在去掉一半以后，代表了数据的“细节”信息。

为了能够恢复原始信息，必须能够从它的N/2个近似分量和N/2个细节分量重建长度为N的原始数据矢量。为了能够实现这样的算法，就必须要求矩阵式(1)为正交的，这样它的逆矩阵就和转置矩阵相同。其矩阵如下，即

为使矩阵式(2)成为矩阵式(1)的逆矩阵，必须满足以下两个等式，即

这样用长度为N（可以看出N为2的幂指数）的数据相乘，就能够分解成近似分量和细节分量，再对近似分量按照上面的方法进行分解，依次类推，从而实现Mallat算法。

 

1.2滤波器角度看小波分解与重构：

考虑如下问题。用N=2的Daubechies小波分解一个有n个取样值（即s0，….，sn-1）的信号。这些样值首先被当作顶级的近似系数aj，然后让其通过与Daubechies小波相关的高通和低通滤波器。这些滤波器分别称为H和L，响应的冲激响应函数序列分别为h和l：

分别用l和h与aj做卷积，然后下取样。那么对于任何aj（不仅是本例中的仅有8个样值的信号）有：

为计算j-1级的各个系数，需要4个相邻的样值，并且要从偶数位置开始。

实际中信号为离散序列，用小波分解信号序列，实际上就是让源信号通过对应的高通和低通滤波器，例如，用N=2的Daubechies小波分解一个8样值点s0，….，s7的信号序列，那么为计算近似系数a₂²，需要s4至s7.，为计算a₃²，需要s6至s9，但实际上并没有s8和s9！由此并不难理解，在此滤波过程中会很快用尽所有样值点，更麻烦的是，该过程仅在k=0,1,2时可进行。这意味着，对有8个样值点的问题来说，仅能得到3个分解系数，而不是原先期望的4个，该现象称为溢出问题，

溢出的根本原因在于不知道源信号前后的准确值，所以必须采取某种方式对源信号进行延拓。常用的延拓方法有零延拓、周期延拓、平滑延拓和对称延拓。

平滑延拓：在源信号两端用线性外插法补充样值，进行延拓。若信号信噪比比较高，或者至少在端点处如此，则该延拓方法是适宜的。

 

在Mallat算法的推导中，假定输入序列是无限长的，而实际应用中常常是分时采样，即输入序列为有限长，此时，滤波器系数与输入序列卷积就会“轮空”现象，因此有必要对原始信号进行边界延拓，减小边界误差，解决办法通常有补零延拓和周期延拓法。

（1）补零延拓是在输入序列的末尾补零。补零延拓的确定是可能认为的造成输入序列边界的不连续，从而使得较高频的小波系数很大。

（2）周期延拓法将原来有限长的输入序列拓展成周期序列。周期延拓可适用于任何小波变换，但可能导致输入序列边缘的不连续，使得高频系数较大，这种方式的延拓卷积后与源信号长度一致

（3）对称延拓（MATLAB默认采取这种方式）可避免输入序列边缘的不连续，但只适用于对称小波变换。

 

对于正交小波变换来说，补零延拓和对称延拓两种方法实现起来比较简单，但重建时会产生边界效应，而且分解的层数越多，产生的边界效应越显著。补零延拓给人一种跳跃的感觉。至于对称延拓，由于正交小波滤波器一般都是非对称性的（Haar小波虽然是正交的，但它是非连续的），重建图像给人一种错位的感觉。相比较而言，只有周期延拓方式可以得到比较精确地结果，它不仅能保证分解和重建正确计算，而且恢复的质量也好。不过，周期延拓虽然是三种方法中比较好的方法，但会导致信号边缘的不连续性，从而使得较高频率（子带）层的小波系数很大，即使信号本身相当平滑，从信号压缩的角度看，大的系数是希望避免的。信号的对称延拓可避免边缘的非连续性问题。然而，对称延拓只能和对称的小波滤波器一起使用。如果降低正交性要求，选择双正交小波变换，对称性延拓不失为一种好的方法，周期延拓可适用于任何小波变换，但可能导致输入序列边缘的不连续，使得高频系数较大。而对称延拓则避免了输入序列边界的不连续，是当前广泛采用的一种延拓方式

 

二、软件实现

离散序列的Mallat算法分解公式：

其中，H(n)、G(n)分别表示所选取的小波函数对应的低通和高通滤波器的抽头系数序列，从Mallat算法的分解原理可知，分解后的序列就是原序列与滤波器序列的卷积再进行隔点抽样而来。

离散序列的Mallat算法重构公式：

其中，h(n)、g(n)分别表示所选取的小波函数对应的低通和高通滤波器的抽头系数序列。

 

设输入信号为f(n)，长度为srcLen，滤波器长度为filterLen，边界延拓后信号长度为srcLen+2*(filter-1)，取值区间[-(filterlen-1), srcLen+filterLen-2]，下面给出信号边界处理几种方法的具体表达式如下：

周期延拓：

例1：以[1 2 3 4 5 6 7 8]长度为8的序列信号为例，当滤波器的长度为4时，其具体的延拓长度为6（单边为3）

[6 7 8 1 2 3 4 5 6 7 8 1 2 3]

 

对称延拓：

例2：以[1 2 3 4 5 6 7 8]长度为8的序列信号为例，当滤波器的长度为4时，其具体的延拓长度为6（单边为3）

[3 2 1 1 2 3 4 5 6 7 8 8 7 6]

 

零值补填：

[0 0 0 1 2 3 4 5 6 7 8 0 0 0]

 

2.1 详细的单级分解过程：

以DB2小波分解为例，通过对称延拓后的详细计算过程：

从上图可知，小波的Mallat算法分解后的序列长度由原序列长srcLen和滤波器长filterLen决定。从Mallat算法（采用对称拓延）的分解原理可知，分解后的序列就是原序列与滤波器序列的卷积再进行隔点抽取而来。即分解抽取的结果长度为(srcLen+filterLen-1)/2。

上述方法也是matlab默认采取的方法，在MATLAB中输入代码，查看结果便知：

对比MATLAB运行结果：

 

C语言实现：

/****************************************************************************************功能：一维小波一层分解(下采样)，提取奇数位置元素*输入参数:@srcData : 源信号*@resData : 分解结果*@srcLen : 源信号长度*@filter:滤波器系数矩阵*@filterLen : 滤波器的长度*@ fpAppCoef:存储近似系数的文件指针*@ fpDetCoef:存储细节系数的文件指针*返回值：无*说明 ： 最终分解结果的长度为延拓后信号长度（srcLen+(filterLen-1)*2）与滤波器卷积长度（filterLen-1）之和，例如信号长度为8，滤波器长度为4，信号经延拓再与滤波器相卷，最终长度为8+(4-1)*2+(4-1)=17***************************************************************************************/static uint DWT_Convlution_Decompose(pElement srcData, uint srcLen,pElement resData, uint resLen,Element (*filter)[DBNLEN], uint filterLen,FILE* fpAppCoef, FILE* fpDetCoef){uint i=0, j=0;uint idx = 0;uint extLen = 0, convLen=0, tmpLen=0;pElement pFilter[filterLen];    extLen = srcLen + (filterLen - 1) * 2;//延拓后信号长度    Element extData[extLen];//延拓后信号    memset(extData, (uchar) 0, sizeof(Element) * extLen);    signalExtend(srcData, srcLen, extData, extLen, METHOD_2, filterLen);//采用对称延拓方式延拓源信号// for(i=0; i<extLen; i++)// {// printf("%lf ", extData[i]);//打印延拓后的信号// }// printf("\n\n");/***************************************求解近似分量********************************************/convLen = extLen + filterLen - 1;  //卷积结果信号长度    Element convData[convLen];memset(convData , 0, sizeof(Element) * convLen);for(i=0; i<filterLen; i++){pFilter[i] = *(filter + 0) + i;}convLen = linear_Convlution(extData, extLen, pFilter, filterLen, convData);//线性卷积结果// for(i=0; i<convLen; i++)// {// printf("%lf ", convData[i]);// }// printf("\n\n");tmpLen = convLen - (filterLen - 1) * 2;//分解结果长度减去延拓部分    Element tmpData[tmpLen];memset(tmpData , (uchar) 0, sizeof(Element) * tmpLen);for(i=(filterLen - 1), j = 0; i<(convLen - filterLen + 1); i++, j++){tmpData[j] = convData[i];//printf("%lf ", tmpData[j]);}//printf("\n\n");    for (i = 1, idx = 0; i < tmpLen; i += 2, idx++)//去掉延拓部分，前后的长度为(filterLen - 1)    {resData[idx] = tmpData[i];  //下采样隔点采样，提取下标为奇数的元素，得到近似系数//printf("%lf ", resData[idx]);fprintf(fpAppCoef, "%lf,", resData[idx]);//将数据流写入文件    }//printf("\n\n");/**************************************求解细节分量*********************************************/convLen = extLen + filterLen - 1;  //卷积结果信号长度memset(convData , (uchar) 0, sizeof(Element) * convLen);for(i=0; i<filterLen; i++){pFilter[i] = *(filter+1) + i;}printf("\n");    convLen = linear_Convlution(extData, extLen, pFilter, filterLen, convData);tmpLen = convLen - (filterLen - 1) * 2;//分解结果长度减去延拓部分memset(tmpData , (uchar) 0, sizeof(Element) * tmpLen);for(i=(filterLen - 1), j = 0; i<=(convLen - filterLen + 1); i++, j++){tmpData[j] = convData[i];//printf("%lf ", tmpData[j]);}//printf("\n\n");    for (i = 1; i < tmpLen; i += 2, idx++)//去掉延拓部分，前后的长度为(filterLen - 1)    {resData[idx] = tmpData[i];  //下采样隔点采样，提取下标为奇数的元素，得到细节系数//printf("%lf ", resData[idx]);fprintf(fpDetCoef, "%lf,", resData[idx]);//将数据流写入文件    }//printf("\n\n");resLen = idx;return  resLen;}

 C语言运算结果：

2.2 详细的单级重构过程

同样以DB2小波为例，通过对称延拓后的详细重构过程：

从上图可知，小波的Mallat算法重构过程是分解过程的逆过程，分解过程得到低频信息（低频系数）和高频信息（高频系数），首先是对分解后的信号隔点插0，并作对称延拓，再与各自的滤波器进行卷积（低频信息与低通滤波器相卷，高频信息与高通滤波器相卷），各自得到的结果再相加，最终去掉延拓部分（filterLen - 1）即可得到原始信号。

注意在此过程中相对于最初的原始信号有过两次延拓，一个是分解过程中，首先是对原始信号进行延拓的，另一个是在重构过程中，对分解结果的信号进行插零后的延拓。

MATLAB运行结果：

为了便于对运算过程的理解，这里对信号的长度进行分析，举个例子：假设原始信号有srcLen = 8点，采用DB2滤波器，filterLen = 4。

分解过程：

（1）延拓：得到信号长度为extLen = srcLen+ (filterLen-1)*2=14

(2)与滤波器卷积：得到信号长度为convLen= extLen + (filterLen-1) = 17

(3)去掉延拓部分：得到信号长度为17 – 6 =11

(4)隔点取样的到分解结果：resLen  = (11-1)/2 = 5（低频或者高频信息）

重构过程：

（1）对低频或者高频信息（分解过程的结果，此时srcLen =5）隔点插零，从第0位置开始，得到信号长度就为11

(2)延拓：得到的信号长度extLen = srcLen+(filterLen-1)*2=17,

(3)与各自的滤波器相卷：得到信号长度convLen= extLen + (filterLen-1) = 20

(4)将上一步各自与高频滤波器和低频滤波器相卷积的结果相加得到结果，此时信号长度不变，还是20，只是多了前后两部分的延拓，因此还原出原始信号长度为20 - (filterLen-1)*2*2=8

 

由上面详细的分解和重构过程对信号长度进行分析，可见，采用Mallat算法对信号进行分解和重构过程，信号长度有很大的冗余，延拓卷积的结果长度是自身的2倍多。

C语言实现：

/****************************************************************************************功能：一维小波重构(上采样)，偶数位置插0*输入参数:@decData : 源信号，已被分解的信号*@decLen : 分解结果长度*@resData : 重构信号*@resLen：重构信号长度*@filter：滤波器系数矩阵*@filterLen : 滤波器的长度*返回值：无***************************************************************************************/static void DWT_Convlution_Rebuild(pElement decData, uint decLen,pElement resData, uint resLen,Element (*filter)[DBNLEN], uint filterLen){uint i=0, j=0;uint convLen = decLen + 1 + (filterLen-1)*2 + (filterLen-1);Element convData[convLen];memset(convData, 0, sizeof(Element) * convLen);uint tmpLen = (convLen - (filterLen - 1));//卷积结果长度减去(滤波器长度-1)Element tmpData[tmpLen];memset(tmpData, 0, sizeof(Element) * tmpLen);convLen = RebuildCoef(decData, decLen, convData, filter, filterLen);for(i = filterLen-1, j = 0; i < tmpLen; i++, j++){tmpData[j] = convData[i];}resLen = (tmpLen - (filterLen - 1) * 2); //再减去减去延拓部分for(i = filterLen-1, j = 0; i < resLen ; i++, j++){resData[j] = tmpData[i];}return ;}

运算结果：

综上，调试程序，运行结果（gcc编译器）：

程序调试成功，对一段音频信号进行分解和重构并且与MATLAB运行对比，最终运行结果良好，精度也很高。附上源程序以及测试数据：

http://download.csdn.net/detail/kuangzuxiaon/9852417

 

参考文献：

[1]Albert Boggess, Francis J.Narcowich. 小波与傅里叶分析基础[M].电子工业出版社.

[2]刘涛, 曾祥利, 曾军. 实用小波分析入门[M]. 国防工业出版社.

[3] http://www.cnblogs.com/IDoIUnderstand/archive/2013/01/30/3280725.html

[4]http://blog.csdn.net/ebowtang/article/details/40433861