Yule-Walker方程-实验报告

评分大理大学实验报告课程名称生物医学信号处理实验名称Yule-Walker方程专业班级姓名羽卒兰cl学号实验日期2016年5月27日实验地点2015—2016学年度第3学期生物医学信号处理第2页共13页一、实验目的学习求解Yule-Walker方程，建立随机信号的AR模型。二、实验环境1、硬件配置：Intel（R）Core(TM)i5-4210UCPU@1.7GHz1.7GHz安装内存（RAM）：4.00GB系统类型：64位操作系统2、软件环境：MATLABR2013b软件三、实验内容(包括本实验要完成的实验问题及需要的相关知识简单概述)编写求解Yule-Walker方程的程序，并对实际生理信号（例如心电、脑电）建立AR模型。对同一数据，使用Matlab信号处理工具箱自带函数aryule计算相同阶数AR模型系数，检验程序是否正确。用伪随机序列（白噪声）驱动AR模型，观察输出是否与真实心电、脑电信号相似，对比真实信号与仿真信号的功率谱。四、实验结果与分析（包括实验原理、数据的准备、运行过程分析、源程序（代码）、图形图象界面等）实验原理随机信号可以看作是由当前激励白噪声w(n)以及若干次以往信号x(n-k)的线性组合产生，即所谓自回归模型（AR模型）pkkkpkkzazHknxanwnx1111模型参数满足Yule-Walker方程01mkmRamRpkxxkxx矩阵形式pRRRaaaRpRpRpRRRpRRRp2102120111021求解Yule-Walker方程，就可以得到AR模型系数当模型阶次较大时，直接用矩阵运算求解的计算量大，不利于实时运算。利用系数矩阵的特性，人们提出了如L-D算法等快速算法。生物医学信号处理第3页共13页源程序：clear;clc;M=1024;loadecgdata;x=ecgdata(1:1024);%导入心电信号的数据%loadeegdata;x=eegdata(1:M);%导入脑电信号的数据%loadicpdata;x=icpdata(1:1024);%导入颅内压信号的数据%loadrespdata;x=respdata(1:1024);%导入个呼吸信号的数据p=1:60;%P的取值范围Sw=zeros(1,length(p));%创建一个1行列长为length(p)的Sw零矩阵或数组E=zeros(1,length(p));%创建一个1行列长为length(p)的E零矩阵或数组FPE=zeros(1,length(p));%创建一个1行列长为length(p)的FPE零矩阵或数组fori=1:60%尝试改变模型阶数，观察效果Rxx=xcorr(x,'biased');%估计随机过程中的互相关序列，biased为有偏的互相关函数估计；Rtemp=zeros(1,i);%创建一个i行1列的Rtemp零矩阵或数组Rl=zeros(i,1);%创建一个i行1列的Rl零矩阵或数组fork=1:length(Rtemp)%for循环从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);%生成自相关系数矩阵(Toeplitz型)A=-inv(Rs)*Rl;%AR模型系数估计Sw(i)=[Rtemp(1),Rl']*[1;A];%白噪声方差估计%采用malab自带函数估计模型系数[a,E(i)]=aryule(x,i);%malab实现L-D算法的AR模型参数估计，a--系数，E--预测误差，k--反射系数%[a,E(i)]=arburg(x,i);%malab实现Burg算法的AR模型参数估计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);%仿真数据endfigure%画图subplot(3,1,1),plot(p',E,'-*'),gridon;%创建窗口画图，并添加网格线title('E随阶数p的变化情况');xlabel('p');ylabel('error');%添加标题和横纵坐标subplot(3,1,2),plot(p',Sw,'-*'),gridon;%创建窗口画图，并添加网格线title('Sw随阶数p的变化情况');xlabel('p');ylabel('白噪声方差估计');subplot(3,1,3),plot(p',FPE,'--*'),gridon;%创建窗口画图，并添加网格线title('FPE随阶数p的变化情况');xlabel('p');ylabel('FEP');%添加标题和横纵坐标心电信号：生物医学信号处理第4页共13页01020304050600123x10-3E随阶数p的变化情况perror01020304050600123x10-3Sw随阶数p的变化情况p白噪声方差估计01020304050600123x10-3FPE随阶数p的变化情况pFEPa）L-D算法01020304050600123x10-3E随阶数p的变化情况perror01020304050600123x10-3Sw随阶数p的变化情况p白噪声方差估计01020304050600123x10-3FPE随阶数p的变化情况pFEPb）Burg算法图1心电信号的不同算法的E、Sw、FPE随阶数p的变化图分析：首先，根据FPE准则找到最佳阶数P；L-D算法;在commandwindow输入find(FPE==min(FPE))，结果：ans=15；在commandwindow输入find(abs(E-1)==min(abs(E-1)))，结果：ans=1；Burg算法：在commandwindow输入find(FPE==min(FPE))，结果：ans=31；在commandwindow输入find(abs(E-1)==min(abs(E-1)))，结果：ans=1。脑电信号：01020304050600.911.11.21.3E随阶数p的变化情况perror01020304050600.911.11.21.3Sw随阶数p的变化情况p白噪声方差估计010203040506011.11.21.31.4FPE随阶数p的变化情况pFEPa）L-D算法生物医学信号处理第5页共13页01020304050600.811.21.41.6E随阶数p的变化情况perror01020304050600.911.11.21.3Sw随阶数p的变化情况p白噪声方差估计01020304050600.911.11.21.3FPE随阶数p的变化情况pFEPb）Burg算法图2脑电信号的不同算法的E、Sw、FPE随阶数p的变化图分析：首先，根据FPE准则找到最佳阶数P；L-D算法：在commandwindow输入find(FPE==min(FPE))，结果：ans=23；在commandwindow输入find(abs(E-1)==min(abs(E-1)))，结果：ans=11；Burg算法：在commandwindow输入find(FPE==min(FPE))，结果：ans=43；在commandwindow输入find(abs(E-1)==min(abs(E-1)))，结果：ans=11。颅内压信号：01020304050600.050.10.150.2E随阶数p的变化情况perror01020304050600.050.10.150.2Sw随阶数p的变化情况p白噪声方差估计01020304050600.050.10.150.2FPE随阶数p的变化情况pFEPa）L-D算法010203040506000.050.10.150.2E随阶数p的变化情况perror01020304050600.050.10.150.2Sw随阶数p的变化情况p白噪声方差估计010203040506000.050.10.150.2FPE随阶数p的变化情况pFEPb）Burg算法图3颅内压信号的不同算法的E、Sw、FPE随阶数p的变化图生物医学信号处理第6页共13页分析：首先，根据FPE准则找到最佳阶数P；L-D算法：在commandwindow输入find(FPE==min(FPE))，结果：ans=43；在commandwindow输入find(abs(E-1)==min(abs(E-1)))，结果：ans=1；Burg算法：在commandwindow输入find(FPE==min(FPE))，结果：ans=57；在commandwindow输入find(abs(E-1)==min(abs(E-1)))，结果：ans=1。呼吸信号：01020304050601200140016001800E随阶数p的变化情况perror01020304050601200140016001800Sw随阶数p的变化情况p白噪声方差估计010203040506014001500160017001800FPE随阶数p的变化情况pFEPa）L-D算法01020304050604006008001000E随阶数p的变化情况perror01020304050601200140016001800Sw随阶数p的变化情况p白噪声方差估计01020304050604006008001000FPE随阶数p的变化情况pFEPb）Burg算法图4呼吸信号的不同算法的E、Sw、FPE随阶数p的变化图分析：首先，根据FPE准则找到最佳阶数P；L-D算法：在commandwindow输入find(FPE==min(FPE))，结果：ans=59；在commandwindow输入find(abs(E-1)==min(abs(E-1)))，结果：ans=60；Burg算法：在commandwindow输入find(FPE==min(FPE))，结果：ans=44；在commandwindow输入find(abs(E-1)==min(abs(E-1)))，结果：ans=60。分析：L-D算法和Burg算法的用find(abs(E-1)==min(abs(E-1)))找到的值是一样的；心电、脑电、颅内压信号的找find(FPE==min(FPE))的值Burg算法比L-D算法大，呼吸信号找到的min(FPE)值，Burg算法比L-D算法小。思考题：生物医学信号处理第7页共13页clear;clc;loadecgdata;x=ecgdata(1:1024);%导入心电信号%loadeegdata;x=eegdata(1:1024);%导入脑电信号%loadicpdata;x=icpdata(1:1024);%导入颅内压信号%loadrespdata;x=respdata(1:1024);%导入个呼吸信号%p=20;%尝试改变模型阶数，观察效果forp=2:2:20%尝试改变模型阶数，观察效果Rxx=xcorr(x,'biased');%估计随机过程中的互相关序列，biased为有偏的互相关函数估计；Rtemp=zeros(1,p);%创建一个i行1列的Rtemp零矩阵或数组Rl=zeros(p,1);%创建一个i行1列的Rl零矩阵或数组fork=1:length(Rtemp)%for循环从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);%生成自相关系