主题:请求懂KALMAN FILTER的给帮忙看看,快绝望了.
总出错误,实在不知道该怎么改了.
请高手们帮帮忙.
是想做一个time-varying system
function [Yt,X,v,REV] = Kalman_ZZ1(Y,X0,L,A,C,P0,Q,R0);
%
% Introduction
% Time-Varying Kalman Filter
% -----------------------------------------------------------
% State Space Model
% X(t+1) = A*X(t) + B*u(t)+G*w(t) Q = E(w(t)*w(t)')
% Yv(t) = C*X(t) + v(t) R = E(v(t)*v(t)')
%
% Input Y = onservations
% X0 = initial state vector
% A = state matrix
% C = observation matrix
% P0 = initial state covariance
% Q = state noise covariance
% R = observation noise covariance
% L =
%
% Output Ye = posteriori estimates vector
% X = parameter vector
% v = error process
% REV = relative error variance MSE/MSY
% ------------------------------------------------------------
T = length(Y);
R = zeros(T,1);
X = X0(ones(T,1),:);
Yt = zeros(T,1);
v = zeros(T,1);
for t=1:T;
if t>1;
Xt = X0;
X_t = X0;
P_t = P0;
Yv = Y + v;
R0 = R(t-1);
R(t) = R0;
% Measurement update
K(t) = P_t(t)*C'/(C*P_t(t)*C'+R(t)); % K(t) = Kalman gain
Xt(t) = X_t(t) + K(t)*(Yv(t)-C*X_t(t)); % X(t) = x[t|t]
P(t) = P_t(t)-K(t)*C*P_t(t); % P(t) = P[t|t]
Yt(t) = C*X(t); % Yt(t) = Y[t|t]
% Time update
X_t(t+1) = A*Xt(t); % x[n+1|n]
P_t(t+1) = A*P(t)*A' + Q; % P[n+1|n]
% Prediction Error
e(t) = Y(t) - C*Xt(t+1);
v(t) = e(t);
if ~isnan(Y);
e2(t) = e(t)*e(t);
Yt(t) = A*X(t);
UC = 0;
ESU(t) = A'*Yt(t);
V0 = V(t-1);
V(t) = V0*(1-UC)+UC*e2;
K(t) = Yt(t) / (ESU(t) + V0); % Kalman Gain
Xt(t) = X_t(t) + K(t)'*e(t);
R(t) = R(t);
end;
end;
end;
请高手们帮帮忙.
是想做一个time-varying system
function [Yt,X,v,REV] = Kalman_ZZ1(Y,X0,L,A,C,P0,Q,R0);
%
% Introduction
% Time-Varying Kalman Filter
% -----------------------------------------------------------
% State Space Model
% X(t+1) = A*X(t) + B*u(t)+G*w(t) Q = E(w(t)*w(t)')
% Yv(t) = C*X(t) + v(t) R = E(v(t)*v(t)')
%
% Input Y = onservations
% X0 = initial state vector
% A = state matrix
% C = observation matrix
% P0 = initial state covariance
% Q = state noise covariance
% R = observation noise covariance
% L =
%
% Output Ye = posteriori estimates vector
% X = parameter vector
% v = error process
% REV = relative error variance MSE/MSY
% ------------------------------------------------------------
T = length(Y);
R = zeros(T,1);
X = X0(ones(T,1),:);
Yt = zeros(T,1);
v = zeros(T,1);
for t=1:T;
if t>1;
Xt = X0;
X_t = X0;
P_t = P0;
Yv = Y + v;
R0 = R(t-1);
R(t) = R0;
% Measurement update
K(t) = P_t(t)*C'/(C*P_t(t)*C'+R(t)); % K(t) = Kalman gain
Xt(t) = X_t(t) + K(t)*(Yv(t)-C*X_t(t)); % X(t) = x[t|t]
P(t) = P_t(t)-K(t)*C*P_t(t); % P(t) = P[t|t]
Yt(t) = C*X(t); % Yt(t) = Y[t|t]
% Time update
X_t(t+1) = A*Xt(t); % x[n+1|n]
P_t(t+1) = A*P(t)*A' + Q; % P[n+1|n]
% Prediction Error
e(t) = Y(t) - C*Xt(t+1);
v(t) = e(t);
if ~isnan(Y);
e2(t) = e(t)*e(t);
Yt(t) = A*X(t);
UC = 0;
ESU(t) = A'*Yt(t);
V0 = V(t-1);
V(t) = V0*(1-UC)+UC*e2;
K(t) = Yt(t) / (ESU(t) + V0); % Kalman Gain
Xt(t) = X_t(t) + K(t)'*e(t);
R(t) = R(t);
end;
end;
end;