最小二乘法

%广义最小二乘递推算法 %z(k+2)=1.5*z(k+1)-0.7*z(k)+u(k+1)+0.5*u(k)+v(k) clear all; clc; %400个产生M序列作为输入 x = [0 1 0 1 1 0 1 1 1]; %初始值 n = 403; M = []; %存放M序列 for i = 1:n temp = xor(x(4),x(9)); M(i) = x(9); for j = 9:-1:2 x(j) = x(j-1); end x(1) = temp; end %产生均值为0，方差为1的高斯白噪声 v = randn(1,400); e = []; e(1) = v(1); e(2) = v(2); for i = 3:400 e(i) = 0*e(i-1)+0*e(i-2)+v(i); end %产生观测序列 z = zeros(400,1); z(1) = -1; z(2) = 0; for i = 3:400 z(i) = 1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-1)+v(i-2); end %变化后的观测序列 zf = []; zf(1) = -1; zf(2) = 0 for i = 3:400 zf(i) = z(i)-0*z(i-1)-0*z(i-2); end %变换后的输入序列 uf = []; uf(1) = M(1); uf(2) = M(2); for i = 3:400 uf(i) = M(i)-0*M(i-1)-0*M(i-2); end %赋初值 p = 100*eye(4); %估计方差 pstore = zeros(4,400); pstore(:,2) = [p(1,1),p(2,2),p(3,3),p(4,4)]; theta = zeros(4,400); %参数的估计值，存放中间过程估计值 theta(:,2) = [3;3;3;3]; k = [10;10;10;10]; %增益矩阵 pe = 10*eye(2); thetae = zeros(2,400); thetae(:,2) = [0.5;0.3]; ke = [10;10]; for i = 3:400 h = [-zf(i-1);-zf(i-2);uf(i-1);uf(i-2)]; k = p*h*inv(h'*p*h+1); theta(:,i) = theta(:,i-1)+k*(z(i)-h'*theta(:,i-1)); p = (eye(4)-k*h')*p; pstore(:,i-1) = [p(1,1),p(2,2),p(3,3),p(4,4)]; he = [-e(i-1);-e(i-2)]; ke = pe*he*inv(1+he'*pe*he); thetae(:,i) = thetae(:,i-1)+ke*(e(i)-he'*thetae(:,i-1)); pe = (eye(2)-ke*he')*pe; end %输出结果及作图 disp('参数a1 a2 b1 b2的估计结果'); theta(:,400); disp('噪声传递函数c1 c2的估计结果为：'); thetae(:,400); i = 1:400; figure; plot(i,theta(1,:),i,theta(2,:),i,theta(3,:),i,theta(4,:)); title('待估参数过渡过程'); figure; plot(i,pstore(1,:),i,pstore(2,:),i,pstore(3,:),i,pstore(4,:)); title('估计方差变化过程'); figure; plot(i,thetae(1,:),i,thetae(1,:));