【数字信号去噪】基于matlab变分贝叶斯卡尔曼滤波器数字信号滤波【含Matlab源码 2256期】

news2024/11/17 21:45:57

⛄一、变分贝叶斯卡尔曼滤波器数字信号滤波

1变分贝叶斯
在参数估计的问题中,在获得了观测样本的数据集Z后,根据贝叶斯准则,核心是参数集 θ 的后验概率密度函数的计算
在这里插入图片描述
而式( 1) 计算的一个难点在于分母,边缘似然概率密度函数p( Z) 的计算。正如引言所述,通常情况下p( Z) 的计算难以得到精确的解析解,引入变分贝叶斯方法,引入一个简单的近似分布函数q( θ) , 并取p( Z) 的对数形式:
在这里插入图片描述

可以发现,每个参数 θi的近似分布可以通过求对数联合概率密度函数关于其他参数分布q( θj≠i) 的期望求得,所以每一个参数分布的计算都依赖于其他参数的分布。这就形成了迭代的机制: 在给定先验知识的情况下,初始化参数值,通过变分贝叶斯算法循环迭代计算,进行参数更新,直至自由能量( 对数边缘似然函数的下界) 达到最大值,判定算法收敛并结束, 此时可以得到系统模型参数的估计值。

⛄二、部分源代码

clear all;
close all;
%%%%%Model parameters%%%%%%%
nxp = 10;
nx = 4; % number of state variables
nz = 2; % number of measures
T = 1;
q = 10;
r = 1;
F = [eye(2) Teye(2);zeros(2) eye(2)];
H = [eye(2) zeros(2)];
Q0 = q * [T^3/3eye(2) T^2/2*eye(2);T2/2eye(2) Teye(2)];
R0 = r * [1 0.5;0.5 1];
L = 500;
Tn = L;
N = 5; %%%%%%The number of variational iteration
times_Of_R = 5;
%%%%%Initial values
x_Ini = [100;100;10;10];
P_Ini = diag([100 100 100 100]);
X = zeros(4,L);
XKF_True = zeros(4,L);
XKF_Const = zeros(4,L);
XKF_VB = zeros(4,L);

mse_Kf_1 = zeros(L,nxp);
mse_Kf_2 = zeros(L,nxp);
mse_Ktf_1 = zeros(L,nxp);
mse_Ktf_2 = zeros(L,nxp);
mse_VB_1 = zeros(L,nxp);
mse_VB_2 = zeros(L,nxp);
for k = 1:nxp
x=x_Ini;
X(:,1) = x;
%%%%Kalman filter with nominal noise covariance matrices (KFNCM)
XKF_Const(:,1) = x;
x1=x;
P_Const = P_Ini;
%%%%Kalman filter with true noise covariance matrices (KFTCM)
XKF_True(:,1) = x;
x2 = x;
P_True = P_Ini;

%%%%Kalman filter of variational Bayesian Approximations
XKF_VB(:,1) = x;
x3 = x;
P_VB = P_Ini;
alfa = [1 1]';
beta = diag(R0);
mk = x;
for t = 2:Tn
    %%%%True noise covariance matrices
    Q = (1+0.5*cos(pi*t/Tn))*Q0;
    R = (1+0.6*cos(pi*t/Tn))*R0;
    %%%%Square-root of noise covariance matrices
    SQ = utchol(Q);
    SR = utchol(R);
    %%%%Simulate true state and measurement
    x = F*x+SQ*randn(nx,1);
    z = H*x+SR*randn(nz,1);
    X(:,t) = x;
    %%%%Filtering
    [x1,P_Const,Ppf] = kf(x1,P_Const,F,H,z,Q0,R0*times_Of_R);
    [x2,P_True,Pptf] = kf(x2,P_True,F,H,z,Q,R);
    [x3,P_VB,alfa,beta,mk] = vbkf(x3,P_VB,alfa,beta,mk,F,H,z,Q0,N);
    %%%%Save data
    XKF_Const(:,t) = x1;
    XKF_True(:,t) = x2;
    XKF_VB(:,t) = x3;
end