阵，根据矩阵求逆定理，矩阵A的逆 矩阵为A-1 = B-BC [ D+CTBC ]-1CtB ·
[0071]
[0072]
[0073]
[0074] =Pf(k)-Pf(k)S-(k)[I+Pf(k)S-(IOr 1Pf(Ic)
[0075] 卡尔曼平滑计算公式如（1-16)所示：
[0076]
(1-16)
[0077] 其中，Kg(k)为卡尔曼平滑增益，其计算公式如（1-17)所示：
[0078] Kg(k)=Pf(k)S-(k)[I+Pf(k)S-(k)]-1 (1-17)
[0079] 以上是传统的扩展卡尔曼平滑计算公式，而一种更为高效的平滑器为RTS平滑器， 其RTS平滑输出值X αμ十算公式如（I -18)所示：
[0080]
.(1-1.8.)
[0081] 其中，RTS 平滑增益 KgahPYkMk+UlOWl^k+Dmk+l)为前向卡尔 曼滤波器预测误差自相关矩阵，F(k+1，k)T为单位矩阵。
[0082] RTS平滑的误差自相关矩阵为：
[0083] P(k)=Pf(k)-Kg(k)[Kf(k,k+l)-P(k+l)]Kg(k) T (1-19)
[0084] 本发明实施例通过前向卡尔曼滤波获取待处理信号的前向卡尔曼滤波数据，根据 所述前向卡尔曼滤波数据和RTS平滑理论，平滑估计出待处理信号中的心电信号，从而实现 心电信号中静电干扰的去除，解决了现有扩展卡尔曼滤波器去除心电信号中静电干扰准确 率不高的问题。
[0085] 实施例二:
[0086] 图2示出了本发明实施例二提供的去除心电信号中静电干扰的方法的实现流程， 所述实现流程详述如下：
[0087] 在步骤S201中，获取待处理信号的相关参数；
[0088] 在本发明实施例中，所述待处理信号为含有静电干扰的心电信号。
[0089] 在本发明实施例中，所述待处理信号的相关参数为过程噪声矢量均值r和前向卡 尔曼滤波观测噪声自相关矩阵R(2，2)。
[0090] 进一步的，所述获取待处理信号的相关参数包括：
[0091] 采集N个心电周期的待处理信号，并存入2 XM矩阵0,矩阵0的第ko个列矢量
其中，kQ=l，2,…，M，M、N为大于零的整数，4。为待处理信号的幅值，k Q为采样点 序号；
[0092]对所述采集的N个心电周期的待处理信号进行R波检测，将所述待处理信号的R峰 的序号ko存入一维列矢量R，计算第i个心电周期的起始点位置
并将其存入一维列矢量Beg，其中，R(i)为第i个待处理信号R峰的采样点序号，INT为取整函 数；
[0093] 在本发明实施例中第i个心电周期信?
[0094] 在本发明实施例中，通过对待处理信号进行R波检测，获得待处理信号的R峰，如3 所示是心电信号R波检测的结果示例图。
[0095] 根据所述第i个心电周期的起始点位置，计算第i个心电周期的采样点数Num1,以 获得所述待处理信号的全局角速房
^其中，fs为采样频率；
[0096]在本发明实施例中，第i个心电周期的采样点数Numi = Beg(i+l)-Beg(i)。 [0097]将第i个心电周期的待处理信号映射成长度为fs的信号序列ECGnewi，并计算所述N
个心电周期的信号平均僅
[0098] 示例性的，若采杆观半Is = ?〇Linz，别付朱IT;L·、PB;可別tfJ付见牲?曰可吹別；3乂叹反 为350的信号序列ECGnewi。
[0099] 如图4所示是多个心电周期的平均ECGmean示例图。
[0100]在本发明实施例中，基于非线性插值法将第i个心电周期的待处理信号映射成长 度为fs的信号序列，具体的：先把第i个心电周期信号ECG1的持续时间均匀地划分匕等分，新 的心电周期中第k e个采样点的值ECGneWl(ke)可以根据其所处的区间，参考原采样点ECG 1 (η)和ECGKn+l)连线后的线段斜率计算得到，具体计算公式为
其中，η为第个新的心电信号所处的区间序号。
[0101] 示例性的，将含有128个采样点的心电周期的待处理信号映射成长度为150个采样 点，映射后的第一个采样点与第一个原采样点一样，区间序号为1;映射后的第二个采样点 处于第一个原采样点与第二个原采样点之间，区间序号为2,依次定义区间信号。
[0102] 根据所述N个心电周期的信号平均值ECGmeanUJ，计算所述待处理信号的模型参
，以获得过程噪声向量均值 C中η为建模不确定性；
[0103] 在本发明实施例中，基于MATLAB平台，利用非线件最小二乘函数Isanonlin函数， 通过ECGmean，对心电信号模型的参数
行非线性拟合
；中，函数func为 在参数为斤条件下，重构的心电信号的第i个数据点与ECGmean中对应的数据点的差值，设 走的参数[仪"，仪"，，仪'%.' K/. ·?"，广 ?-s·，，沒"，，沒"，沒、'V ]的初值为 [45_8，-54.0,429.0,-22.0,76.3，-2.4, -丨'7,0.丨4,0.09,0.67,2.44,2.0,3.丨7,3.32,5_08]。
[0104] 在本发明实施例中，建模不确定性η的均值为0,如图5所示是用计算出的模型参数 重构的心电信号示例图。
[0105] 根据信号序列ECGnewi(ke)和所述N个心电周期的信号平均值ECGmeanUe)，计算所 述待处理信号的观测噪声功率ECGms，以获得前向卡尔曼滤波观测噪声自相关矩阵
其中δ为采样间隔。
[0106] 在步骤S202中，获取所述待处理信号的前向卡尔曼滤波数据；
[0107] 在本发明实施例中，根据待处理信号的相关参数为过程噪声矢量均值承和前向卡 尔曼滤波观测噪声自相关矩阵R(2,2)，计算所述待处理信号的前向卡尔曼滤波数据。
[0108] 进一步的，所述获取待处理信号的前向卡尔曼滤波数据包括：
[0109] 采集所述待处理信号，计算所述待处理信号中第kP个采样点所对应的弧值免卜， 以获得心时刻的观测矩阵
，其中，\为匕时刻的信号幅值，kp为大于零的整数；
[0110] 在本发明实施例中，采集待处理信号，存入2 Xp'矩阵Y'，每一个待处理信号采样 点为
，其中，气为匕时刻的信号幅值，《为一个ECG周期数据点均匀地映射到0 ~2π的范围中，从而计算出第i个心电周期中所包含的采样点所对应的弧值；
[0111] 对所述待处理信号进行R波检测，将得到的R峰的序号ko存入一维列矢量RP，可以得 到矩阵Y '中的第i个心电周期SCGf的起始点位置，并将其存入一维列矢量Begp， Seg
其中，RP(i)为第i个心电R峰的采样点序号，INT为取整函 数；由此可得第i个心电周期中每一个采样数据点的值
ICGf中第kp个采样点的弧
，. NumPi = Begp (i+1) -Begp (i )，从第一个心电周期 的起始位置Begp(I)开始，至Ij最后一个心电周期的结束位置Begp(T) (T为列矢量Begp最后一 个元素的值)结束，将这T-I段的心电周期存入2Xp矩阵Y，矩阵Y的每一列为
矩 阵Y为扩展卡尔曼RTS平滑器的观测矩阵，其中，p = Begp (T)-Begp (1 )，kP = 1，2，…，p。
[0112] 计算kP时刻的前向卡尔曼滤波预测彳]
和kP时 刻的前向卡尔曼滤波预测误差自相关矩阵太％= ，其 中，f为前向卡尔曼滤波的状态转移抽象函数，4^1为kP-l时刻前向卡尔曼滤波的状态转移 抽象函数f对.-i)的偏导数，为kP-l时刻前向卡尔曼滤波的状态转移抽象函数f 对妒的偏导数，Pf (kp-l)为kP-l时刻的前向卡尔曼滤波误差自相关矩阵，Q为前向卡尔曼滤 波观测误差自相关矩阵；
[0113] 根据kP时刻的观测矩阵^，计算kP时刻的新息过程￡\ >以获得前向卡尔曼滤波值
其中，IiijlSkp时刻的前向卡尔曼滤波增益矩阵；
[0114] 根据kP时刻的前向卡尔曼滤波增益矩阵％/p和kP时刻的前向卡尔曼滤波预测误差 自相关矩阵K f( k p - I，k P )，计算k P时刻的前向卡尔曼滤波误差自相关矩阵
，其中，I和C为单位矩阵。
[0115] 在本发明实施例中，定义前向卡尔曼滤波的状态矢量#(〇)= U，前向卡尔曼滤波 误差自相关矩阵的初
以;^iodPPf(O)为输入进行1时刻前向扩展卡尔曼滤波， 保存前向卡尔曼滤波的状态矢量预测值X/丨Oj、前向卡尔曼滤波值^ iij、前向卡尔曼滤波误 差自相关矩阵Pf(I)和前向卡尔曼预测误差自相关矩阵Kf(0，l)，根据扩展卡尔曼滤波理论， 由初始状态计算1时刻前向卡尔曼滤波的具体过程分为以下几个步骤：
[0116] (1)由状态矢量的初值χ/抑和平均参数矢量P，计算1时刻状态矢量的预测值 对⑴叭，其中，前向卡尔曼滤波的状态转移抽象函数为
知+1.为kP+l时刻弧值Θ的预测值，Zy为1^+1时刻心电幅值Z的预测值，Δ6 = ^ - + i e{P,Q,R,S,T};
为 kP+l 时刻 的过程噪声。
[0117] ( 2 )由初始滤波自相关矩阵P f ( 0 )计算1时刻的预测误差自相关矩阵 (0,丨）=4广（〇)4( +6?'，其中喪为状态转移方程组对状态矢量X的偏导数，令为 状态转移方程组对参数矢量W的偏导数，具体计算参见公式(1-13)。
[0118] (3)由1时刻状态矢量的预测值;计算1时刻的新息过程，其 中Yi为1时刻的观测矢量，前向卡尔曼滤波观测抽象函i
[0119