1、 3.1平穩(wěn)隨機信號 3.2平穩(wěn)隨機信號通過線性系統(tǒng) 3.3經典譜估計 3.4現(xiàn)代譜估計醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法 3.1 平穩(wěn)隨機信號 自然界中存在著大量隨機現(xiàn)象如地震波幅、建筑物承受的風載荷、通信系統(tǒng)或控制系統(tǒng)中的各種噪聲和干擾等。這些隨機現(xiàn)象無法用精確的函數(shù)關系來表示 例如電阻內部微觀粒子的隨機熱騷動引起的端電壓稱為熱噪聲電壓，它在任一時刻的值都是隨機變量。對它進行一次觀察，得到一個電壓一時間函數(shù)x1（t）。這個函數(shù)是不可預先確知的，只有通過觀察才能得到，在相同條件下重復觀察，每次都將得到不同的結果。數(shù)字信號處

2、理數(shù)字信號處理 第八章第八章 平穩(wěn)隨機信號平穩(wěn)隨機信號數(shù)字信號處理數(shù)字信號處理 第八章第八章 平穩(wěn)隨機信號平穩(wěn)隨機信號 通常把對隨機現(xiàn)象的一次觀察結果 稱作一個樣本函數(shù)，對在有限時間區(qū)間 上觀察的結果稱為樣本記錄。任一次觀 察結果都僅是無限多個可能產生的結果 中的一個，而所有可能產生的樣本函數(shù) 的集合稱為隨機過程。隨機過程是描述 隨機現(xiàn)象的數(shù)學模型。隨機過程如果是 時間的函數(shù)，則隨機過程和隨機信號這 兩個概念是共通的。數(shù)字信號處理數(shù)字信號處理 第八章第八章 平穩(wěn)隨機信號平穩(wěn)隨機信號數(shù)字信號處理數(shù)字信號處理 第八章第八章 平穩(wěn)隨機信號平穩(wěn)隨機信號數(shù)字信號處理數(shù)字信號處理 第八章第八章 平穩(wěn)隨機信

3、號平穩(wěn)隨機信號 (2) 隨機序列可以用它的統(tǒng)計平均特性表征 隨機變量的取值雖不能先驗確知，卻服從確定的統(tǒng)計規(guī)律，因此可用確定的概率分布特性來統(tǒng)計地描述。我們也可以用統(tǒng)計平均特性來統(tǒng)計地表征隨機變量，因為每種統(tǒng)計平均特性均為隨機變量的某種函數(shù)按概率加權平均的結果，因而可從某個側面反映隨機變量的取值特性。隨機變量的統(tǒng)計平均特性是確定的。醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法 一、隨機序列的描述 (一)概率分布函數(shù)和概率密度函數(shù)醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法 (

4、二)統(tǒng)計平均 用概率函數(shù)可以最充分地在統(tǒng)計意義上描述一個隨機過程。但在工程實際中，要得知一個隨機序列各點隨機變量的概率函數(shù)要通過大量實驗，而許多問題的解決往往只需知道隨機序列概率函數(shù)的某些特征量，即所謂數(shù)字特征。數(shù)字特征是隨機序列的函數(shù)按概率的統(tǒng)計平均，常用的有均值、方差和相關函數(shù)。 1.均值醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法4自相關函數(shù)和自協(xié)方差函數(shù)一個隨機序列中兩個時刻n1和n2上的隨機變量x(n1)和x(n2)之間的

5、自相關函數(shù)定義為:醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻

6、域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法3醫(yī)學信號處理醫(yī)學信號處理 第三章第三

7、章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法4醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法3-1醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法3-2醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法互相關函數(shù)性質:醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的

8、時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法4醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法例3-3：一個零均值的離散時間隨機過程各值xn，xn+m（m不等于零）均互相獨立。求它的自相關函數(shù)和功率譜。自相關函數(shù)：自相關函數(shù)：Rx（m）=Exnxn+mm=0時，時， Rx（0）=Exn2=方差方差m不等于零，不等于零， Rx（m）=Exn Exn+m=0 功率

9、譜：功率譜：Px（w）=DTFTRx（m）= Rx（0） =方差。方差。典型的白噪過程典型的白噪過程醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法例3-4：把上例的xn送入到兩點平均器中，得輸出為，求yn的自相關函數(shù)和功率譜。)xx(2/1y1nnn)xx(41(E)y(E)0(R21nn2ny)xx)(xx(41(E)yy(E) 1 (Rn1n1nn1nny自相關函數(shù)：自相關函數(shù)：Ry（2）以上都為零，平均后數(shù)據(jù)相關性變大）以上都為零，平均后數(shù)據(jù)相關性變大2x21n1nn2n21)x(E41)xx(E21)E(x412x1n1 -n1

10、n1n2n1nn41)xx(E41)xx(E41)x(E41)xE(x41例3-4：把上例的xn送入到兩點平均器中，得輸出為，求yn的自相關函數(shù)和功率譜。)xx(2/1y1nnn功率譜：功率譜：j -yyjyyye ) 1 (R)0(Re ) 1(R)m(R(DTFT)(P)cos1 (22x3-5：已知平穩(wěn)隨機過程x的自相關函數(shù)如下，求其功率譜密度及均方，并根據(jù)所得結果說明該隨機過程是否含有直流分量或周期性分量。cos3cos4e)(Rx16cos25e)(R04x（）（）cos3cos4e)(Rx（）de )(R)(Pjxx)(11)(11 4)3()3(22514)0(R)x(Ex201

11、2 4)0(P2x因為因為所以含有直流分量；所以含有直流分量；因為周期信號的自相關函數(shù)也是周期性的，而因為周期信號的自相關函數(shù)也是周期性的，而R中包含中包含有一個周期性的成分，因此該隨機過程含有周期性分量有一個周期性的成分，因此該隨機過程含有周期性分量。 16cos25e)(R04x（）de )(R)(Pjxx)(161)(161100)(322020411625)0(R)x(Ex2016210032)0(P20 x因為因為所以含有直流分量；所以含有直流分量；因為周期信號的自相關函數(shù)也是周期性的，而因為周期信號的自相關函數(shù)也是周期性的，而R中中沒有包含周期性的成分，因此該隨機過程不含有周沒有包

12、含周期性的成分，因此該隨機過程不含有周期性分量。期性分量。 X(n) h(n) Y(n) H() Fx() Fy()mx myRx(m) Ry(m)Sx() Sy() Rxy(m) Sxy() 僅取決于mx、Rx(m)與Sx()3.2 隨機信號通過線性系統(tǒng)隨機信號通過線性系統(tǒng)醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法y(n)=h(n)x(n)=kknxkh)()(1.輸出均值my(n)=Ey(n)= h(k)Ex(n-k)k= h(k)mx(n-k)k=mx h(k)k=mxH(0)醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信

13、號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法2.互相關函數(shù)Rxy(n,n+m)=Ex(n)y(n+m)=Ex(n) h(k)x(n+m-k)= h(k)Ex(n)x(n+m-k)= h(k)Rx(n,n+m-k)Rxy(m)= h(k)Rx(m-k)=h(m)Rx(m)kkkk醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法3.互功率譜密度Sxy(ej )= Rxy(m)e-j m= h(k)Rx(m-k)e-j m= h(k)e-j k Rx(m-k)e-j (m-k)=H(ej )Sx(ej )mkkmmh(m)Rx(m)醫(yī)

14、學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法同理可得Ryx(m)= h(k)Rx(m+k)=h(-m)Rx(m)Syx(e j )=H(e -j )Sx(e j )k4.輸出的自相關函數(shù)Ry(m)=Ey(n)y(n+m)醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法Ry(m)=E h(k)x(n-k) h(j)x(n+m-j)= h(k)h(j)Ex(n-k)x(n+m-j)= h(k) h(j)Rx(m-j+k)= h(k)h(m+k) Rx(m+k)=h(-m) h(m) R

15、x(m)=h(-m) Rxy(m)=Ryx (m) h(m)kjkjkjkh(m)h(m) mx my my H(0) mx Rxy(m)Rx(m) Ry(m) Ryx(m)h(-m)h(m)h(-m)醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法5.輸出的功率譜密度Sy(e j )= Ry(m)e- j m =Syx(e j )H(e j ) =Sxy(e j )H(e- j ) =Sx(e j )H(e j )H(e -j )mH(e j )H(e- j)Sx(e j )Sy(e j )對應Ry(m)= Rxy(m) h(-m)Ry

16、(m)= Ryx (m) h(m)Ry(m)=Rx(m) h(m) )h(-m)醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法例3-6：已知一線性系統(tǒng)輸入為零均值平穩(wěn)隨機序列x(n)，自相關函數(shù)Rx(m)=x2(m)則：Rxy(m)=Rx(m)h(m) = x2h(m) h(m)=Rxy(m)/ x2給未知系統(tǒng)輸入白噪聲，求出輸入與輸出的互相關函數(shù)，即可得到系統(tǒng)的單位沖激響應，稱為系統(tǒng)辨識。頻域有：Sxy( )=Sx( )H( ) H( )= Sxy( )/ Sx( )= Sxy( )/ x2醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物

17、醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法例3-7：設平穩(wěn)隨機信號x(n)的自相關函數(shù)為x2(m)，線性系統(tǒng)的單位沖激響應是h(k)=rk k 0 |r|r H(e j )=1/(1-r e - j) Sx(e j )=Rx(m) e - j m=x2 Sy(e j )= Sx(e j )H(e j )H(e -j ) = x2/( 1-2rcos +r2)Rxy(m)=Rx(m) h(m)=rk x2(m-k)=rm x2醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法平均功率Ry(m)= Rxy(m) h(-m) =

18、 rk x2r k+m = x2r m/(1-r2) m0 Ry(m)= x2r -m/(1-r2) m0 Ry(m)= x2rm/(1-r2)Dy=Ry(0)= x2/(1-r2) =Ey2(n)=1/(2)Sy( )d功率譜密度醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法例題3-8：線性系統(tǒng)如圖所示，輸入為零均值，方差為1的白噪，求輸出的自相關和功率譜。2acosa11)e (P)e (H)e (P2jx2jjy解：解：Y(z)=aY(z)z-1+X(z)H(z)=1/(1-az-1)ae1/(1)e (Hjjw)e (Pjy+a

19、Z-1x（n） y（n） 2m2jyya11a2acosa11)e (P)m(RIDTFTIDTFT信號與線性系統(tǒng)P168常用傅氏變換對3.3 經典譜估計3.3.1概述 從觀測到的信號數(shù)據(jù)中，定量地預測或推斷被觀測系統(tǒng)的狀態(tài)或參量的問題，就是估計問題。 因為觀測到的信號數(shù)據(jù)是有限的，因此實際估計量與被估計量或真值之間必然存在著估計誤差。如何使誤差最小，使估計量與被估計量近于相等，這就是統(tǒng)計估計從理論上和技術上所要討論和解決的問題。醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法功率譜估計方法經典譜估計現(xiàn)代譜估計周期圖法相關法改進法從方法上

20、分從信號特點分參數(shù)模型法非參數(shù)法ARMAARMA平穩(wěn)非平穩(wěn)短時付氏變換Wigner分布本節(jié)討論要到下節(jié)才會講醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法3.3.2 基本的估計算法理論上：Rx(m) = x(n)x(n+m) =Ex(n)x(n+m)Rx(m) = xk(n)xk(n+m)Sx(e j )= Rx(m) e - j m (1)121limlllln MM1limMk 1求和有無窮多項！m醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法R實際上： x(m)=1/N x

21、(n)x(n+m) -(N-1) m (N-1) x(m)=1/M xk(n)xk(n+m)10NnRMk 1求和為有限項！真正的相關函數(shù)的一個近似。代入（1）式 x(e j )= x(m) e-j m上式為真實功率譜的一個估計。稱作相關法。S1) 1(NNmR類似地用相關法估計互譜為： xy(m)=1/N x(n)y(n+m) (2) xy(e j )= xy(m) e-j m (3)將（2）式代入（3）式 xy(e j )= 1/N xN(n)yN(n+m) e-j m R10NnS1) 1(NNmRS1) 1(NNm10Nn xy(e j )= 1/N xN (n) yN (n+m) e

22、 -j m =1/N xN (n) yN (l) e -j (l-n) m+n=l= 1/N xN (n) e j n yN (l) e -j l= 1/N XN* (e j ) YN (e j ) (4) x(e j )= 1/N XN* (e j ) XN (e j ) (5)S10Nn1) 1(NNm10Nn1)1(NnNnl10Nn10NlS 醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法3.3.3 估計質量的評介設是隨機信號的一個統(tǒng)計特征， 是的估計，它是N個隨機變量x(0)、x(N-1)的函數(shù)，因而也是隨機的。如果 接近于的

23、概率很大，即可以認為是個“好”的估計。其中偏差和方差是評定估計質量時通用的兩個量。估計的偏差定義為：Bia = -E 可以是均值、自相關函數(shù)或功率譜.(p）偏差表示估計的均值與真值的接近程度。若Bia( )為零稱 是 的無偏估計無偏估計。否則稱為有偏估計有偏估計。若用于求均值的樣本記錄趨于無限長時有 Bia =0則稱 是的漸近無偏估計漸近無偏估計。Nlim估計的方差定義為：Var =E( -E )2表示各次估計值相對于估計值均值的分散程度。方差小且無偏方差小但有偏方差大無偏方差為最小的估計稱為最小方差估計。最小方差估計。無偏估計意味著做多次估計后，各次估計值的均值接近于待估計量的真值，并不保證

24、每次估計值都接近于真值，落點可能很分散。最小方差估計的每次估計值都接近于估計的均值，落點很集中。但并不保證這些落點都集中在真值附近。(p）醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法估計的均方誤差均方誤差：E( - )2=(Bia )2+Var 如果隨著樣本數(shù)目的增加，估計的均方誤差趨于零，稱此估計為一致估計為一致估計。偏差與方差都趨于零醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法000) 1(0) 1() 2() 1() 2() 1 () 2() 1 () 0() 1 ()

25、 0(0) 0(000) 1() 1 () 0(NyNyNyNyyyyyyyyyNxxx10)()()(Nnxymnynxmr 3.3.4 相關函數(shù)的估計相關函數(shù)的估計 設序列設序列x,y長度為長度為N點，除區(qū)間點，除區(qū)間0N-1之外皆為零，之外皆為零，用矩陣的形式來表達線性相關：用矩陣的形式來表達線性相關： 計算得到一個計算得到一個2N1點長的行向量，也就是對點長的行向量，也就是對應，應，m(N1)，（，（N1）。如果）。如果x和和y的長度不同，則把短的序列進行補零，使得兩的長度不同，則把短的序列進行補零，使得兩者點長相同，然后者點長相同，然后計算 .醫(yī)學信號處理醫(yī)學信號處理 第三章第三章

26、生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法 設序列設序列x,y長度為長度為N點，除區(qū)間點，除區(qū)間0N-1之外皆為零，之外皆為零，用矩陣的形式來表達循環(huán)相關：用矩陣的形式來表達循環(huán)相關： 計算得到一個計算得到一個N點長的行向量，也就是對應，點長的行向量，也就是對應，m=0，1，（，（N1）。如果）。如果x和和y的長度不同，的長度不同，則把短的序列進行補零，使得兩者點長相同，則把短的序列進行補零，使得兩者點長相同，然后然后計算 .)2()0() 1()0()2() 1 () 1() 1 ()0() 1() 1 ()0(NyyNyyyyNyyyNxxx10)()()()(NnN

27、NxynRmnynxmr醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法一、相關函數(shù)的估計l根據(jù)定義用有限樣本來進行估計：xxN1lim)m(R1N0nmnnNx假設只有假設只有N個數(shù)據(jù)，估計公式為：個數(shù)據(jù)，估計公式為：, 2, 1, 0m,xxN1)m(R1mN0nmnnx1mNn1Nmnxmn不超過數(shù)據(jù)長度為保證醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信

28、號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法偏差：看估計的均值，是有偏估計，但是漸進無偏。方差：估計的方差當N無窮時，趨于零。因此該估計法是一致估計。相關函數(shù)估計的質量：醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法二、FFT法估計相關函數(shù)醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號

29、的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法線性卷積與線性相關的關系：醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法l隨機信號的功率譜反映它的頻率成分以及各成分的相對強弱。l經典估計法：1）先估計相關函數(shù)，然后傅立葉變換；

30、2）對信號傅立葉變換后求模平方再除以N。l估計的方差特性不好，起伏劇烈，數(shù)據(jù)越長越嚴重。醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法3.3.5譜估計1、周期圖法估計功率譜l2、相關法譜估計、相關法譜估計l MNl當M=N時周期圖法和自相關法結果是一致的2jjxj)e (XN1)e (P)n(xDTFT)e (X即可。以計算，再求模平方，除然后點長補零至一致，把為使得結果與自相關法NFFT,1-2Nx, 2, 1, 0m,xxN1)m(R1mN0nmnnx x(e j )= x(m) e-j mS1) 1(MMmR醫(yī)學信號處理醫(yī)學信號處

31、理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法1.周期圖法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法功率譜估計的質量偏差：看估計的均值，是有偏估計，但是漸進無偏。方差：估計的方差當N無窮時，估計值不會趨于零。因此該估計法不是一致估計。醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及

32、頻域分析方法生物醫(yī)學信號的時域及頻域分析方法3、改進法估計功率譜1.平滑：加窗對單一功率譜估計加以平滑。2.平均：對同一隨機過程做多次周期圖法，再加以平均。估計的質量：均值是漸進無偏，方差是趨于零，是一致估計。醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)

33、學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法 3.4現(xiàn)代譜估計醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法 為隨機信號建立參數(shù)模型是研究隨機信

34、號的一種基本方法。 其含義是認為隨機信號是由白噪激勵某一確定系統(tǒng)的響應（圖3.1）。 只要白噪的參數(shù)確定了，研究隨機信號就可以轉化成研究產生隨機信號的系統(tǒng)。)(nw)(nx)(zH)(zW)(zX)(nh圖圖3.1隨機信號的參數(shù)模型隨機信號的參數(shù)模型 經典信號建模法（classical modeling method for signal）：前面已經指出，醫(yī)學信號處理的目的是提取包含于隨機信號中的確定性成分，以便在一定的準確性（最小二乘意義）上進行預測。 這就是建立各種各樣的確定性數(shù)學模型，包括代數(shù)、微分、積分、差分方程模型。這是經典的信號建模方法。 信號的現(xiàn)代建模方法（Modern mode

35、ling method for signal）：是建立在具有最大的不確定性基礎上的預測。 提出了眾多的數(shù)學模型( mathematical models)。根據(jù)Wold的證明：任何平穩(wěn)的ARMA（自回歸移動平均）模型或MA模型均可用無限階或階數(shù)足夠的AR模型去近似。因此本節(jié)著重介紹AR模型的基本原理和方法。隨機信號 由當前的激勵 和若干次過去的激勵 線性組合產生： )(nx)(nw)(knwqkkknwbnx0)()(qkkkzbzWzXzH0)()()(1 三種參數(shù)模型三種參數(shù)模型1) MA模型模型(moving average model)該模型的系統(tǒng)函數(shù)是：該模型的系統(tǒng)函數(shù)是： q表示系

36、統(tǒng)階數(shù)，系統(tǒng)函數(shù)只有零點，沒有極點，所表示系統(tǒng)階數(shù)，系統(tǒng)函數(shù)只有零點，沒有極點，所以該系統(tǒng)一定是穩(wěn)定的系統(tǒng)，也稱為全零點模型，以該系統(tǒng)一定是穩(wěn)定的系統(tǒng)，也稱為全零點模型，用用MA（ q ）來表示。）來表示。隨機信號 由本身的若干次過去值 和當前的激勵值 線性組合產生： )(nx)(knx)(nwpkkknxanwnx1)()()(pkkkzazH111)(2) AR模型模型(auto-regression model)該模型的系統(tǒng)函數(shù)是：該模型的系統(tǒng)函數(shù)是：p是系統(tǒng)階數(shù)，系統(tǒng)函數(shù)中只有極點，無零點，也稱為是系統(tǒng)階數(shù)，系統(tǒng)函數(shù)中只有極點，無零點，也稱為全極點模型，系統(tǒng)由于極點的原因，要考慮到系統(tǒng)

37、全極點模型，系統(tǒng)由于極點的原因，要考慮到系統(tǒng)的穩(wěn)定性，因而要注意極點的分布位置，用的穩(wěn)定性，因而要注意極點的分布位置，用AR（ p ）來表示。）來表示。它既有零點又有極點，所以也稱極零點模型，要考慮極零點的分布位置，保證系統(tǒng)的穩(wěn)定，用ARMR（ p ，q ）表示。pkkqkkknxaknwbnx10)()()(pkkkqkkkzazbzH101)(3 ) ARMA模型模型(auto-regression-moving average model)ARMA是是AR與與MA模型的結合：模型的結合：該模型的系統(tǒng)函數(shù)是：該模型的系統(tǒng)函數(shù)是：2 AR模型參數(shù)的估計隨機信號的建模法最近在生物醫(yī)學信號處理中

38、應用相當普遍，在自發(fā)腦電、誘發(fā)腦電、肌電、心電、胃電等方面都有人嘗試應用模型法進行研究。應用較多的是AR模型，因為建立這種的模型計算工作比較容易。作為數(shù)學逼進，三種模型都可以互相轉換。實際中選用哪一種模型就要考慮到節(jié)約和計算量，選定模型后，剩下的任務就是用適當?shù)乃惴ü烙嬆Ｐ蛥?shù)（ 、 、 、 ），以便用模型對隨機信號進行預測。 kakbpq下面我們就以下面我們就以AR模型為例進行參數(shù)的估計。模型為例進行參數(shù)的估計。 對該式兩邊同時乘以對該式兩邊同時乘以 ，然后求均值：，然后求均值： 1) AR模型參數(shù)和自相關函數(shù)的關系根據(jù)圖3.1：pkkknxanwnx1)()()()(mnx )()()()

39、()()(1pkkmnxknxamnxnwEmnxnxE醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法 系統(tǒng)的單位脈沖響應 是因果的，所以輸出的平穩(wěn)隨機信號和輸入的白噪聲之間的互相關函數(shù)有下列推導： )()()()()()()()(mRmnxnxEkxmkxEmnxnxEmRxxxxpkxxkxwxxkmRamRmR1)()()()(nh)()()(mnwnxEmRxw0)()()()()(kknwkhnhnwnx )()()()()()(1pkkmnxknxamnxnwEmnxnxE 所以 )()()()()()()(00mnwknw

40、EkhmnwknwkhEmRkkxw)()()()()(2200mhkmkhkmRkhwwkkww00)(0)(2mmmhmRwxw0)(0)0()(0)()(121mmRmhkmRamkmRamRxxpkwxxkpkxxkxxpkxxkxwxxkmRamRmR1)()()( 顯然，AR模型輸出信號的自相關函數(shù)具有遞推的性質，即： pkkkzazH111)()()()(1nknhanhpkk1)0(h0)()(1mkmRamRpkxxkxx0)()(01mkmRamRpkxxkxx 由于由于 轉換到時域得到：轉換到時域得到：因而因而上式就是著名的上式就是著名的YuleWalker（Y-W）方程

41、，將上）方程，將上式變換：式變換：把上式的下標簡化并寫成矩陣的形式，可以寫成單一的正規(guī)矩陣方程： 0)()0(12mkRaRpkxxkxxw001) 0() 1()() 1() 0() 1 ()() 1() 0(21wpaaRpRpRpRRRpRRR0)()(01mkmRamRpkxxkxx 【例3-9】已知自回歸信號模型AR（3）為： 式中 是具有方差 =1的平穩(wěn)白噪聲，求 a. 自相關序列 ，m0，1，2，3，4，5。 b. 用a求出的自相關序列來估計AR（3）的參數(shù) ，以及輸入白噪聲的方差 大小。 c. 利用給出的AR模型，用計算機仿真給出32點觀測值 ,用觀測值的自相關序列直接來估計A

42、R(3)的參數(shù) 以及輸入白噪聲的 。)()3(241)2(249)1(2414)(nwnxnxnxnx)(nw2w)(mRxx ka 2w)(nx ka 2w醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法 32點觀測值 0.4282 1.1454 1.5597 1.8994 1.6854 2.3075 2.4679 1.9790 1.6063 1.2804 -0.2083 0.0577 0.0206 0.3572 1.6572 0.7488 1.6666 1.9830 2.6914 1.2521 1.8691 1.6855 0.6242

43、 0.1763 1.3490 0.6955 1.2941 1.0475 0.4319 0.0312 0.5802 -0.6177 解：a. 已知的是模型參數(shù) ， 14/24 9/24， 1/24，來求自相關序列 。 )(nx ka1a2a3a)(mRxx醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法把把 代入，利用自相關函數(shù)的偶對稱，代入，利用自相關函數(shù)的偶對稱，得到一個得到一個44的的矩陣：的的矩陣： 利用式（7-37），001)0() 1()() 1()0() 1 ()() 1()0(21wpaaRpRpRpRRRpRRR ka00

44、0124/124/924/141)0() 1 ()2() 3() 1 ()0() 1 ()2()2() 1 ()0() 1 () 3()2() 1 ()0(RRRRRRRRRRRRRRRR醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法解線性方程組得：R（0）= 4.9377 R（1）4.3287 R（2）4.1964 R（3）3.8654利用式0)()(1mkmRamRpkxxkxx6481. 3)4()4(31kxxkxxkRaR4027. 3)5()5(31kxxkxxkRaR可以求出可以求出R（4），），R（5） 當然還可以求出無

45、窮多的自相關序列值。當然還可以求出無窮多的自相關序列值。醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法b. 已知自相關序列值，來估計已知自相關序列值，來估計3階階AR模型的參數(shù)模型的參數(shù) 以及以及 ,利用式（利用式（7-37）得到矩陣：）得到矩陣： 14/24， 9/24， 1/24， 1 可以發(fā)現(xiàn)對AR模型參數(shù)是無失真的估計，因為已知AR模型，我們可以得到完全的輸出觀測值，因而求得的自相關函數(shù)沒有失真，當然也就可以不失真的估計。ka 2w0001)0() 1 ()2() 3() 1 ()0() 1 ()2()2() 1 ()0() 1

46、 () 3()2() 1 ()0(2321waaaRRRRRRRRRRRRRRRR1 a2 a3 a2w解線性方程組得到：解線性方程組得到：c. 利用給出的32點觀測值，先求自相關序列（按照上節(jié)的樣本自相關定義 計算） 由于偶對稱只給出m0，1，231的 1.9271 1.6618 1.5381 1.3545 1.1349 0.9060 0.8673 0.7520 0.7637 0.8058 0.8497 0.8761 0.9608 0.8859 0.7868 0.7445 0.6830 0.5808 0.5622 0.5134 0.4301 0.3998 0.3050 0.2550 0.19

47、97 0.1282 0.0637 0.0329 -0.0015 -0.0089 -0.0143 -0.0083nimiixxxxnmR11)()(mRxx醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法 把頭4個相關序列值代入矩陣（7-38）求得估計值： 0.6984 , 0.2748 0.0915， 0.4678 1 a2 a3 a2w與真實與真實AR模型參數(shù)誤差為：模型參數(shù)誤差為： e10.1151， e2 0.1002， e3 0.0498，原因在于我們只有，原因在于我們只有一部分的觀測數(shù)據(jù)，使得自相關序列值與理想一部分的觀測數(shù)據(jù)，

48、使得自相關序列值與理想的完全不同的完全不同 醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法式中式中 ，k1，2，m，代表，代表m階預測階預測器的預測系數(shù)。器的預測系數(shù)。 2)YW方程的解法L-D算法把把AR模型和預測系統(tǒng)聯(lián)系起來模型和預測系統(tǒng)聯(lián)系起來若序列的模型已知而用過去觀測的數(shù)據(jù)若序列的模型已知而用過去觀測的數(shù)據(jù)來推求現(xiàn)在和將來的數(shù)據(jù)稱為來推求現(xiàn)在和將來的數(shù)據(jù)稱為前向預測器，表示為：，表示為：mkmknxkanx1)()()( )(kam 假如假如mp，且預測系數(shù)和，且預測系數(shù)和AR模型參數(shù)相同，比模型參數(shù)相同，比較式較式(3-2)

49、和下式，把預測誤差系統(tǒng)框圖和和下式，把預測誤差系統(tǒng)框圖和AR模型模型框圖給出，如圖框圖給出，如圖3.2所示所示 (3-1) 把 e(n) 看成是系統(tǒng)的輸出，x(n) 看成是系統(tǒng)的輸入，得到系統(tǒng)函數(shù)：mkmknxkanxnxnxne1)()()()( )()(mkkmzkazXzE1)(1)()(顯然預測出來的結果與真實的結果存在預測誤差顯然預測出來的結果與真實的結果存在預測誤差或前向預測誤差，設誤差為或前向預測誤差，設誤差為 e(n) ：（3-2）pkkknxanwnx1)()()(pkkkzazH111)(如圖3.2所示，即有 ，即前向預測誤差系統(tǒng)中的輸入為 ，輸出為預測誤差 等于白噪聲。

50、也就是說前向預測誤差系統(tǒng)對觀測信號起了白化的作用。由于AR模型和前向預測誤差系統(tǒng)有著密切的關系，兩者的系統(tǒng)函數(shù)互為倒數(shù)，所以求AR模型參數(shù)就可以通過求預測誤差系統(tǒng)的預測系數(shù)來實現(xiàn)。 )(nw)(nx)(zH)(zW)(zX)(nh)(ne)(nx)(1zH)(zE)(zX)()(nenw)(nx)(ne圖圖3.2預測誤差系統(tǒng)和預測誤差系統(tǒng)和AR模型模型 要使得均方誤差最小，將上式右邊對預測系數(shù)求偏導并且等于零，得到m個等式：) )()()()(212mkmknxkanxEneEmkmlxxmmmkxxmxxklRkalakRkaR111)()()( )()( 2)0(mlklRkalRmkxx

51、mxx, 2 , 1)()()(1mkmknxkanxnxnxne1)()()()( )()(對上式求預測誤差均方值：對上式求預測誤差均方值：(3-3)或或 (3-5)也就是p階預測器的預測系數(shù)等于p階AR模型的參數(shù)，由于 ，所以最小均方預測誤差等于白噪聲方差，即mkxxmxxmkRkaRneE12)()()0()(pkxxkxxpkRaRneE12)()0()(pmkaamk)()(ne)(nw22)(wpneE求得最小均方誤差：求得最小均方誤差： (3-4) 估計AR模型參數(shù)，也即要估計參數(shù) 這里介紹應用廣泛的L-D算法。 L-D算法的基本思想就是根據(jù)Y-W方程式或式(3-3)、(3-4)

52、、(3-5)，自相關序列具有遞推的性質，2,wkpamlklRkalRmkxxmxx, 2 , 1)()()(1mkxxmxxmkRkaRneE12)()()0()(pkxxkxxpkRaRneE12)()0()(22)(wpneE醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法 L-D遞推算法是模型階數(shù)逐漸加大的一種算法， 先計算階次m1時的預測系數(shù) )(kam) 1 (1a21w)(kam) 1 (2a)2(2a22w) 1 (pa)(pap2wp然后計算然后計算m2時的時的一直計算到一直計算到mp階時的階時的這種遞推算法的特點是，每

53、一階次參數(shù)的計算是從這種遞推算法的特點是，每一階次參數(shù)的計算是從低一階次的模型參數(shù)推算出來的，既可減少工作低一階次的模型參數(shù)推算出來的，既可減少工作量又便于尋找最佳的階數(shù)值，滿足精度時就停止量又便于尋找最佳的階數(shù)值，滿足精度時就停止遞推。遞推。 例如按照式(3-3)、(3-5)，取m1，代入，簡化下標，則有： (3-7)2() 1 () 1 () 1 ()0()0() 1 () 1 (111RaRERaR)1 (1)0 (21121aREw)0()2() 1 () 1 ()2() 1 ()2()0() 1 () 1 (2222RaRaRRaRaR) 1 ()0() 1 ()1 ()2() 1

54、()2()0() 1 () 1 ()2()2() 1 ()2() 1 ()0() 1 ()2() 1 () 1 (1212212122RRaaaRRRaRaaaaRRaRa m=2：把（把（37）的（）的（1）代入上式得到：）代入上式得到： (3-8)醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法把(3-8)的（1）代人： (3-9) 2() 1 () 1 () 2() 1 () 1 () 0 () 1 () 1 () 2() 2(11112ERaRRaRRaRa)2()2() 1 () 1 ()0(222RaRaRE)0()2() 1

55、 () 1 ()2() 1 () 1 ()0(2222RaRaaRaR)0()2() 1 ()2(1)1 ()0(2222RaRaaR)0()2() 1 ()2(1)2(1)1 ()0(22221RaRaaaR122122)2(1 )1 () 1 ()0()2(1 EaRaRa根據(jù)根據(jù)(3-5)，估計的方差為：，估計的方差為：醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法 其中 稱為反射系數(shù)，從上式知道整個迭代過程需要已知自相關函數(shù)，給定初始值 以及AR模型的階數(shù) ，就可以按照圖3.7所示流程圖進行估計。 mkkmmwmmmmmmmmk

56、mmkaREmaEkmamakakaEkmRkamRma12122111111)(1 ) 0 ()(1)()()()()()()()()(mam) 0 (0RE 1)0(0ap這樣遞推下去可以得到預測系數(shù)和均方誤差估計的這樣遞推下去可以得到預測系數(shù)和均方誤差估計的通式：通式： (3-10)開始10),(Nnnx輸入)( mRxx計算10000aREmxx），（）：初始化（1m，計算按mmEma)()503(1.1)(mkkampmpkkap.1),(輸出結束是否1mm圖3.7【例3-10】已知自回歸信號模型AR（3）為： 式中 是具有方差 =1的平穩(wěn)白噪聲，利用給出的AR模型，用計算機仿真給出

57、32點觀測值 0.4282 1.1454 1.5597 1.8994 1.6854 2.3075 2.4679 1.9790 1.6063 1.2804 -0.2083 0.0577 0.0206 0.3572 1.6572 0.7488 1.6666 1.9830 2.6914 1.2521 1.8691 1.6855 0.6242 0.1763 1.3490 0.6955 1.2941 1.0475 0.4319 0.0312 0.5802 -0.6177，用L-D算法來估計AR(3)的參數(shù) 以及輸入白噪聲的方差 )()3(241)2(249) 1(2414)(nwnxnxnxnx)(nw2w)(nx ka 2w醫(yī)學信號處理醫(yī)學信號處理 第三章第三章 生物醫(yī)學信號的時域及頻域分析方法生物醫(yī)學信號的時域及頻域分析方法 解：步驟1 利用給出的32點觀測值，先求自相關序列,由于偶對稱，只給出m0，1，231的Rx(m) 1.9271 1.6618 1.5381 1.3545 1.1349 0.9060 0.8673 0.7520 0.7637 0.8058 0.8497 0.8761 0.9608 0.8859 0.7868 0.7445 0.6830