如何用MATLAB对CSI数据进行预处理（卡尔曼滤波篇）

目录

一、前言

二、数据采集

三、卡尔曼滤波算法

1、什么是卡尔曼算法

2、卡尔曼的大致实现过程

3、卡尔曼滤波算法的核心计算式

4、单维度数据处理 

5、多维度数据处理 

一、前言

由于WiFi信号强度受多径效应和噪声的影响会导致定位精度低和性能不稳定等的问题。与RSSI相比，信道状态信息（channel status information，CSI）能有效避免多径效应给定位结果带来的不良影响，因此，采用CSI的值作为定位的特征值，建立Radio Map的位置指纹数据库，通过加权邻近算法匹配k组最近的指纹库数据可以估计出待测点的位置。

二、数据采集

这个阶段不做过多说明，怎么样采集及读取采集的csi数据也是一个难点。基于信道状态信息（channel status information，CSI）的室内定位法，其采用了正交频分复用（Orthogonal Frequency Division Multiplexing，OFDM）技术将通信信道分为多个不同频率的独立子信道，在每个子信道物理层收集CSI值作为指纹测量值。本文利用了配备Intel5300NIC网卡的台式机进行实验，设置64个样本点，每个样本点采集1500个数包，每个数据包有2*3*30条数据（2个发送天线，3个接收天线，30个子载波）。下图为标准采样点（蓝色）和测试点（红色），指纹地图由于涉及到实验保密性，故不展示。

三、卡尔曼滤波算法

1、什么是卡尔曼算法

卡尔曼滤波(Kalman filtering)是一种利用线性系统状态方程，通过系统输入输出观测数据，对系统状态进行最优估计的算法。由于观测数据中包括系统中的噪声和干扰的影响，所以最优估计也可看作是滤波过程。卡尔曼滤波算法在估计值和观测值之间做了一个修正。

2、卡尔曼的大致实现过程

卡尔曼的实现过程：使用上一次的最优结果预测当前的值，同时使用观测值修正当前值，得到最优结果。

3、卡尔曼滤波算法的核心计算式

本文就是根据以下公式，利用MATLAB进行数据处理

4、单维度数据处理 

每个点位采集的CSI数据有2*3*30*1500个幅度信息，我们采用第一根发射天线第一根接收天线的1500个单维幅度值进行卡尔曼滤波，MATLAB代码如下：

close all
clc
r1=xlsread('1.xlsx',1,'A1:A1500' ); %读取采集的csi值
csi_length = length(r1);%获取数据的长度
Z=r1; %获取观测量 
A=1; %状态转移矩阵
X(1)=15; %目标状态向量
H=1; %观测矩阵
P(1)=100;   %第一步预测的状态协方差矩阵
R=10e-4; %观测噪声方差
Q=10e-4; %噪声协方差
%有些变量会随迭代次数发生改变，所以给这些变量预先分配内存
X_=zeros(1,csi_length);
P_=zeros(1,csi_length);
K=zeros(1,csi_length);
 for t=2:csi_length
   %预测
    X_(t)=A*X(t-1);   %状态方程
    P_(t)=A*P(t-1)*A'+Q; %预测的协方差
    %更新校准
    K(t)=P_(t)*H'/(H*P_(t)*H'+R); %卡尔曼增益
    X(t)=X_(t)+K(t)*(Z(t)-H*X_(t)); %状态更新方程
    P(t)=(1-K(t)*H)*P_(t); %协方差更新方程
end
%绘图
x1=1:csi_length;
hold on
plot(x1,Z,'b-',x1,X,'r-');    % 红色线最优化估算结果滤波后的值，%蓝色线表示观测值
legend('观测值','滤波后的值','Location','northwest');
xlabel('子载波(f)');ylabel('幅值(dB)');
title('卡尔曼滤波');
set(gca,'Ylim',[5,25]);
hold off

卡尔曼滤波效果如下图： 

5、多维度数据处理 

其实跟上面类似，只不过现在读取30列*5200行数据而已（换了个数据集)

r2=xlsread('2.xlsx',1,'A1:AD5200' );
csi_length = length(r2);%获取数据的长度
Z=r2; %获取观测量 
A=1; %状态转移矩阵
X(1)=27; %目标状态向量
H=1; %观测矩阵
P(1)=100;   %第一步预测的状态协方差矩阵
R=10e-4; %观测噪声方差
Q=10e-4; %噪声协方差
for t=2:csi_length
   %预测
    X_(t)=A*X(t-1);   %状态方程
    P_(t)=A*P(t-1)*A'+Q; %预测的协方差
    %更新校准
    K(t)=P_(t)*H'/(H*P_(t)*H'+R); %卡尔曼增益
    X(t)=X_(t)+K(t)*(Z(t)-H*X_(t)); %状态更新方程
    P(t)=(1-K(t)*H)*P_(t); %协方差更新方程
end
%绘图
figure('Position', [100 100 1000 500]);
x1=1:csi_length;
hold on
subplot(1,2,1);    plot(x1,Z,'b-');    % 红色线最优化估算结果滤波后的值，%蓝色线表示观测值
legend('观测值','Location','northwest');
xlabel('总子载波数(个)');ylabel('幅值(dB)');
title('原始数据图');
subplot(1,2,2);    plot(x1,X,'r-');
legend('Kalman滤波后的值','Location','northwest');
xlabel('总子载波数(个)');ylabel('幅值(dB)');
title('使用卡尔曼滤波算法之后的数据图');
hold off

 卡尔曼滤波效果如下图： 

注：上述多维度数据处理还是存在一些问题的（滤波前后维度不一致），仅供参考，仅供参考！！！，后续将慢慢优化 。其实我是想做成下图的，横坐标表示第一根发射天线第一根接收天线的子载波数（30个），纵坐标表示每个子载波的幅度值，每条线表示不同的数据包（这里是1500个数据包）