kfdemo.m

kalman filter在MUTILAB下的卡尔曼滤波功能

disp(' ----------------------');
disp('|Bayesian Kalman filter|');
disp('|Stephen Roberts - IMPC|');
disp(' ----------------------');
disp(' ');
disp('to get help about Kalman filter program');
disp('just type help kf');
pause;
load kfdemo;
disp('loaded the demo file');
disp('First example - sine wave that goes to noise');
subplot(4,1,1),plot(x1);
disp('now run the kf over it with model-order of 2');
disp('and initial noise variance set to 1');
[A,y,Ea,Ey,P] = kf(x1,x1,2,1);
disp('now plot the prediction along with error bars');
subplot(4,1,2),errorbar(y,sqrt(Ey));axis([0 200 -5 5]);
disp('Note how the error bars get larger in the noise');
disp('but also note the problem with a *causal* filter-');
disp('we have a settling down time after the noise');
disp('(really we need an *anti-causal* filter as well...)');
pause;
disp(' ');
disp('Now plot one of the model coefficients with error bars');
subplot(4,1,3),errorbar(A(:,1),sqrt(Ea(:,1)));axis([0 200 -2 2]);
disp('Note how the error bars on the model get larger');
disp(' ');
pause;
disp('Now plot the log10 of the innovations probability (likelihood)');
subplot(4,1,4),plot(log10(P+1e-10));axis([0 200 -10 2]);
		% the +1e-10 stops log of zero problems!
disp('Note how the p(innov) goes low at the state changes');
pause;
clf;
disp('Now the EEG example I showed in Marburg');
disp(' ');
disp('approx 2s of EEG activity');
subplot(3,1,1),plot(x_eeg);
pause;
disp('some noise');
subplot(3,1,2),plot(noise);
pause;
disp('added together');
subplot(3,1,3),plot(x_eeg+noise);
disp(' ');
disp('Run the KF over this last sequence, model order = 5');
[A,y,Ea,Ey,P] = kf(x_eeg+noise,x_eeg+noise,5,1);
subplot(4,1,1),plot(x_eeg+noise);
disp('plot the estimated signal noise variance');
subplot(4,1,2),plot(Ey);
pause;
disp('plot the estimated feature noise variance');
subplot(4,1,3),plot(Ea(:,1));
pause;
disp('Now plot the log10 of the innovations probability (likelihood)');
subplot(4,1,4),plot(log10(P));