快速小波變換的定點DSP實現

小波變換具有良好的時——頻局部性，是分析奇異信號的重要方法。定點DSP在工程中的應用十分普遍，具有低成本，高性能的特點。利用DSP實現小波變換可以滿足工程是實時性的要求。文中簡要介紹了小波變換理論及算法，并結合TI公司的16位定點DSP說明算法的實現。
??? 關鍵詞：快速算法? 小波變換? DSP

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

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

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

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

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

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

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

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

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

參考文獻：

［１］彭玉華.《小波變換與工程應用》［Ｍ］．科學出版社
［２］程正興．《小波分析算法與應用》［Ｍ］．西安交大出版社
［３］TMS320C2xx User's Guide［Ｓ］．Texas　lnstruments
［４］楊福生．《小波變換的工程分析與應用》［Ｍ］．科學出版社
［５］李建平．《小波分析與信號處理》［Ｍ］．重慶出版社
［６］秦前清，楊宗凱．《實用小波分析》［Ｍ］．西安電子科技大學出版社