Changeset 9 for matlab

Show
Ignore:
Timestamp:
01/23/08 13:27:35 (16 years ago)
Author:
smidl
Message:
 
Location:
matlab
Files:
1 added
1 modified

Legend:

Unmodified
Added
Removed
  • matlab/testKF.m

    r8 r9  
    1 if 0 
    2 A=[1 -0.5; 1 0]; 
    3 B=[1;0]; 
    4 C=[1 0]; 
    5 D=[0]; 
    6 R=[0.1]; 
    7 Q=[0.2 0 ; 0 0.02];  
     1function testKF(skipgen) 
     2if nargin<1, skipgen=0; 
    83 
    9 sQ = chol(Q)'; 
    10 sR = chol(R)'; 
     4if ~skipgen 
     5        A=[1 -0.5; 1 0]; 
     6        B=[1;0]; 
     7        C=[1 0]; 
     8        D=[0]; 
     9        R=[0.1]; 
     10        Q=[0.2 0 ; 0 0.02];  
    1111 
    12 N =1000; 
    13 mu0 = [0;0]; 
    14 P0 = 200*eye(2); 
     12        sQ = chol(Q)'; 
     13        sR = chol(R)'; 
    1514 
    16 u = ones(1,N); 
    17 x = zeros(2,N); 
    18 y = zeros(1,N); 
     15        N =1000; 
     16        mu0 = [0;0]; 
     17        P0 = 200*eye(2); 
    1918 
    20 x(:,1) = [10;10]; 
    21 for i=2:N; 
    22         et(1:2,i) = sQ*randn(2,1); 
    23         x(:,i) = A*x(:,i-1) + B*u(i) + et(:,i); 
    24         y(:,i) = C*x(:,i) + D*u(i) + sR*randn(1,1); 
     19        u = ones(1,N); 
     20        x = zeros(2,N); 
     21        y = zeros(1,N); 
     22 
     23        x(:,1) = [10;10]; 
     24        for i=2:N; 
     25                et(1:2,i) = sQ*randn(2,1); 
     26                x(:,i) = A*x(:,i-1) + B*u(i) + et(:,i); 
     27                y(:,i) = C*x(:,i) + D*u(i) + sR*randn(1,1); 
     28        end 
     29 
     30        d=[y;u]; 
     31        itsave('testKF.it',d,A,B,C,D,Q,R,P0,mu0) 
     32        save testKF 
     33else 
     34        load testKF 
    2535end 
    26  
    27 d=[y;u]; 
    28 itsave('testKF.it',d,A,B,C,D,Q,R,P0,mu0) 
    29 save testKF 
    30 else 
    31 %!cd ../testKF 
    32  
    33 load testKF 
    3436% init  
    3537mu = mu0; 
     
    4648        K = P*C'*iRy;  
    4749        P = P- K*C*P; % P = P -KCP; 
    48         EP(:,t) = eig(P);  
    49         if any(EP<0), keyboard; end 
    5050        mu = mu + K*(y(t)-C*mu-D*u(t)); 
    5151        Mu(1:2,t)=mu; 
    5252end 
    5353 
    54  
    55  
     54!./testKF 
    5655itload('testKF_res.it'); 
    5756