kalmanMissingData

 wrapper function for the kalman filter toolbox
 code by Fabian Benitez

kalman filter with missing data
input x=the observed data ( dimension by numberSamples), nan where missing
 
      t=the time data indicated wether or not the data is missing
 
 
      A=The system matrix
      H=the observation matrix
      Q=system covariance
      R=observation covariance

0001 % wrapper function for the kalman filter toolbox
0002 % code by Fabian Benitez
0003 %
0004 %kalman filter with missing data
0005 %input x=the observed data ( dimension by numberSamples), nan where missing
0006 %
0007 %      t=the time data indicated wether or not the data is missing
0008 %
0009 %
0010 %      A=The system matrix
0011 %      H=the observation matrix
0012 %      Q=system covariance
0013 %      R=observation covariance
0014 
0015 function [xOut]=kalmanMissingData(x,t,k)
0016 
0017 [numDim tr]=size(x);
0018 A=[eye(numDim) eye(numDim);zeros(numDim) eye(numDim)];
0019 H=[eye(numDim) zeros(numDim)];
0020 Q=eye(numDim*2)*k;
0021 
0022 N=length(t);
0023 idx=1;
0024 idx2=1;
0025 while idx<=N,    
0026     k1=0;
0027     while idx<=N && t(idx)==1 ,
0028         idx=idx+1;
0029         k1=k1+1;
0030     end
0031     if k1>=1,
0032         obs(idx2)=k1;
0033         idx2=idx2+1;
0034     end
0035     k1=0;
0036     while idx<=N && t(idx)<0,
0037         idx=idx+1;
0038         k1=k1+1;
0039     end
0040     if k1>=1,
0041         obs(idx2)=-k1;
0042         idx2=idx2+1;
0043     end
0044 end
0045 
0046 N2=length(obs);
0047 if obs(1)>0,
0048     R=cov(x(:,1:obs(1))')+.0000001*eye(numDim);    
0049     [xfilt Vfilt]=kalman_filter(x(:,1:obs(1)),A,H,Q,R,[x(:,1);zeros(size(x(:,1)))],0);
0050     xPsObs=x(:,1:obs(1));
0051 %     xobs=xfilt(1:numDim,:);
0052     for k1=2:N2,
0053         if obs(k1)<0,
0054             for k2=1:abs(obs(k1)),
0055                 xTemp=A*xfilt(:,end);
0056                 xfilt(:,end+1)=xTemp;
0057                 xPsObs=[xPsObs xTemp(1:numDim)];
0058             end
0059         else
0060             [tr numObs]=size(xPsObs);
0061             xPsObs=[xPsObs x(:,numObs+1:numObs+obs(k1))];
0062             R=cov(xPsObs')+.0000001*eye(numDim);
0063             [xfilt Vfilt]=kalman_filter(xPsObs,A,H,Q,R,[xPsObs(:,1);zeros(size(x(:,1)))],0);            
0064         end        
0065     end
0066 else
0067     x=fliplr(x);
0068     obs=fliplr(obs);
0069     R=cov(x(:,1:obs(1))')+.0000001*eye(numDim);    
0070     [xfilt Vfilt]=kalman_filter(x(:,1:obs(1)),A,H,Q,R,[x(:,1);zeros(size(x(:,1)))],0);
0071     xPsObs=x(:,1:obs(1));
0072 %     xobs=xfilt(1:numDim,:);
0073     for k1=2:N2,
0074         if obs(k1)<0,
0075             for k2=1:abs(obs(k1)),
0076                 xTemp=A*xfilt(:,end);
0077                 xfilt(:,end+1)=xTemp;
0078                 xPsObs=[xPsObs xTemp(1:numDim)];
0079             end
0080         else
0081             [tr numObs]=size(xPsObs);
0082             xPsObs=[xPsObs x(:,numObs+1:numObs+obs(k1))];
0083             R=cov(xPsObs')+.0000001*eye(numDim);
0084             
0085 %             if sum(sum(isnan(R)))
0086 %                 R
0087 %             end
0088             
0089             [xfilt Vfilt]=kalman_filter(xPsObs,A,H,Q,R,[xPsObs(:,1);zeros(size(x(:,1)))],0);            
0090         end        
0091     end
0092     xPsObs=fliplr(xPsObs);
0093 end
0094     
0095 xfilt=kalman_smoother(xPsObs,A,H,Q,R,[x(:,1);zeros(size(x(:,1)))],0);
0096 xOut=xfilt(1:numDim,:);
0097