Révision 455
Ajouté par Henri JOUANNY il y a environ 3 ans
| branch/jouanny2/sp4a3/sp4a3_kalman.c | ||
|---|---|---|
|
double FPFt[4][4] = {{0, 0, 0, 0},
|
||
|
{0, 0, 0, 0},
|
||
|
{0, 0, 0, 0},
|
||
|
{0, 0, 0, 0}}; //création d'une sous-matrice pour la multiplication de FP avec F'
|
||
|
{0, 0, 0, 0}}; //création d'une sous-matrice pour la multiplication de FP avec F'
|
||
|
double Ht[4][2] = {{0, 0,},
|
||
|
{0, 0},
|
||
|
{0, 0},
|
||
|
{0, 0}}; //Transposée de H
|
||
|
double PHt[4][2] = {{0, 0},
|
||
|
{0, 0},
|
||
|
{0, 0},
|
||
|
{0, 0}};
|
||
|
double HPHt[2][2] = {{0, 0,},
|
||
|
{0, 0}};
|
||
|
double HPHtR[2][2] = {{0, 0,},
|
||
|
{0, 0}};
|
||
|
double Kobs[4][1] = {{0},
|
||
|
{0},
|
||
|
{0},
|
||
|
{0}};
|
||
|
double HX [2][1] = {{0},{0}};
|
||
|
double KobsmHX[2][1] = {{0},{0}};
|
||
|
|
||
|
tests_unitaires();
|
||
|
|
||
|
FILE* fichier = fopen("pos_t_x_y.dat","r");
|
||
| ... | ... | |
|
|
||
|
// kalman param
|
||
|
double sigma_etat = 10.0;
|
||
|
double sigma_observation = 2.0;
|
||
|
double sigma_observation = 2.0;
|
||
|
double Obs[2][1] = {{xobs},{yobs}};
|
||
|
double X[4][1] = {{0},{0},{0},{0}};
|
||
|
|
||
|
double P[4][4] = {{sigma_etat*sigma_etat, 0, 0, 0},
|
||
| ... | ... | |
|
// Kalman
|
||
|
|
||
|
// X = F*X
|
||
|
Mul_Mat_Mat(4,4,F,4,1,F,X);
|
||
|
Mul_Mat_Mat(4,4,F,4,1,F,X);
|
||
|
|
||
|
Plot_Mat(X," X(k+1|k) = ");
|
||
|
|
||
|
//P = F*P*F'+Q;
|
||
|
Mul_Mat_Mat(4,4,F,4,4,P,FP);
|
||
|
Mul_Mat_Mat(4,4,FP,4,4,Transpose_Mat(4,4,F,Ft),FPFt);
|
||
|
Add_Mat_Mat(4,4,FPFt,4,4,Q,P);
|
||
|
Transpose_Mat(4,4,F,Ft);
|
||
|
Mul_Mat_Mat(4,4,FP,4,4,Ft,FPFt);
|
||
|
Add_Mat_Mat(4,4,FPFt,4,4,Q,P);
|
||
|
|
||
|
Plot_Mat(P,"P(k+1|k) = F.P(k|k).FT + Q = ");
|
||
|
|
||
|
// K = P*H' / ( H*P*H' + R);
|
||
|
// K = P*H' / ( H*P*H' + R);
|
||
|
Transpose_Mat(2,4,H,Ht);
|
||
|
Mul_Mat_Mat(4,4,P,4,2,Ht,PHt);
|
||
|
Mul_Mat_Mat(2,4,H,4,2,PHt,HPHt);
|
||
|
Add_Mat_Mat(2,2,HPHt,2,2,R,HPHtR);
|
||
|
|
||
|
for (i=0; i<3; i++)
|
||
|
{
|
||
|
for (j=0; j<1; j++)
|
||
|
{
|
||
|
K[i][j]=(PHt[i][j]/HPHtR[i][j]); //assigne les valeurs à K
|
||
|
}
|
||
|
}
|
||
|
Plot_Mat(K,"K = ");
|
||
|
|
||
|
//X = X + K*([xobs(i);yobs(i)]-H*X);
|
||
|
//X = X + K*([xobs(i);yobs(i)]-H*X);
|
||
|
Mul_Mat_Mat(4,2,K,2,1,Obs,Kobs);
|
||
|
Mul_Mat_Mat(2,4,H,4,1,X,HX);
|
||
|
Sub_Mat_Mat(2,1,Kobs,2,1,HX,KobsmHX);
|
||
|
Add_Mat_Mat(2,1,X,2,1,KobsmHX,X);
|
||
|
|
||
|
//Sub_Mat_Mat(2,1,Obs,2,1), pas eu le temps de finir cette ligne
|
||
|
//Plot_Mat(Delta,"DELTA = Obs - H.X(k+1|k)");
|
||
|
Plot_Mat(X," X(k+1|k+1) = X(k+1|k) + K.Delta = ");
|
||
|
|
||
Formats disponibles : Unified diff
-Filtre non finit, arrêt au calcul de Xk+1 en connaissant Xk