bianshi.rar_prbs_prbs system_系统辨识

%PRBS输入信号产生程序 X1=0; X2=1; X3=0; X4=0; m=200; for i=1:m Y4=X4; Y3=X3; Y2=X2; Y1=X1; X4=Y3; X3=Y2; X2=Y1; X1=xor(Y3,Y4); if Y4==0 u(i)=-0.5; else u(i)=0.5; end end figure(1); k=1:m; plot(k,u,k,u,'kx'); xlabel('k'); axis([0 m -0.6 0.6]); %产生输出信号并保存 y=zeros(1,200); for k=4:200 y(k)=0.9631*y(k-1)+4.796*u(k-3) end save('data.txt','y','-ascii'); %产生不同方差的噪声 h=[0 0.5 1 5.0]; for w=1:4 e1=normrnd(0,h(w)^0.5,1,200); y1=zeros(1,200); for k=4:200 y1(k)=0.9631*y1(k-1)+4.796*u(k-3)+e1(k); end %构造正则矩阵H H=zeros(197,4); for i=1:197 H(i,1)=y1(i+2); H(i,2)=u(i+2); H(i,3)=u(i+1); H(i,4)=u(i); end %使用基本最小二乘法进行参数估计 sita=inv(H'*H)*H'*(y1(4:200))' %使用辨识参数得到输出值 y2=zeros(1,200); for k=4:200 y2(k)=sita(1)*y2(k-1)+sita(2)*u(k-1)+sita(3)*u(k-2)+sita(4)*u(k-3); end %输出采样值误差 e1=zeros(1,200); for i=1:200 e1(i)=y(i)-y2(i); end if (w==1) figure(2); subplot(2,1,1); i=1:200; plot(i,y); xlabel('k'); ylabel('y'); title('实际系统输出曲线') hold on grid on %方差为0 subplot(2,1,2); i=1:200; plot(i,y2,'b'); xlabel('k'); ylabel('y2'); title('辨识系统输出曲线') hold on grid on %输出误差曲线 figure(3); i=1:200; plot(i,e1,'b'); xlabel('k'); ylabel('e1'); hold on grid on end if (w==2) %方差为0.5 figure(2) subplot(2,1,2); i=1:200; plot(i,y2,'r'); xlabel('k'); ylabel('y2'); hold on grid on figure(3); i=1:200; plot(i,e1,'r'); xlabel('k'); ylabel('e1'); hold on grid on end if (w==3) %方差为1.0 figure(2) subplot(2,1,2); i=1:200; plot(i,y2,'c'); xlabel('k'); ylabel('y2'); hold on grid on figure(3); i=1:200; plot(i,e1,'c'); xlabel('k'); ylabel('e1'); hold on grid on end if (w==4) %方差为5.0 figure(2) subplot(2,1,2); i=1:200; plot(i,y2,'g'); xlabel('k'); ylabel('y2'); hold on grid on figure(3); i=1:200; plot(i,e1,'g'); xlabel('k'); ylabel('e1'); hold on grid on end end figure(2); subplot(2,1,2); legend('σ=0.0','σ=0.5','σ=1.0','σ=5.0'); figure(3); legend('σ=0.0','σ=0.5','σ=1.0','σ=5.0');