Yule-Walker方程

1、評(píng)分大理大學(xué)實(shí)驗(yàn)報(bào)告課程名稱(chēng)牛物醫(yī)學(xué)信號(hào)處理實(shí)驗(yàn)名稱(chēng)Yule-Walker方程專(zhuān)業(yè)班級(jí)姓名學(xué)號(hào)實(shí)驗(yàn)日期實(shí)驗(yàn)地點(diǎn)20152016學(xué)年度第3學(xué)期實(shí)驗(yàn)?zāi)康膶W(xué)習(xí)求解Yule-Walker方程，建立隨機(jī)信號(hào)的AR模型、實(shí)驗(yàn)環(huán)境1、硬件配置：處理器：AMDA10-5750MAPUwithRadeon(tm)Graphics2.50GHz安裝內(nèi)存：(RAM)4.00GB系統(tǒng)類(lèi)型：win10位操作系統(tǒng)，基于x64位處理器2、軟件環(huán)境：MatlabR2012b三、實(shí)驗(yàn)原理隨機(jī)信號(hào)可以看作是由當(dāng)前激勵(lì)白噪聲w(n)以及若干次以往信號(hào)x(n-k)的線性組合產(chǎn)生，即所謂自回歸模型(AR模型)piakxnkk11pka

2、kZk1模型參數(shù)滿(mǎn)足Yule-Walker方程RxxmpakRxxmk1矩陣形式R0R1R1R0a1a?apR1R2Rp求解Yule-Walker方程，就可以得到AR模型系數(shù)當(dāng)模型階次較大時(shí)，直接用矩陣運(yùn)算求解的計(jì)算量大，不利于實(shí)時(shí)運(yùn)算。利用系數(shù)矩陣的特性，人們提出了如L-D算法等快速算法。四、實(shí)驗(yàn)內(nèi)容編寫(xiě)求解Yule-Walker方程的程序，并對(duì)實(shí)際生理信號(hào)(例如心電、腦電)建立AR模型。對(duì)同一數(shù)據(jù)，使用Matlab信號(hào)處理工具箱自帶函數(shù)aryule計(jì)算相同階數(shù)AR模型系數(shù)，檢驗(yàn)程序是否正確。用偽隨機(jī)序列(白噪聲)驅(qū)動(dòng)AR模型，觀察輸出是否與真實(shí)心電、腦電信號(hào)相似,對(duì)比真實(shí)信號(hào)與仿真信號(hào)的功

3、率譜。五、實(shí)驗(yàn)結(jié)果與分析試驗(yàn)程序：clear;clc;M=1024;loadecgdata;loadeegdata;loadicpdata;loadrespdata;x=ecgdata(1:M);%長(zhǎng)度任意選擇%x=eegdata(1:M);%長(zhǎng)度任意選擇%x=icpdata(1:M);%長(zhǎng)度任意選擇%x=respdata(1:M);%長(zhǎng)度任意選擇i=1:15;Sw=zeros(1,length(i);E=zeros(1,length(i);FPE=zeros(1,length(i);forp=1:15%嘗試改變模型階數(shù)，觀察效果Rxx=xcorr(x,biased);%biased為有偏的互

4、相關(guān)函數(shù)估計(jì)Rtemp=zeros(1,p);Rl=zeros(p,1);fork=1:length(Rtemp)Rtemp(k)=Rxx(length(x)-1+k);Rl(k)=Rxx(length(x)+k);endRs=toeplitz(Rtemp);%生成自相關(guān)系數(shù)矩陣指矩陣中每條自左上至右下的斜線上的元素相同，toeplitz(x)用向量x生成一個(gè)對(duì)稱(chēng)的托普利茲矩陣A=-inv(Rs)*Rl;%AR模型系數(shù)估計(jì)白噪聲方差估計(jì)Sw(p)=Rtemp(1),Rl*1;A;%采用malab自帶函數(shù)估計(jì)模型系數(shù)a,E(p)=aryule(x,p);%a-系數(shù)，E-預(yù)測(cè)誤差，k-反射系數(shù)%a

5、,E(p)=arburg(x,p);%Y=filter(B,A,X),輸入X為濾波前序列，Y為濾波結(jié)果序列,%B/A提供濾波器系數(shù)，B為分子，A為分母da=a(2:end)-A;%自編程序求解是否正確?FPE(p)=E(p)*(M+p+1)/(M-p-1);w=randn(size(x);x2=filter(1,a,w);%仿真數(shù)據(jù)Rx2=xcorr(x2,biased);xx=abs(fft(Rxx);xx2=abs(fft(Rx2);真實(shí)數(shù)據(jù),);仿真數(shù)據(jù),);真實(shí)數(shù)據(jù)功率譜,);仿真數(shù)據(jù)功率譜,);figure;subplot(4,1,1);plot(x);xlim(01024);tit

6、le(subplot(4,1,2);plot(x2);xlim(01024);title(subplot(4,1,3);plot(-1023:1023,xx);title(subplot(4,1,4);plot(-1023:1023,xx2);title(error(p)=mean(x-x2).A2);%求最小均方誤差endPopt=find(FPE=min(FPE)figure,subplot(1,3,1);plot(i,error,-*)title(a)最小均方誤差隨階數(shù)p的變化情況),xlabel(p);ylabel(error);subplot(1,3,2);plot(i,E,-*);

7、gridontitle(b)預(yù)測(cè)誤差隨階數(shù)p的變化情況),xlabel(p);ylabel(E);%figure,stem(i,Sw,-*);gridon%title(白噪聲方差估計(jì)隨階數(shù)p的變化情況,)%xlabel(p);ylabel(白噪聲方差估計(jì),);subplot(1,3,3);plot(i,FPE,-*);title(c)FPE隨階數(shù)p的變化情況);xlabel(p);ylabel(FPE);實(shí)驗(yàn)結(jié)果:1、心電數(shù)據(jù)，階數(shù)p為15,M=1024:(1)最小均方誤差隨p的變化情況、預(yù)測(cè)誤差隨p的變化情況、FPE隨p的變化情況:roe2.51.50.5(a)最小均方誤差隨階數(shù)p白徊頻惻睇

8、確階數(shù)p的變化C)F喇03數(shù)p的變化情況3ppp圖1-1L-D算法圖1-2Burg算法(2)白噪聲方差估計(jì)隨階數(shù)p的變化情況:圖1-4Burg算法真實(shí)數(shù)據(jù)10-101002003004005006007008009001000仿真數(shù)據(jù)圖1-3L-D算法(3)真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)及其功率譜:真實(shí)數(shù)據(jù)111111111jI-0L_-1111111111L01002003004005006007008009001000仿真數(shù)據(jù)0100200300400500600700800900100050i真實(shí)數(shù)據(jù)功率譜2001000050010001500-1500-1000-5004204X10仿真數(shù)據(jù)功率譜

9、050010001500-1500-1000-500200100-1000100015001-1000-5001000-5011hiilj01002003004005006007008009001000真實(shí)數(shù)據(jù)功率譜0-1500x1020-1500-5000仿真數(shù)據(jù)功率譜1500圖1-5L-D算法error=19.234圖1-6Burg算法error=19.9132、腦電數(shù)據(jù)，階數(shù)p為14,M=1024(1)最小均方誤差隨p的變化情況、預(yù)測(cè)誤差隨p的變化情況、FPE隨p的變化情況:(a)最小均方誤差隨階數(shù)p睚顫蒯觸爰隨階數(shù)p的變彳愉聲PE隨階數(shù)p的變化情況(a)ft小均方誤差隨階數(shù)p盟b靂訥囑

10、耀隨階數(shù)p的變化(靖FPE隨階數(shù)p的變化情況圖2-1L-D算法圖2-2Burg算法(2)白噪聲方差估計(jì)隨階數(shù)p的變化情況:圖2-4Burg算法圖2-3L-D算法(3)真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)及其功率譜:真實(shí)數(shù)據(jù)100一，.J。|111L11J11L-l01002003004005006007008009001000真實(shí)數(shù)據(jù)10I.足.ib0-1011111111LJ01002003004005006007008009001000-1500-1000-500050010001500真實(shí)數(shù)據(jù)功率譜-1500-1000-500050010001500仿真數(shù)據(jù)功率譜真實(shí)數(shù)據(jù)功率譜20010005001000

11、1500-1500-1000-5000仿真數(shù)據(jù)50-5010020030040050060070080090010002001000100500仿真數(shù)據(jù)5f0:;小、,*心電一-5I1111fc1;J1101002003004005006007008009001000圖2-5L-D算法error=3.3147圖2-6Burg算法error=3.20433、顱內(nèi)壓數(shù)據(jù)，階數(shù)p為14,M=1024(1)最小均方誤差隨p的變化情況、預(yù)測(cè)誤差隨p的變化情況、FPE隨p的變化情況:45040035030025020015010050010ppppproe0.110.10.090.080.070.060.

12、050.0401020(a)最小均方誤差隨階數(shù)p艙蒯惻B蓬隨階數(shù)p的變化C)FP隨階數(shù)p的變化情況500(a)l小均方誤差隨階數(shù)p白值卿測(cè)I笳隨階數(shù)p的變化礴FP瞄階數(shù)p的變化情況0.12p圖3-1L-D算法圖3-2Burg(2)白噪聲方差估計(jì)隨階數(shù)p的變化情況:0.16白噪聲方差估計(jì)隨階數(shù)p的變化情況0.140.120.10.080.060.040.02計(jì)估差方聲噪白圖3-1L-D算法圖3-2Burg算法(3)真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)及其功率譜:真實(shí)數(shù)據(jù)105001002003004005006007008009001000仿真數(shù)據(jù)50011000-500100200300400500600700

13、x10真實(shí)數(shù)據(jù)功率譜290010000-1500-1000-5000500x104仿真數(shù)據(jù)功率譜1015005050010000-1500-1000-5001500真實(shí)數(shù)據(jù)105500800900100011000-500100200300400500600x1C4真實(shí)數(shù)據(jù)功率譜200100200300400500600仿真數(shù)據(jù)0-1500-1000-5000500x104仿真數(shù)據(jù)功率譜2800900100015001050010000-1500-1000-5001500圖3-3L-D算法error=52.4963圖3-4Burg算法error=103.77014、呼吸壓數(shù)據(jù)，階數(shù)p為7,M=

14、1024(1)最小均方誤差隨p的變化情況、預(yù)測(cè)誤差隨p的變化情況、FPE隨p的變化情況:(a)t小均在屬/隨階數(shù)pl徊戰(zhàn)測(cè)I笳隨階數(shù)p的變傕師E隨階數(shù)p的變化情況3.153.13.0532.952.92.852.82.752.72.650510ppp圖4-1L-D算法圖4-2Burg算法(2)白噪聲方差估計(jì)隨階數(shù)p的變化情況：圖4-3L-D算法圖4-4Burg算法(3)真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)及其功率譜:50真實(shí)數(shù)據(jù)真實(shí)數(shù)據(jù)1000500090010008009001000100200300400500600700仿真數(shù)據(jù)0100200300400500仿真數(shù)據(jù)-501500-1000-500010

15、001500-5000仿真數(shù)據(jù)功率譜x104100200300400500600700真實(shí)數(shù)據(jù)功率譜90010000-1500-10004X101010000-1500500215005-1000-5001000-500100200300400500600x1C8真實(shí)數(shù)據(jù)功率譜480090010000-1500-10004X1010-5000仿真數(shù)據(jù)功率譜10000-15001500圖4-5L-D算法2.7785e05圖4-6Burg算法2.17161e055、心電數(shù)據(jù)，階數(shù)p為5,M=512:(1)最小均方誤差隨p的變化情況、預(yù)測(cè)誤差隨p的變化情況、FPE隨p的變化情況:(a)最小均方誤差隨

16、階數(shù)p的頓I14020012010080E60405pppppp圖5-1L-D算法圖5-2Burg算法(2)白噪聲方差估計(jì)隨階數(shù)p的變化情況:圖5-1L-D算法圖5-2Burg算(3)真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)及其功率譜:真實(shí)數(shù)據(jù)10-1050100150200250300350400450500仿真數(shù)據(jù)真實(shí)數(shù)據(jù)0-1050100150200250300350400450500仿真數(shù)據(jù)500-50050100150200250300350400450500真實(shí)數(shù)據(jù)功率譜500-50050100150200250300350400450500真實(shí)數(shù)據(jù)功率譜50-400400600-2000200仿真數(shù)據(jù)

17、功率譜0-60050100000-6005000-400-20002004000-600-4004X104-2000200仿真數(shù)據(jù)功率譜4006002-400-20002004000-600600600圖5-5L-D算法error=21.1246圖5-6Burg算法error=11.8885腦電數(shù)據(jù)，階數(shù)p為6,M=512:(1)最小均方誤差隨p的變化情況、預(yù)測(cè)誤差隨p的變化情況、FPE隨p的變化情況:0.940.930.920.910.90.890.88pproe0.87.0510p0.940.93,0.92,09.0.90.890.88,pproe(a)l小均方誤差隨階數(shù)p的bEO解灌隨階

18、數(shù)p的變化(情FPE!階數(shù)p的變化情況0.95(a最小均方誤差隨階數(shù)p附I蓑棚般差隨階數(shù)p的變化CFPE隨階數(shù)p的變化情況0.950.87.0510p圖6-1L-D算法圖6-2Burg算法(2)白噪聲方差估計(jì)隨階數(shù)p的變化情況:圖6-3L-D算法圖6-4Burg算法(3)真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)及其功率譜:-100真實(shí)數(shù)據(jù)真實(shí)數(shù)據(jù)100450500504020人力-40020血.10-400-200400600-50150200250300350400450500仿真數(shù)據(jù)0-60001-600-200200150200250300350400450500真實(shí)數(shù)據(jù)功率譜仿真數(shù)據(jù)功率譜50-545050

19、04020-4002010-400-200200100150200250300350仿真數(shù)據(jù)011-600100150200250300350真實(shí)數(shù)據(jù)功率譜-2000200仿真數(shù)據(jù)功率譜0-600圖6-5L-D算法error=2.5564圖6-6Burg算法error=2.71037、顱內(nèi)壓數(shù)據(jù)，階數(shù)p為14,M=512(1)最小均方誤差隨p的變化情況、預(yù)測(cè)誤差隨p的變化情況、FPE隨p的變化情況:(a)最小均方誤差隨階數(shù)p的弱捌幅集隨階數(shù)p的變化(C)FPE隨階數(shù)p的變化情況350*300J-250|.ro200e150110050+蜂洲+0a400E01020p0.2L-*出*,0.190

20、.180.170.160.150.14EPF0.13010200.20.190.180.170.160.1501020p(a)t小均方誤差隨階數(shù)p馳闞則震差隨階數(shù)p的變化C鼾PE隨階數(shù)p的變化情況0.14pproe0.130.120.110.10.090.0801020pp7-1L-D圖7-2Burg(2)白噪聲方差估計(jì)隨階數(shù)p的變化情況:0.2白噪聲方差估計(jì)隨階數(shù)p的變化情況0.180.160.140.120.10.080.060.040.02計(jì)估差方聲噪白圖7-3L-D算法圖7-4Burg算法(3)真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)及其功率譜:真實(shí)數(shù)據(jù)真實(shí)數(shù)據(jù)200250300350400450500仿真

21、數(shù)據(jù)1050200-20150200250300350400450500仿真數(shù)據(jù)150200250300350400450500真實(shí)數(shù)據(jù)功率譜5000-400-200仿真數(shù)據(jù)功率譜0-600500-505000-200仿真數(shù)據(jù)功率譜200250300350400450500真實(shí)數(shù)據(jù)功率譜100000-6005000-400-2000-600-40051-400-20002004000-600X102600圖7-4L-D算法error=17.0841圖7-5Burg算法error=43.72808、呼吸壓數(shù)據(jù)，階數(shù)p為1,M=512(1)最小均方誤差隨p的變化情況、預(yù)測(cè)誤差隨p的變化情況、FPE隨p的變化情況:(a)最小均方嬖楙階數(shù)p眼喇測(cè)靛差隨階數(shù)p的變化C)FP隨階數(shù)p的變化情況2.93122.93122.93122.93122.93121p2.93120圖8-1L-D算法圖8-2Burg算法(2)白噪聲方差估計(jì)隨階數(shù)p的變化情況:圖8-4Burg算法圖8-3L-D算法(3)真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)及其功率譜:50真實(shí)數(shù)據(jù)真實(shí)數(shù)據(jù)10005000-500-40