zoukankan      html  css  js  c++  java
  • 卡尔曼滤波

    原理部分就不过多介绍了,可以参考其他博主的博客。
    这是一个自己写的卡尔曼滤波函数,其中作图部分不是很完善(只能输出二个状态),可以根据自己的需要进行修改,请务必注明出处,谢谢!如有错误,请大家指正。

    %=====================================================================
    function Kalman_By_Vincent(x,x1,A,F,H,w,v,Q,R,P)
    
    % Kalman 滤波的实现
    % x是初始化的状态向量;
    % x1是初始化的预测的状态向量;一般为0向量,或者与x一致
    % A是系统矩阵;
    % H是观测阵;
    % F是状态误差阵;
    % w,v分别是状态噪声和观察噪声;
    % Q和R分别是w和v的方差;
    % P是滤波误差方差阵,一般初始化P{1}=eye(Dim);
    % Author:Vincentqin;
    % Date:2015年12月12日。
    
    Dim  =length(x{1});   %状态维度
    N    =size(w,2);      %噪声的个数
    %% 状态信号和观测信号,用来模拟实际的状态和观测
    for k=1:N-1
        x{k+1}=A*x{k}+F*w(k);   
    end
    for k=1:N
        y(k)  =H*x{k}+v(k);
    end
    
    %% Kalman迭代
    for k=2:N
        %%  预测
        x1_pre    =A*x1{k-1};                     %状态预测值
        y1_pre    =H*x1_pre;                      %测量预测值    
        p1{k-1}   =A*P{k-1}*A'+F*Q*F';            %预报误差方差阵更新  
        %% 更新
        y1_NewMeg =y(k)-y1_pre;                   %新息
        S         =H*p1{k-1}*H'+R;                %新息方差阵
        K{k-1}    =p1{k-1}*H'*inv(S);             %增益矩阵更新
        P{k}      =(eye(Dim)-K{k-1}*H)*p1{k-1};   %滤波误差方差阵更新
        x1{k}     =x1_pre+K{k-1}*y1_NewMeg;       %状态预测更新    
    end
    %cell数据矩阵化
    new_x   =cell2mat(x);
    new_p1  =cell2mat(p1);      %预报方差阵  
    new_K   =cell2mat(K);      %增益矩阵
    new_P   =cell2mat(P);      %滤波方差阵
    new_x1  =cell2mat(x1);     %预测
    
    %% 预测数据
    y1=H*new_x1;
    %% 显示输出图像
    figure(3)
    set(gcf,'position',[400 150 800 600]);
    plot(y,'b','linewidth',2);
    grid on;hold on;
    plot(y1,'r');
    legend('原始观测信号','Kalman输出观测信号');
    % title('原始信号与预测信号示意图','fontsize',15);
    hold off;
    
    figure(4)%状态量1
    set(gcf,'position',[400 150 800 600]);
    plot(new_x(1,:),'b','linewidth',2);grid on;hold on;
    plot(new_x1(1,:),'r','linestyle','-');
    legend('状态一','Kalman滤波输出状态一');
    hold off;
    
    figure(5)%状态量2
    set(gcf,'position',[400 150 800 600]);
    plot(new_x(2,:),'b','linewidth',2);grid on;hold on;
    plot(new_x1(2,:),'r','linestyle','-');
    legend('状态二','Kalman滤波输出状态二');
    hold off;
    
    figure(6)%增益
    set(gcf,'position',[400 150 800 600]);
    plot(new_K(1,:),'b','linewidth',2);grid on;hold on;
    plot(new_K(2,:),'g','linewidth',2);
    legend('状态一增益','状态二增益');
    hold off;
    
    %---------------------------------------------------------------------
  • 相关阅读:
    Android Service启动原理分析
    线程池原理分析
    仿EventBus做一个简单的基于订阅发布的事件总线
    EventBus原理以及源代码分析
    Android从点击Launcher图标开始到App打开流程分析
    使用LruCache和DiskLruCache手写一个ImageLoader
    OkHttp2连接池复用原理分析
    OkHttp执行流程源码分析
    Android使用动态代理模仿Retrofit的create方法,使其可以返回任意的接口类型
    Android模仿Retrofit的建造者模式
  • 原文地址:https://www.cnblogs.com/vincentQin/p/5406012.html
Copyright © 2011-2022 走看看