離散小波變換

1、1.離散小波變換長期以來，離散小波變換(Discrete Wavelet Transform)在數(shù)字信號(hào)處理、石油勘探、地 震預(yù)報(bào)、醫(yī)學(xué)斷層診斷、編碼理論、量子物理及概率論等領(lǐng)域中都得到了廣泛的應(yīng)用。各種 快速傅氏變換(FFT)和離散小波變換(DWT)算法不斷出現(xiàn)，成為數(shù)值代數(shù)方面最活躍的一個(gè) 研究領(lǐng)域，而其意義遠(yuǎn)遠(yuǎn)超過了算法研究的范圍，進(jìn)而為諸多科技領(lǐng)域的研究打開了一個(gè)嶄新的局面。本章分別對(duì)FFT和DWT的基本算法作了簡(jiǎn)單介紹，若需在此方面做進(jìn)一步研究， 可參考文獻(xiàn)2。1.1離散小波變換DWT1.1.1離散小波變換DWT及其串行算法先對(duì)一維小波變換作一簡(jiǎn)單介紹。設(shè)f(x)為一維輸入信號(hào)，記-

2、jk (x) =2一"2 "2X _ k)，'-：jk(x) T-i/G&Tx-k)，這里(x)與(x)分別稱為定標(biāo)函數(shù)與子波函數(shù)，jk(X)與p jk (x)為二個(gè)正交基函數(shù)的集合。記Pof=f,在第j級(jí)上的一維 離散小波變換DWT(DiscreteWavelet Transform)通過正交投影 Pjf與Qf將Pj-f分解為：Pj 丄f =Pj f Qj fck k 二 dkjjkkk:P： .j Pif -其中：ck= S h(n)c2“，dkj= S g(n曲竊(j=1,2,L,k =0,1,N/2j1)，這里,h(n)n =0n =0與g(n)分別

3、為低通與高通權(quán)系數(shù)，它們由基函數(shù) jk(x)與f jk(x)來確定，p為權(quán)系數(shù)的長度。cn0為信號(hào)的輸入數(shù)據(jù)，n為輸入信號(hào)的長度，l為所需的級(jí)數(shù)。由上式可見，每級(jí)一維DWT與一維卷積計(jì)算很相似。所不同的是：在DWT中，輸出數(shù)據(jù)下標(biāo)增加1時(shí)，權(quán)系數(shù)在輸入數(shù)據(jù)的對(duì)應(yīng)點(diǎn)下標(biāo)增加2,這稱為“間隔取樣”。算法22.3 維離散小波變換串行算法輸入：co=do(co0, c10,cn-10)h=(h 0, h1,，hL-1) g=(g 0, g1,，g)輸出：cj , dij (i=0, 1,，N/2j-1, j > 0)Begin(1) j=0, n=N(2) While (n> 1) do(

4、2.1)for i=0 to n-1 docij+1 =0, d ij+1 =0(2.1.2)for k=0 to L-1 do電 i) mod n，d/ 1gkdj(k -2 i) mod nend forend for(2.2)j=j+1, n=n/2end whileEnd顯然，算法22.3的時(shí)間復(fù)雜度為O(N*L)。在實(shí)際應(yīng)用中，很多情況下采用緊支集小波(Compactly Supported Wavelets),這時(shí)相應(yīng)的尺度系數(shù)和小波系數(shù)都是有限長度的，不失一般性設(shè)尺度系數(shù)只有有限個(gè)非零值： h1, ,hN，N為偶數(shù)，同樣取小波使其只有有限個(gè)非零值：g!, ,gN。為簡(jiǎn)單起見，設(shè)尺

5、度系數(shù)與小波函數(shù)都是實(shí)數(shù)。對(duì)有限長度的輸入數(shù)據(jù)序列：C0二xn,n =1,2，M (其余點(diǎn)的值都看成0)，它的離散小波變換為：Ck = ' Cn hn _2kn Z -dk + =Z Crj gn/kn.Zj =0,1，J -1其中J為實(shí)際中要求分解的步數(shù)，最多不超過logzM，其逆變換為C*二 ' Ck hn _2k V Ck hn-2kk目k目j = J, ,1注意到尺度系數(shù)和輸入系列都是有限長度的序列，上述和實(shí)際上都只有有限項(xiàng)。若完全按照上述公式計(jì)算，在經(jīng)過J步分解后，所得到的 J+1個(gè)序列dkj, j = 0,1，J -1和ck的非零項(xiàng)的個(gè)數(shù)之和一般要大于M,究竟這個(gè)項(xiàng)

6、目增加到了多少？下面來分析一下上述計(jì)算過程。j=0時(shí)計(jì)算過程為1 MCk 二 xnhn_2kn=1Mdk 二' xng n_2kn土不難看出，C：的非零值范圍為：k =N+11,0件11,即有 k2|2 Ik =-； T + |M =J T 個(gè)非零值。dk的非零值范圍相同。繼續(xù)往下分解時(shí)，非零項(xiàng)出 現(xiàn)的規(guī)律相似。分解多步后非零項(xiàng)的個(gè)數(shù)可能比輸入序列的長度增加較多。例如，若輸入序列長度為100，N=4，則d：有51項(xiàng)非零，d2有27項(xiàng)非零，d3有15項(xiàng)非零，d：有9項(xiàng)非 零，dk有6項(xiàng)非零，d6有4項(xiàng)非零，c6有4項(xiàng)非零。這樣分解到 6步后得到的序列的非 零項(xiàng)個(gè)數(shù)的總和為116，超過了輸

7、入序列的長度。在數(shù)據(jù)壓縮等應(yīng)用中，希望總的長度基本 不增加，這樣可以提高壓縮比、減少存儲(chǔ)量并減少實(shí)現(xiàn)的難度?？梢圆捎蒙晕⒏淖冇?jì)算公式的方法，使輸出序列的非零項(xiàng)總和基本上和輸入序列的非零項(xiàng)數(shù)相等，并且可以完全重構(gòu)。這種方法也相當(dāng)于把輸入序列進(jìn)行延長（增加非零項(xiàng)），因而稱為延拓法。只需考慮一步分解的情形，下面考慮第一步分解（j=1）。將輸入序列作延拓，若 M為偶數(shù)，直接將其按 M為周期延拓；若 M為奇數(shù)，首先令xM t =0。然后按M + 1為周期延拓。作了這種延拓后再按前述公式計(jì)算，相應(yīng)的變換矩陣已不再是H和G,事實(shí)上這時(shí)的變換矩陣類似于循環(huán)矩陣。例如，當(dāng)M=8，N=4時(shí)矩陣H變?yōu)?h3h400

8、00h2hh2h3h4000000h1h2h3h4000000h1h2h3h4h3h40000h2當(dāng)M=7, N=4時(shí)矩陣H變?yōu)椋篽3h40000hihh2h3h400000hih2h3h400000hih2h3h3h40000hi從上述的矩陣表示可以看出， 兩種情況下的矩陣內(nèi)都有完全相同的行，這說明作了重復(fù)計(jì)算，因而從矩陣中去掉重復(fù)的那一行不會(huì)減少任何信息量，也就是說，這時(shí)我們可以對(duì)矩陣進(jìn)行截短（即去掉一行），使得所得計(jì)算結(jié)果仍然可以完全恢復(fù)原輸入信號(hào)。當(dāng)M=8，N=4時(shí)截短后的矩陣為：-h3h40000hih2hh2h3h40000H00hh2h3h4000000hih2h3h4當(dāng)M=7，

9、N=4時(shí)截短后的矩陣為h3h40000hhh2h3h4000H=00h1h2h3h400000hih2h3這時(shí)的矩陣都只有葺行。分解過程成為:c1 = HC0D1 二GC0向量C1和D1都只有 M 個(gè)元素。重構(gòu)過程為：|2C -H *C1 G* D1可以完全重構(gòu)。矩陣H , G有等式H*H + G*G = I一般情況下，按上述方式保留矩陣的M I行，可以完全恢復(fù)原信號(hào)。|2 I這種方法的優(yōu)點(diǎn)是最后的序列的非0元素的個(gè)數(shù)基本上和輸入序列的非0元素個(gè)數(shù)相同，特別是若輸入序列長度為2的幕，則完全相同，而且可以完全重構(gòu)輸入信號(hào)。其代價(jià)是得到的變換系數(shù)Dj中的一些元素已不再是輸入序列的離散小波變換系數(shù)，對(duì)某些應(yīng)用可能是不適合的，但在數(shù)據(jù)壓縮等應(yīng)用領(lǐng)域，這種方法是可行的。Begin對(duì)所有處理器my_rank(my_rank=0, , , p-1)同時(shí)執(zhí)行如下的算法：(1) j=0；(2) while (j<r) do(2.1) 將數(shù)據(jù)c*(n =0,1,，N/2j _1)按塊分配給p臺(tái)處理器(2.2) 將處理器i+1中前L-1個(gè)數(shù)據(jù)發(fā)送給處理器i(2.3) 處理器i負(fù)責(zé)cn*和dn出(n =iN帀,，(i +1)N帀一1)的計(jì)算(2.4) j=j+1end whileEnd這里每一步分解后數(shù)據(jù)c 1和dnj 1已經(jīng)是按塊存儲(chǔ)在P臺(tái)處理器上，因此算法第一步中的數(shù)據(jù)分配除了