最小二乘法进行系统辨识

%% %最小二乘法辨识系统 clear clc %% %生成4阶的M序列，幅值为1 %%移位寄存器产生M序列 Np=63;%周期 deltaT=1;%节拍 a=1;%幅值 %逻辑“0”为a,逻辑“1”为-a; %初始序列（010101） x1=1; x2=0; x3=1; x4=0; x5=1; x6=0; %进行循环移位，生成M序列 for i=1:deltaT:Np y4=x4;y3=x3;y2=x2;y1=x1; x4=y3;x3=y2;x2=y1;x1=xor(y3,y4); if y4==0 %u(i)=-1 ; %逻辑“0”为a u(i)=a; else %u(i)=y4; %逻辑“1”为-a; u(i)=-a; end end %绘图 i1=i; k=1:deltaT:i1; plot(k,u,'b'); xlabel('k'); ylabel('M序列'); title('移位寄存器产生的M序列') %% %生成服从正态分布的白噪声N(0,1) %初始化 A=65539; N=1200; x0=1; M=2147483647; C=0;%混合同余法直接变成乘同余法 %乘法递推N次 for k=1:N x2=A*x0; x1=mod(x2,M); v1=x1/2147483647;%将x1中的数除以M得到小于1的随机数 v(:,k)=v1; x0=x1; end v2=v;%保存0-1随机数到v2 ave=mean(v) var=var(v) k1=k; save v; %绘图程序 %k=1:k1; %plot(k,v,'r'); %xlabel('k');ylabel('v');title('(0-1)'); [num,val]=hist(v,10); %生成服从正态分布的随机数 for n=1:100 sum1(:,n)=sum(v(:,(12*(n-1)+1):12*n));%将v中的数分为10组，每组12个求和 vy(:,n)=sum1(:,n)-6;%vy近似服从N（0,1）分布 end figure n=1:100; plot(n,vy,'b'); xlabel('n');ylabel('vy');title('服从N(0,1)的白噪声') %% %最小二乘系统辨识,输入信号取值是从k=1到k=16的M序列 z=zeros(17,1);%定义观测输出值的长度17行*1列 ZL=zeros(15,1);%定义观测输出值的长度15行*1列 for k=3:17 z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+vy(k);%用理想输出作为观测值 ZL(k-2)=z(k); end %构建观测矩阵HL for m=2:16 HL=[-z(m) -z(m-1) u(m) u(m-1)]; for n=1:4 H((m-1),n)=HL(n); end end figure m=3:17; plot(m,z(m),'b'); xlabel('m');ylabel('z');title('观测输出值'); %求待辨识矩阵c=[a1 a2 b1 b2]',c=((HL'*HL)^(-1))*HL'*ZL c1=H'*H; c2=inv(c1); c3=H'*ZL; c=c2*c3; c,a1=c(1),a2=c(2),b1=c(3),b2=c(4)