| 23 | | int main() { |
| 24 | | // Pocet dat |
| 25 | | int Ndat = 900; |
| 26 | | |
| 27 | | // Objekt pro ukladani vysledku |
| 28 | | memlog L(Ndat); |
| 29 | | |
| | 23 | int main( int argc, char* argv[] ) { |
| | 24 | /* const char *fname; |
| | 25 | if ( argc>1 ) {fname = argv[1]; } |
| | 26 | else { fname = "k1.cfg"; } |
| | 27 | UIFile F ( fname ); //protected by exceptions, should complain if not found*/ |
| | 28 | |
| | 29 | FileDS DS("data.it","Data"); // Data Source |
| | 30 | int ndat = DS.ndat(); // number of data |
| | 31 | memlog L(ndat); // Logger |
| | 32 | string resfile="k1_results.it"; // name of output file |
| | 33 | vec dQ="0.1 0.1"; |
| | 34 | vec dR="0.2 0.2"; // TODO: read from config file <== broken on windows |
| | 35 | |
| 41 | | vec Qdiag ( "1.0 10. 10 10" ); |
| 42 | | vec Rdiag ( "1 1" ); //var(diff(xth)) = "0.034 0.034" |
| 43 | | mat Q =diag( Qdiag ); |
| 44 | | mat R =diag ( Rdiag ); |
| 45 | | EKFfull Efix ( RVstav,RVpozor,RVut ); |
| 46 | | // pocatecni podminky |
| 47 | | Efix.set_est ( mu0, 1*eye ( 4 ) ); // nulova |
| 48 | | // nastaveni modelu pro EKF |
| 49 | | Efix.set_parameters ( &fxu,&hxu,Q,R); |
| | 46 | EKFCh Efix; //Extended KF s fix. variancemi |
| | 47 | Efix.set_parameters ( &fxu,&hxu,diag(dQ),diag(dR)); |
| | 48 | Efix.set_statistics ( zeros(fxu.dimension()), 100*eye ( fxu.dimension() ) ); // nulova |
| 51 | | int L_xt = L.add(RVstav, "xt"); // Tady se rika jak velky vektor (pomoci obj. RV) se bude logovat |
| 52 | | // A jak se bude jsmenovat vysledek |
| 53 | | int L_ut = L.add(RVut, "ut"); // Tady se rika jak velky vektor (pomoci obj. RV) se bude logovat |
| 54 | | // A jak se bude jsmenovat vysledek |
| 55 | | int L_mean = L.add(RVstav, "odh_xt"); // Tady se rika jak velky vektor (pomoci obj. RV) se bude logovat |
| 56 | | // A jak se bude jsmenovat vysledek |
| 57 | | |
| | 50 | // Definovat co se bude logovat |
| | 51 | Efix.log_add(&L,"E"); |
| 59 | | // Priprava poli pro simulaci |
| 60 | | vec ut(RVut.count()); |
| 61 | | vec xt(RVstav.count()); |
| 62 | | vec dt(RVpozor.count()); |
| 63 | | // minuly stav |
| 64 | | vec xtm=zeros(RVstav.count()); // nulovy pocatecni stav |
| 65 | | for ( int t=1;t<Ndat;t++ ) { |
| 66 | | // Nastaveni vstupu |
| 67 | | ut(0) = 1+sin((double)t/10); // V ut jsou same jednicky |
| 68 | | ut(1) = 1+cos((double)t/10); // V ut jsou same jednicky |
| 69 | | |
| 70 | | // Generovani DAT modelem |
| 71 | | xt = fxu.eval(xtm,ut); |
| 72 | | dt = hxu.eval(xt,ut); |
| 73 | | xtm = xt; //save xt for the next step |
| 74 | | |
| 75 | | //ESTIMATE |
| 76 | | Efix.bayes(concat(dt,ut)); |
| | 53 | |
| | 54 | vec dt; |
| | 55 | for ( int t=1;t<ndat;t++ ) { |
| | 56 | DS.getdata(dt); // dt is allocated |
| | 57 | Efix.bayes(dt); //ESTIMATE |