1 Mallat演算法

離散序列的Mallat演算法分解公式如下：

其中，H(n)、G(n)分別表示所選取的小波函式對應的低通和高通濾波器的抽頭係數序列。

從Mallat演算法的分解原理可知，分解後的序列就是原序列與濾波器序列的卷積再進行隔點抽取而來。

離散序列的Mallat演算法重構公式如下：

其中，h(n)、g(n)分別表示所選取的小波函式對應的低通和高通濾波器的抽頭係數序列。

2 小波變換實現過程(C/C++)

2.1       小波變換結果序列長度

      小波的Mallat演算法分解後的序列長度由原序列長SoureLen和濾波器長FilterLen決定。從Mallat演算法的分解原理可知，分解後的序列就是原序列與濾波器序列的卷積再進行隔點抽取而來。即分解抽取的結果長度為(SoureLen+FilterLen-1)/2。

2.2       獲取濾波器組

對於一些通用的小波函式，簡單起見，可以通過Matlab的wfilters(‘wavename’)獲取4個濾波器；特殊的小波函式需要自行構造獲得。

下面以db1小波函式(Haar小波)為例，其變換與重構濾波器組的結果如下：

//matlab輸入獲取命令
>> [Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('db1')

//獲取的結果
Lo_D =
    0.7071    0.7071
Hi_D =
   -0.7071    0.7071
Lo_R =
    0.7071    0.7071
Hi_R =
    0.7071   -0.7071

2.3       訊號邊界延拓

        在Mallat演算法中，假定輸入序列是無限長的，而實際應用中輸入的訊號是有限的取樣序列，這就會出現訊號邊界處理問題。對於邊界訊號的延拓一般有3種方法，即零延拓、對稱延拓和週期延拓。

3種延拓方法比較情況如下：

        對於正交小波變換來說，前兩種延拓方法實現起來比較簡單，但重建時會產生邊界效應，而且分解的層數越多，產生的邊界效應越顯著。零延拓方法給人一種跳躍的感覺。至於對稱性延拓，由於正交小波濾波器一般都是非對稱性的（Haar小波基雖然是正交的，但它是非連續的），重建圖象給人一種錯位的感覺。相比較而言，只有最後一種延拓方式可以得到比較精確的重建結果，它不僅能保證分解與重建正確計算，而且恢復的質量也好。不過，週期性延拓方法雖然是常用的三種方法中比較好的方法，但會導致訊號邊緣的非連續性，從而會使得較高頻率（子帶）層的小波係數很大，即使訊號本身相當平滑。從訊號壓縮的角度看，大的係數是希望避免的。訊號的對稱延拓可避免邊緣的非連續性問題。然而，對稱延拓只能和對稱的小波濾波器一起適用。如果降低正交性要求，選擇雙正交小波變換，對稱性延拓不失為一種好的方法。週期延拓可適用於任何小波變換，但可能導致輸入序列邊緣的不連續，使得高頻係數較大。而對稱延拓則避免了輸入序列邊界的不連續，是當前廣泛採用的一種延拓方法。

下式中給出了最常用的對稱延拓表示式。

當原序列長sLEN為偶數時延拓後的序列長為sLEN+2*(filterLEN)，而原序列長為奇數時則需要在右端再延拓一個元素。注：在Matlab中預設使用了對稱延拓。

2.4       小波分解

在db1小波函式下，離散序列的Mallat演算法分解公式展開如下：

其它的db小波，不再詳述。小波分解C++原始碼如下。

/**
 * @brief 小波變換之分解
 * @param sourceData 源資料
 * @param dataLen 源資料長度
 * @param db 過濾器型別
 * @param cA 分解後的近似部分序列-低頻部分
 * @param cD 分解後的細節部分序列-高頻部分
 * @return
 *         正常則返回分解後序列的資料長度，錯誤則返回-1
 */
int Wavelet::Dwt(double *sourceData, int dataLen, Filter db, double *cA, double *cD)
{
    if(dataLen < 2)
        return -1;
    if((NULL == sourceData)||(NULL == cA)||(NULL == cD))
        return -1;

    m_db = db;
    int filterLen = m_db.length;
    int n,k,p;
    int decLen = (dataLen+filterLen-1)/2;
    double tmp=0;
    cout<<"decLen="<<decLen<<endl;

    for(n=0; n<decLen; n++)
    {
        cA[n] = 0;
        cD[n] = 0;
        for(k=0; k<filterLen; k++)
        {
            p = 2*n-k;  // 感謝網友onetrain指正，此處由p = 2*n-k+1改為p = 2*n-k，解決小波重構後邊界異常問題--2013.3.29 by hmm

            // 訊號邊沿對稱延拓
            if((p<0)&&(p>=-filterLen+1))
                tmp = sourceData[-p-1];
            else if((p>dataLen-1)&&(p<=dataLen+filterLen-2))
                tmp = sourceData[2*dataLen-p-1];
            else if((p>=0)&&(p<dataLen-1))
                tmp = sourceData[p];
            else
                tmp = 0;

            // 分解後的近似部分序列-低頻部分
            cA[n] += m_db.lowFilterDec[k]*tmp;

            // 分解後的細節部分序列-高頻部分
            cD[n] += m_db.highFilterDec[k]*tmp;
        }

    }

    return decLen;
}

2.5      小波重構

/**
 * @brief 小波變換之重構
 * @param cA 分解後的近似部分序列-低頻部分
 * @param cD 分解後的細節部分序列-高頻部分
 * @param cALength 輸入資料長度
 * @param db 過濾器型別
 * @param recData 重構後輸出的資料
 */
void  Wavelet::Idwt(double *cA, double *cD,  int cALength, Filter db, double *recData)
{
    if((NULL == cA)||(NULL == cD)||(NULL == recData))
        return;

    m_db = db;
    int filterLen = m_db.length;

    int n,k,p;
    int recLen = 2*cALength-filterLen+1;
    cout<<"recLen="<<recLen<<endl;

    for(n=0; n<recLen; n++)
    {
        recData[n] = 0;
        for(k=0; k<cALength; k++)
        {
            p = n-2*k+filterLen-1;

            // 訊號重構
            if((p>=0)&&(p<filterLen))
            {
                recData[n] += m_db.lowFilterRec[p]*cA[k] + m_db.highFilterRec[p]*cD[k];
                //cout<<"recData["<<n<<"]="<<recData[n]<<endl;
            }

        }
    }
}

2.6      c++實現結果

3 小波變換實現(MATLAB)

使用matlab小波工具箱實現db4的情況如下。

1、MatlabDB4.m檔案內容。

%載入txt資料示例
s=importdata('data2.txt'); %load data2.txt;
subplot(521);plot(s); %畫出原始訊號的波形圖
title('原始訊號');

[cA,cD]=dwt(s,'db4'); %採用db4小波並對訊號進行一維離散小波分解。
y=idwt(cA,cD,'db4'); %一維離散小波反變換
subplot(522);
plot(cA); %畫出波形圖
title('MATLAB低頻部分dwt-cA');

subplot(523);
plot(cD); %畫出波形圖
title('MATLAB高頻部分dwt-cD');

subplot(524);
plot(y); %畫出波形圖
title('MATLAB重構idwt');

2、波形圖如下。

4 小結

        在此，採用C++實現了dbN系列小波的單層一維離散小波變換，通過對比db4小波變換髮現（筆者也同樣對比驗證了db1-db3），C++實現效果和Matlab效果完全一樣。通過這一過程，可以很好地幫助理解小波變換的Mallat演算法原理，並將C++程式碼快速應用到工程實踐，同時對MATLAB小波工具箱的實現細節也有更深入的理解。相關文獻表明，C++程式碼實現的DWT演算法比Matlab小波工具箱dwt方法效率更高。

注：對小波變換見識還不深，有問題者，歡迎提出討論交流，對原始碼細節感興趣的tx，可以到如下連結下載。