Kalman Filter

Simulated GPS time series in position(Top), the long tern trend (middle), and periodic annual signal (bottom). Set y=10+5.5t+1.9sin(2πt)+3sin(4πt)+err , initial states were set to be [0,0,0,0].

Kalman Filter could automatically fit the real state (red number in above equation) in several iterations, even though initial guesses were wrong.

METHOD

The figure shows the conception of Kalman filter. In real case, we have our observations and a physical model to explain the observations. For example, a free-fall object could following the equation of (S=Vt+0.5gt^2). However, the measuring position of an object could be varied by measuring error, also, the equation cannot fully explain the motion of free-fall (e.g. the drag force). Thus, we need Kalman filter to correct these variation.

The core conception of Kalman filter is "optimize" and "update":

Where x is the position, the minus means previous state, and will be update to the next state (i.e. the position x at time t is related to the position at time t-1 ). P and Q are the model noise and data noise coveriance matrix, respectively.(e.g. drag force and measure error). Phi denotes the transition state matrix, which transits the state from previous to "ready" (i.e. x1=x0+0.5g^t2) but not yet optimizing.

The optimizing process is depended on the factor K, which also called Kalman gain. Kalman gain is similar to our conception of averaging. For example for 2 measurements {9 , 10} (with analogize to observation and model), the optimize value will be the average of 9 and 10 and will be 9.5. The detail calculation is 9*(1/2)+10*(1/2), where (1/2) is the blending factor. And the reason why (1/2) is because the measurements {9 , 10} are equally important.

But in reality, observation error and model error might not equal. That is why we need the Kalman gain.

The figure shows flow chart of Kalman filter

t=linspace(1,length(y),length(y));

r=1;q=0.01; %error of Obs&model

p0=10;x0=-10; %initial value(guess)

x=x0;p=p0;sav_x=[];

for i=1:length(y)

x_proc=x;

p_proc=p+q;

kg=p_proc*(1/(p_proc+r));

x=x_proc+kg*(y(i)-x_proc);

p=(1-kg)*p_proc;

sav_x=[sav_x;x]

end

plot(t,y,'k.',t,x_sav,'r-')