維納濾波的使用

第五章維納濾波

〔WienerFiltering〕王靜生物醫(yī)學(xué)工程系1a本章內(nèi)容：5.1維納濾波器的時(shí)域解5.2維納預(yù)測(cè)器5.3維納濾波器的應(yīng)用2a引言〔1〕維納維納于1894年生于美國(guó)密蘇里州哥倫比亞市的一個(gè)猶太人的家庭中。他18歲就獲得哈佛大學(xué)數(shù)學(xué)和哲學(xué)兩個(gè)博士學(xué)位，是信息論的前驅(qū)和控制論的奠基人。

3a維納：開創(chuàng)維納信息論他從帶直流電流或者至少可看作直流電流的電路出發(fā)來研究信息論，將統(tǒng)計(jì)方法引入通訊工程，奠定了信息論的理論根底。創(chuàng)立控制論1947年10月，維納寫出了劃時(shí)代的著作?控制論?，1948年出版后，立即風(fēng)行世界。它揭示了機(jī)器中的通信和控制機(jī)能與人的神經(jīng)、感覺機(jī)能的共同規(guī)律；為現(xiàn)代科學(xué)技術(shù)研究提供了嶄新的科學(xué)方法。4a引言在前面章節(jié)中已經(jīng)介紹了噪聲中確定性信號(hào)的檢測(cè)和估計(jì)；隨機(jī)性是生物醫(yī)學(xué)信號(hào)的特點(diǎn)之一，因此，討論噪聲中隨機(jī)信號(hào)的估計(jì)具有現(xiàn)實(shí)意義。維納濾波器正是解決該問題。維納濾波器的限制：要求被估計(jì)的隨機(jī)信號(hào)是平穩(wěn)的。5a〔2〕維納濾波技術(shù)：從噪聲中提取有用的平穩(wěn)隨機(jī)信號(hào)。h(n)x(n)=s(n)+w(n)y(n)=6a7a8a根底知識(shí)點(diǎn)回憶1.卷積運(yùn)算2.相關(guān)運(yùn)算3.信號(hào)與系統(tǒng)9a設(shè)有一個(gè)線性系統(tǒng)，它的單位脈沖響應(yīng)是h(n)，當(dāng)輸入一個(gè)觀測(cè)到的隨機(jī)信號(hào)x(n)，簡(jiǎn)稱觀測(cè)值，且該信號(hào)包含噪聲w(n)和有用信號(hào)s(n)，也即x(n)=s(n)+w(n)那么輸出y(n)為y(n)=x(n)*h(n)=

10a我們希望輸出得到的y(n)與有用信號(hào)s(n)盡量接近，因此稱y(n)為s(n)的估計(jì)值，用來表示y(n)，我們就有了維納濾波器的系統(tǒng)框圖．這個(gè)系統(tǒng)的單位脈沖響應(yīng)也稱為對(duì)于s(n)的一種估計(jì)器。h(n)x(n)=s(n)+w(n)y(n)=11a維納濾波技術(shù)可應(yīng)用于以下三個(gè)方面：濾波：用當(dāng)前的和過去的觀測(cè)值來估計(jì)當(dāng)前的信號(hào)，稱為濾波；預(yù)測(cè)：用過去的觀測(cè)值來估計(jì)當(dāng)前的或?qū)淼男盘?hào)，稱為預(yù)測(cè)；平滑或內(nèi)插：用過去的觀測(cè)值來估計(jì)過去的信號(hào)，稱為平滑或者內(nèi)插。12a系統(tǒng)框圖中估計(jì)到的信號(hào)和我們期望得到的有用信號(hào)s(n)不可能完全相同，這里用e(n)來表示真值和估計(jì)值之間的誤差

h(n)x(n)=s(n)+w(n)y(n)=顯然e(n)是隨機(jī)變量，維納濾波和卡爾曼濾波的誤差準(zhǔn)那么就是最小均方誤差準(zhǔn)那么

13a５.１維納濾波器的時(shí)域解

〔TimedomainsolutionoftheWienerfilter〕設(shè)計(jì)維納濾波器的過程就是尋求在最小均方誤差下濾波器的單位脈沖響應(yīng)h(n)或傳遞函數(shù)H(z)的表達(dá)式，其實(shí)質(zhì)就是解維納－霍夫〔Wiener－Hopf〕方程。

14a5.1.1因果的維納濾波器，即維納-霍夫方程

5.1.2有限脈沖響應(yīng)法求解維納-霍夫方程

5.1.3預(yù)白化法求解維納-霍夫方程15a因果的維納濾波器設(shè)h(n)是物理可實(shí)現(xiàn)的，也即是因果序列：h(n)=0,當(dāng)n<0

16a要使得均方誤差最小，那么將上式對(duì)各h(m),m＝0，1，…，求偏導(dǎo)，并且等于零，得： (5-7)即 (5-8) (5-9)17a從維納－霍夫方程中解出的h就是最小均方誤差下的最正確h，hopt(n)。求到hopt(n)，這時(shí)的均方誤差為最?。?8a有限脈沖響應(yīng)法求解維納－霍夫方程

設(shè)h(n)是一個(gè)因果序列且可以用有限長(zhǎng)〔N點(diǎn)長(zhǎng)〕的序列去逼進(jìn)它，那么式(5-5)－(5-10)分別發(fā)生變化： (5-11) (5-12)

(5-13)19a

(5-14) (5-15)于是得到N個(gè)線性方程：

20a寫成矩陣形式有：

(5-16)簡(jiǎn)化形式：

RxxH=Rxs(5-17)

式中，H＝[h(0)h(1)…h(huán)(N-1)]′是待求的單位脈沖響應(yīng):21aRxxH=Rxs

只要Rxx是非奇異的，就可以求到H：

H=Rxx－1Rxs

(5-18)求得hopt(n)后，這時(shí)的均方誤差為最?。?/p>

22a23a假設(shè)信號(hào)與噪聲互不相關(guān)，即，Rsw(m)=Rws(m)=0那么有Rxs(m)=E[x(n)s(n+m)]=E[(s(n)+w(n))s(n+m)]=E[s(n)s(n+m)+w(n)s(n+m)]=Rss(m)Rxx(m)=E[(s(n)+w(n))(s(n+m)+w(n+m))]=Rss(m)+Rww(m)24a那么式〔5－15〕和式〔5－19〕化為： (5-20) (5-21)【例5-1】如圖，x(n)=s(n)+w(n),且s(n)與w(n)統(tǒng)計(jì)獨(dú)立，其中s(n)的自相關(guān)序列 ,w(n)是方差為1的單位白噪聲，試設(shè)計(jì)一個(gè)N＝2的維納濾波器來估計(jì)s(n)，并求最小均方誤差。25a解：依題意，信號(hào)的自相關(guān)和噪聲的自相關(guān)為：代入式〔5－20〕得:解得：h(0)＝，h(1)＝。將上述結(jié)果代入式〔5－21〕，求得最小均方誤差：h(n)x(n)=s(n)+w(n)y(n)=26a5.1.1因果的維納濾波器，即維納-霍夫方程

5.1.2有限脈沖響應(yīng)法求解維納-霍夫方程

5.1.3預(yù)白化法求解維納-霍夫方程27ah(n)x(n)=s(n)+w(n)y(n)=G(z)x(n)y(n)=1/B(z)w1(n)28a預(yù)白化法求解維納－霍夫方程

隨機(jī)信號(hào)都可以看成是由一白色噪聲w1(n)鼓勵(lì)一個(gè)物理可實(shí)現(xiàn)的系統(tǒng)或模型的響應(yīng)，如下圖．

圖5.2s(n)的信號(hào)模型A(z)w1(n)s(n)29a由于x(n)=s(n)+w(n)，在圖的根底上給出x(n)的信號(hào)模型，圖所示。把這兩個(gè)模型合并最后得到維納濾波器的信號(hào)模型，圖所示，其中傳遞函數(shù)用B〔z〕表示。

圖5.3x(n)的信號(hào)模型s(n)A(z)w1(n)w(n)x(n)30a維納濾波器的輸入信號(hào)模型B(z)w1(n)x(n)31a白噪聲的自相關(guān)函數(shù)為它的z變換就等于。圖中輸出信號(hào)的自相關(guān)函數(shù)為，根據(jù)卷積性質(zhì)有令l=r-k，上式32a令代入上式得

(5-22)對(duì)式〔5－22〕進(jìn)行Z變換得到系統(tǒng)函數(shù)和相關(guān)函數(shù)的z變換之間的關(guān)系： (5-23)同樣，對(duì)圖進(jìn)行z變換得 (5-24)33a圖中利用卷積性質(zhì)還可以找到互相關(guān)函數(shù)之間的關(guān)系：

兩邊z變換得到

(5-25)

34a如果觀測(cè)信號(hào)的自相關(guān)函數(shù)，求它的z變換，然后找到該函數(shù)的成對(duì)零點(diǎn)、極點(diǎn)，取其中在單位圓內(nèi)的那一半零點(diǎn)、極點(diǎn)構(gòu)成B(z),另外在單位圓外的零、極點(diǎn)構(gòu)成B(z-1)，這樣就保證了B(z)是因果的，并且是最小相位系統(tǒng)。從圖可得(5-26)35a由于系統(tǒng)函數(shù)B(z)的零點(diǎn)和極點(diǎn)都在單位圓內(nèi)，即是一個(gè)物理可實(shí)現(xiàn)的最小相位系統(tǒng)，那么1/B(z)也是一個(gè)物理可實(shí)現(xiàn)的最小相移網(wǎng)絡(luò)函數(shù)。我們就可以利用式〔5－25〕對(duì)x(n)進(jìn)行白化，即把x(n)當(dāng)作輸入，w1(n)當(dāng)作輸出，1/B(z)系統(tǒng)傳遞函數(shù)。36ah(n)x(n)=s(n)+w(n)y(n)=G(z)x(n)y(n)=1/B(z)w1(n)H(z)=G(z)/B(z)37a白化法求解維納－霍夫方程步驟如下：１〕對(duì)觀測(cè)信號(hào)x(n)的自相關(guān)函數(shù)Rxx(m)求z變換得到Rxx(z)２〕利用等式找到最小相位系統(tǒng)B(z)３〕利用均方誤差最小原那么求解因果的G〔z〕４〕H(z)=G(z)/B(z),即得到維納－霍夫方程的系統(tǒng)函數(shù)解38a步驟3的求解過程按圖〔b〕有 (5-28)均方誤差為39a由于代入上式，并且進(jìn)行配方得

(5-29)

均方誤差最小也就是上式的中間一項(xiàng)最小，所以

(5-30)40a注意，這里的g(m)是因果的。對(duì)該式求z變換,得到

(5-31)

表示對(duì)求單邊z變換。所以維納－霍夫方程的系統(tǒng)函數(shù)解表示為41a由式

(5-32)因果的維納濾波器的最小均方誤差為：

(5-33)42a利用帕塞伐爾定理，上式可用z域來表示

(5-34)例5－2。。。見書43a5.2.1因果的維納預(yù)測(cè)器5.2.2純預(yù)測(cè)器〔N步〕5.2.3一步線性預(yù)測(cè)器5.2維納預(yù)測(cè)器44a因果的維納預(yù)測(cè)器圖就是維納預(yù)測(cè)器的模型，N>0，yd(n)是希望得到的輸出，而y(n)表示實(shí)際的估計(jì)值。

圖維納預(yù)測(cè)器y(n)=yd(n)=h(n)x(n)=s(n)+w(n)45a設(shè)h(n)是物理可實(shí)現(xiàn)的，也即是因果序列：h(n)=0,當(dāng)n<0，那么有

(5-35)

(5-36)要使得均方誤差最小，那么將上式對(duì)各h(m),m＝0，1，…，求偏導(dǎo)，并且等于零，得：46a

(5-37)即用相關(guān)函數(shù)R來表達(dá)上式： (5-39)由yd(n)=s(n+N)，那么47az變換得

(5-40)因果的預(yù)測(cè)器的傳遞函數(shù)為：

(5-41)最小均方誤差為

(5-42)48a利用帕塞伐爾定理，上式可用z域來表示

(5-43)49a例5-3圖中，x(n)=s(n)+w(n),且s(n)與w(n)統(tǒng)計(jì)獨(dú)立，其中s(n)的自相關(guān)序列為，w(n)是方差為1的單位白噪聲，試設(shè)計(jì)一個(gè)物理可實(shí)現(xiàn)的維納預(yù)測(cè)器估s(n+1)，并求最小均方誤差。解：依題意50a求z變換

由于，容易找到最小相位系統(tǒng)和白噪聲方差

51a由式〔5－41〕，N＝1，

對(duì)括號(hào)里面求z反變換，注意括號(hào)內(nèi)的收斂域?yàn)?/p>

52a取因果局部，也就是第一項(xiàng)，所以

把上式寫成差分方程形式有：53a最小均方誤差為：54a純預(yù)測(cè)器〔N步〕純預(yù)測(cè)器指的是w(n)＝0的情況下，對(duì)s(n+N)的預(yù)測(cè)。如下圖。

圖5.7N步純預(yù)測(cè)器這時(shí)，用白化法來求解預(yù)測(cè)器的系統(tǒng)函數(shù)。y(n)=yd(n)=h(n)x(n)=s(n)55a因?yàn)?，從而? (5-44)將上式代入式〔5－41〕、〔5－43〕得

(5-45)56a假設(shè)B〔z〕是b〔n〕的z變換，且b〔n〕是實(shí)序列，那么上式可以利用帕塞伐爾定理進(jìn)一步化簡(jiǎn)：上式說明最小均方誤差隨著N的增加而增加，也即預(yù)測(cè)距離越遠(yuǎn)誤差越大。57a例5-4圖中，x(n)=s(n)其中s(n)的自相關(guān)序列為，試設(shè)計(jì)一個(gè)物理可實(shí)現(xiàn)的維納預(yù)測(cè)器來估計(jì)，并求最小均方誤差。解：依題意，，那么

58a因?yàn)槿菀渍业阶钚∠辔幌到y(tǒng)和白噪聲方差利用式〔5－45〕59a因?yàn)橹蝗〉木植?，有代入?最小均方誤差為：60a它說明當(dāng)N越大，誤差越大，當(dāng)N＝0時(shí)，沒有誤差

圖5.8例題5-3的純預(yù)測(cè)模型y(n)==x(n)=s(n)61a一步線性預(yù)測(cè)器對(duì)于純預(yù)測(cè)問題，有

然而預(yù)測(cè)的問題常常是要求在過去的p個(gè)觀測(cè)值的根底上來預(yù)測(cè)當(dāng)前值，也就是這就是一步線性預(yù)測(cè)公式62a一步線性預(yù)測(cè)公式，常常用以下符合表示(5-48)式中p為階數(shù)，。預(yù)測(cè)的均方誤差為

(5-49)63a要使得均方誤差最小，將上式右邊對(duì)求偏導(dǎo)并且等于零，得到p個(gè)等式： (5-50)最小均方誤差： (5-51)式〔5－50〕就是Yule－Walker〔Y-W〕方程.64a例5-5圖中，x(n)=s(n),其中s(n)的自相關(guān)序列為，試設(shè)計(jì)一個(gè)p＝2的可實(shí)現(xiàn)的一步線性預(yù)測(cè)器，并求最小均方誤差。65a解：

利用Y-W方程

可以列出2個(gè)方程式66a解得：結(jié)果和例〔5-4〕N=1時(shí)一致67a維納濾波器的應(yīng)用〔ApplicationofWienerfilter〕要設(shè)計(jì)維納濾波器必須知道觀測(cè)信號(hào)和估計(jì)信號(hào)之間的相關(guān)函數(shù)，即先驗(yàn)知識(shí)。如果我們不知道它們之間的相關(guān)函數(shù)，就必須先對(duì)它們的統(tǒng)計(jì)特性做估計(jì)，然后才能設(shè)計(jì)出維納濾波器，這樣設(shè)計(jì)出的濾波器被稱為“后驗(yàn)維納濾波器〞。68a在生物醫(yī)學(xué)信號(hào)處理中比較典型的應(yīng)用就是關(guān)于誘發(fā)腦電信號(hào)的提取。大腦誘發(fā)電位〔EvokedPotential，EP〕指在外界刺激下，從頭皮上記錄到的特異電位，它反映了外周感覺神經(jīng)、感覺通路及中樞神經(jīng)系統(tǒng)中相關(guān)結(jié)構(gòu)在特定刺激情況下的狀態(tài)反響。在神經(jīng)學(xué)研究以及臨床診斷、手術(shù)監(jiān)護(hù)中有重要意義。EP信號(hào)十分微弱，一般都淹沒在自發(fā)腦電〔EEG〕之中，從EEG背景中提取誘發(fā)電位一直是個(gè)難題：EP的幅度比自發(fā)腦電低一個(gè)數(shù)量級(jí)，無法從一次觀察中直接得到；EP的頻譜與自發(fā)腦電頻譜完全重迭，使得頻率濾波失效；在統(tǒng)計(jì)上EP是非平穩(wěn)的、時(shí)變的腦誘發(fā)電位。通過屢次刺激得到的腦電信號(hào)進(jìn)行疊加來提取EP，這是現(xiàn)今最為廣泛使用的EP提取方法69a為了解決誘發(fā)電位提取問題，研究者利用維納濾波來提高信噪比，先后有Walter、Doyle、Weerd等對(duì)維納濾波方法