Yule-Walker方程-实验报告.doc

生物医学信号处理评分大理大学实验报告 课程名称 生物医学信号处理 实验名称 Yule-Walker方程 专业班级 姓 名 羽卒兰cl 学 号 实验日期 2016年5月27日 实验地点 20152016学年度第 3 学期一、 实验目的学习求解Yule-Walker方程，建立随机信号的AR模型。二、实验环境 1、硬件配置：Intel（R） Core(TM) i5-4210U CPU 1.7GHz 1.7GHz 安装内存（RAM）：4.00GB 系统类型：64位操作系统 2、软件环境：MATLAB R2013b软件三、实验内容(包括本实验要完成的实验问题及需要的相关知识简单概述)编写求解Yule-Walker方程的程序，并对实际生理信号（例如心电、脑电）建立AR模型。对同一数据，使用Matlab信号处理工具箱自带函数aryule计算相同阶数AR模型系数，检验程序是否正确。用伪随机序列（白噪声）驱动AR模型，观察输出是否与真实心电、脑电信号相似，对比真实信号与仿真信号的功率谱。四、实验结果与分析 （包括实验原理、数据的准备、运行过程分析、源程序（代码）、图形图象界面等） 实验原理随机信号可以看作是由当前激励白噪声w(n)以及若干次以往信号x(n-k)的线性组合产生，即所谓自回归模型（AR模型）模型参数满足Yule-Walker方程矩阵形式求解Yule-Walker方程，就可以得到AR模型系数当模型阶次较大时，直接用矩阵运算求解的计算量大，不利于实时运算。利用系数矩阵的特性，人们提出了如L-D算法等快速算法。 源程序：clear; clc;M=1024;load ecgdata;x = ecgdata(1:1024); %导入心电信号的数据%load eegdata;x = eegdata(1:M); %导入脑电信号的数据%load icpdata;x = icpdata(1:1024); %导入颅内压信号的数据%load respdata;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零矩阵或数组for i=1:60 %尝试改变模型阶数，观察效果Rxx = xcorr(x,biased);%估计随机过程中的互相关序列，biased为有偏的互相关函数估计；Rtemp = zeros(1,i); %创建一个i行1列的Rtemp零矩阵或数组Rl = zeros(i,1); %创建一个i行1列的Rl零矩阵或数组for k = 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,-*),grid on; %创建窗口画图，并添加网格线title(E随阶数p的变化情况);xlabel(p);ylabel(error);%添加标题和横纵坐标subplot(3,1,2),plot(p,Sw,-*),grid on; %创建窗口画图，并添加网格线title(Sw随阶数p的变化情况);xlabel(p);ylabel(白噪声方差估计);subplot(3,1,3),plot(p,FPE,-*),grid on; %创建窗口画图，并添加网格线title(FPE随阶数p的变化情况);xlabel(p);ylabel(FEP);%添加标题和横纵坐标心电信号：a）L-D算法 b）Burg算法 图1 心电信号的不同算法的E、Sw、FPE随阶数p的变化图分析：首先，根据FPE准则找到最佳阶数P； L-D算法; 在command window 输入find(FPE=min(FPE)，结果：ans=15； 在command window 输入find(abs(E-1)=min(abs(E-1)，结果：ans=1； Burg算法：在command window 输入find(FPE=min(FPE)，结果：ans=31； 在command window 输入find(abs(E-1)=min(abs(E-1)，结果：ans=1。脑电信号： a）L-D算法b）Burg算法图2 脑电信号的不同算法的E、Sw、FPE随阶数p的变化图分析：首先，根据FPE准则找到最佳阶数P； L-D算法： 在command window 输入find(FPE=min(FPE)，结果：ans=23； 在command window 输入find(abs(E-1)=min(abs(E-1)，结果：ans=11；Burg算法：在command window 输入find(FPE=min(FPE)，结果：ans=43； 在command window 输入find(abs(E-1)=min(abs(E-1)，结果：ans=11。颅内压信号：a）L-D算法b）Burg算法 图3颅内压信号的不同算法的E、Sw、FPE随阶数p的变化图分析：首先，根据FPE准则找到最佳阶数P；L-D算法： 在command window 输入find(FPE=min(FPE)，结果：ans=43； 在command window 输入find(abs(E-1)=min(abs(E-1)，结果：ans=1；Burg算法：在command window 输入find(FPE=min(FPE)，结果：ans=57； 在command window 输入find(abs(E-1)=min(abs(E-1)，结果：ans=1。 呼吸信号：a）L-D算法b）Burg算法图4 呼吸信号的不同算法的E、Sw、FPE随阶数p的变化图分析：首先，根据FPE准则找到最佳阶数P； L-D算法： 在command window 输入find(FPE=min(FPE)，结果：ans=59； 在command window 输入find(abs(E-1)=min(abs(E-1)，结果：ans=60；Burg算法：在command window 输入find(FPE=min(FPE)，结果：ans=44； 在command window 输入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算法小。 思考题：clear; clc;load ecgdata;x = ecgdata (1:1024); %导入心电信号%load eegdata;x = eegdata (1:1024); %导入脑电信号%load icpdata;x = icpdata(1:1024); %导入颅内压信号%load respdata;x = respdata (1:1024); %导入个呼吸信号%p =20; %尝试改变模型阶数，观察效果for p=2:2:20 %尝试改变模型阶数，观察效果Rxx = xcorr(x,biased);%估计随机过程中的互相关序列，biased为有偏的互相关函数估计；Rtemp = zeros(1,p); %创建一个i行1列的Rtemp零矩阵或数组Rl = zeros(p,1); %创建一个i行1列的Rl零矩阵或数组for k = 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(p/2) = Rtemp(1),Rl*1;A; %白噪声方差估计% 采用malab自带函数估计模型系数a,E = aryule(x,p); %a-系数，E-预测误差，k-反射系数%a,E = arburg(x,p); %malab实现Burg算法的AR模型参数估计da = a(2:end)-A %自编程序求解是否正确？stem(da);title(参数估计偏差) %画图w = randn(size(x); %产生随机噪声信号x2 = filter(1,a,w); %仿真数据figure;subplot(2,1,1);plot(x);title(真实数据); %绘制真实数据的图像subplot(2,1,2);plot(x2);title(仿真数据);%绘制仿真数据的图像error(p/2)=mean(x-x2).2); %显示最小均方误差的计算结果endfiguresubplot(1,2,1),plot(2:2:20,error,-*) %画图title(最小均方误差随阶数p的变化情况);xlabel(p);ylabel(error);%绘制最小均方误差随阶数p的变化情况图subplot(1,2,2),stem(2:2:20,Sw,-*),grid ontitle(白噪声方差估计随阶数p的变化情况);xlabel(p);ylabel(白噪声方差估计);%绘制白噪声方差估计随阶数p的变化情况图Rxx2=xcorr(x2,biased); %求仿真数据的自相关函数 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变化散点图，以下是统计得到的结果 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.比较四种信号的参数估计偏差 a）心电 b）脑电 c）颅内压 d）呼吸图5 心电、脑电、颅内压、呼吸信号的的真实数据和参数估计偏差分析：由图5可以看出脑电信号的参数估计偏差比其它三个信号的参数估计偏差要小许多，脑电信号的自编程序跟 MATLAB 信号处理工具箱自带函数aryule 的处理更接近。由此可知，L-D 算法与自编程序相比较，自编程序对估计的参数比较精确一些。2.比较四种信号的真实数据与仿真数据a）心电 b）脑电 c）颅内压 d）呼吸图6 心电、脑电、颅内压、呼吸信号的的真实数据与仿真数据分析：从图6可以看出，对心电、脑电、颅内压、呼吸信号建模后，心电、颅内压、呼吸信号的真实数据和仿真数据相差很大，脑电模型产生的信号更能真实的反映脑电信号的特征。不同噪声的激励得到的信号时不同的。3.比较四种信号L-D 算法与 Burg 算法的功率谱L-D 算法： a）心电 b）脑电 c）颅内压 d）呼吸Burg 算法：a）心电 b）脑电 c）颅内压 d）呼吸图7 心电、脑电、颅内压、呼吸信号的L-D 算法与 Burg 算法功率谱分析：由图7看出，在比较心电、脑电、颅内压、呼吸信号的L-D 算法与 Burg 算法功率谱发现，对信号进行功率谱估计时， 脑电和颅内压信号的功率谱的纵坐标相差不大，但是相对而说心电和呼吸信号的相差是极大的。4.比较四种信号的L-D 算法与 Burg 算法L-D 算法 a）心电 b）脑电 c）颅内压 d）呼吸图8四种信号L-D算法，阶数p与均方误差error和噪声方差估计值Sw之间的关系Burg 算法 a）心电 b）脑电 c）颅内压 d）呼吸图9四种信号Burg算法，阶数p与均方误差error和噪声方差估计值Sw之间的关系分析：由图8和图9可以看出，虽然四种信号的L-D 算法与 Burg 算法阶数p与均方误差error的变化是不稳定的，但是它们的阶数p和噪声方差估计值Sw之间的变化是大体一致的。4.1心电信号的L-D 算法与 Burg 算法比较 a）L-D算法 b）Burg算法 a）L-D算法 b）Burg算法图10心电信号不同算法的真实数据与参数估计 图11心电信号不同算法的真实与仿真数据 表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脑电信号不同算法的真实数据与参数估计 图13脑电信号不同算法的真实与仿真数据表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颅内压信号的L-D 算法与 Burg 算法比较a）L-D算法 b）Burg算法 a）L-D算法 b）Burg算法图14颅内压信号不同算法的真实数据与参数估计 图15颅内压信号不同算法真实与仿真数据表3 颅内压信号的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呼吸信号不同算法的真实数据与参数估计 图17呼吸信号不同算法的真实与仿真数据表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算法的优点是：求得的AR模型是稳定的，较高的计算效率，但递推还是用的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和