Yule-Walker方程-實驗報告.doc

生物醫(yī)學(xué)信號處理評分大理大學(xué)實驗報告 課程名稱 生物醫(yī)學(xué)信號處理 實驗名稱 Yule-Walker方程 專業(yè)班級 姓 名 羽卒蘭cl 學(xué) 號 實驗日期 2016年5月27日 實驗地點 20152016學(xué)年度第 3 學(xué)期一、 實驗?zāi)康膶W(xué)習(xí)求解Yule-Walker方程，建立隨機信號的AR模型。二、實驗環(huán)境 1、硬件配置：Intel（R） Core(TM) i5-4210U CPU 1.7GHz 1.7GHz 安裝內(nèi)存（RAM）：4.00GB 系統(tǒng)類型：64位操作系統(tǒng) 2、軟件環(huán)境：MATLAB R2013b軟件三、實驗內(nèi)容(包括本實驗要完成的實驗問題及需要的相關(guān)知識簡單概述)編寫求解Yule-Walker方程的程序，并對實際生理信號（例如心電、腦電）建立AR模型。對同一數(shù)據(jù)，使用Matlab信號處理工具箱自帶函數(shù)aryule計算相同階數(shù)AR模型系數(shù)，檢驗程序是否正確。用偽隨機序列（白噪聲）驅(qū)動AR模型，觀察輸出是否與真實心電、腦電信號相似，對比真實信號與仿真信號的功率譜。四、實驗結(jié)果與分析 （包括實驗原理、數(shù)據(jù)的準備、運行過程分析、源程序（代碼）、圖形圖象界面等） 實驗原理隨機信號可以看作是由當(dāng)前激勵白噪聲w(n)以及若干次以往信號x(n-k)的線性組合產(chǎn)生，即所謂自回歸模型（AR模型）模型參數(shù)滿足Yule-Walker方程矩陣形式求解Yule-Walker方程，就可以得到AR模型系數(shù)當(dāng)模型階次較大時，直接用矩陣運算求解的計算量大，不利于實時運算。利用系數(shù)矩陣的特性，人們提出了如L-D算法等快速算法。 源程序：clear; clc;M=1024;load ecgdata;x = ecgdata(1:1024); %導(dǎo)入心電信號的數(shù)據(jù)%load eegdata;x = eegdata(1:M); %導(dǎo)入腦電信號的數(shù)據(jù)%load icpdata;x = icpdata(1:1024); %導(dǎo)入顱內(nèi)壓信號的數(shù)據(jù)%load respdata;x = respdata (1:1024);%導(dǎo)入個呼吸信號的數(shù)據(jù)p=1:60; %P的取值范圍Sw=zeros(1,length(p); %創(chuàng)建一個1行列長為length(p)的Sw零矩陣或數(shù)組E=zeros(1,length(p); %創(chuàng)建一個1行列長為length(p)的E零矩陣或數(shù)組FPE=zeros(1,length(p); %創(chuàng)建一個1行列長為length(p)的FPE零矩陣或數(shù)組for i=1:60 %嘗試改變模型階數(shù)，觀察效果Rxx = xcorr(x,biased);%估計隨機過程中的互相關(guān)序列，biased為有偏的互相關(guān)函數(shù)估計；Rtemp = zeros(1,i); %創(chuàng)建一個i行1列的Rtemp零矩陣或數(shù)組Rl = zeros(i,1); %創(chuàng)建一個i行1列的Rl零矩陣或數(shù)組for k = 1:length(Rtemp) %for循環(huán)從1到length(Rtemp)Rtemp(k) = Rxx(length(x)-1+k); %取Rxx（0）到Rxx（p-1）賦值Rl(k) = Rxx(length(x)+k); %取Rxx（1）到Rxx（p）賦值endRs = toeplitz(Rtemp); %生成自相關(guān)系數(shù)矩陣(Toeplitz型)A = -inv(Rs)*Rl; %AR模型系數(shù)估計Sw(i)= Rtemp(1),Rl*1;A; %白噪聲方差估計% 采用malab自帶函數(shù)估計模型系數(shù)a,E(i) = aryule(x,i); %malab實現(xiàn)L-D算法的AR模型參數(shù)估計，a-系數(shù)，E-預(yù)測誤差，k-反射系數(shù)%a,E(i) = arburg(x,i); %malab實現(xiàn)Burg算法的AR模型參數(shù)估計da = a(2:end)-A; %自編程序求解是否正確？FPE(i)=E(i)*(M+i+1)/(M-i-1); %FPE算法w = randn(size(x); %生成隨機噪聲x2 = filter(1,a,w); %仿真數(shù)據(jù)endfigure %畫圖subplot(3,1,1),plot(p,E,-*),grid on; %創(chuàng)建窗口畫圖，并添加網(wǎng)格線title(E隨階數(shù)p的變化情況);xlabel(p);ylabel(error);%添加標題和橫縱坐標subplot(3,1,2),plot(p,Sw,-*),grid on; %創(chuàng)建窗口畫圖，并添加網(wǎng)格線title(Sw隨階數(shù)p的變化情況);xlabel(p);ylabel(白噪聲方差估計);subplot(3,1,3),plot(p,FPE,-*),grid on; %創(chuàng)建窗口畫圖，并添加網(wǎng)格線title(FPE隨階數(shù)p的變化情況);xlabel(p);ylabel(FEP);%添加標題和橫縱坐標心電信號：a）L-D算法 b）Burg算法 圖1 心電信號的不同算法的E、Sw、FPE隨階數(shù)p的變化圖分析：首先，根據(jù)FPE準則找到最佳階數(shù)P； L-D算法; 在command window 輸入find(FPE=min(FPE)，結(jié)果：ans=15； 在command window 輸入find(abs(E-1)=min(abs(E-1)，結(jié)果：ans=1； Burg算法：在command window 輸入find(FPE=min(FPE)，結(jié)果：ans=31； 在command window 輸入find(abs(E-1)=min(abs(E-1)，結(jié)果：ans=1。腦電信號： a）L-D算法b）Burg算法圖2 腦電信號的不同算法的E、Sw、FPE隨階數(shù)p的變化圖分析：首先，根據(jù)FPE準則找到最佳階數(shù)P； L-D算法： 在command window 輸入find(FPE=min(FPE)，結(jié)果：ans=23； 在command window 輸入find(abs(E-1)=min(abs(E-1)，結(jié)果：ans=11；Burg算法：在command window 輸入find(FPE=min(FPE)，結(jié)果：ans=43； 在command window 輸入find(abs(E-1)=min(abs(E-1)，結(jié)果：ans=11。顱內(nèi)壓信號：a）L-D算法b）Burg算法 圖3顱內(nèi)壓信號的不同算法的E、Sw、FPE隨階數(shù)p的變化圖分析：首先，根據(jù)FPE準則找到最佳階數(shù)P；L-D算法： 在command window 輸入find(FPE=min(FPE)，結(jié)果：ans=43； 在command window 輸入find(abs(E-1)=min(abs(E-1)，結(jié)果：ans=1；Burg算法：在command window 輸入find(FPE=min(FPE)，結(jié)果：ans=57； 在command window 輸入find(abs(E-1)=min(abs(E-1)，結(jié)果：ans=1。 呼吸信號：a）L-D算法b）Burg算法圖4 呼吸信號的不同算法的E、Sw、FPE隨階數(shù)p的變化圖分析：首先，根據(jù)FPE準則找到最佳階數(shù)P； L-D算法： 在command window 輸入find(FPE=min(FPE)，結(jié)果：ans=59； 在command window 輸入find(abs(E-1)=min(abs(E-1)，結(jié)果：ans=60；Burg算法：在command window 輸入find(FPE=min(FPE)，結(jié)果：ans=44； 在command window 輸入find(abs(E-1)=min(abs(E-1)，結(jié)果：ans=60。分析：L-D算法和Burg算法的用find(abs(E-1)=min(abs(E-1)找到的值是一樣的； 心電、腦電、顱內(nèi)壓信號的找find(FPE=min(FPE)的值Burg算法比L-D算法大，呼吸 信號找到的min(FPE)值，Burg算法比L-D算法小。 思考題：clear; clc;load ecgdata;x = ecgdata (1:1024); %導(dǎo)入心電信號%load eegdata;x = eegdata (1:1024); %導(dǎo)入腦電信號%load icpdata;x = icpdata(1:1024); %導(dǎo)入顱內(nèi)壓信號%load respdata;x = respdata (1:1024); %導(dǎo)入個呼吸信號%p =20; %嘗試改變模型階數(shù)，觀察效果for p=2:2:20 %嘗試改變模型階數(shù)，觀察效果Rxx = xcorr(x,biased);%估計隨機過程中的互相關(guān)序列，biased為有偏的互相關(guān)函數(shù)估計；Rtemp = zeros(1,p); %創(chuàng)建一個i行1列的Rtemp零矩陣或數(shù)組Rl = zeros(p,1); %創(chuàng)建一個i行1列的Rl零矩陣或數(shù)組for k = 1:length(Rtemp) %for循環(huán)從1到length(Rtemp)Rtemp(k) = Rxx(length(x)-1+k); %取Rxx（0）到Rxx（p-1）賦值Rl(k) = Rxx(length(x)+k); %取Rxx（1）到Rxx（p）賦值endRs = toeplitz(Rtemp); %生成自相關(guān)系數(shù)矩陣(Toeplitz型)A = -inv(Rs)*Rl; %AR模型系數(shù)估計Sw(p/2) = Rtemp(1),Rl*1;A; %白噪聲方差估計% 采用malab自帶函數(shù)估計模型系數(shù)a,E = aryule(x,p); %a-系數(shù)，E-預(yù)測誤差，k-反射系數(shù)%a,E = arburg(x,p); %malab實現(xiàn)Burg算法的AR模型參數(shù)估計da = a(2:end)-A %自編程序求解是否正確？stem(da);title(參數(shù)估計偏差) %畫圖w = randn(size(x); %產(chǎn)生隨機噪聲信號x2 = filter(1,a,w); %仿真數(shù)據(jù)figure;subplot(2,1,1);plot(x);title(真實數(shù)據(jù)); %繪制真實數(shù)據(jù)的圖像subplot(2,1,2);plot(x2);title(仿真數(shù)據(jù));%繪制仿真數(shù)據(jù)的圖像error(p/2)=mean(x-x2).2); %顯示最小均方誤差的計算結(jié)果endfiguresubplot(1,2,1),plot(2:2:20,error,-*) %畫圖title(最小均方誤差隨階數(shù)p的變化情況);xlabel(p);ylabel(error);%繪制最小均方誤差隨階數(shù)p的變化情況圖subplot(1,2,2),stem(2:2:20,Sw,-*),grid ontitle(白噪聲方差估計隨階數(shù)p的變化情況);xlabel(p);ylabel(白噪聲方差估計);%繪制白噪聲方差估計隨階數(shù)p的變化情況圖Rxx2=xcorr(x2,biased); %求仿真數(shù)據(jù)的自相關(guān)函數(shù) Px=abs(fft(Rxx); %計算真實信號自功率譜 Px2=abs(fft(Rxx2); %計算仿真信號自功率譜figuresubplot(2,1,1);stem(-1023:1023,Px);title(真實信號功率譜); %繪制火柴桿圖subplot(2,1,2);stem(-1023:1023,Px2);title(仿真信號功率譜); %繪制火柴桿圖%繪制Sw、E隨p變化散點圖，以下是統(tǒng)計得到的結(jié)果 p=8 9 10 11 12 13 14 15 16;Sw=1.0527 1.0301 1.0150 0.9961 0.9943 0.9942 0.9896 0.9877 0.9872; E=3.8741 3.4337 3.3200 3.7314 3.3683 3.4809 3.5826 3.6253 3.4176; figuresubplot(2,1,1);stem(p,Sw); %繪制火柴桿圖xlabel(p);ylabel(Sw);title(Sw隨p變化散點圖);%添加標題和橫縱坐標subplot(2,1,2);plot(p,E,-o); %繪制散點圖xlabel(p);ylabel(E);title(E隨p變化散點圖); %添加標題和橫縱坐標1.比較四種信號的參數(shù)估計偏差 a）心電 b）腦電 c）顱內(nèi)壓 d）呼吸圖5 心電、腦電、顱內(nèi)壓、呼吸信號的的真實數(shù)據(jù)和參數(shù)估計偏差分析：由圖5可以看出腦電信號的參數(shù)估計偏差比其它三個信號的參數(shù)估計偏差要小許多，腦電信號的自編程序跟 MATLAB 信號處理工具箱自帶函數(shù)aryule 的處理更接近。由此可知，L-D 算法與自編程序相比較，自編程序?qū)烙嫷膮?shù)比較精確一些。2.比較四種信號的真實數(shù)據(jù)與仿真數(shù)據(jù)a）心電 b）腦電 c）顱內(nèi)壓 d）呼吸圖6 心電、腦電、顱內(nèi)壓、呼吸信號的的真實數(shù)據(jù)與仿真數(shù)據(jù)分析：從圖6可以看出，對心電、腦電、顱內(nèi)壓、呼吸信號建模后，心電、顱內(nèi)壓、呼吸信號的真實數(shù)據(jù)和仿真數(shù)據(jù)相差很大，腦電模型產(chǎn)生的信號更能真實的反映腦電信號的特征。不同噪聲的激勵得到的信號時不同的。3.比較四種信號L-D 算法與 Burg 算法的功率譜L-D 算法： a）心電 b）腦電 c）顱內(nèi)壓 d）呼吸Burg 算法：a）心電 b）腦電 c）顱內(nèi)壓 d）呼吸圖7 心電、腦電、顱內(nèi)壓、呼吸信號的L-D 算法與 Burg 算法功率譜分析：由圖7看出，在比較心電、腦電、顱內(nèi)壓、呼吸信號的L-D 算法與 Burg 算法功率譜發(fā)現(xiàn)，對信號進行功率譜估計時， 腦電和顱內(nèi)壓信號的功率譜的縱坐標相差不大，但是相對而說心電和呼吸信號的相差是極大的。4.比較四種信號的L-D 算法與 Burg 算法L-D 算法 a）心電 b）腦電 c）顱內(nèi)壓 d）呼吸圖8四種信號L-D算法，階數(shù)p與均方誤差error和噪聲方差估計值Sw之間的關(guān)系Burg 算法 a）心電 b）腦電 c）顱內(nèi)壓 d）呼吸圖9四種信號Burg算法，階數(shù)p與均方誤差error和噪聲方差估計值Sw之間的關(guān)系分析：由圖8和圖9可以看出，雖然四種信號的L-D 算法與 Burg 算法階數(shù)p與均方誤差error的變化是不穩(wěn)定的，但是它們的階數(shù)p和噪聲方差估計值Sw之間的變化是大體一致的。4.1心電信號的L-D 算法與 Burg 算法比較 a）L-D算法 b）Burg算法 a）L-D算法 b）Burg算法圖10心電信號不同算法的真實數(shù)據(jù)與參數(shù)估計 圖11心電信號不同算法的真實與仿真數(shù)據(jù) 表1 心電信號的L-D 算法與 Burg 算法的Da值變化表Da值L-D算法-8.990e-13-5.223e-12-1.340e-11-8.120e-12-2.144e-124.546e-122.100e-113.596e-115.641e-113.163e-11-1.130e-113.8190e-123.5763e-113.2489e-113.8141e-113.1735e-111.5867e-11-5.012e-12-3.589e-123.7491e-14Burg算法-0.2810.499-0.092-0.2750.1030.0610.026-0.0970.0530.027-0.0435-0.01710.01170.1069-0.14150.04950.0247-0.05730.0670-0.02844.2腦電信號的L-D 算法與 Burg 算法比較a）L-D算法 b）Burg算法 a）L-D算法 b）Burg算法圖12腦電信號不同算法的真實數(shù)據(jù)與參數(shù)估計 圖13腦電信號不同算法的真實與仿真數(shù)據(jù)表2 腦電信號的L-D 算法與 Burg 算法的Da值變化表Da值L-D算法5.5511e-161.3877e-17-6.939e-180-1.110e-16-1.665e-169.7144e-174.8572e-17-8.3266-170-8.326e-172.1510e-16-2.082e-171.3877e-17-2.168e-161.2490e-16-1.145e-16-2.081e-171.1102e-165.5511e-17Burg算法0.0009-0.0003-0.0006-0.0004-0.00140.00040.0013-0.0003-0.0006-0.00050.0005-0.00170.00030.00090.00050.00150.0009-0.0001-0.00170.00024.3顱內(nèi)壓信號的L-D 算法與 Burg 算法比較a）L-D算法 b）Burg算法 a）L-D算法 b）Burg算法圖14顱內(nèi)壓信號不同算法的真實數(shù)據(jù)與參數(shù)估計 圖15顱內(nèi)壓信號不同算法真實與仿真數(shù)據(jù)表3 顱內(nèi)壓信號的L-D 算法與 Burg 算法的Da值變化表Da值L-D算法3.1286e-131.4018e-121.1497e-124.9960e-15-1.567e-121.4885e-134.3890e-123.5907e-12-3.099e-124.4496e-128.5008e-121.0558e-12-7.490e-12-6.810e-123.4846e-121.3522e-12-1.628e-116.4571e-124.0563e-122.9531e-14Burg算法-0.16210.08530.10570.0371-0.0042-0.0803-0.0145-0.03080.00220.02500.03560.2582-0.39560.04790.07670.03840.0087-0.04940.0169-0.00254.4呼吸信號的L-D 算法與 Burg 算法比較a）L-D算法 b）Burg算法 a）L-D算法 b）Burg算法圖16呼吸信號不同算法的真實數(shù)據(jù)與參數(shù)估計 圖17呼吸信號不同算法的真實與仿真數(shù)據(jù)表4 呼吸信號的L-D 算法與 Burg 算法的Da值變化表Da值L-D算法-2.735e-133.042e-138.882e-124.431e-12-1.890e-12-7.866e-13-2.557e-12-9.929e-13-1.172e-11-1.041e-112.6015e-118.3792e-128.1914e-121.0824e-114.7703e-121.5861e-116.6561e-12-5.402e-12-2.922e-12-6.889e-14Burg算法-0.0117-0.0006-0.00080.00010.0006-0.00060.00110.00120.00060.00130.00190.00070.00160.00050.00110.0009-0.00070.0007-0.00040.0007分析：Burg算法的優(yōu)點是：求得的AR模型是穩(wěn)定的，較高的計算效率，但遞推還是用的L-D算法，因此仍然存在明顯的缺點。5.腦電信號，改變P值，比較L-D算法與Burg算法 a）L-D算法 b）Burg算法圖18 腦電信號改變P值，比較L-D算法與Burg算法Sw和E值變化圖表5 腦電信號改變P值，比較L-D算法與Burg算法Sw和E值表P值8910111213141516L-D算法Sw1.05271.03011.01500.99610.99430.99420.98960.98770.9872E3.46873.20103.85303.38883.28743.39453.93154.15733.3942Burg算法Sw1.0521.0301.0150.9960.9940.9940.9890.9870.987E3.8743.4333.3203.7313.3683.4803.5823.6253.417分析：觀察圖18和表5可以看出：不管是 Burg 算法還是 L-D 算法在改變了P的大小之后， 噪聲方差估計Sw和