马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
- N=456;
- B1=[1 0.3544 0.3508 0.1736 0.2401];
- A1=[1 -1.3817 1.5632 -0.8843 0.4096];
- 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)).*[zeros(1,200),ones(1,256)];
- [Py11,F]=pcov(y1,4,512,1);%AR(4)的估计%
- [Py12,F]=pcov(y1,8,512,1);%AR(8)的估计%
- [Py13,F]=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=[0:255];
- z=fliplr(y);nz=-fliplr(ny);
- nb=ny(1)+nz(1);ne=ny(length(y))+nz(length(z));
- n=[nb:ne];
- 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=[1,a4']%AR的参数a(1)-a(4)的估计值
- B14=fliplr(conv(fliplr(B1),fliplr(A14)));%MA模型的分子
- y24=filter(B14,A1,randn(1,N));%.*[zeros(1,200),ones(1,256)];%由估计出的MA模型产生数据
- %---因为(q=4)<<L<<N=256,所以选取L=32---%
- [Ama4,Ema4]=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=[1,a8']%AR的参数a(1)-a(8)的估计值
- B18=fliplr(conv(fliplr(B1),fliplr(A18)))%MA模型的分子
- y28=filter(B18,A1,randn(1,N));%.*[zeros(1,200),ones(1,256)];%由估计出的MA模型产生数据
- %---因为(q=8)<<L<<N=256,所以选取L=48---%
- [Ama8,Ema8]=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([0 0.5 0 2000]);
[ 本帖最后由 eight 于 2007-2-16 16:38 编辑 ] |
发表于 2006-5-9 21:38
提升卡
变色卡
显身卡
