💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势：🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。
⛳️座右铭：行百里者，半于九十。
📋📋📋本文目录如下：🎁🎁🎁
⛳️赠与读者
👨💻做科研，涉及到一个深在的思想系统，需要科研者逻辑缜密，踏实认真，但是不能只是努力，很多时候借力比努力更重要，然后还要有仰望星空的创新点和启发点。建议读者按目录次序逐一浏览，免得骤然跌入幽暗的迷宫找不到来时的路，它不足为你揭示全部问题的答案，但若能解答你胸中升起的一朵朵疑云，也未尝不会酿成晚霞斑斓的别一番景致，万一它给你带来了一场精神世界的苦雨，那就借机洗刷一下原来存放在那儿的“躺平”上的尘埃吧。
或许，雨过云收，神驰的天地更清朗.......🔎🔎🔎
💥1 概述
基于扩展卡尔曼滤波器和无迹卡尔曼滤波器的窄带信号时变频率估计研究
摘要：窄带信号在通信、雷达、声呐等领域广泛应用，其时变频率的准确估计对信号解调、参数估计及系统控制至关重要。本文聚焦于扩展卡尔曼滤波器（EKF）和无迹卡尔曼滤波器（UKF）在窄带信号时变频率估计中的应用。通过建立窄带信号的非线性状态空间模型，分别运用EKF和UKF进行递归估计，对比分析两者在估计精度、鲁棒性及计算复杂度等方面的性能特点。实验结果表明，UKF在非线性程度较高的情况下具有更高的估计精度和鲁棒性，而EKF计算量相对较小，适用于非线性程度较低的场景。研究为窄带信号时变频率估计提供了有效的理论支持和实践方法。
关键词：窄带信号；时变频率估计；扩展卡尔曼滤波器；无迹卡尔曼滤波器；非线性状态空间模型
一、引言
本文讲解和比较了基于卡尔曼滤波器的频率跟踪方法的能力，例如扩展卡尔曼滤波器 (EKF) 和无味卡尔曼滤波器 (UKF)，以跟踪窄带谐波信号的时变频率。这些结果与 Savaresi 等人的结果进行了比较。在 [3] 中。为了评估算法实现的估计质量，使用了两个标准：性能指数 (PI) 和鲁棒性指数 (RI)，如 [3] 中所述。引入了辅助收敛比以更好地解释 PI 图。考虑两个信号：从现实世界记录的单频正弦信号和生成的噪声信号，该信号具有随时间变化的频率，遵循阶跃分布。该算法能够以~10-7 的精度估计真实世界的信号，并实现与 [3] 一致的性能。在最后一节中显示和讨论了获得的结果。
本文讨论了估计嵌入噪声的窄带谐波信号的频率问题；特别是基于卡尔曼滤波器的方法，例如扩展卡尔曼滤波器 (EKF) 和无迹卡尔曼滤波器 (UKF)。为了评估算法实现的估计质量，引入了两个标准：性能指数（PI）和鲁棒性指数（RI），以及辅助收敛比。已对记录的真实信号和生成的噪声信号进行了测试。详细文章见第4部分下载。
窄带信号作为一种常见的信号类型，其带宽远小于中心频率，广泛存在于通信、雷达、声呐以及生物医学信号处理等众多领域。在无线通信中，精确的载波频率跟踪是实现相干解调和高效数据传输的关键；在雷达和声呐系统中，目标的多普勒频移提供了关于目标速度和运动状态的重要信息，时变频率的准确估计是实现目标跟踪和识别的前提；在生物医学领域，心电图（ECG）和脑电图（EEG）信号的频率成分会随生理状态的变化而改变，对这些时变频率的精确分析有助于疾病的早期诊断和状态监测。
然而，现实环境中的窄带信号往往受到噪声的干扰，且其频率可能随时间发生变化，这使得准确的时变频率估计面临巨大挑战。传统的基于傅里叶变换的方法，如短时傅里叶变换（STFT），在时间和频率分辨率之间存在固有矛盾，对于快速变化的频率跟踪性能较差；小波变换虽然在一定程度上缓解了这一问题，但选择合适的小波基函数和分解层数较为困难，且计算复杂度较高。因此，寻找一种鲁棒且高效的时变频率估计方法具有重要的理论意义和实际应用价值。
卡尔曼滤波器作为一种优化的递归数据处理算法，能够有效地从含有噪声的测量数据中估计系统的状态，在处理时变系统和噪声方面展现出优越性。但传统的卡尔曼滤波器仅适用于线性系统和高斯噪声环境，对于非线性系统，其性能会显著下降。为了应对非线性问题，扩展卡尔曼滤波器和无迹卡尔曼滤波器应运而生，并在窄带信号时变频率估计中得到了广泛应用。本文将深入探讨基于EKF和UKF的窄带信号时变频率估计方法，对比分析两者的性能特点，并通过实验验证其有效性。
二、窄带信号模型与状态空间描述
2.1 窄带信号模型
假设窄带信号可以表示为如下形式：
其中，T为采样间隔。
2.2 状态空间模型
为了利用卡尔曼滤波器进行估计，需要建立系统的状态空间模型。将瞬时相位θ(n)和频率f(n)作为状态变量，定义状态向量为X(n)=[θ(n),f(n)]T。
状态方程描述了状态变量随时间的变化规律，假设过程噪声反映了模型的不确定性和频率变化的随机性，状态方程可以表示为：
测量方程则描述了状态变量与观测信号之间的关系：
至此，建立了窄带信号时变频率估计的状态空间模型，该模型包含一个线性状态方程和一个非线性测量方程。
三、扩展卡尔曼滤波器（EKF）原理及应用
3.1 EKF原理
扩展卡尔曼滤波器通过对非线性测量方程进行一阶泰勒展开，将其线性化，然后应用标准卡尔曼滤波器的框架进行状态估计。其基本步骤包括预测步骤和更新步骤。
预测步骤：
3.2 EKF在窄带信号时变频率估计中的应用
将上述EKF算法应用于窄带信号时变频率估计的状态空间模型中。在预测步骤中，根据状态方程预测下一时刻的状态和协方差；在更新步骤中，利用观测信号和预测值计算卡尔曼增益，修正状态预测值和协方差估计值。通过不断地递归预测和更新，实现对窄带信号时变频率的实时估计。
EKF的核心在于对非线性测量方程进行线性化，但这种线性化方法仅在一阶泰勒展开点附近有效。当非线性程度较高时，线性化误差会显著影响滤波器的性能，甚至导致滤波器发散。此外，EKF需要计算非线性函数的雅可比矩阵，这在某些复杂的系统中可能需要进行符号运算或数值求导，增加了算法的实现难度和计算复杂度。
四、无迹卡尔曼滤波器（UKF）原理及应用
4.1 UKF原理
无迹卡尔曼滤波器采用了一种不同的策略来处理非线性问题。它不是直接对非线性函数进行线性化，而是通过一组称为Sigma点的样本点来近似状态的概率分布。这些Sigma点经过非线性函数变换后，其统计特性（均值和协方差）可以更准确地反映状态的真实分布。UKF的基本步骤如下：
生成Sigma点：
计算权重：
预测步骤：
更新步骤：
4.2 UKF在窄带信号时变频率估计中的应用
将UKF算法应用于窄带信号时变频率估计的状态空间模型中。首先生成Sigma点，然后通过状态方程和测量方程对Sigma点进行传播和变换，计算预测的均值、协方差、测量协方差和互协方差，进而得到卡尔曼增益，最后更新状态估计值和协方差估计值。通过不断地递归计算，实现对窄带信号时变频率的实时估计。
UKF通过无迹变换逼近非线性变换对分布的影响，避免了EKF中对非线性函数进行线性化带来的误差，理论上精度更高，鲁棒性更好，尤其适用于非线性程度较高的情况。此外，UKF不需要计算非线性函数的雅可比矩阵，避免了复杂的求导过程，在实现上相对更便捷。然而，UKF需要生成和传播多个Sigma点，计算复杂度相对较高。
五、实验结果与分析
5.1 实验设置
为了评估EKF和UKF在窄带信号时变频率估计中的性能，进行了一系列实验。实验中生成具有时变频率的窄带信号，信号幅度A=1，采样间隔T=0.01s，信号频率随时间呈阶跃变化。过程噪声协方差矩阵Q和观测噪声协方差矩阵R根据实际情况进行设置。分别使用EKF和UKF对信号的时变频率进行估计，并记录估计结果。
5.2 性能评估指标
为了定量评估EKF和UKF的估计性能，引入性能指数（PI）和鲁棒性指数（RI）作为评估指标。性能指数反映了估计值与真实值之间的平均误差，鲁棒性指数反映了算法在不同噪声环境下的稳定性。
5.3 实验结果
实验结果表明，在非线性程度较低的情况下，EKF和UKF都能较好地估计窄带信号的时变频率，估计精度相差不大。但随着非线性程度的增加，UKF的估计精度明显高于EKF，且鲁棒性更好。例如，当信号频率变化较快或信噪比较低时，UKF能够更准确地跟踪频率的变化，而EKF的估计误差较大，甚至可能出现发散现象。
此外，从计算复杂度方面来看，EKF的计算量相对较小，因为其不需要生成和传播多个Sigma点，也不需要计算复杂的权重。而UKF由于需要处理多个Sigma点，计算复杂度相对较高，但随着计算机性能的不断提高，计算复杂度的问题在一定程度上可以得到缓解。
六、结论
本文深入探讨了基于扩展卡尔曼滤波器和无迹卡尔曼滤波器的窄带信号时变频率估计方法。通过建立窄带信号的非线性状态空间模型，分别运用EKF和UKF进行递归估计，并对比分析了两者在估计精度、鲁棒性及计算复杂度等方面的性能特点。实验结果表明，UKF在非线性程度较高的情况下具有更高的估计精度和鲁棒性，能够更准确地跟踪窄带信号的时变频率；而EKF计算量相对较小，适用于非线性程度较低的场景。
在实际应用中，应根据具体的需求和场景选择合适的滤波器。如果对估计精度和鲁棒性要求较高，且计算资源允许，UKF是更好的选择；如果对实时性要求较高且非线性程度可控，EKF可能更合适。未来的研究可以进一步探索如何自适应地调整过程噪声和观测噪声的协方差矩阵，以提高滤波器在未知噪声环境下的性能，以及将频率估计与信号幅度、相位等其他参数的估计相结合，实现对窄带信号的全面状态跟踪。
📚2 运行结果
%% RI
clear all
close all
clc
fclose('all');
global logpath;
addpath('./ukf')
addpath('./ekf')
addpath('./generate')
addpath('./ri')
PLOT_PREDICTION = false;
PLOT_PROFILE = false;
COMPUTE = true;
% Parameters
initial_omega = pi/2;
step_profile = [pi/11 -pi/5.5 pi/4 -pi/3.4 pi/2.8 -pi/2.5 pi/2 -pi/1.6 pi/1.4]; % Approximation of Fig.2 on the paper
step_length = 400;
sigma = 1e+1;
q = 1e-5;
r = sigma*q;
n_simulations = 500;
initialization_noise_sigma = 0.001;
threshold = 120;
% sigma_noise = 2e-2;
% sigma_omega_noise = 1e-3;
sigma_noise = 4e-1;
sigma_omega_noise = 3e-1;
if (PLOT_PREDICTION)
ri_figure = figure('visible','off');
set(ri_figure, 'Position', [0, 0, 720, 720],'Resize','off')
end
namestring = sprintf('q%1.2e_s%1.2e_sn%1.2e_sno%1.2e',q,sigma,sigma_noise,sigma_omega_noise);
logpath = strcat(pwd,'\data\ri\',namestring,'\')
if ~exist(logpath,'dir'), mkdir(logpath), end
if (COMPUTE)
ri_ekf = zeros(1,n_simulations);
ri_ukf = ri_ekf;
for ii = 1:n_simulations
fprintf('***** Iteration %d *****\n',ii)
% Generate signal
[signal, omega]=generate_signal_step(step_length,initial_omega,step_profile,sigma_noise,sigma_omega_noise);
%% Track
x_pred_0 = [0, 0, normrnd(omega(1),omega(1)*initialization_noise_sigma)]; % Initialization
pred_vec_ekf = ekf( ...
signal, ...%signal
x_pred_0, ...%x_pred_0
initialization_noise_sigma, ...%sigma_init
r, ... %r,
q ... %q
);
pred_vec_ukf = ukf( ...
signal, ...%signal
x_pred_0, ...%x_pred_0
initialization_noise_sigma, ...%sigma_init
r, ... %r,
q ... %q
);
% Take only the state we are interested in
pred_omega_ekf = pred_vec_ekf(3,:);
pred_omega_ukf = pred_vec_ukf(3,:);
disp('************ EKF *************');
ri_ekf(ii) = compute_ri(pred_omega_ekf,omega,step_length,threshold);
disp('************ UKF *************');
ri_ukf(ii) = compute_ri(pred_omega_ukf,omega,step_length,threshold);
% Plot ground truth and predictions
if PLOT_PREDICTION
set(0, 'currentfigure', ri_figure);
clf;
title('Predictions for RI');
t = 1:length(omega);
stairs(t, omega, 'k');
xlabel('Samples')
ylabel('Omega')
set(gca,'YLim',[0,pi]);
set(gca,'YTick',0:pi/4:pi);
set(gca,'YTickLabel',{'0' 'pi/4','pi/2','3/4pi','pi'});
hold on
plot(t,pred_omega_ekf,'ro');
plot(t,pred_omega_ukf,'bx');
legend('True','EKF','UKF')
RI(ii) = getframe;
end
end
if (PLOT_PREDICTION)
generate_video(strcat('ri_',namestring),RI);
end
save(strcat(logpath,'ri_ekf.mat'),'ri_ekf')
save(strcat(logpath,'ri_ukf.mat'),'ri_ukf')
else
if exist('ri_ekf.mat','file') && exist('ri_ukf.mat','file')
load(strcat(logpath,'ri_ekf.mat'),'ri_ekf')
load(strcat(logpath,'ri_ukf.mat'),'ri_ukf')
else
error('No saved curves. Recompute.')
end
end
n_steps = length(step_profile);
steps_axis = 1:n_steps;
psi_ekf = hist(ri_ekf, steps_axis);
psi_ukf = hist(ri_ukf, steps_axis);
if PLOT_PROFILE
figure(1)
hold off
stairs(omega)
xlabel('Samples')
ylabel('Omega')
set(gca,'YLim',[initial_omega-pi/2,initial_omega+pi/2]);
set(gca,'YTick',0:pi/4:pi);
set(gca,'YTickLabel',{'0' 'pi/4','pi/2','3/4pi','pi'});
end
ri_figure = figure(2);
title('RI')
bar(steps_axis,[psi_ekf', psi_ukf'])
legend('EKF','UKF')
saveas(ri_figure, strcat(logpath,'ri_figure_',namestring,'.png'), 'png');
🎉3参考文献
文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。(文章内容仅供参考，具体效果以运行结果为准)
🌈4Matlab代码、文章下载
资料获取，更多粉丝福利，MATLAB|Simulink|Python资源获取
