%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');