Yule-Walker方程-實(shí)驗(yàn)報(bào)告.doc

生物醫(yī)學(xué)信號(hào)處理評(píng)分大理大學(xué)實(shí)驗(yàn)報(bào)告 課程名稱 生物醫(yī)學(xué)信號(hào)處理 實(shí)驗(yàn)名稱 Yule-Walker方程 專業(yè)班級(jí) 姓 名 羽卒蘭cl 學(xué) 號(hào) 實(shí)驗(yàn)日期 2016年5月27日 實(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、硬件配置：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軟件三、實(shí)驗(yàn)內(nèi)容(包括本實(shí)驗(yàn)要完成的實(shí)驗(yàn)問(wèn)題及需要的相關(guān)知識(shí)簡(jiǎn)單概述)編寫求解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)的功率譜。四、實(shí)驗(yàn)結(jié)果與分析 （包括實(shí)驗(yàn)原理、數(shù)據(jù)的準(zhǔn)備、運(yùn)行過(guò)程分析、源程序（代碼）、圖形圖象界面等） 實(shí)驗(yàn)原理隨機(jī)信號(hào)可以看作是由當(dāng)前激勵(lì)白噪聲w(n)以及若干次以往信號(hào)x(n-k)的線性組合產(chǎn)生，即所謂自回歸模型（AR模型）模型參數(shù)滿足Yule-Walker方程矩陣形式求解Yule-Walker方程，就可以得到AR模型系數(shù)當(dāng)模型階次較大時(shí)，直接用矩陣運(yùn)算求解的計(jì)算量大，不利于實(shí)時(shí)運(yùn)算。利用系數(shù)矩陣的特性，人們提出了如L-D算法等快速算法。 源程序：clear; clc;M=1024;load ecgdata;x = ecgdata(1:1024); %導(dǎo)入心電信號(hào)的數(shù)據(jù)%load eegdata;x = eegdata(1:M); %導(dǎo)入腦電信號(hào)的數(shù)據(jù)%load icpdata;x = icpdata(1:1024); %導(dǎo)入顱內(nèi)壓信號(hào)的數(shù)據(jù)%load respdata;x = respdata (1:1024);%導(dǎo)入個(gè)呼吸信號(hào)的數(shù)據(jù)p=1:60; %P的取值范圍Sw=zeros(1,length(p); %創(chuàng)建一個(gè)1行列長(zhǎng)為length(p)的Sw零矩陣或數(shù)組E=zeros(1,length(p); %創(chuàng)建一個(gè)1行列長(zhǎng)為length(p)的E零矩陣或數(shù)組FPE=zeros(1,length(p); %創(chuàng)建一個(gè)1行列長(zhǎng)為length(p)的FPE零矩陣或數(shù)組for i=1:60 %嘗試改變模型階數(shù)，觀察效果Rxx = xcorr(x,biased);%估計(jì)隨機(jī)過(guò)程中的互相關(guān)序列，biased為有偏的互相關(guān)函數(shù)估計(jì)；Rtemp = zeros(1,i); %創(chuàng)建一個(gè)i行1列的Rtemp零矩陣或數(shù)組Rl = zeros(i,1); %創(chuàng)建一個(gè)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ù)估計(jì)Sw(i)= Rtemp(1),Rl*1;A; %白噪聲方差估計(jì)% 采用malab自帶函數(shù)估計(jì)模型系數(shù)a,E(i) = aryule(x,i); %malab實(shí)現(xiàn)L-D算法的AR模型參數(shù)估計(jì)，a-系數(shù)，E-預(yù)測(cè)誤差，k-反射系數(shù)%a,E(i) = arburg(x,i); %malab實(shí)現(xiàn)Burg算法的AR模型參數(shù)估計(jì)da = a(2:end)-A; %自編程序求解是否正確？FPE(i)=E(i)*(M+i+1)/(M-i-1); %FPE算法w = randn(size(x); %生成隨機(jī)噪聲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);%添加標(biāo)題和橫縱坐標(biāo)subplot(3,1,2),plot(p,Sw,-*),grid on; %創(chuàng)建窗口畫圖，并添加網(wǎng)格線title(Sw隨階數(shù)p的變化情況);xlabel(p);ylabel(白噪聲方差估計(jì));subplot(3,1,3),plot(p,FPE,-*),grid on; %創(chuàng)建窗口畫圖，并添加網(wǎng)格線title(FPE隨階數(shù)p的變化情況);xlabel(p);ylabel(FEP);%添加標(biāo)題和橫縱坐標(biāo)心電信號(hào)：a）L-D算法 b）Burg算法 圖1 心電信號(hào)的不同算法的E、Sw、FPE隨階數(shù)p的變化圖分析：首先，根據(jù)FPE準(zhǔn)則找到最佳階數(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。腦電信號(hào)： a）L-D算法b）Burg算法圖2 腦電信號(hào)的不同算法的E、Sw、FPE隨階數(shù)p的變化圖分析：首先，根據(jù)FPE準(zhǔn)則找到最佳階數(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)壓信號(hào)：a）L-D算法b）Burg算法 圖3顱內(nèi)壓信號(hào)的不同算法的E、Sw、FPE隨階數(shù)p的變化圖分析：首先，根據(jù)FPE準(zhǔn)則找到最佳階數(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。 呼吸信號(hào)：a）L-D算法b）Burg算法圖4 呼吸信號(hào)的不同算法的E、Sw、FPE隨階數(shù)p的變化圖分析：首先，根據(jù)FPE準(zhǔn)則找到最佳階數(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)壓信號(hào)的找find(FPE=min(FPE)的值Burg算法比L-D算法大，呼吸 信號(hào)找到的min(FPE)值，Burg算法比L-D算法小。 思考題：clear; clc;load ecgdata;x = ecgdata (1:1024); %導(dǎo)入心電信號(hào)%load eegdata;x = eegdata (1:1024); %導(dǎo)入腦電信號(hào)%load icpdata;x = icpdata(1:1024); %導(dǎo)入顱內(nèi)壓信號(hào)%load respdata;x = respdata (1:1024); %導(dǎo)入個(gè)呼吸信號(hào)%p =20; %嘗試改變模型階數(shù)，觀察效果for p=2:2:20 %嘗試改變模型階數(shù)，觀察效果Rxx = xcorr(x,biased);%估計(jì)隨機(jī)過(guò)程中的互相關(guān)序列，biased為有偏的互相關(guān)函數(shù)估計(jì)；Rtemp = zeros(1,p); %創(chuàng)建一個(gè)i行1列的Rtemp零矩陣或數(shù)組Rl = zeros(p,1); %創(chuàng)建一個(gè)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ù)估計(jì)Sw(p/2) = Rtemp(1),Rl*1;A; %白噪聲方差估計(jì)% 采用malab自帶函數(shù)估計(jì)模型系數(shù)a,E = aryule(x,p); %a-系數(shù)，E-預(yù)測(cè)誤差，k-反射系數(shù)%a,E = arburg(x,p); %malab實(shí)現(xiàn)Burg算法的AR模型參數(shù)估計(jì)da = a(2:end)-A %自編程序求解是否正確？stem(da);title(參數(shù)估計(jì)偏差) %畫圖w = randn(size(x); %產(chǎn)生隨機(jī)噪聲信號(hào)x2 = filter(1,a,w); %仿真數(shù)據(jù)figure;subplot(2,1,1);plot(x);title(真實(shí)數(shù)據(jù)); %繪制真實(shí)數(shù)據(jù)的圖像subplot(2,1,2);plot(x2);title(仿真數(shù)據(jù));%繪制仿真數(shù)據(jù)的圖像error(p/2)=mean(x-x2).2); %顯示最小均方誤差的計(jì)算結(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(白噪聲方差估計(jì)隨階數(shù)p的變化情況);xlabel(p);ylabel(白噪聲方差估計(jì));%繪制白噪聲方差估計(jì)隨階數(shù)p的變化情況圖Rxx2=xcorr(x2,biased); %求仿真數(shù)據(jù)的自相關(guān)函數(shù) Px=abs(fft(Rxx); %計(jì)算真實(shí)信號(hào)自功率譜 Px2=abs(fft(Rxx2); %計(jì)算仿真信號(hào)自功率譜figuresubplot(2,1,1);stem(-1023:1023,Px);title(真實(shí)信號(hào)功率譜); %繪制火柴桿圖subplot(2,1,2);stem(-1023:1023,Px2);title(仿真信號(hào)功率譜); %繪制火柴桿圖%繪制Sw、E隨p變化散點(diǎn)圖，以下是統(tǒng)計(jì)得到的結(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變化散點(diǎn)圖);%添加標(biāo)題和橫縱坐標(biāo)subplot(2,1,2);plot(p,E,-o); %繪制散點(diǎn)圖xlabel(p);ylabel(E);title(E隨p變化散點(diǎn)圖); %添加標(biāo)題和橫縱坐標(biāo)1.比較四種信號(hào)的參數(shù)估計(jì)偏差 a）心電 b）腦電 c）顱內(nèi)壓 d）呼吸圖5 心電、腦電、顱內(nèi)壓、呼吸信號(hào)的的真實(shí)數(shù)據(jù)和參數(shù)估計(jì)偏差分析：由圖5可以看出腦電信號(hào)的參數(shù)估計(jì)偏差比其它三個(gè)信號(hào)的參數(shù)估計(jì)偏差要小許多，腦電信號(hào)的自編程序跟 MATLAB 信號(hào)處理工具箱自帶函數(shù)aryule 的處理更接近。由此可知，L-D 算法與自編程序相比較，自編程序?qū)烙?jì)的參數(shù)比較精確一些。2.比較四種信號(hào)的真實(shí)數(shù)據(jù)與仿真數(shù)據(jù)a）心電 b）腦電 c）顱內(nèi)壓 d）呼吸圖6 心電、腦電、顱內(nèi)壓、呼吸信號(hào)的的真實(shí)數(shù)據(jù)與仿真數(shù)據(jù)分析：從圖6可以看出，對(duì)心電、腦電、顱內(nèi)壓、呼吸信號(hào)建模后，心電、顱內(nèi)壓、呼吸信號(hào)的真實(shí)數(shù)據(jù)和仿真數(shù)據(jù)相差很大，腦電模型產(chǎn)生的信號(hào)更能真實(shí)的反映腦電信號(hào)的特征。不同噪聲的激勵(lì)得到的信號(hào)時(shí)不同的。3.比較四種信號(hào)L-D 算法與 Burg 算法的功率譜L-D 算法： a）心電 b）腦電 c）顱內(nèi)壓 d）呼吸Burg 算法：a）心電 b）腦電 c）顱內(nèi)壓 d）呼吸圖7 心電、腦電、顱內(nèi)壓、呼吸信號(hào)的L-D 算法與 Burg 算法功率譜分析：由圖7看出，在比較心電、腦電、顱內(nèi)壓、呼吸信號(hào)的L-D 算法與 Burg 算法功率譜發(fā)現(xiàn)，對(duì)信號(hào)進(jìn)行功率譜估計(jì)時(shí)， 腦電和顱內(nèi)壓信號(hào)的功率譜的縱坐標(biāo)相差不大，但是相對(duì)而說(shuō)心電和呼吸信號(hào)的相差是極大的。4.比較四種信號(hào)的L-D 算法與 Burg 算法L-D 算法 a）心電 b）腦電 c）顱內(nèi)壓 d）呼吸圖8四種信號(hào)L-D算法，階數(shù)p與均方誤差error和噪聲方差估計(jì)值Sw之間的關(guān)系Burg 算法 a）心電 b）腦電 c）顱內(nèi)壓 d）呼吸圖9四種信號(hào)Burg算法，階數(shù)p與均方誤差error和噪聲方差估計(jì)值Sw之間的關(guān)系分析：由圖8和圖9可以看出，雖然四種信號(hào)的L-D 算法與 Burg 算法階數(shù)p與均方誤差error的變化是不穩(wěn)定的，但是它們的階數(shù)p和噪聲方差估計(jì)值Sw之間的變化是大體一致的。4.1心電信號(hào)的L-D 算法與 Burg 算法比較 a）L-D算法 b）Burg算法 a）L-D算法 b）Burg算法圖10心電信號(hào)不同算法的真實(shí)數(shù)據(jù)與參數(shù)估計(jì) 圖11心電信號(hào)不同算法的真實(shí)與仿真數(shù)據(jù) 表1 心電信號(hào)的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腦電信號(hào)的L-D 算法與 Burg 算法比較a）L-D算法 b）Burg算法 a）L-D算法 b）Burg算法圖12腦電信號(hào)不同算法的真實(shí)數(shù)據(jù)與參數(shù)估計(jì) 圖13腦電信號(hào)不同算法的真實(shí)與仿真數(shù)據(jù)表2 腦電信號(hào)的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)壓信號(hào)的L-D 算法與 Burg 算法比較a）L-D算法 b）Burg算法 a）L-D算法 b）Burg算法圖14顱內(nèi)壓信號(hào)不同算法的真實(shí)數(shù)據(jù)與參數(shù)估計(jì) 圖15顱內(nèi)壓信號(hào)不同算法真實(shí)與仿真數(shù)據(jù)表3 顱內(nèi)壓信號(hào)的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呼吸信號(hào)的L-D 算法與 Burg 算法比較a）L-D算法 b）Burg算法 a）L-D算法 b）Burg算法圖16呼吸信號(hào)不同算法的真實(shí)數(shù)據(jù)與參數(shù)估計(jì) 圖17呼吸信號(hào)不同算法的真實(shí)與仿真數(shù)據(jù)表4 呼吸信號(hào)的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)點(diǎn)是：求得的AR模型是穩(wěn)定的，較高的計(jì)算效率，但遞推還是用的L-D算法，因此仍然存在明顯的缺點(diǎn)。5.腦電信號(hào)，改變P值，比較L-D算法與Burg算法 a）L-D算法 b）Burg算法圖18 腦電信號(hào)改變P值，比較L-D算法與Burg算法Sw和E值變化圖表5 腦電信號(hào)改變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的大小之后， 噪聲方差估計(jì)Sw和