Yule-Walker方程實驗報告

1、評分大理大學實驗報告課程名稱生物醫(yī)學信號處理實驗名稱Yu le W al ke r方程專業(yè)班級姓 名羽卒蘭cl學號實驗日期2016年5月27日實驗地點2015- 2 01 6學年度第學期、實驗目得學習求解Y ul e -Walke r方程,建立隨機信號得A R模型。二、實驗環(huán)境1、 硬件配置:In tel（R） Core（TM）i5 -42 1 0U CPU 1、7 G Hz 1、7GHz 安裝內存（RAM） : 4、00G B系統(tǒng)類型:64位操作系統(tǒng)2、軟件環(huán)境：M ATLAB R2013b軟件三、實驗內容（包括本實驗要完成得實驗問題及需要得相關知識簡單概述）編寫求解Yul e Walke

2、r方程得程序，并對實際生理信號（例如心電、腦電）建立A R模型。對同一數(shù)據(jù)，使用M atla b信號處理工具箱自帶函數(shù)aryul e計算相同階數(shù) AR模型系數(shù)，檢驗程序就是否正確用偽隨機序列（白噪聲）驅動AR模型，觀察輸出就是否與真實心電、腦電信號相似，對比真實信號與仿真信號得功率譜。四、實驗結果與分析（包括實驗原理、數(shù)據(jù)得準備、運行過程分析、源程序（代碼）、圖形圖象界面等）實驗原理隨機信號可以瞧作就是由當前激勵白噪聲w（n ）以及若干次以往信號 x （nk）得線性組合產(chǎn)生，即所謂自回歸模型（AR模型）模型參數(shù)滿足 Yule W alker方程矩陣形式求解Yul e Walke r方程,就可以

3、得到A R模型系數(shù) 當模型階次較大時，直接用矩陣運算求解得計算量大， 人們提出了如LD算法等快速算法。源程序：c lear ; cl c;M =1024;loa d ecgda t a; x不利于實時運算禾U用系數(shù)矩陣得特性%1 oa d ee g d a t a; x% load ic p d a ta;x = i c%loa d r e spd atp=1:60;Sw =zeros(1, 1 e n得Sw零矩陣或數(shù)組E=ze ro s (1 ,l e零矩陣或數(shù)組FPE=zeros(1 , l e ngth (p);矩陣或數(shù)組for i=1: 6 0=ec g data( 1:1 02 4

4、);=ee g dat a (1:M ); pdata (1:10 2 4);r espd a ta (1:102a; xg th(p)ngth(p )；%導入心電信號得數(shù)據(jù)%導入腦電信號得數(shù)據(jù)%導入顱內壓信號得數(shù)據(jù)4) ; %導入個呼吸信號得數(shù)據(jù) %P得取值范圍%創(chuàng)建一個1行列長為le%創(chuàng)建一個1行列長為len%創(chuàng)建一個1行列長為1 en g thngth ( p )%嘗試改變模型階數(shù)，觀察效果gth(p )得 E(p)得FPE零Rx x = x c orr (x，/ b ias 關函數(shù)估計；R temp = z eros(1, i );R1 = ze ros(i , 1);組for k =

5、 1:1 engt h( Rtemp)Rt emp( k ) = Rxx (length(x) 1+k);Rl( k) = Rxx (lengt h (x)+k) ;%e n dRs = toep 1 it z (Rtem p );A = inv(R s ) * Rl;Sw( i) =: Rtemp (1), R l 1; A;% 采用malab自帶函數(shù)估計模型系數(shù)a,E (i ) = aryule(x , i);估計,ed );%估計隨機過程中得互相關序列,參數(shù)估計d a = a(2 : e n d) A;FP E (i ) =E( i )* ( M+ i+1 ) / ( M-i 1w =

6、r andn(si z e (x);x 2= fi 1t er (1, a ,w);%a, E (i)= ar bu rg (x,i );en d” biased為有偏得互相%創(chuàng)建一個i行1列得Rtemp零矩陣或數(shù)組 %創(chuàng)建一個i行1列得Rl零矩陣或數(shù)%f o r 循環(huán)從 1 至U length( R temp)%取 Rxx (0)至U Rxx(p-1 )賦值取 R xx( 1)至U Rx x (p)賦值%生成自相關系數(shù)矩陣(Toeplitz型)%AR模型系數(shù)估計% 白噪聲方差估計% malab實現(xiàn)L -D算法得AR模型參數(shù)a-系數(shù)，E-預測誤差，k -反射系數(shù)%malab實現(xiàn)Bur g算法得

7、A R模型%自編程序求解就是否正確？% FPE算法%生成隨機噪聲%仿真數(shù)據(jù)fi g u re%畫圖sub p l o t(3, 1, 1),plot (p/ ,E, * ) , gr i d on;%創(chuàng)建窗口畫圖，并添加網(wǎng)格線title( / E 隨階數(shù) p 得變化情況);xlabe 1 ( p ) ; y l abe l ( e rr or ) ; %添加標題與 橫縱坐標s u bp 1 ot(3,1 , 2) , p lot( p ,Sw, / -*), g r i d on;%創(chuàng)建窗口畫圖，并添加網(wǎng)格線ti tle( / Sw隨階數(shù)p得變化情況)；xlab el ( pz) ;ylab

8、e l (白噪聲方差估計);subplo t ( 3 , 1, 3),plot(p ,F P E, ) , g rid on；%創(chuàng)建窗口畫圖,并添加網(wǎng)格線t it le (FPE隨階數(shù)p得變化情況);xlabe 1 (/ p); ylabe 1( F EPZ );%添加標題與橫 縱坐標心電信號：x 10-3E隨階數(shù)p的變化情況0010503pSw隨階數(shù)p的變化情況X 10-3210210x 10-3pFPE 隨階數(shù)p的變化情況p-332a) L D算法0仃 4- + 4- 4- 4 4- -14 4- 4- * * *1=牛4*4*4*x 10pSw隨階數(shù)p的變化情況10102030pFPE

9、隨階數(shù)p的變化情況4050603*1-3x 10210 10203040p50060b)Burg 算法 圖1心電信號得不同算法得E、Sw FPE隨階數(shù)p得變化圖分析：首先,根據(jù)FPE準則找到最佳階數(shù)P;LD算法；在 ma nd wi ndow 輸入 fi n d (FPE= mi n ( FP E)，結果：ans= 1 5;在 mand wi ndo w 輸入f i nd (abs ( E1)= = mi n(abs(E 1),結果：a ns=1; B u rg算法：在 mand wi n d o w 輸入 f ind(F P E= =mi n (FPE)，結果：an s =31;在m and

10、 wi n dow 輸入 fi n d (abs (E 1)=min( a bs (E -1),結果：ans=1。腦電信號0.9W- -a1(E隨階數(shù)p的變化情況PSw隨階數(shù)p的變化情況0.9h.一 一舟PFPE 隨階數(shù)p的變化情況1.4a) LD算法1.6 *+ + *p * * * *r* H-4- + *斗P*+- +*+*卜*蒼* * * *秸廿*E隨階數(shù)p的變化情況況 情 化 變 的 數(shù) 階 隨計肇、書噪白況 情 化 變 的 數(shù) 階 隨4- 4-1-亠士 rr * .丄 J, JL 厶 a A20b) B u r g算法圖2腦電信號得不同算法得E、Sw FPE隨階數(shù)p得變化圖分析：首

11、先,根據(jù)F PE準則找到最佳階數(shù) P；L D算法：在 mand wi nd o w 輸入f in d(F PE = = min( FP E),結果：a n s =2 3;在m and w i nd ow 輸入 fin d(a b s (E-1) = min( a b s (E- 1 ),結果：a ns= 1 1; Burg算法:在 ma nd wi n d o w 輸入 fin d(FPE= = mi n(FPE),結果:a n s=4 3 ;在 man d wi ndow 輸入 fi n d(ab s (E 1) =min (abs(E 1)，結果：an s=11. 顱內壓信號：-計#近囁白

12、 a) L- D算法0.15P 0.10.050Pb) Bu r g算法圖3顱內壓信號得不同算法得E、Sw、FPE隨階數(shù)p得變化圖分析：首先，根據(jù)FP E準則找到最佳階數(shù)P；L-D算法：在 ma nd w in dow 輸入 f i n d(F PE= =mi n( F PE),結果：an s=43;在m and w in dow 輸入 f i n d (abs(E-1 )= = mi n (abs( E 1),結果：ans=1;Bur g算法：在 ma nd wi n dow 輸入 fin d(FPE= = mi n(FPE),結果：a ns=57;在 mand windo w 輸入 f i

13、nd ( abs(E- 1)=m i n(a bs(E 1),結果:ans =1呼吸信號：0102030405060Pb) Burg算法圖4呼吸信號得不同算法得E、Sw FPE隨階數(shù)p得變化圖分析：首先，根據(jù)FPE準則找到最佳階數(shù) P;L -D算法：在 ma nd wi ndow輸入 f i n d(FP E =mi n( FP E)，結果：ans =59;在 ma nd w i ndow 輸入 find(ab s( E-1 ) =mi n (abs(E- 1 ),結果：ans = 60;Bu r g算法：在 mand wi n dow 輸入 find(FP E =min (F PE),結果：

14、a n s= 44;在m and wi ndow 輸入 fi nd (abs (E 1) = =mi n ( abs ( E 1),結果：an s =60。 分析：L- D算法與E urg算法得用find (abs (E -1 ) =mi n(abs( E -1 )找到得值就是一樣 得；心電、腦電、顱內壓信號得找find(F P E=min(FPE)得值Burg算法比LD算法大，呼吸信號找到得 mi n(F PE)值，B u r g算法比L D算法小。思考題：c lear; c1 c;1 oad ec gdata;x =ecgdata (1 :1024 );%導入心電信號%lo ade egd

15、ata ; x= eegdata(1:102 4 );%導入腦電信號%1 oadicpdata; x =icpd a t a(1:1024 );%導入顱內壓信號%1 oad res p da t a;x=respd ata (1:1024 );%導入個呼吸信號%p =20;%嘗試改變模型階數(shù)，觀察效果f o r p=2:2:20%嘗試改變模型階數(shù)，觀察效果R xx = xcorr(x,biased ;%估計隨機過程中得互相關序列，” bi a s e d為有偏得互相關函數(shù)估計；Rt emp = zeros (1, p);%創(chuàng)建一個i行1列得Rtem p零矩陣或數(shù)組R1 = zero s( p,

16、 1);%創(chuàng)建一個i行1列得Rl零矩陣或數(shù)組fork = 1: 1e ng t h(Rtem p )%f o r 循環(huán)從 1 至U l e ngth( R temp )R temp (k) = Rxx( le n g t h(x) 1+k);%取 Rxx (0 )至 Rxx(p 1)賦值R1 (k )= Rxx (l en g t h(x ) +k);%取 Rxx( 1)至U Rx x(p)賦值endRs = toepli t z(Rtemp);%生成自相關系數(shù)矩陣(Toeplitz 型A = -inv ( Rs)*R l ;% AR模型系數(shù)估計Sw( p/2 )= Rtemp(1) , R

17、l * 1;A ;%白噪聲方差估計%米用malab自帶函數(shù)估計模型系數(shù)a,E = ary u le(x,p);%a 系數(shù)，E預測誤差，k 反射系數(shù)%a,E = arburg( x , p );% malab實現(xiàn)Bu r g算法得AR模型參數(shù)估計da = a (2:end) A%自編程序求解就是否正確 ？stem (d a) ; ti t le(參數(shù)估計偏差)%畫圖w = r an d n (s ize(x );%產(chǎn)生隨機噪聲信號x2 = fi lt e r (1,a , w);%仿真數(shù)據(jù)fig u re ;subp 1 ot (2,1,1);plo t( x) ;tit 1 e (真實數(shù)據(jù));

18、%繪制真實數(shù)據(jù)得圖像subplot (2,1 , 2);pl o t (x2) ;t i t le(仿真數(shù)據(jù));繪制仿真數(shù)據(jù)得圖像error (p/ 2) =m ean (x -x2)、A 2 ); end%顯示最小均方誤差得計算結果f ig u r esub plot (1,2 , 1) , pl o t(2:2 : 2 0, e rror ti t I eC最小均方誤差隨階數(shù)p得變化情況-* )%畫圖);x label(p ) ;y lab el(error );%繪制最小均方誤差隨階數(shù)p得變化情況圖subplot (1, 2,2 ),stem(2:2: 20,Sw, * ) , grid

19、 o nti tl e (白噪聲方差估計隨階數(shù)p得變化情況);x la b el(p ); y label(白噪聲方差Rxx 2=x co r r(x2, / b iased );Px = a bs( f ft(Rxx);Px2 = a bs(ff fi gur e subplo t (2,1 火柴桿圖 subp l ot(2, 桿圖%繪制S w、E隨p變化散點圖，以下就是統(tǒng)計得到得結果Sw= 1、0527 1、0 301 1、015098 7 2;E= 3、874 14 176:; figu r e s u bp l ot(2,1, 1);st e m( p, Sw);xlabe l ( p

20、 );ylab el ( Sw );ti su bp lot (2 ,1 , 2 );pl ot (p, E, xlabel( p );ylt(Rx x2);9繪制白噪聲方差估計隨階數(shù)p得變化情況圖%求仿真數(shù)據(jù)得自相關函數(shù)計算真實信號自功率譜%計算仿真信號自功率譜,1) ; s t em(1023:1 023, Px ) ; titl1 ,2 ) ; s tem ( 1 0 23:1023,Px2 ) ; t i3、4 3 37 3、32000、 9 9610、3、7 3 143、e （真實信號功率譜）；繪制le （仿真信號功率譜）；p= 8 9 1 01 1 12 13 14994 3 0、

21、99 4 2 0、 9 896 0、%繪制火柴151 6;9 877 0、3683 3、4809 3、5826 3、62 5 33、tl e(Sw隨p變化散點圖 /o)；a bel(E ) ; t i t le( E 隨 p 變化散點圖);1、比較四種信號得參數(shù)估計偏差10-5-1001(真實數(shù)據(jù)參數(shù)估十偏差02468101214161820-16x 10 1520040060080010001200真實數(shù)據(jù)%繪制火柴桿圖）；%添加標題與橫縱坐標%繪制散點圖%添加標題與橫縱坐標真實數(shù)據(jù)x 101參數(shù)估計偏差, 一一 丨 l L , r II.J 11-11 1x 1011參數(shù)估計偏差02468

22、101214161820h 11II102468101214161820估計);a）心電b）腦電c）顱內壓d ）呼吸圖5心電、腦電、顱內壓、呼吸信號得得真實數(shù)據(jù)與參數(shù)估計偏差分析：由圖5可以瞧出腦電信號得參數(shù)估計偏差比其它三個信號得參數(shù)估計偏差要小許多，腦電信號得自編程序跟MATLAB信號處理工具箱自帶函數(shù)ary u le得處理更接近。由此可知丄-D算法與自編程序相比較，自編程序對估計得參數(shù)比較精確一些 2、比較四種信號得真實數(shù)據(jù)與仿真數(shù)據(jù)真實數(shù)據(jù)0 撫加11 “-5-100屮嚴八帥啊壯也仿真數(shù)據(jù)1050-5a）心電105真實數(shù)據(jù)真實數(shù)據(jù))0200400600HoIlli, Pill-1. i

23、f：|H 贋曲曲刪200400600800100000c） 顱內壓真實數(shù)據(jù)d） 呼吸b）腦電圖6心電、腦電、顱內壓、呼吸信號得得真實數(shù)據(jù)與仿真數(shù)據(jù)分析：從圖6可以瞧出，對心電、腦電、顱內壓、呼吸信號建模后，心電、顱內壓、呼吸信號得真實數(shù)據(jù)與仿真數(shù)據(jù)相差很大，腦電模型產(chǎn)生得信號更能真實得反映腦電信號得特征。不同噪聲得激勵得到得信號時不同得。3、比較四種信號L- DL D算法：真實信號功率譜算法與Bur g 算法得功率譜150真實信號功率譜150100仿真言號功率譜40真實信號功率譜a）心電b）腦電c）顱內壓d）呼吸B ur g算法真實信號功率譜仿真信號功率譜真實信號功率譜8X 10仿真信號功率譜

24、a）心電b） 腦電c） 顱內壓圖7心電、腦電、顱內壓、呼吸信號得LD算法與B u r g分析：由圖7瞧出，在比較心電、腦電、顱內壓、呼吸信號得L-D 算法與d）呼吸算法功率譜Bu r g算法功率譜發(fā)現(xiàn),對信號進行功率譜估計時，腦電與顱內壓信號得功率譜得縱坐標相差不大，但就是相對而說心電與呼吸信號得相差就是極大得。4、比較四種信號得L- D 算法與B u rg算法L D算法240最小均方誤差隨階數(shù)|的變化情況220200180160140120100p1X，0白噪聲方差估計隨階數(shù)p的變化情況0.90.80.70.6計估 差方5聲噪白0.40.30.20.10p最小均萬誤差 隨階數(shù)p的變化情況白

25、噪聲萬 差估計 隨階數(shù)P的變化情況1.48 6計估差方聲噪白rpea）心電b ）腦電況rDFe況計善方士陽白況 情c）顱內壓圖8四種信號LD算法，階數(shù)Burg 算法d ）呼吸p與均方誤差erro r與噪聲方差估計值Sw之間得關系最小均方誤差隨階數(shù)p的變化450戀差萬聲噪白7 6 5 a a a噪聲10方差估計隨階數(shù)p的1變化情況10p20計估差方聲噪白況數(shù) 階 隨 差 誤 方 均 最況 情 化 變 的數(shù) 階 隨 計 估 差 方 聲 噪齊白況 情 化 變 的最小X均0方誤差隨階數(shù)p的變化情況rare58275227淸白80噪聲方差估計隨階數(shù)p的變化情況16001400啦1000眈60040020

26、020c）顱內壓d）呼吸圖9四種信號Bur g算法，階數(shù)p與均方誤差er r o r與噪聲方差估計值 Sw之間得關系分析：由圖8與圖9可以瞧出，雖然四種信號得L D算法與 B urg算法階數(shù)p與均方誤差error得變化就是不穩(wěn)定得，但就是它們得階數(shù)p與噪聲方差估計值 Sw之間得變化就是大體一致得。4、1心電信號得L-D 算法與 Bu r g 算法比較真實數(shù)據(jù)真實數(shù)據(jù)仿真數(shù)據(jù)b)Bu r g算法a)L D 算法 b)B urg算法真實數(shù)據(jù)真實數(shù)據(jù)仿真數(shù)據(jù)/ 1 *1 I宀rJI,-I1 f1 riiififif1618a)L D算法Da值L-D-8 、-5 、-1 、-8 、-2 、4、2 、3

27、、5、3、算法990e-223e-340e-120e-144e-546e-100e-596e-641e-163e-13121112121211111111-1 、3、3、3、3、3、1 、-5、-3 、3、130e-8190e5763e2489e8141e1735e5867e012e-589e-7491e11-12-11-11-11-11-111212-14Burg算法-0 、2810、499-0 、092-0 、2750、1030、0610、026-0 、0970、0530、027-0 、-0 、0 、0 、-0 、0 、0 、-0 、0 、-0 、043501710117106914150

28、4950247057306700284圖10心電信號不同算法得真實數(shù)據(jù)與參數(shù)估計圖11心電信號不同算法得真實與仿真數(shù)據(jù)表1心電信號得L D算法與Burg 算法得D a值變化表4、2腦電信號得LD算法與B urg算法比較真實數(shù)據(jù)1050-5-10真實數(shù)據(jù)1(210，-1-21 /I i J,11 1 1 1B urg算法020040060080010001200x 10-3參數(shù)估計偏差1(0200400600 800 1000 1200020040060080010001200仿真數(shù)據(jù)5(汕神斛-11IIII020040060080010001200仿真數(shù)據(jù)真實數(shù)據(jù)a)L-D算法 法b)a)L

29、D算法b ) Bur g 算圖1 2腦電信號不同算法得真實數(shù)據(jù)與參數(shù)估計圖1 3腦電信號不同算法得真實與仿真數(shù)據(jù)表2腦電信號得L -D算法與B u rg算法得Da值變化表Da值L-D5、1 、-6 、0-1 、-1 、9、4、-8 、0算法5511e3877e939e-110e-665e-7144e8572e3266-16-17181616-17-1717-8 、2 、-2 、1 、-2 、1 、-1 、-2 、1 、5、326e-1510e082e-3877e168e-2490e145e-081e-1102e5511e17-1617-1716-161617-16-17Burg0 、-0 、-

30、0 、-0 、-0 、0 、0 、-0 、-0 、-0 、算法00090003000600040014000400130003000600050 、-0 、0 、0 、0 、0 、0 、-0 、-0 、0 、00050017000300090005001500090001001700024、3顱內壓信號得L D算法與 B ur g算法比較真實數(shù)據(jù)600o1 1 1 1 1 0 _ III, 1 . Il1 1 |2x 10-11參數(shù)估計偏差1真實數(shù)據(jù)6真實數(shù)據(jù)，-iliH： 1r.|：,II f ； i：:?：. ；i； - K2004T0600 800 1000200 0仿真數(shù)據(jù)20仿真數(shù)據(jù)

31、j .Mrr1i |-2-4200200 .00 600 800 “00 200a)L D算法b) Burg算法a)L D算法b) Burg算Da值L-D3、1 、1 、4、-1 、1 、4、3、-3 、4、算法1286e4018e1497e9960e567e-4885e3890e5907e099e-4496e-13-12-12-1512-13-12-1212-128 、1 、-7 、-6 、3、1 、-1 、6 、4、2 、5008e0558e490e-810e-4846e3522e628e-4571e0563e9531e-12-121212-12-1211-12-12-14Burg-0 、

32、0 、0 、0 、-0 、-0 、-0 、-0 、0 、0 、算法16210853105703710042080301450308002202500 、0 、-0 、0 、0 、0 、0 、-0 、0 、-0 、0356258239560479076703840087049401690025法圖14顱內壓信號不同算法得真實數(shù)據(jù)與參數(shù)估計圖15顱內壓信號不同算法真實與仿真數(shù)據(jù)表3顱內壓信號得L D算法與Bur g 算法得Da值變化表4、4呼吸信號得L-D 算法與Bur g算法比較100050000A h A1 l 1 V I X hn- y yV i真實數(shù)據(jù)真實數(shù)據(jù)0111 i 11iJx 10

33、-11參數(shù)估計偏差402468101214161820x 10-3參數(shù)估計偏差真實數(shù)據(jù)真實數(shù)據(jù)Da值L-D-2 、3、8 、4、-1 、-7、-2 、-9 、-1 、-1 、算法735e-042e-882e-431e-890e-866e-557e-929e-172e-041e-131312121213121311112 、8 、8 、1 、4、1 、6 、-5 、-2 、-6 、6015e3792e1914e0824e7703e5861e6561e402e-922e-889e-11-12-12-11-12-11-12121214Burg-0 、-0 、-0 、0 、0 、-0 、0 、0 、0

34、 、0 、算法01170006000800010006000600110012000600130 、0 、0 、0 、0 、0 、-0 、0 、-0 、0 、0019000700160005001100090007000700040007a)L D算法b )Burg算法圖16呼吸信號不同算法得真實數(shù)據(jù)與參數(shù)估計a) L - D算法b)B u rg算法圖17呼吸信號不同算法得真實與仿真數(shù)據(jù)表4呼吸信號得L-D算法與Burg算法得Da值變化表分析:Burg算法得優(yōu)點就是：求得得AR模型就是穩(wěn)定得，較高得計算效率，但遞推還就是用得LD算法，因此仍然存在明顯得缺點5、腦電信號,改變P值，比較L D算法與Bu rg算法Sw 隨p變化散點圖P1.5Sw 隨 p變化散點圖0EpE隨 p變化散點 圖Pa) L- D算法圖1 8腦電信號改變P值，比較b) Bu rg算法 L D算法與Burg算法Sw與E值變化圖P值8910111213141516L-D算法Sw1、05271、03011、0150099610、99430、99420、98960、98770、9872E3、46873、20103、8530338883、28743、39453、93154、15733、3942