Yule-Walker方程

评分大理大学实验报告课程名称生物医学信号处理实验名称Yule-Walker方程专业班级姓名学号实验日期实验地点2015—2016学年度第3学期《生物医学信号处理》实验报告第2页共13页一、实验目的学习求解Yule-Walker方程，建立随机信号的AR模型。二、实验环境1、硬件配置：处理器：AMDA10-5750MAPUwithRadeon(tm)Graphics2.50GHz安装内存：（RAM）4.00GB系统类型：win10位操作系统，基于x64位处理器2、软件环境：MatlabR2012b三、实验原理随机信号可以看作是由当前激励白噪声w(n)以及若干次以往信号x(n-k)的线性组合产生，即所谓自回归模型（AR模型）pkkkpkkzazHknxanwnx1111模型参数满足Yule-Walker方程01mkmRamRpkxxkxx矩阵形式pRRRaaaRpRpRpRRRpRRRp2102120111021求解Yule-Walker方程，就可以得到AR模型系数当模型阶次较大时，直接用矩阵运算求解的计算量大，不利于实时运算。利用系数矩阵的特性，人们提出了如L-D算法等快速算法。四、实验内容编写求解Yule-Walker方程的程序，并对实际生理信号（例如心电、脑电）建立AR模型。对同一数据，使用Matlab信号处理工具箱自带函数aryule计算相同阶数AR模型《生物医学信号处理》实验报告第3页共13页系数，检验程序是否正确。用伪随机序列（白噪声）驱动AR模型，观察输出是否与真实心电、脑电信号相似，对比真实信号与仿真信号的功率谱。五、实验结果与分析试验程序：clear;clc;M=1024;loadecgdata;loadeegdata;loadicpdata;loadrespdata;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));forp=1:15%尝试改变模型阶数，观察效果Rxx=xcorr(x,'biased');%biased为有偏的互相关函数估计Rtemp=zeros(1,p);Rl=zeros(p,1);fork=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];%白噪声方差估计%采用malab自带函数估计模型系数[a,E(p)]=aryule(x,p);%a--系数，E--预测误差，k--反射系数%[a,E(p)]=arburg(x,p);%Y=filter(B,A,X)，输入X为滤波前序列，Y为滤波结果序列，《生物医学信号处理》实验报告第4页共13页%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));xx2=abs(fft(Rx2));figure;subplot(4,1,1);plot(x);xlim([01024]);title('真实数据');subplot(4,1,2);plot(x2);xlim([01024]);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);%求最小均方误差endPopt=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,'-*');gridontitle('(b)预测误差随阶数p的变化情况'),xlabel('p');ylabel('E');%figure,stem(i,Sw,'-*');gridon%title('白噪声方差估计随阶数p的变化情况')%xlabel('p');ylabel('白噪声方差估计');subplot(1,3,3);plot(i,FPE,'-*');title('(c)FPE随阶数p的变化情况');xlabel('p');ylabel('FPE');《生物医学信号处理》实验报告第5页共13页实验结果：1、心电数据，阶数p为15,M=1024:（1）最小均方误差随p的变化情况、预测误差随p的变化情况、FPE随p的变化情况：05101520406080100120140160180200(a)最小均方误差随阶数p的变化情况perror0510150.511.522.53x10-3(b)预测误差随阶数p的变化情况pE0510150.511.522.53x10-3(c)FPE随阶数p的变化情况pFPE01020050100150200250300350400(a)最小均方误差随阶数p的变化情况perror0102000.511.522.53x10-3(b)预测误差随阶数p的变化情况pE0102000.511.522.53x10-3(c)FPE随阶数p的变化情况pFPE图1-1L-D算法图1-2Burg算法（2）白噪声方差估计随阶数p的变化情况：05101500.511.522.53x10-3白噪声方差估计随阶数p的变化情况p白噪声方差估计05101500.511.522.53x10-3白噪声方差估计随阶数p的变化情况p白噪声方差估计图1-3L-D算法图1-4Burg算法（3）真实数据和仿真数据及其功率谱：01002003004005006007008009001000-101真实数据01002003004005006007008009001000-50050仿真数据-1500-1000-5000500100015000100200真实数据功率谱-1500-1000-500050010001500024x104仿真数据功率谱01002003004005006007008009001000-101真实数据01002003004005006007008009001000-50050仿真数据-1500-1000-5000500100015000100200真实数据功率谱-1500-1000-500050010001500012x104仿真数据功率谱图1-5L-D算法error=19.234图1-6Burg算法error=19.913《生物医学信号处理》实验报告第6页共13页2、脑电数据，阶数p为14，M=1024：（1）最小均方误差随p的变化情况、预测误差随p的变化情况、FPE随p的变化情况：010203.13.23.33.43.53.63.73.83.94(a)最小均方误差随阶数p的变化情况perror010200.9511.051.11.151.21.25(b)预测误差随阶数p的变化情况pE0102011.051.11.151.21.25(c)FPE随阶数p的变化情况pFPE010203.23.43.63.844.24.44.64.8(a)最小均方误差随阶数p的变化情况perror010200.9511.051.11.151.21.25(b)预测误差随阶数p的变化情况pE0102011.051.11.151.21.25(c)FPE随阶数p的变化情况pFPE图2-1L-D算法图2-2Burg算法（2）白噪声方差估计随阶数p的变化情况：0246810121400.20.40.60.811.21.4白噪声方差估计随阶数p的变化情况p白噪声方差估计0246810121400.20.40.60.811.21.4白噪声方差估计随阶数p的变化情况p白噪声方差估计图2-3L-D算法图2-4Burg算法（3）真实数据和仿真数据及其功率谱：01002003004005006007008009001000-10010真实数据01002003004005006007008009001000-505仿真数据-1500-1000-5000500100015000100200真实数据功率谱-1500-1000-500050010001500050100仿真数据功率谱01002003004005006007008009001000-10010真实数据01002003004005006007008009001000-505仿真数据-1500-1000-5000500100015000100200真实数据功率谱-1500-1000-50005001000150002040仿真数据功率谱图2-5L-D算法error=3.3147图2-6Burg算法error=3.2043《生物医学信号处理》实验报告第7页共13页3、颅内压数据，阶数p为14，M=1024：(1)最小均方误差随p的变化情况、预测误差随p的变化情况、FPE随p的变化情况：01020050100150200250300350400450500(a)最小均方误差随阶数p的变化情况perror010200.080.090.10.110.120.130.140.150.16(b)预测误差随阶数p的变化情况pE010200.080.090.10.110.120.130.140.150.16(c)FPE随阶数p的变化情况pFPE0102050100150200250300350400450(a)最小均方误差随阶数p的变化情况perror010200.040.050.060.070.080.090.10.110.12(b)预测误差随阶数p的变化情况pE010200.040.050.060.070.080.090.10.110.12(c)FPE随阶数p的变化情况pFPE图3-1L-D算法图3-2Burg算法(2)白噪声方差估计随阶数p的变化情况：0246810121400.020.040.060.080.10.120.140.16白噪声方差估计随阶数p的变化情况p白噪声方差估计0246810121400.020.040.060.080.10.120.140.16白噪声方差估计随阶数p的变化情况p白噪声方差估计图3-1L-D算法图3-2Burg算法(3)真实数据和仿真数据及其功率谱：010020030040050060070080090010000510真实数据01002003004005006007008009001000-50050仿真数据-1500-1000-500050010001500012x104真实数据功率谱-1500-1000-5000500100015000510x104仿真数据功率谱010020030040050060070080090010000510真实数据01002003004005006007008009001000-50050仿真数据-1500-1000-5000500100