一維小波變換的C++實現

　　將小波展開系數當成離散信號，尺度函數和小波函數的MRA方程系數看成數字濾波器組，根據Mallat快速算法的原理，小波變換對數據的處理方法可簡化成對信號逐級采樣和濾波的過程。

圖1 小波變換的濾波器實現

(a)分解算法 (b)重構算法

　　一層小波分解算法流程如圖2所示，信號將先經過小波分解低通濾波器和高通濾波器，隨后被降采樣，實現數據重構。而濾波算法可簡化為待處理信號與濾波器數組卷積的過程，為了保證卷積前和卷積后數組的長度相同，結合小波變換中數組延拓的思想，在實際編程過程中，可以將超過信號長度的那段數據以前端對齊的方式與前面一段數據相加。將卷積后的數組每2個點采樣一次，即可獲得小波分解后的尺度系數和小波系數。

圖2 一層小波分解算法

（X：待分解數組；H，G：小波分解濾波器；C，D：小波重構后數組）

　　小波重構算法是小波分解算法的逆運算，其流程為升采樣和濾波，最后數據相加實現重構。小波重構算法中濾波可視為系數與小波重構濾波器的卷積，與小波正變換類似，在重構算法中，需要將卷積后的數組末位對齊相加，獲得與原數組長度相同的卷積結果。將小波系數和尺度系數以2為步長進行升采樣，將獲得的新數組分別經過小波重構低通濾波器和高通濾波器，再將濾波后的兩組數據相加，即實現了一層小波重構。

#define LENGTH 512
#define LEVEL 4
#define L_core 6

static void Covlution(double data[], double core[], double cov[], int LEN)
{
        double temp[LENGTH + L_core - 1] = {0};
        int i = 0;
        int j = 0;

    for(i = 0; i < LEN; i++)
        {
                for(j = 0; j  < L_core; j++)
                {
            temp[i + j] += data[i] * core[j];
        }
    }

        for(i = 0; i < LEN; i++)
        {
                if(i < L_core - 1)
                        cov[i] = temp[i] + temp[LEN + i];
                else
                        cov[i] = temp[i];
        }

}

static void Covlution2(double data[], double core[], double cov[], int LEN)
{
        double temp[LENGTH + L_core - 1] = {0};
        int i = 0;
        int j = 0;

    for(i = 0; i < LEN; i++)
        {
                for(j = 0; j  < L_core; j++)
                {
            temp[i + j] += data[i] * core[j];
        }
    }

        for(i = 0; i < LEN; i++)
        {
                if(i < L_core - 1)
                        cov[i + LEN - L_core + 1] = temp[i] + temp[LEN + i];
                else
                        cov[i - L_core + 1] = temp[i];
        }

}

static void DWT1D(double input[], double output[], double LF[], double HF[], int l)
{
    int i = 0;
        double temp[LENGTH] = {0};
        int LEN = LENGTH / pow(2, l - 1);

    Covlution(input, LF, temp, LEN);
    for(i = 1; i < LEN; i += 2)
    {
        output[i/2] = temp[i];
    }

    Covlution(input, HF, temp, LEN);
    for(i = 1; i < LEN; i += 2)
    {
        output[LEN/2 + i/2] = temp[i];
    }
}

static void DWT(double input[], double output[], double LF[], double HF[], int len[])
{
        int i;
        int j;

        len[0] = len[1] = LENGTH / pow(2, LEVEL);
        for(i = 2; i <= LEVEL; i++) len[i] = len[i - 1] * 2;

        DWT1D(input, output, LF, HF, 1);
        for(i = 2; i <= LEVEL; i++)
    {
                for(j = 0; j < len[LEVEL + 2 - i]; j++) input[j] = output[j];
        DWT1D(input, output, LF, HF, i);
    }
}

static void IDWT1D(double input[], double output[], double LF[], double HF[], int l, int flag)
{
    int i = 0;
        double temp[LENGTH] = {0};
        int LEN = l * 2;

        if(flag) Covlution2(input, HF, temp, LEN);
        else Covlution2(input, LF, temp, LEN);

        for(i = 0; i < LEN; i++)
    {
        output[i] = temp[i];
    }
}

static void IDWT(double input[], double output[], double LF[], double HF[], int len[], int level)
{
        int i;
        int j;
        for(j = 0; j < len[LEVEL + 1 - level]; j++)
        {
                output[2 * j] = 0;
                output[2 * j + 1] = input[j];
        }
        for(j = 0; j < 2 * len[LEVEL + 1 - level]; j++)
        {
                input[j] = output[j];
        }
    IDWT1D(input, output, LF, HF, len[LEVEL + 1 - level], 1);

        for(i = level - 1; i > 0; i--)
    {
                for(j = 0; j < len[LEVEL + 1 - i]; j++)
                {
                        input[2 * j] = 0;
                        input[2 * j + 1] = output[j];
                }
        IDWT1D(input, output, LF, HF, len[LEVEL + 1 - i], 0);
    }
}

 用C++算法實現的小波變換結果與MATLAB實現的小波變換結果對比（心電信號，db5小波，5層分解）。計算的結果幾乎相差無異，大多數位於中間的數據計算結果完全相同，兩端的數據則存在一定誤差，這可能與對信號采取的延拓方式不同有關。

 

 

 

(a)小波變換結果（紅線：MATLAB計算結果；藍線：C++算法計算結果）

(b)小波變換誤差（兩個計算結果相減）

參考文獻：網上各類資料，時間有點久，當時忘記記下了。