Yule-Walker方程

1、生物醫(yī)學(xué)信號(hào)處理實(shí)驗(yàn)報(bào)告評(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模型。2、 實(shí)驗(yàn)環(huán)境 1、硬件配置： 處理器：AMD A10-5750M APU with Radeon(tm) Graphics 2.50GHz 安裝內(nèi)存：（RAM）4.00GB 系統(tǒng)類(lèi)型：win10位操作系統(tǒng)，基于x64位處理器 2、軟件環(huán)境：Matlab R2012b3、 實(shí)驗(yàn)原理隨機(jī)信號(hào)可以看作是由當(dāng)前激勵(lì)白噪聲w(n)以及若干次

2、以往信號(hào)x(n-k)的線性組合產(chǎn)生，即所謂自回歸模型（AR模型）模型參數(shù)滿(mǎn)足Yule-Walker方程矩陣形式求解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; load ecgdata;load eegdata;load icpdata;load respdata;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);for p=1:15 % 嘗試改變模型階數(shù)，觀察效果Rxx = xcorr(x

4、,'biased');%biased為有偏的互相關(guān)函數(shù)估計(jì)Rtemp = zeros(1,p);Rl = zeros(p,1);for k = 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ì)Sw(p)= Rtemp(1),Rl'*1;A; %

5、白噪聲方差估計(jì)% 采用malab自帶函數(shù)估計(jì)模型系數(shù)a,E(p) = aryule(x,p); % a-系數(shù)，E-預(yù)測(cè)誤差，k-反射系數(shù) %a,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 = ab

6、s(fft(Rxx);xx2 = abs(fft(Rx2);figure;subplot(4,1,1);plot(x); xlim(0 1024);title('真實(shí)數(shù)據(jù)');subplot(4,1,2);plot(x2);xlim(0 1024);title('仿真數(shù)據(jù)');subplot(4,1,3); plot(-1023:1023,xx);title('真實(shí)數(shù)據(jù)功率譜'); subplot(4,1,4); plot(-1023:1023,xx2);title('仿真數(shù)據(jù)功率譜');error(p)=mean(x-x2).2

7、);%求最小均方誤差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,'-*');grid ontitle('(b)預(yù)測(cè)誤差隨階數(shù)p的變化情況'),xlabel('p');ylabel('E');%figure,stem(i,Sw,

8、'-*');grid on%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的變化情況： 圖1-1 L-D算法 圖1-

9、2 Burg算法（2）白噪聲方差估計(jì)隨階數(shù)p的變化情況： 圖1-3 L-D算法 圖1-4 Burg算法（3）真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)及其功率譜： 圖1-5 L-D算法error=19.234 圖1-6 Burg算法 error=19.9132、腦電數(shù)據(jù)，階數(shù)p為14，M=1024：（1）最小均方誤差隨p的變化情況、預(yù)測(cè)誤差隨p的變化情況、FPE隨p的變化情況： 圖2-1 L-D算法 圖2-2 Burg算法（2） 白噪聲方差估計(jì)隨階數(shù)p的變化情況：圖 2-3 L-D算法 圖2-4 Burg算法（3）真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)及其功率譜： 圖2-5 L-D算法 error=3.3147 圖2-6 Burg算法

10、error=3.20433、顱內(nèi)壓數(shù)據(jù)，階數(shù)p為14，M=1024：(1)最小均方誤差隨p的變化情況、預(yù)測(cè)誤差隨p的變化情況、FPE隨p的變化情況： 圖3-1 L-D算法 圖3-2 Burg算法(2) 白噪聲方差估計(jì)隨階數(shù)p的變化情況： 圖3-1 L-D算法 圖3-2 Burg算法(3) 真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)及其功率譜： 圖3-3 L-D算法error=52.4963 圖3-4 Burg算法 error=103.77014、 呼吸壓數(shù)據(jù)，階數(shù)p為7，M=1024：(1)最小均方誤差隨p的變化情況、預(yù)測(cè)誤差隨p的變化情況、FPE隨p的變化情況：圖4-1 L-D算法 圖4-2 Burg算法(2)白噪

11、聲方差估計(jì)隨階數(shù)p的變化情況：圖4-3 L-D算法 圖4-4 Burg算法（3）真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)及其功率譜： 圖4-5 L-D算法 2.7785e05 圖4-6 Burg算法2.17161e05 5、心電數(shù)據(jù)，階數(shù)p為5,M=512: （1）最小均方誤差隨p的變化情況、預(yù)測(cè)誤差隨p的變化情況、FPE隨p的變化情況：圖5-1 L-D算法 圖5-2 Burg算法（2）白噪聲方差估計(jì)隨階數(shù)p的變化情況：圖5-1 L-D算法 圖5-2 Burg算（3）真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)及其功率譜：圖5-5 L-D算法 error=21.1246 圖5-6 Burg算法error=11.8885 6、腦電數(shù)據(jù)，階數(shù)p

12、為6,M=512: （1）最小均方誤差隨p的變化情況、預(yù)測(cè)誤差隨p的變化情況、FPE隨p的變化情況：圖6-1 L-D算法 圖6-2 Burg算法（2）白噪聲方差估計(jì)隨階數(shù)p的變化情況：圖6-3 L-D算法 圖6-4 Burg算法（3）真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)及其功率譜： 圖6-5 L-D算法 error=2.5564 圖6-6 Burg算法 error=2.71037、顱內(nèi)壓數(shù)據(jù)，階數(shù)p為14，M=512：(1)最小均方誤差隨p的變化情況、預(yù)測(cè)誤差隨p的變化情況、FPE隨p的變化情況：圖7-1 L-D算法 圖7-2 Burg算法（2） 白噪聲方差估計(jì)隨階數(shù)p的變化情況：圖7-3 L-D算法 圖7-4

13、 Burg算法（3）真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)及其功率譜： 圖7-4 L-D算法 error=17.0841 圖7-5 Burg算法error=43.72808、呼吸壓數(shù)據(jù)，階數(shù)p為1，M=512：(1)最小均方誤差隨p的變化情況、預(yù)測(cè)誤差隨p的變化情況、FPE隨p的變化情況：圖8-1 L-D算法 圖8-2 Burg算法（2）白噪聲方差估計(jì)隨階數(shù)p的變化情況：圖8-3 L-D算法 圖8-4 Burg算法（3）真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)及其功率譜： 圖8-5 L-D算法 error=2.9724e05 圖8-6 Burg算法 error=2.9917e05結(jié)果分析：1、 對(duì)信號(hào)進(jìn)行功率譜估計(jì)是，腦電信號(hào)的縱坐標(biāo)相差不大，而心電信號(hào)、顱內(nèi)壓信號(hào)、呼吸信號(hào)的縱坐標(biāo)相差較大。2、 當(dāng)AR模型結(jié)束p確定時(shí)，L-D算法的均方誤差比Brug算法的均方誤差要小。六、實(shí)驗(yàn)小結(jié)：1、L-D遞推算法是模型階數(shù)逐漸加大的一種算法，為了克服其誤差，提出了Burg算法，其基本思想是對(duì)觀測(cè)數(shù)據(jù)進(jìn)行向前或向后預(yù)測(cè)，然后讓兩者的均方誤差之和最小做為估計(jì)準(zhǔn)則估計(jì)反射系數(shù)進(jìn)而根據(jù)L-D算法的遞推公式求出AR模型參數(shù)。2、Burg算法的優(yōu)點(diǎn)是：求得出