设计卡尔曼滤波器

问题描述 投票:0回答:1

所以,我被要求设计一个卡尔曼滤波器,并运行一些模拟来看看结果.给定的方程。

x(k+1) = -0.8*x(k) + sin(0.1k) + w(k) y(k) = x(k) + v(k) 哪儿 w(k)v(k) 是处理和测量的噪音 w(k)~N(0,0.01) v(k)~N(0,1)

我想确保我的代码中一切都能正常工作,并且有意义。

A = -0.8;
B = 1;
C = 1;

% Covariance matrices
% Processing noise
W = 0.01;
Q=W;

% Measurement noise
V = 1;
R=V;

% Initial conditions
x0 = 0;
P0 = 0;
xpri = x0;
Ppri = P0;
xpost = x0;
Ppost = P0;

% State
Y = zeros(1, size(t,2));
X = Y;
X(1) = x0;

% xpri - x priori
% xpost - x posteriori
% Ppri - P priori
% Ppost - P posteriori

for i = 1:size(t,2)
    %Simulation
    Y(i) = C*sin((i-1)*0.1)+normrnd(0,sqrt(V));

    if i > 1

        % Prediction
        xpri = A*xpost+B*sin((i-1)*0.1);
        Ppri = A*Ppost*A' + Q;

        eps = Y(i) - C*xpri;
        S = C*Ppri*C' + R;
        K = Ppri*C'*S^(-1);
        xpost = xpri + K*eps;
        Ppost = Ppri - K*S*K';
        X(i) = xpost;
    end
end

plot(t, Y, 'b', t, X, 'r')
title('Kalman's Filter')
xlabel('Time')
ylabel('Value')
legend('Measured value', 'Estimated value')

这个KF能正常工作吗?如果不是,是什么问题?

enter image description here

python matlab filter kalman-filter estimation
1个回答
1
投票

代码看起来没问题,结果也是如此。是什么让你怀疑?

这里有一个通用的卡尔曼滤波实现的函数,如果你想仔细检查,你可以检查出来。但是,对于滤波器,它的所有内容都是关于协方差矩阵的调整 P, Q, R...

function [v_f,xyz] = KF(u,z,A,B,C,D,x0,P0,Q,R)
%% initialization
persistent x P
if isempty(x)||isempty(P)
   % state vector
   x = reshape(x0,numel(x0),1);  % ensure vector size

   if nargin < 8 || isempty(P0)
       P0 = eye(length(x)); %default initialization
   emd
   P = P0;    % covariance matrix
end


%% covariance matrices 
if nargin < 9 || isempty(Q)
    Q = diag( ones(length(x) )*1e1;  % proess-noise-covariance
end
% % Q = 0 -> perfect model. Q =/= 0 -> pay more attention to the measurements.
%
if nargin < 10
     R = diag( ones(length(z)) ) *1e1;     % measurement-noise-covariance
end
% The higher R, the less trusting the measurements


%% prediction
x_P = A*x + B*u;
P = A*P*A.' + Q;   % covariance matrix

%% update
H = C;

K = P*H.'*(H*P*H.' +R)^-1;
x = x_P + K*(z - H*x_P);
P = (eye(length(x)) - K*H)*P;

%% Output
y = C*x;
end

.无论如何,我建议把这个问题转到 "国际刑事法院"。编码审查页 因为你显然没有问题,只是想审查你的代码,对吗?

© www.soinside.com 2019 - 2024. All rights reserved.