声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 13367|回复: 11

[共享资源] 利用ARMA、AR、MA模型,以及周期图等进行系统参数估计

[复制链接]
发表于 2006-5-9 21:38 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
  1. N=456;
  2. B1=[1 0.3544 0.3508 0.1736 0.2401];
  3. A1=[1 -1.3817 1.5632 -0.8843 0.4096];
  4. w=linspace(0,pi,512);
  5. H1=freqz(B1,A1,w);%产生信号的频域响应
  6. Ps1=abs(H1).^2;
  7. SPy11=0;%20次AR(4)
  8. SPy12=0;%20次AR(8)
  9. SPy13=0;%20次AR周期图
  10. SPy14=0;%20次ARMA(4,4)
  11. SPy15=0;%20次ARMA(8,8)
  12. VSPy11=0;%20次AR(4)
  13. VSPy12=0;%20次AR(8)
  14. VSPy13=0;%20次AR周期图
  15. VSPy14=0;%20次ARMA(4,4)
  16. VSPy15=0;%20次ARMA(8,8)
  17. for k=1:20
  18. %采用自协方差法对AR模型参数进行估计%
  19. %gA1:AR模型的参数;gE1:激励白噪声的方差%
  20. y1=filter(B1,A1,randn(1,N)).*[zeros(1,200),ones(1,256)];
  21. [Py11,F]=pcov(y1,4,512,1);%AR(4)的估计%
  22. [Py12,F]=pcov(y1,8,512,1);%AR(8)的估计%
  23. [Py13,F]=periodogram(y1,[],512,1);
  24. SPy11=SPy11+Py11;
  25. SPy12=SPy12+Py12;
  26. SPy13=SPy13+Py13;
  27. VSPy11=VSPy11+abs(Py11).^2;
  28. VSPy12=VSPy12+abs(Py12).^2;
  29. VSPy13=VSPy13+abs(Py13).^2;
  30. figure(1)
  31. plot(w./(2*pi),Ps1,F,Py11);
  32. legend('真实功率谱','20次AR(4)估计图');
  33. hold on;
  34. figure(2)
  35. plot(w./(2*pi),Ps1,F,Py12);
  36. legend('真实功率谱','20次AR(8)估计图');
  37. hold on;
  38. figure(3)
  39. plot(w./(2*pi),Ps1,F,Py13);
  40. legend('真实功率谱','20次周期图法估计图');
  41. hold on;
  42. %------------ARMA模型---------------%
  43. y=zeros(1,256);
  44. for i=1:256
  45. y(i)=y1(200+i);
  46. end
  47. ny=[0:255];
  48. z=fliplr(y);nz=-fliplr(ny);
  49. nb=ny(1)+nz(1);ne=ny(length(y))+nz(length(z));
  50. n=[nb:ne];
  51. Ry=conv(y,z);
  52. R4=zeros(8,4);%ARMA(4,4)的R
  53. r4=zeros(8,1);%ARMA(4,4)的r
  54. for i=1:8
  55. r4(i,1)=-Ry(260+i);
  56. for j=1:4
  57. R4(i,j)=Ry(260+i-j);
  58. end
  59. end
  60. R4%R矩阵
  61. r4%r矩阵
  62. a4=inv(R4'*R4)*R4'*r4%利用最小二乘法得到的a的估计参数
  63. %----------------ARMA(4,4)对MA的参数b(1)-b(4)进行估计----------------------%
  64. A1
  65. A14=[1,a4']%AR的参数a(1)-a(4)的估计值
  66. B14=fliplr(conv(fliplr(B1),fliplr(A14)));%MA模型的分子
  67. y24=filter(B14,A1,randn(1,N));%.*[zeros(1,200),ones(1,256)];%由估计出的MA模型产生数据
  68. %---因为(q=4)<<L<<N=256,所以选取L=32---%
  69. [Ama4,Ema4]=arburg(y24,32),%利用数据y2估计AR(32)的参数
  70. B1
  71. b4=arburg(Ama4,4)%求出MA模型的参数
  72. %---求功率谱---%
  73. w=linspace(0,pi,512);
  74. %H1=freqz(B1,A1,w)
  75. H14=freqz(b4,A14,w);%产生信号的频域响应
  76. %Ps1=abs(H1).^2;%真实谱
  77. Py14=abs(H14).^2;%估计谱
  78. %if Py14>200
  79. % PPy14=200;
  80. %elseif Py14<200
  81. % PPy14=Py14;
  82. %end
  83. SPy14=SPy14+Py14;
  84. VSPy14=VSPy14+abs(Py14).^2;
  85. figure(4)
  86. plot(w./(2*pi),Ps1,w./(2*pi),Py14);
  87. legend('真实功率谱','20次ARMA(4,4)的估计图');
  88. hold on;
  89. R8=zeros(16,8);%ARMA(8,8)的R
  90. r8=zeros(16,1);%ARMA(8,8)的r
  91. for i=1:16
  92. r8(i,1)=-Ry(264+i);
  93. for j=1:8
  94. R8(i,j)=Ry(264+i-j);
  95. end
  96. end
  97. R8%R矩阵
  98. r8%r矩阵
  99. a8=inv(R8'*R8)*R8'*r8%利用最小二乘法得到的a的估计参数
  100. %----------------ARMA(8,8)对MA的参数b(1)-b(8)进行估计----------------------%
  101. A1
  102. A18=[1,a8']%AR的参数a(1)-a(8)的估计值
  103. B18=fliplr(conv(fliplr(B1),fliplr(A18)))%MA模型的分子
  104. y28=filter(B18,A1,randn(1,N));%.*[zeros(1,200),ones(1,256)];%由估计出的MA模型产生数据
  105. %---因为(q=8)<<L<<N=256,所以选取L=48---%
  106. [Ama8,Ema8]=arburg(y28,48),%利用数据y2估计AR(32)的参数
  107. B1
  108. b8=arburg(Ama8,8)%求出MA模型的参数
  109. %---求功率谱---%
  110. %w=linspace(0,pi,512);
  111. %H1=freqz(B1,A1,w)
  112. H18=freqz(b8,A14,w);%产生信号的频域响应
  113. %Ps1=abs(H1).^2;%真实谱
  114. Py15=abs(H18).^2;%估计谱
  115. %if Py15>200
  116. % PPy15=200;
  117. %elseif Py15<200
  118. % PPy15=Py15;
  119. %end
  120. SPy15=SPy15+Py15;
  121. VSPy15=VSPy15+abs(Py15).^2;
  122. figure(5)
  123. plot(w./(2*pi),Ps1,w./(2*pi),Py15);
  124. legend('真实功率谱','20次ARMA(8,8)的估计图');
  125. hold on;
  126. %-----------------------------------%
  127. end
  128. V1=VSPy11/20-abs(SPy11/20).^2;
  129. V2=VSPy12/20-abs(SPy12/20).^2;
  130. V3=VSPy13/20-abs(SPy13/20).^2;
  131. V4=VSPy14/20-abs(SPy14/20).^2;
  132. V5=VSPy15/20-abs(SPy15/20).^2;
  133. figure(6)
  134. plot(w./(2*pi),Ps1,F,SPy11/20);
  135. legend('真实功率谱','20次AR(4)估计的均值');
  136. figure(7)
  137. plot(w./(2*pi),Ps1,F,SPy12/20);
  138. legend('真实功率谱','20次AR(8)估计的均值');
  139. figure(8)
  140. plot(w./(2*pi),Ps1,F,SPy13/20);
  141. legend('真实功率谱','20次周期图估计的均值');
  142. figure(9)
  143. plot(w./(2*pi),Ps1,w./(2*pi),SPy14/20);
  144. legend('真实功率谱','20次ARMA(4,4)估计的平均值');
  145. figure(10)
  146. plot(w./(2*pi),Ps1,w./(2*pi),SPy15/20);
  147. legend('真实功率谱','20次ARMA(8,8)估计的平均值');
  148. figure(12)
  149. plot(F,V1,F,V2,F,V3,w./(2*pi),V4,w./(2*pi),V5);
  150. legend('AR(4)方差','AR(8)方差','周期图方差','ARMA(4,4)方差','ARMA(8,8)方差');
  151. axis([0 0.5 0 2000]);
复制代码

[ 本帖最后由 eight 于 2007-2-16 16:38 编辑 ]

评分

1

查看全部评分

回复
分享到:

使用道具 举报

发表于 2006-5-11 18:57 | 显示全部楼层
suffer主任好:
请问,arma是不是自回归滑动平均模型,我现在在搞卡尔曼滤波,需要用arma来估计系统噪声,请问,有没有比较基础的定义,基础的资料,这个对我很重要,期待您的帮助!
谢谢您共享的资料!!

[ 本帖最后由 suffer 于 2006-10-9 19:24 编辑 ]
发表于 2006-5-13 20:35 | 显示全部楼层
找一本时间序分析列的术看一看
发表于 2006-5-23 15:50 | 显示全部楼层
<P>看看现代信号处理,谱估计里面就有关于arma的详细的讲述。</P>
发表于 2006-6-14 10:42 | 显示全部楼层
请问:输入的已知条件是什么??
发表于 2006-8-21 00:35 | 显示全部楼层
非常感谢,我也是刚了解了很少关于这几个函数的应用,这下可以借鉴一下了
发表于 2006-8-25 16:04 | 显示全部楼层
有没有可靠性评估的参数估计程序啊?通过已知的数据对可靠性模型(G-O模型,Musa模型等)进行参数估计,有哪位高人用过VC++中运用Matlab编写过此类程序,给小妹指点一下吧。
发表于 2008-2-22 20:48 | 显示全部楼层
GANXIE
发表于 2011-2-9 18:24 | 显示全部楼层
发表于 2011-4-11 16:38 | 显示全部楼层
谢谢了。。。会好好看的。。。
发表于 2012-5-9 22:43 | 显示全部楼层
请问这个例子中的B1 A1的物理含义是什么的?谢谢!
发表于 2017-3-21 16:23 | 显示全部楼层
还没看懂,不过人会只是时间的长短问题
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-17 07:10 , Processed in 0.066035 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表