快速小波變換在DSP中的實(shí)現(xiàn)方法

1、快速小波變換在DSP中的實(shí)現(xiàn)方法2008-03-25 09:21姜新華范征宇（上海交通大學(xué)信息與控制工程系上海，200240）摘要小波分析是分析非穩(wěn)定信號的一種非常有效的方法。Mallat快速算法使得小波分析的廣泛應(yīng)用成為現(xiàn)實(shí)。在實(shí)時信號分析中，離散小波變換在DSP上的有效應(yīng)用受到了特別的關(guān)注。本文簡單介紹了FWT，詳細(xì)闡述了在DSP上FWT的周期性擴(kuò)展的實(shí)現(xiàn)。其中，特別介紹了DSP的循環(huán)尋址。最后，給出了針對TI公司的TMS320C3X系列DSP的相應(yīng)的匯編代碼。關(guān)鍵詞：小波變換；快速算法；數(shù)字信號處理器引言小波變換是從信號理論領(lǐng)域發(fā)展形成的。小波變換涵蓋許多領(lǐng)域，它在這些領(lǐng)域得到了（或?qū)⒌?/p>

2、到）廣泛的應(yīng)用。作為信號處理的經(jīng)典工具，傅里葉變換的缺點(diǎn)是它沒有時域局部化性質(zhì)。它僅僅給出信號的頻率成分表達(dá)，但一般不能從變換后的信號中得出頻率與時域局部有關(guān)的信息，其原因就在于分解函數(shù)的周期性（三角函數(shù)）。這些函數(shù)可在頻域定位，但在時域則不行。因此，傅立葉變換是表示穩(wěn)定信號的理想工具。小波變換的分解函數(shù)，即所謂的“小波”在時域和頻域都具有局部性。因此，小波變換是表示包含不連續(xù)或不穩(wěn)定成分的信號的理想工具。在實(shí)際的信號處理應(yīng)用中，常采用的是小波變換的離散形式，也即所謂的離散小波變換（DWT）。1988年，Mallat在Burt和Adelson的圖像分解和重構(gòu)的塔形算法啟示下，基于多分辨分析框架

3、，給出了信號與圖像（或函數(shù)）分解為不同頻率通道的算法及其重構(gòu)算法Mallat算法1，大大提高了小波變換的計(jì)算速度，其重要性相當(dāng)于傅里葉變換中的FFT算法。從而使小波變換具有了明顯的工程應(yīng)用價值。通用的微處理器在運(yùn)算速度上難以適應(yīng)信號實(shí)時處理的要求。DSP處理器就是為了適應(yīng)這種需求而開發(fā)的。DSP處理器中集成有高速的乘法器硬件，能快速地進(jìn)行大量數(shù)據(jù)的乘法和加法運(yùn)算。小波變換和DSP處理器的結(jié)合使得小波變換具有很強(qiáng)的可實(shí)現(xiàn)性。1快速小波變換（FWT）基于多分辨分析的理論，Mallat給出了快速小波變換（FWT）的算法Mallat算法。對于一個多分辨分析VjjZ，和分別為相應(yīng)的尺度函數(shù)和小波函數(shù)，函

4、數(shù)f屬于由尺度函數(shù)生成的MRA的基本子空間V0。它可以表述成5對上述表達(dá)式的幾點(diǎn)說明如下：（1）hk是多分辨分析的尺度系數(shù)，gk是多分辨分析的小波系數(shù)，兩者的關(guān)系是gn（1）nh2k1n，kZ；（2）dmn和cmn分別是函數(shù)f在第m層尺度上的小波分解系數(shù)和尺度逼近系數(shù)；（3）實(shí)際應(yīng)用中，函數(shù)f不是連續(xù)函數(shù)，而是所要分析的信號的離散采樣序列fk。當(dāng)采樣率相當(dāng)高時，信號的采樣就非常近似于展開系數(shù)c0。這樣，計(jì)算離散小波變換可以描述成遞歸計(jì)算式（1）和式（2）兩個和式。因此，分析算法有如下形式：輸入c0（輸入序列）M（分解層次）計(jì)算for m1，M2有限序列的FWT以上算法推導(dǎo)時，假定序列是無限的。

5、然而，實(shí)際應(yīng)用時序列總是有限的。為了滿足多分辨分析的要求，可以用合適的數(shù)來擴(kuò)展有限序列，常采用的有零插值法和周期性擴(kuò)展法。零插值法是將序列支撐以外的元素都設(shè)為零。如果序列不是以零開頭或結(jié)尾，這就在序列的邊界處人為地產(chǎn)生了不連續(xù)性，從而使得在dm相應(yīng)位置出現(xiàn)較大的值。另外，零插值法中，輸入序列與系數(shù)序列的支撐的長度不是相互獨(dú)立的，這增加了應(yīng)用時的計(jì)算量。周期性擴(kuò)展法是把序列當(dāng)作一個周期信號進(jìn)行擴(kuò)展，序列支撐的長度就是周期的長度。周期性擴(kuò)展克服了零插值法的不足，而且這種周期性能沿著尺度傳遞。另外，可以利用DSP處理器的循環(huán)尋址功能實(shí)現(xiàn)序列的周期性擴(kuò)展，從而更加提高了FWT的計(jì)算速度。因此，本文采用

6、周期性擴(kuò)展法在DSP處理器上實(shí)現(xiàn)FWT。對式（1，2）變換下標(biāo)可得 輸入序列定義為x（n），尺度展開系數(shù)為z（k），而系數(shù)為f（l）（包括尺度系數(shù)h（l）和小波系數(shù)g（l）。當(dāng)算法應(yīng)用于DSPs上時，這些序列的下標(biāo)必須從0開始。系數(shù)序列的下標(biāo)的范圍為：0，lupper，lupper表示下標(biāo)的最大值 如果輸入序列x（n）的長度N2u，uZ，則z（k）的長度沿尺度不斷地被二等分；因此，下標(biāo)k的范圍是：0，N/2-1。     這樣，就可以通過Mallat算法計(jì)算某個尺度值：將系數(shù)序列放至輸入序列的開頭，各個系數(shù)與相應(yīng)的輸入值相乘再累加。最后所得就是尺度展開序列上的第一

7、個值。然后，系數(shù)序列向前跳過輸入序列的兩個位置（基于2的向下抽樣），再計(jì)算各個系數(shù)與相應(yīng)的輸入值的乘積的和，這樣就得到第二個值。依次類推，直到系數(shù)序列跳出輸入序列的范圍。此時為了使每個系數(shù)都與輸入序列中的值相乘，把輸入序列的開始部分的一些值加到序列的結(jié)尾。整個計(jì)算完成后，得到的序列就是下一個尺度計(jì)算的輸入序列。算法的原理如圖1所示。圖中長度為4的系數(shù)序列f（l）（l0，3）平移過長度為N輸入序列x（n）（n0，N1；N2u，uZ）。當(dāng)k0，計(jì)算尺度“m1”上的第一個值 xm1（0）xm1（0）f（0）·xm（0）f（1）·xm（1）   

8、60;       f（2）·xm（2）f（3）·xm（3）系數(shù)序列然后平移至“k1”，計(jì)算第二個值xm1（1）xm1（1）f（0）·xm（2）f（1）·xm（3）       f（2）·xm（4）f（3）·xm（5）如此反復(fù)，直到移至“kN21”。這時，只有頭兩個系數(shù)可以與輸入序列相乘。為此，將xm（0）和xm（1）加到序列的結(jié)尾處。這樣，尺度“m1”的最后一個值xm1（N21）就可以由下式計(jì)算xm+1（N/21）

9、=f（0）·xm(N-2)+              f（1）·xm（N1）              f（2）·xm（0）          f（3）·xm（1）然后，把xm1（n）（n0，N21）作為輸入

10、值計(jì)算下一個尺度“m2”。3循環(huán)尋址大多數(shù)的DSPs的指令系統(tǒng)中有“循環(huán)尋址”（Circular addressing）的指令。周期性擴(kuò)展可以通過循環(huán)尋址很容易地實(shí)現(xiàn)。循環(huán)尋址是間接編址的特殊情況，它用到一個循環(huán)緩存區(qū)。這種循環(huán)緩存區(qū)通常用于卷積、相關(guān)、數(shù)字濾波等算法中。它可以看成一個滑動的窗口，窗口內(nèi)是最新的待處理的數(shù)據(jù)，如果有新的數(shù)據(jù)輸入，就覆蓋掉原有的數(shù)據(jù)，從而形成循環(huán)緩存。    數(shù)據(jù)x（n）的開頭由一個指針定位，x（n）之前的數(shù)據(jù)從這個指針開始按順時針方向連續(xù)寫入。當(dāng)?shù)玫揭粋€新數(shù)據(jù)后，其將被放置到x（n）。在有關(guān)的算法計(jì)算中（比如FIR濾波器的應(yīng)用），通

11、過指針平移來連續(xù)定位每一個數(shù)據(jù)，直到最后一個輸出數(shù)據(jù)的得出。這時，指針又回到了x（n）。然后，指針再逆時針方向移動一個位置，指向x（nN）。因?yàn)檫@個值是“最舊的”，不再有用，下一個數(shù)據(jù)就寫到x（nN）的位置。周而復(fù)始，直至整個計(jì)算完成。在文2中討論了循環(huán)緩存區(qū)在數(shù)字濾波器中的應(yīng)用。所有的微處理器（包括DSP處理器）都提供線性存儲空間，而循環(huán)尋址使得線性存儲空間可以用作循環(huán)緩存空間。本文中，開始地址和循環(huán)緩存區(qū)的大小必須預(yù)先設(shè)定。指針指向開始地址，并且每次循環(huán)執(zhí)行完成后隨即增加“1”以指向下一個較高的存儲位置，直到指針指向緩存尾部。循環(huán)緩存的這種結(jié)構(gòu)很好地滿足了輸入信號的周期性擴(kuò)展的要求，進(jìn)而可

12、用于實(shí)現(xiàn)FWT算法。上面提到為了在序列的邊界上，每個系數(shù)都有輸入值與之相乘，將序列開始部分的一些值添加到序列的末尾。如果輸入序列應(yīng)用循環(huán)緩存區(qū)，那么在到達(dá)數(shù)據(jù)緩沖末尾時，指針又一次跳到開始部分。也就是說，添加過程是自動完成的。循環(huán)尋址的原理在文3中闡述。 4小波變換在DSP中的應(yīng)用對DSP編寫代碼主要有兩種方法：（1）用高級語言編寫，常用的如C語言；（2）用匯編語言編寫。兩者有各自的優(yōu)缺點(diǎn)，見表1。大多數(shù)的DSP程序在一小段的代碼上往往占去極大多數(shù)的執(zhí)行時間。根據(jù)經(jīng)驗(yàn)法則，90的執(zhí)行時間將被10的代碼占據(jù)。為了充分利用高級語言和匯編語言的優(yōu)點(diǎn)，最好的方法是編寫一個高效的匯編語言的子程序，在C語

13、言程序中象普通函數(shù)一樣調(diào)用它4。本文以TMS320C3X系列為例來闡述小波變換在DSP中的應(yīng)用。Texas Instruments的TMS320C3X系列在相應(yīng)的用戶手冊上有詳細(xì)的說明。匯編語言子程序的主要部分是計(jì)算一個尺度上的尺度展開值，以下是針對TMS320C3X DSP的子程序的主要部分： LDIXARRAY，AR0LDI      ZARRAY，AR1LDI      WLETTAB，AR2LDI      0，IR0LDI &

14、#160;    DX，BKLDI      DZ，RCSUBI     1，RCRPTB     WTBLOCKLDF      00，R0LDF      00，R2PUSH     STPUSH RSPUSH     REPUSH  

15、   RCRPTS     LUPPERMPYF     AR0，AR2，R0ADDF   R0，R2ADDF     R0，R2POP      RCPOP      REPOP      RSPOP      STSTF &#

16、160;    R2，AR1ADDI     2，IR0LDI      XARRAY，AR0NOP      AR0（IR0）WTBLOCKLDIWLETTAB，AR2尺度計(jì)算是在兩個嵌套的循環(huán)中完成的。在內(nèi)循環(huán)中，每次通過系數(shù)與輸入值相乘累加求得一個尺度值（見式（6）。外循環(huán)保證每次內(nèi)循環(huán)執(zhí)行完畢后系數(shù)序列能平移到輸入序列的正確位置（基于2向下抽樣）。外循環(huán)運(yùn)行的次數(shù)對應(yīng)于該尺度下所要計(jì)算的尺度值的個數(shù)。比如輸入序列的長度

17、是N（偶數(shù)），那么尺度值個數(shù)就是N2。輸入值、尺度值和小波系數(shù)都由指針來定位。這些指針必須在程序的開頭初始化 LDIXARRAY，AR0LDIZARRAY，AR1LDIWLETTAB，AR2指針AR0指向輸入序列，AR1指向要計(jì)算的尺度值的數(shù)據(jù)緩存區(qū)，AR2指向所用的小波系數(shù)表。外循環(huán)通過循環(huán)WTBLOCK代碼塊來計(jì)算某個尺度所有的尺度值（RPTB指令，Repeat block of code）。決定循環(huán)執(zhí)行次數(shù)的循環(huán)次數(shù)寄存器RC必須賦予所要分析的尺度的長度（DZ），最小是1。LDIDZ，RCSUBI 1，RC內(nèi)循環(huán)通過指令RPTS（Repeat single instruction）來計(jì)算

18、其中一個尺度值。RPTSsrc指令循環(huán)執(zhí)行下一條指令src1次。在程序中，循環(huán)次數(shù)與所用小波系數(shù)序列長度相等（RPTSLUPPER）。尺度的任何一個值正是通過循環(huán)執(zhí)行下一條指令計(jì)算得出：把輸入序列與小波系數(shù)相乘，再累加，這些乘法和加法由指令的平行操作結(jié)構(gòu)來實(shí)現(xiàn)。在這種結(jié)構(gòu)中，兩個指令是同時執(zhí)行的。這種結(jié)構(gòu)具有高度的并發(fā)性。平行操作的標(biāo)識符是“”，放在作平行操作的兩條指令的第二條指令前。關(guān)于平行操作指令的詳細(xì)說明可參考文3。相應(yīng)的代碼是RPTSLUPPERMPYF    AR0，AR2，R0ADDF   R0，R2 ADDF  

19、;   R0，R2在執(zhí)行指令RPTS前，寄存器R0和R2必須置零（LDF00，R0和LDF00，R2）。指針AR0指向的輸入值與指針AR2指向的小波系數(shù)相乘，結(jié)果存入寄存器R0。此外，R0中的值加到寄存器R2中。緊接著，指針AR0和AR2各增加1以指向下一個更高的存儲位置。指針AR0是按循環(huán)尋址遞增的（以“”標(biāo)識），這樣，當(dāng)?shù)竭_(dá)數(shù)據(jù)緩存區(qū)的末尾時，AR0在下一次操作中將重新指向緩存的開頭。為了保證循環(huán)緩存區(qū)的正常工作，其大小必須在第一次執(zhí)行外循環(huán)前定義。為此，寄存器BK被賦予輸入序列的長度值（DZ）LDIDX，BK如果平行操作指令循環(huán)lupper次，寄存器R0中值反復(fù)地被加

20、到寄存器R2中，最終在R2中得到該尺度下的一個尺度值。TMS320C3X執(zhí)行指令采用流水線技術(shù)3。一條指令的執(zhí)行經(jīng)過四個獨(dú)立的步驟：取指、譯碼、訪問數(shù)據(jù)、執(zhí)行。流水線技術(shù)是將各指令的執(zhí)行時間重疊起來，使得每一條指令的最終執(zhí)行時間是在單個指令周期內(nèi)完成的，大大提高了速度。所有的寄存器都在執(zhí)行段開始時讀出，結(jié)束時寫入。也就是說在上述例程第一次執(zhí)行平行操作指令時，在把輸入值與小波系數(shù)的乘積寫入寄存器R0前，R0中的值（0）先被讀出并加到R2中。第一次執(zhí)行完成后，R2中的值仍是0，并不是第一次乘法操作的結(jié)果。因此，RPTS指令執(zhí)行完成后，R0中的值必須再加入R2中一次。第二條ADDF指令累加了最后一次

21、的乘積。緊接著，計(jì)算得出的尺度值（R2中的內(nèi)容）存儲到由指針AR1指向的存儲位置中。AR1指向存儲尺度值的存儲空間。AR1在執(zhí)行下面的指令后增加1STFR2，AR1在循環(huán)塊WTBLOCK的結(jié)尾，指針AR0重新指向輸入序列的開頭，其內(nèi)容增加IR0中的值LDIXARRAY，AR0NOPAR0（IR0）寄存器IR0在第一次執(zhí)行循環(huán)模塊WTBLOCK前初始化為0，在模塊的結(jié)尾其值增加2。第一次執(zhí)行時，AR0指向X（0），第二次指向X（2），第三次指向X（4），等。這樣就實(shí)現(xiàn)了Mallat算法中的基于2的向下抽樣。 指向小波系數(shù)的指針AR2在每次執(zhí)行新的循環(huán)前必須指向小波的第一個系數(shù)LDIWLETTAB，AR2當(dāng)執(zhí)行嵌套循環(huán)操作時（在指令RPTB中執(zhí)行指令RPTS），必須注意這樣一個事實(shí)：指