利用ARMA、AR、MA模型，以及周期图等进行系统参数估计
N=456;B1=;
A1=;
w=linspace(0,pi,512);
H1=freqz(B1,A1,w);%产生信号的频域响应
Ps1=abs(H1).^2;
SPy11=0;%20次AR(4)
SPy12=0;%20次AR(8)
SPy13=0;%20次AR周期图
SPy14=0;%20次ARMA(4,4)
SPy15=0;%20次ARMA(8,8)
VSPy11=0;%20次AR(4)
VSPy12=0;%20次AR(8)
VSPy13=0;%20次AR周期图
VSPy14=0;%20次ARMA(4,4)
VSPy15=0;%20次ARMA(8,8)
for k=1:20
%采用自协方差法对AR模型参数进行估计%
%gA1：AR模型的参数；gE1：激励白噪声的方差%
y1=filter(B1,A1,randn(1,N)).*;
=pcov(y1,4,512,1);%AR(4)的估计%
=pcov(y1,8,512,1);%AR(8)的估计%
=periodogram(y1,[],512,1);
SPy11=SPy11+Py11;
SPy12=SPy12+Py12;
SPy13=SPy13+Py13;
VSPy11=VSPy11+abs(Py11).^2;
VSPy12=VSPy12+abs(Py12).^2;
VSPy13=VSPy13+abs(Py13).^2;
figure(1)
plot(w./(2*pi),Ps1,F,Py11);
legend('真实功率谱','20次AR(4)估计图');
hold on;
figure(2)
plot(w./(2*pi),Ps1,F,Py12);
legend('真实功率谱','20次AR(8)估计图');
hold on;
figure(3)
plot(w./(2*pi),Ps1,F,Py13);
legend('真实功率谱','20次周期图法估计图');
hold on;
%------------ARMA模型---------------%
y=zeros(1,256);
for i=1:256
y(i)=y1(200+i);
end
ny=;
z=fliplr(y);nz=-fliplr(ny);
nb=ny(1)+nz(1);ne=ny(length(y))+nz(length(z));
n=;
Ry=conv(y,z);
R4=zeros(8,4);%ARMA(4,4)的R
r4=zeros(8,1);%ARMA(4,4)的r
for i=1:8
r4(i,1)=-Ry(260+i);
for j=1:4
R4(i,j)=Ry(260+i-j);
end
end
R4%R矩阵
r4%r矩阵
a4=inv(R4'*R4)*R4'*r4%利用最小二乘法得到的a的估计参数
%----------------ARMA(4,4)对MA的参数b(1)-b(4)进行估计----------------------%
A1
A14=%AR的参数a(1)-a(4)的估计值
B14=fliplr(conv(fliplr(B1),fliplr(A14)));%MA模型的分子
y24=filter(B14,A1,randn(1,N));%.*;%由估计出的MA模型产生数据
%---因为(q=4)<<L<<N=256,所以选取L=32---%
=arburg(y24,32),%利用数据y2估计AR(32)的参数
B1
b4=arburg(Ama4,4)%求出MA模型的参数
%---求功率谱---%
w=linspace(0,pi,512);
%H1=freqz(B1,A1,w)
H14=freqz(b4,A14,w);%产生信号的频域响应
%Ps1=abs(H1).^2;%真实谱
Py14=abs(H14).^2;%估计谱
%if Py14>200
% PPy14=200;
%elseif Py14<200
% PPy14=Py14;
%end
SPy14=SPy14+Py14;
VSPy14=VSPy14+abs(Py14).^2;
figure(4)
plot(w./(2*pi),Ps1,w./(2*pi),Py14);
legend('真实功率谱','20次ARMA(4,4)的估计图');
hold on;
R8=zeros(16,8);%ARMA(8,8)的R
r8=zeros(16,1);%ARMA(8,8)的r
for i=1:16
r8(i,1)=-Ry(264+i);
for j=1:8
R8(i,j)=Ry(264+i-j);
end
end
R8%R矩阵
r8%r矩阵
a8=inv(R8'*R8)*R8'*r8%利用最小二乘法得到的a的估计参数
%----------------ARMA(8,8)对MA的参数b(1)-b(8)进行估计----------------------%
A1
A18=%AR的参数a(1)-a(8)的估计值
B18=fliplr(conv(fliplr(B1),fliplr(A18)))%MA模型的分子
y28=filter(B18,A1,randn(1,N));%.*;%由估计出的MA模型产生数据
%---因为(q=8)<<L<<N=256,所以选取L=48---%
=arburg(y28,48),%利用数据y2估计AR(32)的参数
B1
b8=arburg(Ama8,8)%求出MA模型的参数
%---求功率谱---%
%w=linspace(0,pi,512);
%H1=freqz(B1,A1,w)
H18=freqz(b8,A14,w);%产生信号的频域响应
%Ps1=abs(H1).^2;%真实谱
Py15=abs(H18).^2;%估计谱
%if Py15>200
% PPy15=200;
%elseif Py15<200
% PPy15=Py15;
%end
SPy15=SPy15+Py15;
VSPy15=VSPy15+abs(Py15).^2;
figure(5)
plot(w./(2*pi),Ps1,w./(2*pi),Py15);
legend('真实功率谱','20次ARMA(8,8)的估计图');
hold on;
%-----------------------------------%
end
V1=VSPy11/20-abs(SPy11/20).^2;
V2=VSPy12/20-abs(SPy12/20).^2;
V3=VSPy13/20-abs(SPy13/20).^2;
V4=VSPy14/20-abs(SPy14/20).^2;
V5=VSPy15/20-abs(SPy15/20).^2;
figure(6)
plot(w./(2*pi),Ps1,F,SPy11/20);
legend('真实功率谱','20次AR(4)估计的均值');
figure(7)
plot(w./(2*pi),Ps1,F,SPy12/20);
legend('真实功率谱','20次AR(8)估计的均值');
figure(8)
plot(w./(2*pi),Ps1,F,SPy13/20);
legend('真实功率谱','20次周期图估计的均值');
figure(9)
plot(w./(2*pi),Ps1,w./(2*pi),SPy14/20);
legend('真实功率谱','20次ARMA(4,4)估计的平均值');
figure(10)
plot(w./(2*pi),Ps1,w./(2*pi),SPy15/20);
legend('真实功率谱','20次ARMA(8,8)估计的平均值');
figure(12)
plot(F,V1,F,V2,F,V3,w./(2*pi),V4,w./(2*pi),V5);
legend('AR(4)方差','AR(8)方差','周期图方差','ARMA(4,4)方差','ARMA(8,8)方差');
axis();
[ 本帖最后由 eight 于 2007-2-16 16:38 编辑 ] suffer主任好：
请问，arma是不是自回归滑动平均模型，我现在在搞卡尔曼滤波，需要用arma来估计系统噪声，请问，有没有比较基础的定义，基础的资料，这个对我很重要，期待您的帮助！
谢谢您共享的资料！！
[ 本帖最后由 suffer 于 2006-10-9 19:24 编辑 ] 找一本时间序分析列的术看一看 <P>看看现代信号处理，谱估计里面就有关于arma的详细的讲述。</P> 请问：输入的已知条件是什么？？ 非常感谢，我也是刚了解了很少关于这几个函数的应用，这下可以借鉴一下了 有没有可靠性评估的参数估计程序啊？通过已知的数据对可靠性模型（G-O模型，Musa模型等）进行参数估计，有哪位高人用过VC＋＋中运用Matlab编写过此类程序，给小妹指点一下吧。 GANXIE {:{44}:} 谢谢了。。。会好好看的。。。 请问这个例子中的B1 A1的物理含义是什么的？谢谢！ 还没看懂，不过人会只是时间的长短问题
页:
[1]