Yule-Walker方程实验报告

1、评分大理大学实验报告课程名称生物医学信号处理实验名称Yu le W al ke r方程专业班级姓 名羽卒兰cl学号实验日期2016年5月27日实验地点2015- 2 01 6学年度第学期、实验目得学习求解Y ul e -Walke r方程,建立随机信号得A R模型。二、实验环境1、 硬件配置:In tel（R） Core（TM）i5 -42 1 0U CPU 1、7 G Hz 1、7GHz 安装内存（RAM） : 4、00G B系统类型:64位操作系统2、软件环境：M ATLAB R2013b软件三、实验内容（包括本实验要完成得实验问题及需要得相关知识简单概述）编写求解Yul e Walke

2、r方程得程序，并对实际生理信号（例如心电、脑电）建立A R模型。对同一数据，使用M atla b信号处理工具箱自带函数aryul e计算相同阶数 AR模型系数，检验程序就是否正确用伪随机序列（白噪声）驱动AR模型，观察输出就是否与真实心电、脑电信号相似，对比真实信号与仿真信号得功率谱。四、实验结果与分析（包括实验原理、数据得准备、运行过程分析、源程序（代码）、图形图象界面等）实验原理随机信号可以瞧作就是由当前激励白噪声w（n ）以及若干次以往信号 x （nk）得线性组合产生，即所谓自回归模型（AR模型）模型参数满足 Yule W alker方程矩阵形式求解Yul e Walke r方程,就可以

3、得到A R模型系数 当模型阶次较大时，直接用矩阵运算求解得计算量大， 人们提出了如LD算法等快速算法。源程序：c lear ; cl c;M =1024;loa d ecgda t a; x不利于实时运算禾U用系数矩阵得特性%1 oa d ee g d a t a; x% load ic p d a ta;x = i c%loa d r e spd atp=1:60;Sw =zeros(1, 1 e n得Sw零矩阵或数组E=ze ro s (1 ,l e零矩阵或数组FPE=zeros(1 , l e ngth (p);矩阵或数组for i=1: 6 0=ec g data( 1:1 02 4

4、);=ee g dat a (1:M ); pdata (1:10 2 4);r espd a ta (1:102a; xg th(p)ngth(p )；%导入心电信号得数据%导入脑电信号得数据%导入颅内压信号得数据4) ; %导入个呼吸信号得数据 %P得取值范围%创建一个1行列长为le%创建一个1行列长为len%创建一个1行列长为1 en g thngth ( p )%尝试改变模型阶数，观察效果gth(p )得 E(p)得FPE零Rx x = x c orr (x，/ b ias 关函数估计；R temp = z eros(1, i );R1 = ze ros(i , 1);组for k =

5、 1:1 engt h( Rtemp)Rt emp( k ) = Rxx (length(x) 1+k);Rl( k) = Rxx (lengt h (x)+k) ;%e n dRs = toep 1 it z (Rtem p );A = inv(R s ) * Rl;Sw( i) =: Rtemp (1), R l 1; A;% 采用malab自带函数估计模型系数a,E (i ) = aryule(x , i);估计,ed );%估计随机过程中得互相关序列,参数估计d a = a(2 : e n d) A;FP E (i ) =E( i )* ( M+ i+1 ) / ( M-i 1w =

6、r andn(si z e (x);x 2= fi 1t er (1, a ,w);%a, E (i)= ar bu rg (x,i );en d” biased为有偏得互相%创建一个i行1列得Rtemp零矩阵或数组 %创建一个i行1列得Rl零矩阵或数%f o r 循环从 1 至U length( R temp)%取 Rxx (0)至U Rxx(p-1 )赋值取 R xx( 1)至U Rx x (p)赋值%生成自相关系数矩阵(Toeplitz型)%AR模型系数估计% 白噪声方差估计% malab实现L -D算法得AR模型参数a-系数，E-预测误差，k -反射系数%malab实现Bur g算法得

7、A R模型%自编程序求解就是否正确？% FPE算法%生成随机噪声%仿真数据fi g u re%画图sub p l o t(3, 1, 1),plot (p/ ,E, * ) , gr i d on;%创建窗口画图，并添加网格线title( / E 随阶数 p 得变化情况);xlabe 1 ( p ) ; y l abe l ( e rr or ) ; %添加标题与 横纵坐标s u bp 1 ot(3,1 , 2) , p lot( p ,Sw, / -*), g r i d on;%创建窗口画图，并添加网格线ti tle( / Sw随阶数p得变化情况)；xlab el ( pz) ;ylab

8、e l (白噪声方差估计);subplo t ( 3 , 1, 3),plot(p ,F P E, ) , g rid on；%创建窗口画图,并添加网格线t it le (FPE随阶数p得变化情况);xlabe 1 (/ p); ylabe 1( F EPZ );%添加标题与横 纵坐标心电信号：x 10-3E随阶数p的变化情况0010503pSw随阶数p的变化情况X 10-3210210x 10-3pFPE 随阶数p的变化情况p-332a) L D算法0仃 4- + 4- 4- 4 4- -14 4- 4- * * *1=牛4*4*4*x 10pSw随阶数p的变化情况10102030pFPE

9、随阶数p的变化情况4050603*1-3x 10210 10203040p50060b)Burg 算法 图1心电信号得不同算法得E、Sw FPE随阶数p得变化图分析：首先,根据FPE准则找到最佳阶数P;LD算法；在 ma nd wi ndow 输入 fi n d (FPE= mi n ( FP E)，结果：ans= 1 5;在 mand wi ndo w 输入f i nd (abs ( E1)= = mi n(abs(E 1),结果：a ns=1; B u rg算法：在 mand wi n d o w 输入 f ind(F P E= =mi n (FPE)，结果：an s =31;在m and

10、 wi n dow 输入 fi n d (abs (E 1)=min( a bs (E -1),结果：ans=1。脑电信号0.9W- -a1(E随阶数p的变化情况PSw随阶数p的变化情况0.9h.一 一舟PFPE 随阶数p的变化情况1.4a) LD算法1.6 *+ + *p * * * *r* H-4- + *斗P*+- +*+*卜*蒼* * * *秸廿*E随阶数p的变化情况况 情 化 变 的 数 阶 随计肇、書噪白况 情 化 变 的 数 阶 随4- 4-1-亠士 rr * .丄 J, JL 厶 a A20b) B u r g算法图2脑电信号得不同算法得E、Sw FPE随阶数p得变化图分析：首

11、先,根据F PE准则找到最佳阶数 P；L D算法：在 mand wi nd o w 输入f in d(F PE = = min( FP E),结果：a n s =2 3;在m and w i nd ow 输入 fin d(a b s (E-1) = min( a b s (E- 1 ),结果：a ns= 1 1; Burg算法:在 ma nd wi n d o w 输入 fin d(FPE= = mi n(FPE),结果:a n s=4 3 ;在 man d wi ndow 输入 fi n d(ab s (E 1) =min (abs(E 1)，结果：an s=11. 颅内压信号：-计#近囁白

12、 a) L- D算法0.15P 0.10.050Pb) Bu r g算法图3颅内压信号得不同算法得E、Sw、FPE随阶数p得变化图分析：首先，根据FP E准则找到最佳阶数P；L-D算法：在 ma nd w in dow 输入 f i n d(F PE= =mi n( F PE),结果：an s=43;在m and w in dow 输入 f i n d (abs(E-1 )= = mi n (abs( E 1),结果：ans=1;Bur g算法：在 ma nd wi n dow 输入 fin d(FPE= = mi n(FPE),结果：a ns=57;在 mand windo w 输入 f i

13、nd ( abs(E- 1)=m i n(a bs(E 1),结果:ans =1呼吸信号：0102030405060Pb) Burg算法图4呼吸信号得不同算法得E、Sw FPE随阶数p得变化图分析：首先，根据FPE准则找到最佳阶数 P;L -D算法：在 ma nd wi ndow输入 f i n d(FP E =mi n( FP E)，结果：ans =59;在 ma nd w i ndow 输入 find(ab s( E-1 ) =mi n (abs(E- 1 ),结果：ans = 60;Bu r g算法：在 mand wi n dow 输入 find(FP E =min (F PE),结果：

14、a n s= 44;在m and wi ndow 输入 fi nd (abs (E 1) = =mi n ( abs ( E 1),结果：an s =60。 分析：L- D算法与E urg算法得用find (abs (E -1 ) =mi n(abs( E -1 )找到得值就是一样 得；心电、脑电、颅内压信号得找find(F P E=min(FPE)得值Burg算法比LD算法大，呼吸信号找到得 mi n(F PE)值，B u r g算法比L D算法小。思考题：c lear; c1 c;1 oad ec gdata;x =ecgdata (1 :1024 );%导入心电信号%lo ade egd

15、ata ; x= eegdata(1:102 4 );%导入脑电信号%1 oadicpdata; x =icpd a t a(1:1024 );%导入颅内压信号%1 oad res p da t a;x=respd ata (1:1024 );%导入个呼吸信号%p =20;%尝试改变模型阶数，观察效果f o r p=2:2:20%尝试改变模型阶数，观察效果R xx = xcorr(x,biased ;%估计随机过程中得互相关序列，” bi a s e d为有偏得互相关函数估计；Rt emp = zeros (1, p);%创建一个i行1列得Rtem p零矩阵或数组R1 = zero s( p,

16、 1);%创建一个i行1列得Rl零矩阵或数组fork = 1: 1e ng t h(Rtem p )%f o r 循环从 1 至U l e ngth( R temp )R temp (k) = Rxx( le n g t h(x) 1+k);%取 Rxx (0 )至 Rxx(p 1)赋值R1 (k )= Rxx (l en g t h(x ) +k);%取 Rxx( 1)至U Rx x(p)赋值endRs = toepli t z(Rtemp);%生成自相关系数矩阵(Toeplitz 型A = -inv ( Rs)*R l ;% AR模型系数估计Sw( p/2 )= Rtemp(1) , R

17、l * 1;A ;%白噪声方差估计%米用malab自带函数估计模型系数a,E = ary u le(x,p);%a 系数，E预测误差，k 反射系数%a,E = arburg( x , p );% malab实现Bu r g算法得AR模型参数估计da = a (2:end) A%自编程序求解就是否正确 ？stem (d a) ; ti t le(参数估计偏差)%画图w = r an d n (s ize(x );%产生随机噪声信号x2 = fi lt e r (1,a , w);%仿真数据fig u re ;subp 1 ot (2,1,1);plo t( x) ;tit 1 e (真实数据);

18、%绘制真实数据得图像subplot (2,1 , 2);pl o t (x2) ;t i t le(仿真数据);绘制仿真数据得图像error (p/ 2) =m ean (x -x2)、A 2 ); end%显示最小均方误差得计算结果f ig u r esub plot (1,2 , 1) , pl o t(2:2 : 2 0, e rror ti t I eC最小均方误差随阶数p得变化情况-* )%画图);x label(p ) ;y lab el(error );%绘制最小均方误差随阶数p得变化情况图subplot (1, 2,2 ),stem(2:2: 20,Sw, * ) , grid

19、 o nti tl e (白噪声方差估计随阶数p得变化情况);x la b el(p ); y label(白噪声方差Rxx 2=x co r r(x2, / b iased );Px = a bs( f ft(Rxx);Px2 = a bs(ff fi gur e subplo t (2,1 火柴杆图 subp l ot(2, 杆图%绘制S w、E随p变化散点图，以下就是统计得到得结果Sw= 1、0527 1、0 301 1、015098 7 2;E= 3、874 14 176:; figu r e s u bp l ot(2,1, 1);st e m( p, Sw);xlabe l ( p

20、 );ylab el ( Sw );ti su bp lot (2 ,1 , 2 );pl ot (p, E, xlabel( p );ylt(Rx x2);9绘制白噪声方差估计随阶数p得变化情况图%求仿真数据得自相关函数计算真实信号自功率谱%计算仿真信号自功率谱,1) ; s t em(1023:1 023, Px ) ; titl1 ,2 ) ; s tem ( 1 0 23:1023,Px2 ) ; t i3、4 3 37 3、32000、 9 9610、3、7 3 143、e （真实信号功率谱）；绘制le （仿真信号功率谱）；p= 8 9 1 01 1 12 13 14994 3 0、

21、99 4 2 0、 9 896 0、%绘制火柴151 6;9 877 0、3683 3、4809 3、5826 3、62 5 33、tl e(Sw随p变化散点图 /o)；a bel(E ) ; t i t le( E 随 p 变化散点图);1、比较四种信号得参数估计偏差10-5-1001(真实数据参数估十偏差02468101214161820-16x 10 1520040060080010001200真实数据%绘制火柴杆图）；%添加标题与横纵坐标%绘制散点图%添加标题与横纵坐标真实数据x 101参数估计偏差, 一一 丨 l L , r II.J 11-11 1x 1011参数估计偏差02468

22、101214161820h 11II102468101214161820估计);a）心电b）脑电c）颅内压d ）呼吸图5心电、脑电、颅内压、呼吸信号得得真实数据与参数估计偏差分析：由图5可以瞧出脑电信号得参数估计偏差比其它三个信号得参数估计偏差要小许多，脑电信号得自编程序跟MATLAB信号处理工具箱自带函数ary u le得处理更接近。由此可知丄-D算法与自编程序相比较，自编程序对估计得参数比较精确一些 2、比较四种信号得真实数据与仿真数据真实数据0 撫加11 “-5-100屮严八帥啊壮也仿真数据1050-5a）心电105真实数据真实数据)0200400600HoIlli, Pill-1. i

23、f：|H 贋曲曲删200400600800100000c） 颅内压真实数据d） 呼吸b）脑电图6心电、脑电、颅内压、呼吸信号得得真实数据与仿真数据分析：从图6可以瞧出，对心电、脑电、颅内压、呼吸信号建模后，心电、颅内压、呼吸信号得真实数据与仿真数据相差很大，脑电模型产生得信号更能真实得反映脑电信号得特征。不同噪声得激励得到得信号时不同得。3、比较四种信号L- DL D算法：真实信号功率谱算法与Bur g 算法得功率谱150真实信号功率谱150100仿真言号功率谱40真实信号功率谱a）心电b）脑电c）颅内压d）呼吸B ur g算法真实信号功率谱仿真信号功率谱真实信号功率谱8X 10仿真信号功率谱

24、a）心电b） 脑电c） 颅内压图7心电、脑电、颅内压、呼吸信号得LD算法与B u r g分析：由图7瞧出，在比较心电、脑电、颅内压、呼吸信号得L-D 算法与d）呼吸算法功率谱Bu r g算法功率谱发现,对信号进行功率谱估计时，脑电与颅内压信号得功率谱得纵坐标相差不大，但就是相对而说心电与呼吸信号得相差就是极大得。4、比较四种信号得L- D 算法与B u rg算法L D算法240最小均方误差随阶数|的变化情况220200180160140120100p1X，0白噪声方差估计随阶数p的变化情况0.90.80.70.6计估 差方5声噪白0.40.30.20.10p最小均万误差 随阶数p的变化情况白

25、噪声万 差估计 随阶数P的变化情况1.48 6计估差方声噪白rpea）心电b ）脑电况rDFe况计善方士陽白况 情c）颅内压图8四种信号LD算法，阶数Burg 算法d ）呼吸p与均方误差erro r与噪声方差估计值Sw之间得关系最小均方误差随阶数p的变化450恋差万声噪白7 6 5 a a a噪声10方差估计随阶数p的1变化情况10p20计估差方声噪白况数 阶 随 差 误 方 均 最况 情 化 变 的数 阶 随 计 估 差 方 声 噪齐白况 情 化 变 的最小X均0方误差随阶数p的变化情况rare58275227淸白80噪声方差估计随阶数p的变化情况16001400啦1000眈60040020

26、020c）颅内压d）呼吸图9四种信号Bur g算法，阶数p与均方误差er r o r与噪声方差估计值 Sw之间得关系分析：由图8与图9可以瞧出，虽然四种信号得L D算法与 B urg算法阶数p与均方误差error得变化就是不稳定得，但就是它们得阶数p与噪声方差估计值 Sw之间得变化就是大体一致得。4、1心电信号得L-D 算法与 Bu r g 算法比较真实数据真实数据仿真数据b)Bu r g算法a)L D 算法 b)B urg算法真实数据真实数据仿真数据/ 1 *1 I宀rJI,-I1 f1 riiififif1618a)L D算法Da值L-D-8 、-5 、-1 、-8 、-2 、4、2 、3

27、、5、3、算法990e-223e-340e-120e-144e-546e-100e-596e-641e-163e-13121112121211111111-1 、3、3、3、3、3、1 、-5、-3 、3、130e-8190e5763e2489e8141e1735e5867e012e-589e-7491e11-12-11-11-11-11-111212-14Burg算法-0 、2810、499-0 、092-0 、2750、1030、0610、026-0 、0970、0530、027-0 、-0 、0 、0 、-0 、0 、0 、-0 、0 、-0 、043501710117106914150

28、4950247057306700284图10心电信号不同算法得真实数据与参数估计图11心电信号不同算法得真实与仿真数据表1心电信号得L D算法与Burg 算法得D a值变化表4、2脑电信号得LD算法与B urg算法比较真实数据1050-5-10真实数据1(210，-1-21 /I i J,11 1 1 1B urg算法020040060080010001200x 10-3参数估计偏差1(0200400600 800 1000 1200020040060080010001200仿真数据5(汕神斛-11IIII020040060080010001200仿真数据真实数据a)L-D算法 法b)a)L

29、D算法b ) Bur g 算图1 2脑电信号不同算法得真实数据与参数估计图1 3脑电信号不同算法得真实与仿真数据表2脑电信号得L -D算法与B u rg算法得Da值变化表Da值L-D5、1 、-6 、0-1 、-1 、9、4、-8 、0算法5511e3877e939e-110e-665e-7144e8572e3266-16-17181616-17-1717-8 、2 、-2 、1 、-2 、1 、-1 、-2 、1 、5、326e-1510e082e-3877e168e-2490e145e-081e-1102e5511e17-1617-1716-161617-16-17Burg0 、-0 、-

30、0 、-0 、-0 、0 、0 、-0 、-0 、-0 、算法00090003000600040014000400130003000600050 、-0 、0 、0 、0 、0 、0 、-0 、-0 、0 、00050017000300090005001500090001001700024、3颅内压信号得L D算法与 B ur g算法比较真实数据600o1 1 1 1 1 0 _ III, 1 . Il1 1 |2x 10-11参数估计偏差1真实数据6真实数据，-iliH： 1r.|：,II f ； i：:?：. ；i； - K2004T0600 800 1000200 0仿真数据20仿真数据

31、j .Mrr1i |-2-4200200 .00 600 800 “00 200a)L D算法b) Burg算法a)L D算法b) Burg算Da值L-D3、1 、1 、4、-1 、1 、4、3、-3 、4、算法1286e4018e1497e9960e567e-4885e3890e5907e099e-4496e-13-12-12-1512-13-12-1212-128 、1 、-7 、-6 、3、1 、-1 、6 、4、2 、5008e0558e490e-810e-4846e3522e628e-4571e0563e9531e-12-121212-12-1211-12-12-14Burg-0 、

32、0 、0 、0 、-0 、-0 、-0 、-0 、0 、0 、算法16210853105703710042080301450308002202500 、0 、-0 、0 、0 、0 、0 、-0 、0 、-0 、0356258239560479076703840087049401690025法图14颅内压信号不同算法得真实数据与参数估计图15颅内压信号不同算法真实与仿真数据表3颅内压信号得L D算法与Bur g 算法得Da值变化表4、4呼吸信号得L-D 算法与Bur g算法比较100050000A h A1 l 1 V I X hn- y yV i真实数据真实数据0111 i 11iJx 10

33、-11参数估计偏差402468101214161820x 10-3参数估计偏差真实数据真实数据Da值L-D-2 、3、8 、4、-1 、-7、-2 、-9 、-1 、-1 、算法735e-042e-882e-431e-890e-866e-557e-929e-172e-041e-131312121213121311112 、8 、8 、1 、4、1 、6 、-5 、-2 、-6 、6015e3792e1914e0824e7703e5861e6561e402e-922e-889e-11-12-12-11-12-11-12121214Burg-0 、-0 、-0 、0 、0 、-0 、0 、0 、0

34、 、0 、算法01170006000800010006000600110012000600130 、0 、0 、0 、0 、0 、-0 、0 、-0 、0 、0019000700160005001100090007000700040007a)L D算法b )Burg算法图16呼吸信号不同算法得真实数据与参数估计a) L - D算法b)B u rg算法图17呼吸信号不同算法得真实与仿真数据表4呼吸信号得L-D算法与Burg算法得Da值变化表分析:Burg算法得优点就是：求得得AR模型就是稳定得，较高得计算效率，但递推还就是用得LD算法，因此仍然存在明显得缺点5、脑电信号,改变P值，比较L D算法与Bu rg算法Sw 随p变化散点图P1.5Sw 随 p变化散点图0EpE随 p变化散点 图Pa) L- D算法图1 8脑电信号改变P值，比较b) Bu rg算法 L D算法与Burg算法Sw与E值变化图P值8910111213141516L-D算法Sw1、05271、03011、0150099610、99430、99420、98960、98770、9872E3、46873、20103、8530338883、28743、39453、93154、15733、3942