卡尔曼滤波实现一阶马尔可夫形式的滤波|价格滤波|MATLAB代码|无需下载，复制后即可运行【付费专栏试读】

本文链接：https://blog.csdn.net/callmeup/article/details/141816946

在这里插入图片描述

一阶马尔可夫

一阶马尔可夫噪声是一种具有马尔可夫性质的随机过程。在这种噪声中，当前时刻的状态只与前一时刻的状态有关，与更早的状态无关。

一阶马尔可夫噪声可以用一个状态转移矩阵表示，矩阵的每个元素表示从一个状态转移到另一个状态的概率。

滤波模型

状态量的迭代模型如下：
在这里插入图片描述
观测量为X的第一维，所以观测方程也就是取X的第一维。

运行结果

应用背景为价格滤波，所以对比X真值和滤波值的第一维（对应的是价格），曲线如下：
在这里插入图片描述

对比滤波前和滤波后的价格误差，曲线如下：

在这里插入图片描述
误差的CDF图像显示了各种误差下的概率，如下：

在这里插入图片描述
蓝线KF为滤波后的曲线，红线为滤波前的，能看出来滤波后几乎所有时刻的误差都在1以内，而滤波前仅由10%（纵坐标0.1）的情况对应了1的误差。很明显：

滤波明显改善了价格估计的精度

源代码

源代码如下：

% 简单的马尔可夫价格（一阶马尔可夫噪声）KF预测，留EKF接口（三维状态量，留一维备用）
% author:matlabfilter
% date: 2024-01-19/Ver1
% 2024-09-02/Ver2：添加逐行注释
clear;clc;close all;
rng(0); %固定随机数种子（保证每次滤波得到的结果一样）
%% initial
T = 0.1; %采样率
t = [T:T:100]; %构建时间序列
Q = 0.1*diag([1,1,1]);w=sqrt(Q)*randn(size(Q,1),length(t)); %系统噪声
a = 1*randn(1,length(t)); %构建马尔可夫型数据的系数，服从方差为1、均值为0的正态分布
R = 0.01*diag([1]);v=sqrt(R)*randn(size(R,1),length(t)); %状态噪声
P = 0.1*eye(3); %状态协方差
P_num = zeros(length(t),size(P,1),size(P,2)); %存放每次迭代的P
P_num(1,:,:) = P; %记录协方差
X=zeros(3,length(t)); %给真实值预留空间
X_kf=zeros(3,length(t)); %准备储存滤波后的值
X_kf(:,1)=X(:,1); %给滤波值赋初值
Z=zeros(1,length(t)); %定义观测值形式
Z(:,1)=X(1,1)+v(:,1); %观测量
%% model
X_=zeros(3,length(t)); %给未滤波值预留空间
X_(:,1)=X(:,1); %给未滤波值赋初值
for i1 = 2:length(t) %构造for循环
    X(:,i1) = [X(1,i1-1) + T*X(2,i1-1);
        X(2,i1-1)+T*a(i1);
        0]; %真实值迭代
    X_(:,i1) = [X_(1,i1-1) + T*X(2,i1-1);
        X_(2,i1-1)+T*a(i1);
        0] + w(:,i1-1);%未滤波的值迭代
    Z(:,i1) = X(1,i1) + v(i1); %观测值计算
end
%% EKF
w=sqrt(Q)*randn(size(Q,1),length(t)); %生成状态误差omega
for k = 2 : length(t) %使用for构造滤波循环
    Xpre = [X_kf(1,k-1)+ T*X_kf(2,k-1);
        X_kf(2,k-1)+T*a(i1);
        0]+w(:,k-1); %状态一步转移（状态预测）
    Z_hat =Xpre(1); %预测的观测量
    F = [1,T,0;
        0,T,0;
        0,0,0]; %状态转移矩阵
    H = [1,0,0]; %观测矩阵
    PP=F*P*F'+Q; %协方差计算
    Kk=PP*H'/(H*PP*H'+R); %计算增益
    X_kf(:,k)=Xpre+Kk*(Z(:,k)-Z_hat); %状态预测，此状态为滤波输出状态
    P=PP-Kk*H*PP; %协方差更新
    P_num(k,:,:) = P; %存储协方差
end


%% display
figure; %创建绘图窗口
plot(t,X(1,:),t,X_kf(1,:)); %绘制真值和KF滤波后的值
title('价格对比'); %添加图片标题
legend('real','KF'); %给图片添加图例

figure; %创建绘图窗口
plot(t,X_kf(1,:)-X(1,:),t,X_(1,:)-X(1,:));
title('价格估计误差对比'); %添加图片标题
legend('KF','without KF'); %给图片添加图例
% [~,p1] = max(abs(X_kf(1,:)-X(1,:)));
% hh = gca();
% h = hh.Children(1);
% datatip(h,p1*T-1,X_kf(1,p1)-X(1,p1));
% [~,p2] = max(abs(X_(1,:)-X(1,:)));
% datatip(h,p2*T-1,X_(1,p2)-X(1,p2));

figure; %创建绘图窗口
cdfplot(abs(X_kf(1,:)-X(1,:))); %绘制cdf图像
hold on %固定窗口内容，后续绘制的曲线在原有图像上覆盖
cdfplot(abs(X_(1,:)-X(1,:))); %绘制未滤波时误差的CDF图像
legend('KF','without KF'); %给图片添加图例
title('累计概率密度函数'); %添加图片标题