function [xt,Rx,yp]=Kalman(xt,yt,ut,M,N,F,A,B,G,Rw,Rv,Rx)
  // Kalman filter for state estimtion with the model
  //                             xt = M*xt + N*ut + F + w
  //                             yt = A*xt + B*ut + G + v
  // xt    state
  // Rx    state estimate covariance matrix
  // yp    output prediction
  // yt    output
  // ut    input
  // M,N,F state model parameters
  // A,B,G output model parameters
  // Rw    state model covariance
  // Rv    output model covariance

  nx=size(M,1);
  ny=size(A,1);
  if isempty(F), F=zeros(nx,1); end
  if isempty(G), G=zeros(ny,1); end

  xt=M*xt+N*ut+F;                      // time update of the state
  Rx=Rw+M*Rx*M';                       // time updt. of state covariance

  yp=A*xt+B*ut+G;                      // output prediction
  Ry=Rv+A*Rx*A';                       // noise covariance update
  Rx=Rx-Rx*A'*inv(Ry)*A*Rx;            // state est. coariance update
  ey=yt-yp;                            // prediction error 
  KG=Rx*A'*inv(Rv);                    // Kalman gain
  xt=xt+KG*ey;                         // data update of the state
endfunction
