roboforum.ruТехнический форум по робототехнике. |
|
|
# Kalman filter example demo in Python
# A Python implementation of the example given in pages 11-15 of "An
# Introduction to the Kalman Filter" by Greg Welch and Gary Bishop,
# University of North Carolina at Chapel Hill, Department of Computer
# Science, TR 95-041,
# http://www.cs.unc.edu/~welch/kalman/kalmanIntro.html
# by Andrew D. Straw
import numpy
import pylab
# intial parameters
n_iter = 50
sz = (n_iter,) # size of array
x = -0.37727 # truth value (typo in example at top of p. 13 calls this z)
z = numpy.random.normal(x,0.1,size=sz) # observations (normal about x, sigma=0.1)
Q = 1e-5 # process variance
# allocate space for arrays
xhat=numpy.zeros(sz) # a posteri estimate of x
P=numpy.zeros(sz) # a posteri error estimate
xhatminus=numpy.zeros(sz) # a priori estimate of x
Pminus=numpy.zeros(sz) # a priori error estimate
K=numpy.zeros(sz) # gain or blending factor
R = 0.1**2 # estimate of measurement variance, change to see effect
# intial guesses
xhat[0] = 0.0
P[0] = 1.0
for k in range(1,n_iter):
# time update
xhatminus[k] = xhat[k-1]
Pminus[k] = P[k-1]+Q
# measurement update
K[k] = Pminus[k]/( Pminus[k]+R )
xhat[k] = xhatminus[k]+K[k]*(z[k]-xhatminus[k])
P[k] = (1-K[k])*Pminus[k]
pylab.figure()
pylab.plot(z,'k+',label='noisy measurements')
pylab.plot(xhat,'b-',label='a posteri estimate')
pylab.axhline(x,color='g',label='truth value')
pylab.legend()
pylab.xlabel('Iteration')
pylab.ylabel('Voltage')
pylab.figure()
valid_iter = range(1,n_iter) # Pminus not valid at step 0
pylab.plot(valid_iter,Pminus[valid_iter],label='a priori error estimate')
pylab.xlabel('Iteration')
pylab.ylabel('$(Voltage)^2$')
pylab.setp(pylab.gca(),'ylim',[0,.01])
pylab.show()
/** A simple kalman filter example by Adrian Boeing
www.adrianboeing.com
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double frand() {
return 2*((rand()/(double)RAND_MAX) - 0.5);
}
int main() {
//initial values for the kalman filter
float x_est_last = 0;
float P_last = 0;
//the noise in the system
float Q = 0.022;
float R = 0.617;
float K;
float P;
float P_temp;
float x_temp_est;
float x_est;
float z_measured; //the 'noisy' value we measured
float z_real = 0.5; //the ideal value we wish to measure
srand(0);
//initialize with a measurement
x_est_last = z_real + frand()*0.09;
float sum_error_kalman = 0;
float sum_error_measure = 0;
for (int i=0;i<30;i++) {
//do a prediction
x_temp_est = x_est_last;
P_temp = P_last + Q;
//calculate the Kalman gain
K = P_temp * (1.0/(P_temp + R));
//measure
z_measured = z_real + frand()*0.09; //the real measurement plus noise
//correct
x_est = x_temp_est + K * (z_measured - x_temp_est);
P = (1- K) * P_temp;
//we have our new system
printf("Ideal position: %6.3f \n",z_real);
printf("Mesaured position: %6.3f [diff:%.3f]\n",z_measured,fabs(z_real-z_measured));
printf("Kalman position: %6.3f [diff:%.3f]\n",x_est,fabs(z_real - x_est));
sum_error_kalman += fabs(z_real - x_est);
sum_error_measure += fabs(z_real-z_measured);
//update our last's
P_last = P;
x_est_last = x_est;
}
printf("Total error if using raw measured: %f\n",sum_error_measure);
printf("Total error if using kalman filter: %f\n",sum_error_kalman);
printf("Reduction in error: %d%% \n",100-(int)((sum_error_kalman/sum_error_measure)*100));
return 0;
}
Мне понятна в принципе работа ФК
все это калман нам отфильтровывает и выдает уже нормальное значение
Сейчас этот форум просматривают: нет зарегистрированных пользователей и гости: 2