% 方差為0.1的白噪聲  最小二乘辨識
% 清除變量y
  clear all;
% 輸入u(k)
  u=[1.147, 0.201, -0.787, -1.589, -1.052, 0.866, 1.152, 1.573, 0.626, 0.433, -0.958, ...
     0.81, -0.044, 0.947, -1.474, -0.719, -0.086, -1.099, 1.45, 1.151, 0.485, 1.633, ...
     0.043, 1.326, 1.706, -0.340, 0.890, 0.144, 1.172, -0.390];    
  randn('state',1);
  eibuqizh=randn(1,30);
  eibuqizh=eibuqizh*sqrt(0.1);                    
  kesi=eibuqizh;
% 計算y(k)
  for(k=1:30)
    if k==1
       y(k)=kesi(k);
    elseif  k==2
       y(k)=-1.642*y(k-1)+0.39*u(k-1)+kesi(k);
    else
       y(k)=-1.642*y(k-1)-0.715*y(k-2)+0.39*u(k-1)+0.35*u(k-2)+kesi(k);
    end
  end
  y
  y=y';
  save y_ls y;
% 計算測量矩陣Fai_ls,n=2,N=28
  for (i=1:28)
    for(j=1:2)
        A(i,j)=-y(2+i-j);    
    end
    for(l=1:3)
        B(i,l)=u(3+i-l);
    end
  end
  Fai_ls=[A,B];
% 計算估計參數guji_ls
  e=y(3:30);
  guji_ls=inv(Fai_ls'*Fai_ls)*Fai_ls'*e
  %guji_ls(3)=0  % 第三項為０，即b0=０，
  save guji_ls guji_ls;
  save Fai_ls Fai_ls;