快速小波變換的定點(diǎn)DSP實(shí)現(xiàn)

1 引言
    小波變換是近年來發(fā)展起來的一種數(shù)學(xué)理論和方法。作為一種新興的理論，小波分析是數(shù)學(xué)發(fā)展史上的重要成果，對工程應(yīng)用產(chǎn)生了深遠(yuǎn)的影響。廣泛應(yīng)用于語音信號處理、圖像信號處理、信號檢測、語音與圖像編碼、多尺度邊緣提取與重建等領(lǐng)域。近年來，在電力系統(tǒng)中也開始應(yīng)用小波分析進(jìn)行故障檢測及故障定位，并取得了有效的成果。
    計算機(jī)只能處理數(shù)字信號，所以在實(shí)際信號處理中，常采用離散形式的小波變換(Discrete Wavelet Transform，DWT)。由于小波變換算法的復(fù)雜性，盡管當(dāng)今處理器芯片運(yùn)算速度得到了大幅度的提高，仍然在實(shí)時性上不能滿足要求。為了簡化計算過程，人們發(fā)展了一些快速算法，如Mallat塔式算法，及利用調(diào)頻Z變換(chirped Z Transform，CZT)，梅林變換(Mellia Transform)進(jìn)行快速計算等算法。其中，尤其以Mallat塔式算法在實(shí)際應(yīng)用比較廣泛。
    在數(shù)字信號處理領(lǐng)域，通常使用專用的數(shù)字信號處理器芯片(DSP)以完成特定的運(yùn)算要求。美國TI公司是全球最大的DSP供應(yīng)商，其生產(chǎn)的TMS320C2xx系列16位定點(diǎn)DSP芯片具有高性能、低價格等特點(diǎn)，具有廣泛的應(yīng)用領(lǐng)域。本文中用該系列DSP芯片實(shí)現(xiàn)小波變換的快速算法。
    本文將小波變換快速算法用DSP加以實(shí)現(xiàn)，既可利用小波變換實(shí)現(xiàn)應(yīng)用要求，又可降低成本，增強(qiáng)市場競爭力。尤其在當(dāng)今，隨著電力系統(tǒng)的不斷發(fā)展，及用戶對電能質(zhì)量的要求越來越高，對電力系統(tǒng)運(yùn)行監(jiān)控及保護(hù)的采樣點(diǎn)數(shù)越來越多的情況下，利用此方法可以解決運(yùn)算量大、運(yùn)算精度高的問題。
2 小波變換及算法
2.1 小波變換
    小波函數(shù)的確切定義為：
    設(shè)ψ(t)為—平方可積函數(shù)，若其傅里葉變換ψ(ω)滿足條件

    則稱ψ(t)為一個基本小波或小波母函數(shù)。將其進(jìn)行伸縮和平移后，得到小波基函數(shù)

    所謂小波變換就是把信號在上述小波基下進(jìn)行展開。當(dāng)然，此變換必須存在逆變換，否則，不能恢復(fù)原信號，該變換就沒有什么意義了。
2.2 多分辨率分析
    多分辨率分析在正交小波變換理論中具有非常重要的地位，在多分辨率分析理論產(chǎn)生之前，人們構(gòu)造正交小波基函數(shù)要憑借技巧，具有一定的難度。自從有了多分辨率分析理論，這項工作變得容易的多。當(dāng)然，要尋找合適的基函數(shù)還是需要一定的經(jīng)驗(yàn)的。當(dāng)找到了合適的濾波器系數(shù)后，就可以利用Mallat給出的快速小波算法來計算小波變換了。
    通俗的講，多分辨率分析就是把空間V0上的函數(shù)f(t)分解為細(xì)節(jié)部分W1(小波空間)和大尺度逼近部分V1(尺度空間)，然后將大尺度逼近部分V1進(jìn)一步分解，如此重復(fù)就可得到任意尺度(或分辨率)上的逼近部分和細(xì)節(jié)部分。
2.3 濾波器系數(shù)
    根據(jù)多分辨率分析理論，如果φ(t)，ψ(t)分別為尺度空間V0及小波空間W0的一個標(biāo)準(zhǔn)正交基函數(shù)，則在任意相鄰尺度j，j-1之間，都有二尺度空間基函數(shù)關(guān)系

其中的h(n)，g(n)即為濾波器系數(shù)，由尺度函數(shù)φ(t)和小波系數(shù)ψ(t)決定。
2.4?Mallat塔式算法
    當(dāng)有了一組小波基函數(shù)后，剩下的事就是計算分解了，即把信號用小波基函數(shù)表示出來，從而關(guān)鍵問題是求出表示式中的系數(shù)。根據(jù)多分辨率分析，將信號f(t)□Vj-1分解一次(即分別投影到Vj、Ｗj空間)，此時cj,k和dj,k為j尺度上的展開系數(shù)，經(jīng)過不算復(fù)雜的推導(dǎo)，可得

[!--empirenews.page--]

其中cj,k和dj,k分別稱為j尺度空間的剩余系數(shù)和小波系數(shù)，上式說明它們可由j-1尺度空間的剩余系數(shù)cj-1，k經(jīng)濾波器系數(shù)進(jìn)行加權(quán)求和得到。實(shí)際中的濾波器h，g的長度都是有限長的或近似有限長的，因此分解運(yùn)算非常簡單。將cj,k進(jìn)一步分解下去，可分別得到Vj+1、Wj+1空間的剩余系數(shù)Cj+1，k和小波系數(shù)dj+1,k

    從而得到任意尺度空間上的分解。分解過程如圖所示

    在上述算法中必須要有一個初始輸入序列Cj-1，k，分解才能順利進(jìn)行，這是一個問題。在大多數(shù)應(yīng)用中，為了簡便，常用輸入信號的采樣序列來近似作為C0，k。在一些文獻(xiàn)里也給出了其它幾種確定C0，k的方法。
3?算法在DSP上的實(shí)現(xiàn)
    假設(shè)輸入信號x(t)，采樣頻率N(=2n)，得采樣序列x(k)，k=0，…，N-1，作為初始輸入序列C0，k。濾波器系數(shù)h(m),g(m)，m=0，…，L-1。為了應(yīng)用簡便，(1)、(2)式可變?yōu)?/p>

    從第一個尺度j=1開始，求出濾波器系數(shù)與剩余系數(shù)的加權(quán)和，分別得到cj,k與dj,k，并且它們的長度均為N／2。依次求出j=2，3，…各尺度值，cj,k與dj,k的長度也將變?yōu)镹／4,N／8,…。應(yīng)注意到，濾波器系數(shù)序列與輸入信號序列相乘時，各個系數(shù)依次相乘然后累加即為cj,k或dj,k值。計算完一個后，要將濾波器系數(shù)序列向后移兩個位置，再與輸入信號相乘。最后，只剩下兩個值時，再從第一個位置繼續(xù)，從而構(gòu)成一個圓周形式，得到最后一個cj,k與dj,k。以j=1為例來說明這一點(diǎn)。

    同樣，可得d1,N/2-1。?
    在TMS320C2xx系列定點(diǎn)DSP中，并未提供直接實(shí)現(xiàn)上述算法的尋址方式，可以利用循環(huán)指令實(shí)現(xiàn)。其中，要用到的兩條重要指令是MAC(乘累加指令)和RPT(重復(fù)指令)。MAC指令是DSP指令中最有特色的指令之一，當(dāng)RPT流水線啟動后，通過MAC指令可以在單指令周期內(nèi)實(shí)現(xiàn)乘加操作。算法的關(guān)鍵是一個卷積計算，其過程用如下幾條語句即可實(shí)現(xiàn)。假設(shè)程序存儲器地址0xFF00h開始，存放了小波濾波器系數(shù)h(k)，k=0，1，2，…,L-1的值，從0xFF80h開始，存放了小波濾波器系數(shù)g(k)，k=0，1，2，…，L-1的值。數(shù)據(jù)寄存器地?址0x1000h(用cc表示)開始，存放輸入信號。計算c1，k程序如下
        RPTL-1
        MAC0FF00h，cc+2*k
    下面簡要介紹一下單個尺度上的計算過程。仍假設(shè)輸入信號為N點(diǎn)采樣值，小波濾波器長度為L。
    由于有上述圓周形式的循環(huán)算法，直接計算所有值會有很多不便，因此將整個過程分為兩部分。第一部分，保存c0，k(k=N-2,N-1，0，1，2，…，L-3)到一連續(xù)存儲空間，然后計算k＜N／2-L時的d1，k的值，并保存到臨時存儲空間datad(需要N個單元)中。第二部分計算k=N／2-L,…,N／2-1時的d1，k的值，保存到datad中(從第N／2-L單元開始)。與前面相同方法計算c１，0，c1，1，…，的值，保存時覆蓋掉原來的c0,0，c0,1，…，即可。如果為了系數(shù)重構(gòu)方便，可以把datad中前N／2個單元內(nèi)容移到原來的c0,k(k=N／2,…，N-1)位置。如此方法便可得到分解系數(shù)。?
4 結(jié)論
    小波分析具有良好的時頻局部性，是當(dāng)今很受歡迎的分析方法之一，利用定點(diǎn)DSP進(jìn)行小波變換，滿足了實(shí)時性，具有良好的精度和低成本，是工控領(lǐng)域的理想選擇。