Daubechies小波系數(shù)新解法

張旭俊，張宇

(國網(wǎng)江西省電力有限公司電力科學(xué)研究院， 南昌 330096)

0 引 言

在大數(shù)據(jù)時代，小波是大量信號處理的有效時頻分析工具[1-2]，在電氣工程領(lǐng)域的諧波檢測[3-4]、設(shè)備故障診斷[5]、故障選線[6]和繼電保護[7]等方面有很高的實用價值。在經(jīng)典小波中Daubechies小波(以下簡稱DB小波)是應(yīng)用最廣的，它充分體現(xiàn)了Daubechies思想的深邃和數(shù)學(xué)才華。但是經(jīng)典的DB小波推演過程十分深奧，需要較多的數(shù)學(xué)知識，不便普及推廣，文章的目的就是用普通的高等數(shù)學(xué)給出各種DB小波系數(shù)的計算方法，透徹的理解才能更好地靈活應(yīng)用。為了深入認(rèn)識，還把離散小波正交的概念引向正交的小波矩陣，即完全在時域中分析小波的分解與重構(gòu)，可避免傅里葉積分變換的困擾。

1 雙尺度分解

經(jīng)典小波理論認(rèn)為，尺度函數(shù)φ(t)和小波函數(shù)ψ(t)只在某一區(qū)間內(nèi)有值，其余區(qū)間都是零，它和傅里葉級數(shù)分解不同，以分析128點采樣數(shù)據(jù)為例，小波分解是用分段位移的逐次雙尺度分解來分析函數(shù)的構(gòu)成。對比傅里葉級數(shù)，小波是在選定的區(qū)間，用尺度函數(shù)對應(yīng)濾出 0～32 Hz低次諧波部分，用小波函數(shù)濾出33 Hz ～64 Hz高次諧波部分。接著又對0～32 Hz低頻諧波部分，進行雙尺度分解為0～16 Hz更低的低頻部分和 17 Hz～32 Hz的次高頻部分。繼續(xù)雙尺度分解以達(dá)到對采樣函數(shù)的多尺度分解。

尺度函數(shù)Φ(t)或小波函數(shù)ψ(t)可以由高一倍頻率的尺度函數(shù)Φ(2t)及其位移Φ(2t-n)的線性組合構(gòu)成，即：

(1)

(2)

它們的傅里葉變化函數(shù)Φ(ω)和Ψ(ω)可推出如下：

(3)

同理有：

(4)

且有如下關(guān)系：

(5)

從物理意義上說，尺度函數(shù)Φ(ω/2)的頻率比Φ(ω)頻率高一倍，而H(ω/2)相當(dāng)于濾波函數(shù)。

令z=e-jω，濾波函數(shù)H(ω)可表達(dá)如下：

(6)

(7)

尺度函數(shù)和另一平移后的尺度函數(shù)應(yīng)當(dāng)是正交的，用如下公式表示：

ejωmdω=δ0,m

(8)

為了達(dá)到正交的目的，只當(dāng)m=0時，才會有δ0,0=1，這需要滿足式(9)的恒等式關(guān)系。其中應(yīng)用了雙尺度關(guān)系式(3)。

(9)

根據(jù)H(ω)=H(ω+2mπ)的周期性，在式(10)推導(dǎo)中，并將其中整數(shù)k依奇、偶數(shù)分割成2m和2m+1兩部分，推導(dǎo)得出：

(10)

依尺度正交的恒等式(9)，也應(yīng)有如下關(guān)系：

(11)

將式(11)代入式(10)可得：

(12)

或?qū)懗扇缦滦问剑?/p>

|H(ω)|2+|H(ω+π)|2≡1

(13)

2 Daubechies濾波函數(shù)H(ω)系數(shù)的計算

依據(jù)三角函數(shù)立即可看出如下的一組解，設(shè)H1(ω)=cos(ω/2)，則H1(ω+π)=cos[(ω+π)/2]=-sin(ω/2)，它符合式(13)，這就是Haar小波。依式(6)可得：

(14)

(15)

Daubechies思路的巧妙在于，依據(jù)三角函數(shù)的二項式高次冪展開，而得出的更多層次的DB小波。為了進一步說明Daubechies各種層次下的DB小波系數(shù)的計算方法，可依據(jù)恒等式(16)展開，并將它分成兩部分：

(16)

為了書寫簡化，令C=cos(ω/2)；S=sin(ω/2)，取二項式展開的前半部分為|H(ω)|2，如公式(17)所示，二項式的后半部分為|H(ω+π)|2，如公式(18)所示：

(17)

(18)

式(17)中直接將S?C和C?S即可得式(18)。上面的系數(shù)就是二項式展開的楊輝系數(shù)，計算簡單。關(guān)鍵是要求出濾波函數(shù)H(ω)的系數(shù)。

(19)

下面舉例如下： 令N=5，依式(17)可得：

(20)

令：

A(ω)=C8+9C6S2+36C4S4+84C2S6+126S8

(21)

令：

(22)

可得高次多項式，通過解高次方程，將其分解成多個二項式的乘積，即：

(23)

將T還原到三角函數(shù)：

A(ω)=|B(ω)|2=|D1(ω)·D2(ω)|2=

(24)

為了求出|B(ω)|，可將每個二次多項式配方，轉(zhuǎn)成|Dk(ω)|2=|u2+v2|模式，求出共軛復(fù)數(shù)，在式(25)中任取一個復(fù)數(shù)就是它的解。

Dk(ω)=uk±jvk

(25)

(26)

對其中每一個因子分式代入C2=(1+cosω)/2；

S2=(1-cosω)/2；CS=sinω/2。并令z=ejω，并|z|=1可得：

(27)

其中：

(28)

注意C5=z-5/2[(1+z)/2]5，于是一般形式可得：

|H(ω)|=|C5·B(ω)|=

(29)

可取：

(30)

當(dāng)式(28)中取正符號時有M2k>M0k，而取負(fù)符號時有M2k

表1 實例計算比較Daubechies的DB小波和SYM對稱小波尺度波形系數(shù)

Daubechies小波DB10和SYM10小波波形的對比圖見圖1。可見，SYM10的小波相比于DB10小波波形更具有對稱性。

圖1 Daubechies小波DB10和SYM10小波波形的對比

3 多尺度DB小波系數(shù)的計算

依據(jù)式(3)連續(xù)變換Φ(ω/2k)乘積可得：

(31)

采用連乘積的辦法得到式(31)，似乎不能求解出Φ(2ω)，其實從遞推公式的角度看Daubechies尺度函數(shù)的解已經(jīng)求出。首先注意對于等時間隔采樣的函數(shù)而言，最高次頻率是有限的，ω≤ωmax所以當(dāng)N→∞時有ξ=ω/2N?0。因而有Φ(ω/2N)=Φ(ξ)?Φ(0)，而Φ(0)≠0，否則會有Φ(2ω)≡0，這不合理。可取Φ(ξ)=Φ(0)=1，于是有遞推關(guān)系式：

(32)

由于傅里葉反變換，有：

zn=ejnξ?δ(t-n)

(33)

求出尺度波形是一串離散的數(shù)據(jù)，即：

(34)

(35)

(36)

對于多尺度小波的計算，可從遞推公式著手，并注意傅里葉變換的對等關(guān)系：

(37)

(38)

其中uk(k=0,1,2,...,M-1)共有M項，一般說對于DB2N小波在第k尺度下，尺度波形的點數(shù)M可通過式(39)計算，且uk滿足式(40)的關(guān)系。

M=(2N-1)·(2k-1)

(39)

(40)

第3尺度的尺度波形共有M個點，即：

(41)

圖2、圖3中分別示出DB10小波和SYM10近似對稱小波在3種尺度波形和小波波形的對照圖。可見，不論是DB10小波還是SYM10小波，其尺度1～3下的小波波形和尺度波形的特征相似，時間跨度依次以2倍遞增。

圖2 DB10小波在3種尺度下的尺度波形Φ(t)和小波函數(shù)ψ(t)

圖3 SYM10小波在3種尺度下的尺度波形Φ(t)和小波函數(shù)ψ(t)

4 小波矩陣與小波的分解與重構(gòu)

為了清楚地表達(dá)小波的正交性和小波分析的分解和重構(gòu)，特引出小波矩陣[9-10]，現(xiàn)以DB4小波為例，其尺度波形的數(shù)據(jù)如下[11]：

(42)

對應(yīng)的也有小波波形4個數(shù)據(jù)gk=(1)kh3k，限于矩陣篇幅，設(shè)定采樣數(shù)據(jù)只有8個點(一般可擴展到2N個數(shù)據(jù))，小波矩陣B[8×8]，見式(43)，從第0行開始排列尺度波形系數(shù)，而第1行參數(shù)相對第0行位移2個列，其余都是0。到第3行有h2、h3超出矩陣B[8×8]范圍，需要循環(huán)移位到第0列、第1列中。同理從第4行開始排列小波波形系數(shù)。B矩陣中上面4行構(gòu)成尺度波形矩陣H[4,8]，下面4行構(gòu)成小波波形矩陣。B矩陣的優(yōu)點是它具有正交性，即它和自身的轉(zhuǎn)置矩陣互為逆矩陣，即BTB=E，其中E是單位矩陣,符合公式(44)。

(43)

(44)

B矩陣可以按2次冪擴展，如對128點采樣數(shù)據(jù)F[128]作小波的分解與重構(gòu)，B矩陣就是128階的方陣B[128×128]，小波分解計算如下：

(45)

其中C[64]是反映采樣波形沿分段的時間軸與尺度波形的相關(guān)系數(shù)。而D[64]是反映采樣波形沿分段的時間軸與小波波形的相關(guān)系數(shù)，數(shù)據(jù)減少一半，稱為二倍壓縮。雙尺度分解可以通過小波重構(gòu)得到低頻的主波F1[128]和高頻的細(xì)波X1[128]，稱為二倍膨脹。計算公式如下：

F1[128]=HT[128,64]·C[64]

(46)

X1[128]=GT[128,64]·D[64]

(47)

依據(jù)小波矩陣正交的式(43)可知，低頻主波F1[128]和高頻細(xì)波X1[128]的合成還原了原采樣波形F[128]，如下：

F[128]=F1[128]+X1[128]=HTH·F+GTG·F=F

(48)

5 更多的正交小波引入

依照文中推證的方法，可以很方便地引出一系列另類的正交小波，如把原來式(20)的|H(ω)|2經(jīng)過式(49)改寫成|Hm(ω)|2：

(49)

當(dāng)其中參變數(shù)m若滿足126>m>0的關(guān)系時，則|Hm(ω)|2也滿足如下的恒等式關(guān)系：

|Hm(ω)|2+|Hm(ω+π)|2≡1

(50)

6 工程應(yīng)用案例

小波分析應(yīng)用的范圍很廣，涉及各個領(lǐng)域。如圖4所示,對三種128點采樣波形，經(jīng)過采用DB10小波的雙尺度分解后，得出的“低頻主波”和“高頻細(xì)波”[12]，依據(jù)波形頻次不同，三種波形分解和重構(gòu)所采用尺度等級各不相同，三種波形的特征如圖4所示。

圖4 典型電氣工程波形的小波分解Fig.4 Wavelet decomposition of typical electrical engineering waveform

故障時刻：小波分析主要指出故障發(fā)生時刻，和暫態(tài)電流的高次衰減波形，后者有時可用于接地選線的判斷分析。其尺度時譜和低頻主波的趨勢是一致的。

勵磁涌流：即變壓器空載合閘電流，可通過低頻主波的分析，可用二次諧波制動來避免合閘勵磁涌流的誤動，也可分析剩余磁通衰減過程。

瓷瓶(絕緣子)放電：這是多個電容串聯(lián)分壓的模型，放電瓷絕緣子相當(dāng)于放電電容并聯(lián)一放電間隙。當(dāng)正常充電達(dá)到該瓷絕緣子間隙的放電電壓時，瓷絕緣子放電，相應(yīng)電容電壓降為零，絕緣恢復(fù)。接著充電又開始，反復(fù)放電，其平均電壓很低，它的低頻主波也不是正弦波。

目前柔性輸電、諧波和無功治理等專業(yè)領(lǐng)域都涉及開關(guān)函數(shù)波形，圖5所示例子為采用雙正交小波對開關(guān)函數(shù)波形進行3級尺度的小波分解，得出低頻主波為變化的三角波。

圖5 用雙正交小波分解開關(guān)函數(shù)波形

在小波分解和重構(gòu)中還存在邊緣失真的問題，如圖6中128點的采樣波形，開始點的數(shù)值比末尾點幅值高很多，即首尾數(shù)值存在變化劇烈的差值。由于正交小波矩陣是循環(huán)移位構(gòu)成的，會發(fā)生首尾邊緣失真問題，如圖6中E、F兩點所示。

圖6 波形首尾數(shù)值變化劇烈而產(chǎn)生的邊緣失真Fig.6 Edge distortion caused by drastic changes in the

解決邊緣失真問題的方法，可使用鏡像法，將原始128點數(shù)據(jù)擴展成256點，此時首尾點的數(shù)值相等，再進行256點擴展數(shù)據(jù)的小波分析。如此就能完全消除邊緣失真效應(yīng)，如圖7所示。為了和圖6比較，其中只顯示前128點的小波分析結(jié)果。

圖7 通過鏡像擴展一倍點數(shù)后進行小波分析的結(jié)果Fig.7 Result of wavelet analysis after doubling the number of points by mirroring

7 結(jié)束語

文章完整地把經(jīng)典小波理論用一般的高等數(shù)學(xué)加以演繹推導(dǎo)，通俗易懂。按照文中的思路，提出了更多新的正交小波尺度波形構(gòu)造方法。工程應(yīng)用中，對于波形首尾數(shù)值變化劇烈而產(chǎn)生的邊緣失真問題，采用鏡像擴展一倍點數(shù)后進行小波分析可有效消除失真效應(yīng)。