• 卡尔曼滤波


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

    %=====================================================================
    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;
    
    %---------------------------------------------------------------------
  • 相关阅读:
    php mysql增删查改
    php文件上传代码解析
    php验证码
    html加javascript和canvas类似超级玛丽游戏
    javascript图片轮播
    php分页
    《Hadoop权威指南》笔记 第二章 Hadoop Streaming
    《Hadoop权威指南》笔记 第一章&第二章 MapReduce初探
    《深入理解Java虚拟机》笔记 第七章 虚拟机加载机制及双亲委派模型
    《深入理解Java虚拟机》笔记 第三章 内存分配与回收策略
  • 原文地址:https://www.cnblogs.com/vincentQin/p/5406012.html
Copyright © 2020-2023  润新知