LS1.rar_辅助变量法

function A=LS1 T=15000; %% 生成M序列 M=zeros(T,4); M(1,:)=ones(1,4); for i=2:T temp=bitxor(M(i-1,3),M(i-1,4));%4,3位亦或 for j=4:-1:2 M(i,j)=M(i-1,j-1);%向右移位 end M(i,1)=temp;%首位=temp end M_sequence=M(:,4); M_sequence(find(M(:,4)==0))=-1; u=M_sequence;%输入信号 %% 生成仿真数据 a=[1.642 0.715]; b=[0.39 0.35];%输入真实值 ture=[a b]'; y=zeros(1,T); n=y; NoiseData=normrnd(0,0.1,1,T);%高斯噪声 NoiseData=NoiseData-mean(NoiseData);%输入的噪声 for k=3:T n(k)=NoiseData(k)+a(1)*NoiseData(k-1)+a(2)*NoiseData(k-2); y(k)=-a(1)*y(k-1)-a(2)*y(k-2)+b(1)*u(k-1)+b(2)*u(k-2); end %观测值 figure(8) plot(y);%绘出观测的数据 title('仿真数据及噪声') xlabel('仿真次数') ylabel('信号及噪声大小') hold on plot(n,'r')%绘出噪声 hold off y=y+n; %带有噪声的测量值 %% 递推最小二乘 aa2=zeros(4,T);%递推最小二乘估计值阵 K=zeros(4,T);%K阵 P=cell(T,1);%P P{2}=5000*eye(4);%p的初值 aa2(:,2)=zeros(4,1);%列2恒为0 for k=3:T q=[-y(k-1) -y(k-2) u(k-1) u(k-2)]';%pusai K(:,k)=P{k-1}*q*inv(1+q'*P{k-1}*q);%K=P*pusai*（1+pusai'*p*pusai) aa2(:,k)=aa2(:,k-1)+K(:,k)*(y(k)-q'*aa2(:,k-1));%theta=theta+K*(y-pusai'*theta) P{k}=P{k-1}-P{k-1}*q*inv(1+q'*P{k-1}*q)*q'*P{k-1};%p=p-p*pusai*(1+pusai'*p*pusai)*pusai'*p end figure(1) h=plot(aa2');%绘出得到的观测值 title('递推最小二乘收敛过程') xlabel('仿真次数') ylabel('参数估计') text([T-100 ;T-100; T-100; T-100],[1.7; 0.8; 0.5 ;0.2],['a1'; 'a2'; 'b1'; 'b2' ]) figure(2) plot(sum((aa2-repmat(ture,1,T)).^2));%绘出方差 v1=(aa2-repmat(ture,1,T)).^2; var=v1(:,T) title('递推最小二乘误差指标函数曲线') xlabel('仿真次数') ylabel('递推最小二乘误差指标') %% 递推辅助变量法 aa4=zeros(4,T); K=zeros(4,T); P=cell(T,1); P{2}=5000*eye(4); aa4(:,2)=zeros(4,1); Q=zeros(4,T); Z=Q;%辅助变量阵 % 前15步 用最小二乘递推 for k=3:60 q=[-y(k-1) -y(k-2) u(k-1) u(k-2)]'; Q(:,k-2)=q; %构建测量阵 K(:,k)=P{k-1}*q*inv(1+q'*P{k-1}*q); aa4(:,k)=aa4(:,k-1)+K(:,k)*(y(k)-q'*aa4(:,k-1)); P{k}=P{k-1}-P{k-1}*q*inv(1+q'*P{k-1}*q)*q'*P{k-1}; end %递推最小二乘求出粗略的参数 % 15步以后转换到辅助变量法 for k=61:T q=[-y(k-1) -y(k-2) u(k-1) u(k-2)]';%pusai Q(:,k-2)=q;%测量阵 yy(k-1)=Q(:,k-3)'*aa4(:,k-1); yy(k-2)=Q(:,k-4)'*aa4(:,k-2);%求辅助变量 z=[-yy(k-1) -yy(k-2) u(k-1) u(k-2)]'; Z(:,k-2)=z;%构造辅助阵 K(:,k)=P{k-1}*z*inv(1+q'*P{k-1}*z);%thetaN+1=thetaN+Kn+1(yn+1-pusaiN+1'*thetaN) aa4(:,k)=aa4(:,k-1)+K(:,k)*(y(k)-q'*aa4(:,k-1));%Pn+1=pn-Kn+1*pusai'*pn P{k}=P{k-1}-P{k-1}*q*inv(1+q'*P{k-1}*q)*q'*P{k-1};%Kn+1=Pn*zn+1*(pusai'*Pn*zn+1)^-1 end figure(3) h=plot(aa4'); title('递推辅助变量法收敛过程') xlabel('仿真次数') ylabel('参数估计') text([T-100 ;T-100; T-100; T-100],[1.7; 0.8; 0.5 ;0.2],['a1'; 'a2'; 'b1'; 'b2' ]) figure(4) plot(sum((aa4-repmat(ture,1,T)).^2)); v2=(aa4-repmat(ture,1,T)).^2; title('递推辅助变量法误差指标函数曲线') xlabel('仿真次数') ylabel('递推辅助变量法误差指标') %% 广义最小二乘 c0=zeros(4,1); %被辨识参数的初始值直接给取，取一个充分小的实向量 ce0=zeros(1,1); %滤波器C(q-1)的参数初值 p0=5000*eye(4,4); %初始状态P0也采用直接取方式，直接给出初始状态P0，取一个充分大的实数单位矩阵 pe0=5000*eye(1,1); %计算滤波器的P矩阵初值 yf=zeros(1,T); uf=zeros(1,T); c=zeros(4,T); %被辨识参数矩阵的初始值及大小 e=zeros(4,T); %相对误差-残差的初始值及大小 ce=zeros(1,T); %滤波器参数每次辨识的结果 ee=zeros(1,T); for k=3:T; yf(k)=y(k)+ce0(1)*y(k-1); %滤波后的输出序列 uf(k)=u(k)+ce0(1)*u(k-1); %滤波后的输入序列 h1=[-yf(k-1),-yf(k-2),uf(k-1),uf(k-2)]'; %求h(k)%滤波后的输出输出组成的辨识数据 k1=p0*h1*inv(h1'*p0*h1+1);%求出K的值 c1=c0+k1*(yf(k)-h1'*c0);%求被辨识参数c e1=c1-c0;%求参数当前值与上一次的值的差值 e2=e1./c0;%求参数的相对变化 e(:,k)=e2; %把当前相对变化的列向量加入误差矩阵的最后一列 c0=c1; %新获得的参数作为下一次递推的旧参数 c(:,k)=c1; %把当前所辨识参数的c1列向量加入辨识参数矩阵的最后一列 p1=p0-k1*k1'*[h1'*p0*h1+1]; %求p(k)值 p0=p1; %把当前的p(k)值给下次用 ee(k)=y(k)-[-y(k-1),-y(k-2),u(k-1),u(k-2)]*c1; %计算残差e(k),注意这里的h的组成 he1=[-ee(k-1)]; ke1=pe0*he1*inv(1+he1'*pe0*he1); %求当前的K矩阵 ce1=ce0+ke1*(ee(k)-he1'*ce0); ce(k)=ce1; %保存滤波器参数 ce0=ce1; %为下一时刻做准备 pe1=pe0-ke1*ke1'*[he1'*pe0*he1+1]; %为下一时刻赋值 pe0=pe1; %为下一时刻P矩阵赋值 end figure(5) h=plot(c'); title('广义最小二乘系统参数收敛过程') xlabel('仿真次数') ylabel('参数估计') text([T-100 ;T-100; T-100; T-100],[1.7; 0.8; 0.5 ;0.2],['a1'; 'a2'; 'b1'; 'b2' ]); figure(6) plot(sum((c-repmat(ture,1,T)).^2)); v3=(c-repmat(ture,1,T)).^2; title('广义最小二乘误差指标函数曲线') xlabel('仿真次数') ylabel('广义最小二乘误差指标') figure(7) plot(ce') title('广义最小二乘成形滤波器参数收敛过程') xlabel('仿真次数') ylabel('参数估计') %% 结果 disp('第一列：真值') disp('第二列：递推最小二乘') disp('第三列：递推辅助变量') disp('第四列：广义最小二乘') A=[[a b]' aa2(:,T) aa4(:,T),c(:,T)];