cor.m

关于系统辨识的一些仿真程序

%系统辨识--二步法一次完成
%产生M序列
clc
clear
close all
J=16;
y1=1;y2=1;y3=1;y4=1;y5=0;
for i=1:J
    x1=xor(y3,y5);
    x2=y1;
    x3=y2;
    x4=y3;
    x5=y4;
    y(i)=y5;
    if y(i)>0.5,u(i)=-1;
    else u(i)=1;
    end
    y1=x1;y2=x2;y3=x3;y4=x4;y5=x5;
end
z=zeros(1,16);
N=16;
var_v=1;
v=0.02*sqrt(var_v)*randn(1,N);%均值为0,方差为var_v的白噪声信号
for k=5:16
    z(k)=3.808*z(k-1)-5.434*z(k-2)+3.445*z(k-3)-0.8187*z(k-4)+u(k)-3.935*u(k-1)+5.806*u(k-2)-3.807*u(k-3)+0.9362*u(k-4)+v(k);
end
ruz=zeros(1,35);
ruz=xcorr(u,z,17);
plot(-17:17,ruz);
HL=[-ruz(16) -ruz(15) -ruz(14) -ruz(13) 1 -1/31 -1/31 -1/31 -1/31;
    -ruz(17) -ruz(16) -ruz(15) -ruz(14) -1/31 1 -1/31 -1/31 -1/31;
    -ruz(18) -ruz(17) -ruz(16) -ruz(15) -1/31 -1/31 1 -1/31 -1/31;
    -ruz(19) -ruz(18) -ruz(17) -ruz(16) -1/31 -1/31 -1/31 1 -1/31;
    -ruz(20) -ruz(19) -ruz(18) -ruz(17) -1/31 -1/31 -1/31 -1/31 1;
    -ruz(21) -ruz(20) -ruz(19) -ruz(18) -1/31 -1/31 -1/31 -1/31 -1/31;
    -ruz(22) -ruz(21) -ruz(20) -ruz(19) -1/31 -1/31 -1/31 -1/31 -1/31;
    -ruz(23) -ruz(22) -ruz(21) -ruz(20) -1/31 -1/31 -1/31 -1/31 -1/31;
    -ruz(24) -ruz(23) -ruz(22) -ruz(21) -1/31 -1/31 -1/31 -1/31 -1/31;
    -ruz(25) -ruz(24) -ruz(23) -ruz(22) -1/31 -1/31 -1/31 -1/31 -1/31;
    -ruz(26) -ruz(25) -ruz(24) -ruz(23) -1/31 -1/31 -1/31 -1/31 -1/31;
    -ruz(27) -ruz(26) -ruz(25) -ruz(24) -1/31 -1/31 -1/31 -1/31 -1/31;
    -ruz(28) -ruz(27) -ruz(26) -ruz(25) -1/31 -1/31 -1/31 -1/31 -1/31;
    -ruz(29) -ruz(28) -ruz(27) -ruz(26) -1/31 -1/31 -1/31 -1/31 -1/31;
    -ruz(30) -ruz(29) -ruz(28) -ruz(27) -1/31 -1/31 -1/31 -1/31 -1/31;
    -ruz(31) -ruz(30) -ruz(29) -ruz(28) -1/31 -1/31 -1/31 -1/31 -1/31];
ZL=[ruz(17);ruz(18);ruz(19);ruz(20);ruz(21);ruz(22);ruz(23);ruz(24);ruz(25);ruz(26);ruz(27);ruz(28);ruz(29);ruz(30);ruz(31);ruz(32)];
c1=HL'*HL; c2=inv(c1); c3=HL'*ZL; c=c2*c3 ;
a1=c(1), a2=c(2), a3=c(3), a4=c(4),b1=c(5), b2=c(6), b3=c(7), b4=c(8),b5=c(9)
ze=zeros(1,16);
for k=5:16
    z(k)=3.808*z(k-1)-5.434*z(k-2)+3.445*z(k-3)-0.8187*z(k-4)+u(k)-3.935*u(k-1)+5.806*u(k-2)-3.807*u(k-3)+0.9362*u(k-4)+v(k);
    ze(k)=-a1*ze(k-1)-a2*ze(k-2)-a3*ze(k-3)-a4*ze(k-4)+b1*u(k)+b2*u(k-1)+b3*u(k-2)+b4*u(k-3)+b5*u(k-4)+v(k);
    ee(k)=abs((ze(k)-ruz(k))/ruz(k));
end
figure(3)
i=1:16; 
plot(i,ee,'r')
title('COR-LS算法模型误差')