YuleWalker方程

1、评分大理大学实验报告 课程名称 生物医学信号处理 实验名称 Yule-Walker方程 专业班级 姓 名 学 号 实验日期 实验地点 20152016学年度第 3 学期一、 实验目的学习求解Yule-Walker方程，建立随机信号的AR模型。2、 实验环境 1、硬件配置： 处理器：AMD A10-5750M APU with Radeon(tm) Graphics 2.50GHz 安装内存：（RAM）4.00GB 系统类型：win10位操作系统，基于x64位处理器 2、软件环境：Matlab R2012b3、 实验原理随机信号可以看作是由当前激励白噪声w(n)以及若干次以往信号x(n-k)的线

2、性组合产生，即所谓自回归模型（AR模型）模型参数满足Yule-Walker方程矩阵形式求解Yule-Walker方程，就可以得到AR模型系数当模型阶次较大时，直接用矩阵运算求解的计算量大，不利于实时运算。利用系数矩阵的特性，人们提出了如L-D算法等快速算法。四、实验内容编写求解Yule-Walker方程的程序，并对实际生理信号（例如心电、脑电）建立AR模型。对同一数据，使用Matlab信号处理工具箱自带函数aryule计算相同阶数AR模型系数，检验程序是否正确。用伪随机序列（白噪声）驱动AR模型，观察输出是否与真实心电、脑电信号相似，对比真实信号与仿真信号的功率谱。 五、 实验结果与分析试验程

3、序：clear; clc;M=1024; load ecgdata;load eegdata;load icpdata;load respdata;x= ecgdata(1:M); % 长度任意选择%x= eegdata(1:M); % 长度任意选择%x= icpdata(1:M); % 长度任意选择%x= respdata(1:M); % 长度任意选择 i=1:15;Sw=zeros(1,length(i);E=zeros(1,length(i);FPE=zeros(1,length(i);for p=1:15 % 尝试改变模型阶数，观察效果Rxx = xcorr(x,'biased

4、');%biased为有偏的互相关函数估计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); % 生成自相关系数矩阵指矩阵中每条自左上至右下的 % 斜线上的元素相同，toeplitz(x)用向量x生成一 % 个对称的托普利兹矩阵A = -inv(Rs)*Rl; % AR模型系数估计Sw(p)= Rtemp(1),Rl'*1;A; % 白噪声方差估计% 采用m

5、alab自带函数估计模型系数a,E(p) = aryule(x,p); % a-系数，E-预测误差，k-反射系数 %a,E(p) = arburg(x,p);%Y = filter(B,A,X) ，输入X为滤波前序列，Y为滤波结果序列，%B/A 提供滤波器系数，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); % 仿真数据Rx2 = xcorr(x2,'biased');xx = abs(fft(Rxx);x

6、x2 = abs(fft(Rx2);figure;subplot(4,1,1);plot(x); xlim(0 1024);title('真实数据');subplot(4,1,2);plot(x2);xlim(0 1024);title('仿真数据');subplot(4,1,3); plot(-1023:1023,xx);title('真实数据功率谱'); subplot(4,1,4); plot(-1023:1023,xx2);title('仿真数据功率谱');error(p)=mean(x-x2).2);%求最小均方误差en

7、dPopt=find(FPE=min(FPE)figure,subplot(1,3,1);plot(i,error,'-*')title('(a)最小均方误差随阶数p的变化情况'),xlabel('p');ylabel('error');subplot(1,3,2);plot(i,E,'-*');grid ontitle('(b)预测误差随阶数p的变化情况'),xlabel('p');ylabel('E');%figure,stem(i,Sw,'-*'

8、);grid on%title('白噪声方差估计随阶数p的变化情况')%xlabel('p');ylabel('白噪声方差估计');subplot(1,3,3);plot(i,FPE,'-*');title('(c)FPE随阶数p的变化情况');xlabel('p');ylabel('FPE');实验结果：1、 心电数据，阶数p为15,M=1024: （1）最小均方误差随p的变化情况、预测误差随p的变化情况、FPE随p的变化情况： 图1-1 L-D算法 图1-2 Burg算法（2）白

9、噪声方差估计随阶数p的变化情况： 图1-3 L-D算法 图1-4 Burg算法（3）真实数据和仿真数据及其功率谱： 图1-5 L-D算法error=19.234 图1-6 Burg算法 error=19.9132、脑电数据，阶数p为14，M=1024：（1）最小均方误差随p的变化情况、预测误差随p的变化情况、FPE随p的变化情况： 图2-1 L-D算法 图2-2 Burg算法（2） 白噪声方差估计随阶数p的变化情况：图 2-3 L-D算法 图2-4 Burg算法（3）真实数据和仿真数据及其功率谱： 图2-5 L-D算法 error=3.3147 图2-6 Burg算法 error=3.2043

10、3、颅内压数据，阶数p为14，M=1024：(1)最小均方误差随p的变化情况、预测误差随p的变化情况、FPE随p的变化情况： 图3-1 L-D算法 图3-2 Burg算法(2) 白噪声方差估计随阶数p的变化情况： 图3-1 L-D算法 图3-2 Burg算法(3) 真实数据和仿真数据及其功率谱： 图3-3 L-D算法error=52.4963 图3-4 Burg算法 error=103.77014、 呼吸压数据，阶数p为7，M=1024：(1)最小均方误差随p的变化情况、预测误差随p的变化情况、FPE随p的变化情况：图4-1 L-D算法 图4-2 Burg算法(2)白噪声方差估计随阶数p的变化

11、情况：图4-3 L-D算法 图4-4 Burg算法（3）真实数据和仿真数据及其功率谱： 图4-5 L-D算法 2.7785e05 图4-6 Burg算法2.17161e05 5、心电数据，阶数p为5,M=512: （1）最小均方误差随p的变化情况、预测误差随p的变化情况、FPE随p的变化情况：图5-1 L-D算法 图5-2 Burg算法（2）白噪声方差估计随阶数p的变化情况：图5-1 L-D算法 图5-2 Burg算（3）真实数据和仿真数据及其功率谱：图5-5 L-D算法 error=21.1246 图5-6 Burg算法error=11.8885 6、脑电数据，阶数p为6,M=512: （1

12、）最小均方误差随p的变化情况、预测误差随p的变化情况、FPE随p的变化情况：图6-1 L-D算法 图6-2 Burg算法（2）白噪声方差估计随阶数p的变化情况：图6-3 L-D算法 图6-4 Burg算法（3）真实数据和仿真数据及其功率谱： 图6-5 L-D算法 error=2.5564 图6-6 Burg算法 error=2.71037、颅内压数据，阶数p为14，M=512：(1)最小均方误差随p的变化情况、预测误差随p的变化情况、FPE随p的变化情况：图7-1 L-D算法 图7-2 Burg算法（2） 白噪声方差估计随阶数p的变化情况：图7-3 L-D算法 图7-4 Burg算法（3）真实

13、数据和仿真数据及其功率谱： 图7-4 L-D算法 error=17.0841 图7-5 Burg算法error=43.72808、呼吸压数据，阶数p为1，M=512：(1)最小均方误差随p的变化情况、预测误差随p的变化情况、FPE随p的变化情况：图8-1 L-D算法 图8-2 Burg算法（2）白噪声方差估计随阶数p的变化情况：图8-3 L-D算法 图8-4 Burg算法（3）真实数据和仿真数据及其功率谱： 图8-5 L-D算法 error=2.9724e05 图8-6 Burg算法 error=2.9917e05结果分析：1、 对信号进行功率谱估计是，脑电信号的纵坐标相差不大，而心电信号、颅内压信号、呼吸信号的纵坐标相差较大。2、 当AR模型结束p确定时，L-D算法的均方误差比Brug算法的均方误差要小。六、实验小结：1、L-D递推算法是模型阶数逐渐加大的一种算法，为了克服其误差，提出了Burg算法，其基本思想是对观测数据进行向前或向后预测，然后让两者的均方误差之和最小做为估计准则估计反射系数进而根据L-D算法的递推公式求出AR模型参数。2、Burg算法