終於成功仿了一次Kalman濾波器

終於成功仿了一次Kalman濾波器

首先是測試了從網上down的一段代碼

% KALMANF - updates a system state vector estimate based upon an
% observation, using a discrete Kalman filter.
%
% Version 1.0, June 30, 2004
%
% This tutorial function was written by Michael C. Kleder
% (Comments are appreciated at: public@kleder.com)
%
% INTRODUCTION
%
% Many people have heard of Kalman filtering, but regard the topic
% as mysterious. While it's true that deriving the Kalman filter and
% proving mathematically that it is "optimal" under a variety of
% circumstances can be rather intense, applying the filter to
% a basic linear system is actually very easy. This Matlab file is
% intended to demonstrate that.
%
% An excellent paper on Kalman filtering at the introductory level,
% without detailing the mathematical underpinnings, is:
% "An Introduction to the Kalman Filter"
% Greg Welch and Gary Bishop, University of North Carolina
% http://www.cs.unc.edu/~welch/kalman/kalmanIntro.html
%
% PURPOSE:
%
% The purpose of each iteration of a Kalman filter is to update
% the estimate of the state vector of a system (and the covariance
% of that vector) based upon the information in a new observation.
% The version of the Kalman filter in this function assumes that
% observations occur at fixed discrete time intervals. Also, this
% function assumes a linear system, meaning that the time evolution
% of the state vector can be calculated by means of a state transition
% matrix.
%
% USAGE:
%
% s = kalmanf(s)
%
% "s" is a "system" struct containing various fields used as input
% and output. The state estimate "x" and its covariance "P" are
% updated by the function. The other fields describe the mechanics
% of the system and are left unchanged. A calling routine may change
% these other fields as needed if state dynamics are time-dependent;
% otherwise, they should be left alone after initial values are set.
% The exceptions are the observation vectro "z" and the input control
% (or forcing function) "u." If there is an input function, then
% "u" should be set to some nonzero value by the calling routine.
%
% SYSTEM DYNAMICS:
%
% The system evolves according to the following difference equations,
% where quantities are further defined below:
%
% x = Ax + Bu + w meaning the state vector x evolves during one time
% step by premultiplying by the "state transition
% matrix" A. There is optionally (if nonzero) an input
% vector u which affects the state linearly, and this
% linear effect on the state is represented by
% premultiplying by the "input matrix" B. There is also
% gaussian process noise w.
% z = Hx + v meaning the observation vector z is a linear function
% of the state vector, and this linear relationship is
% represented by premultiplication by "observation
% matrix" H. There is also gaussian measurement
% noise v.
% where w ~ N(0,Q) meaning w is gaussian noise with covariance Q
% v ~ N(0,R) meaning v is gaussian noise with covariance R
%
% VECTOR VARIABLES:
%
% s.x = state vector estimate. In the input struct, this is the
% "a priori" state estimate (prior to the addition of the
% information from the new observation). In the output struct,
% this is the "a posteriori" state estimate (after the new
% measurement information is included).
% s.z = observation vector
% s.u = input control vector, optional (defaults to zero).
%
% MATRIX VARIABLES:
%
% s.A = state transition matrix (defaults to identity).
% s.P = covariance of the state vector estimate. In the input struct,
% this is "a priori," and in the output it is "a posteriori."
% (required unless autoinitializing as described below).
% s.B = input matrix, optional (defaults to zero).
% s.Q = process noise covariance (defaults to zero).
% s.R = measurement noise covariance (required).
% s.H = observation matrix (defaults to identity).
%
% NORMAL OPERATION:
%
% (1) define all state definition fields: A,B,H,Q,R
% (2) define intial state estimate: x,P
% (3) obtain observation and control vectors: z,u
% (4) call the filter to obtain updated state estimate: x,P
% (5) return to step (3) and repeat
%
% INITIALIZATION:
%
% If an initial state estimate is unavailable, it can be obtained
% from the first observation as follows, provided that there are the
% same number of observable variables as state variables. This "auto-
% intitialization" is done automatically if s.x is absent or NaN.
%
%x = inv(H)*z
%P = inv(H)*R*inv(H')
%
% This is mathematically equivalent to setting the initial state estimate
% covariance to infinity.
%
% SCALAR EXAMPLE (Automobile Voltimeter):
%
% % Define the system as a constant of 12 volts:
function T
clear s
s.x = 12;
s.A = 1;
% % Define a process noise (stdev) of 2 volts as the car operates:
s.Q = 2^2; % variance, hence stdev^2
% Define the voltimeter to measure the voltage itself:
s.H = 1;
% % Define a measurement error (stdev) of 2 volts:
s.R = 2^2; % variance, hence stdev^2
%Do not define any system input (control) functions:
s.B = 0;
s.u = 0;
% % Do not specify an initial state:
s.x = nan;
s.P = nan;
% % Generate random voltages and watch the filter operate.
tru=[]; % truth voltage
for t=1:20
    tru(end+1) = randn*2+12;
    s(end).z = tru(end) + randn*2; % create a measurement
    s(end+1)=kalmanf(s(end)); % perform a Kalman filter iteration
    % end
    % figure
    % hold on
    % grid on
    % % plot measurement data:
    hz=plot([s(1:end-1).z],'r');hold on
    % % plot a-posteriori state estimates:
    hk=plot([s(2:end).x],'b-');hold on
    ht=plot(tru,'g-');hold on
    legend('observations','Kalman output','true voltage',0)
    title('Automobile Voltimeter Example')
    % hold off
end    
    
function s = kalmanf(s)

% set defaults for absent fields:
if ~isfield(s,'x'); s.x=nan*z; end
if ~isfield(s,'P'); s.P=nan; end
if ~isfield(s,'z'); error('Observation vector missing'); end
if ~isfield(s,'u'); s.u=0; end
if ~isfield(s,'A'); s.A=eye(length(x)); end
if ~isfield(s,'B'); s.B=0; end
if ~isfield(s,'Q'); s.Q=zeros(length(x)); end
if ~isfield(s,'R'); error('Observation covariance missing'); end
if ~isfield(s,'H'); s.H=eye(length(x)); end

if isnan(s.x)
    % initialize state estimate from first observation
    if diff(size(s.H))
        error('Observation matrix must be square and invertible for state autointialization.');
    end
    s.x = inv(s.H)*s.z;
    s.P = inv(s.H)*s.R*inv(s.H'); 
else
    
    % This is the code which implements the discrete Kalman filter:
    
    % Prediction for state vector and covariance:
    s.x = s.A*s.x + s.B*s.u;
    s.P = s.A * s.P * s.A' + s.Q;
    
    % Compute Kalman gain factor:
    K = s.P*s.H'*inv(s.H*s.P*s.H'+s.R);
    
    % Correction based on observation:
    s.x = s.x + K*(s.z-s.H*s.x);
    s.P = s.P - K*s.H*s.P;
    
    % Note that the desired result, which is an improved estimate
    % of the sytem state vector x and its covariance P, was obtained
    % in only five lines of code, once the system was defined. (That's
    % how simple the discrete Kalman filter is to use.) Later,
    % we'll discuss how to deal with nonlinear systems.
    
end

後來不過癮，本身寫了一個，沒想到稍微改了一改居然成功了，效果還不錯

% 狀態
% xk=A•xk-1+B•uk+wk
% zk=H•xk+vk,
% p(w) ~ N(0,Q)
% p(v) ~ N(0,R),

% 預測
% x'k=A•xk+B•uk
% P'k=A•P(k-1)*AT + Q
% 修正
% Kk=P'k•HT•(H•P'k•HT+R)-1
% xk=x'k+Kk•(zk-H•x'k)
% Pk=(I-Kk•H)•P'k

%要注意的是：必須把系統狀態和kalman濾波器內部預測的狀態分開
function Test
A=[1 0.1;0 1];
B=0;
Xp=rand(2,1)*0.1; 
X=[0 0]';
H=[1 0];
Q=eye(2)*1e-5;
R=eye(1)*0.1; 
P=eye(2);% P'(k)
angle=[];
angle_m=[];
angle_real=[];
for i=1:500
    angle_real=[angle_real X(1)]; %實際角度
    
    [Xp,P]=Predict(A,Xp,P,Q); 
    
    X=A*X+rand(2,1)*1e-5;
    z_m=H*X+rand(1,1)*0.1-0.05;  
    angle_m=[angle_m z_m(1)];   %測量的角度
        
    [Xp,P]=Correct(P,H,R,X,z_m);
    angle=[angle Xp(1)];     %預測的角度    
end
t=1:500;
plot(t,angle,'r',t,angle_m,'g',t,angle_real,'b')
legend('預測值','測量值','實際值')

figure
plot(t,angle-angle_real,'r',t,angle_m-angle_real,'g')
legend('濾波後的偏差','測量的偏差')
title('偏差分析')
xlabel('time');
ylabel('error');

function [Xk,Pk]=Predict(A,Xk,Pk_1,Q)
Xk=A*Xk;
Pk=A*Pk_1*A'+Q;

function [Xk,Pk]=Correct(Pk,H,R,Xk,zk)
Kk=Pk * H' * inv(H * Pk * H' + R);
Xk=Xk+ Kk*(zk-H*Xk);
Pk=(eye(size(Pk,1)) - Kk*H)*Pk;

 

 

順便附上幾張仿真圖
這是狀態圖

這是偏差分析