xitongbianshi.rar_matlab参数辨识_参数辨识_系统辨识_辅助变量法

clear all close all % 生成输入的M序列[-2 2] M = [1,0,1,0,1,0,0,1,1,1]; % 设置寄存器的初始状态 N = length(M); T =2^N-1; % M 序列的周期 u = zeros(1,T); for i = 1:T; m = xor(M(10),M(5)); % 反馈逻辑运算 for j = N:-1:2 M(j) = M(j-1); % 移位运算 end M(1) = m; u(i) = M(1); if u(i)==0 u(i)=-2; else u(i)=2; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % z的前三个初始值置零 z(1)=0;z(2)=0;z(3)=0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% v=randn(1,T); e(1)=0.6; r=0.9; %控制信噪比 for i=2:T e(i)=-0.1*e(i-1)+r*v(i); end % 构造输出数据(或含噪声数据) v = normrnd(0,1,T,1); for k=4:T; z(k)=1.5*z(k-1)-0.9*z(k-2)+u(k-1)+0.5*u(k-2)+e(k); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 准备阶次及参数辨识 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 置参数theta的初始值(很小的实向量) theta0=10^(-5)*diag(eye(6,6)); % 置协方差阵P的初始值(a^2*I,a充分大) p0=10^6*eye(6,6); % 相对误差 epsilon=10^(-8); cc=[zeros(5,T-1)]; ee=zeros(6,T); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %AIC阶次辨识(根据书P353及P245编写） %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% a=4;b=4;length=400; L=max(size(u)); aic=zeros(a,b); for n=1:a for m=1:b theta=10^(-9)*ones(2*n+m,1); P=eye(2*n+m); aa=max(m,n); zf=zeros(1,L); uf=zeros(1,L); vf(1:aa)=0; for k=(aa+1):L h=[-z((k-1):-1:(k-n)) u((k-1):-1:(k-n)) v((k-1):-1:(k-m))]; v(k)=z(k)-h*theta; vf(k)=v(k)-vf((k-1):-1:(k-m))*theta((2*n+1):(2*n+m)); zf(k)=z(k)-zf((k-1):-1:(k-m))*theta((2*n+1):(2*n+m)); uf(k)=u(k)-uf((k-1):-1:(k-m))*theta((2*n+1):(2*n+m)); hf=[-zf((k-1):-1:(k-n)) uf((k-1):-1:(k-n)) vf((k-1):-1:(k-m))]; K=P*hf'*inv(hf*P*hf'+1); P=(eye(2*n+m)-K*hf)*P; theta=theta+K*(z(k)-h*theta);, end v=zeros(1,L); for k=(aa+1):L v(k)=z(k)+[z((k-1):-1:(k-n)) -u((k-1):-1:(k-n)) -v((k-1):-1:(k-m))]*theta; end sigma=1/L*v(length+1:L)*v(length+1:L)'; AIC(n,m)=(L-length)*log(sigma)+2*(2*n+m); end end AIC AA=min(min(AIC)); for i=1:n for j=1:m if AIC(i,j)==AA c=[i i j]; end end end sprintf('分母，分子，噪声的阶次为: %d %d %d',c(1:3)) close all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 基于自适应滤波的辅助变量法对参数进行辨识 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% theta=10^(-5)*diag(eye(2*c(1),2*c(1))); P=10^5*eye(2*c(1)); aa=max(c(1),c(2)); l=max(size(u)); x(1:aa)=0; cc=2; y(:,1)=theta; v(:,aa)=0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=(aa+1):l x(i)=[-x((i-1):-1:(i-c(1))) u((i-1):-1:(i-c(2)))]*theta; h=[-z((i-1):-1:(i-c(1))) u((i-1):-1:(i-c(2)))]'; h1=[-x((i-1):-1:(i-c(1))) u((i-1):-1:(i-c(2)))]'; K=P*h1*inv(h'*P*h1+1); P=(eye(2*c(1))-K*h')*P; t1=z(i)-h'*theta; theta=theta+K*t1; y(:,cc)=theta; cc=cc+1; v(:,i)=t1; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% thetazhen=[-1.5 0.9 1 0.5]'; theta=[theta thetazhen] figure(1) plot(y') Xlabel('K') Ylabel('参数值') title('参数变化曲线') figure(2) plot(v) Xlabel('K') Ylabel('残差值') title('残差变化曲线')