[发明专利]基于线性预测的膈肌肌电实时处理方法及存储介质

权利要求书

1.一种基于线性预测的膈肌肌电实时处理方法，其特征在于：

设置线性预测模型，该线性预测模型作为利用历史时刻心电干扰的测量值来预测当前及未来时刻心电干扰值的膈肌肌电滤波器；

采用膈肌肌电滤波器对被心电干扰的膈肌肌电信号进行卷积运算，对卷积后信号作超阈值置零处理，得到信号y(k)；

求取含心电的膈肌肌电片段，计算该片段的线性预测模型系数，利用该模型系数自适应调整膈肌肌电滤波器系数；

对信号y(k)依次进行二阶高通滤波和M_conv阶低通滤波，得到降噪后膈肌肌电信号z(k)；对降噪后膈肌肌电信号z(k)进行包络计算，得到降噪后膈肌肌电信号的包络信号z_e(k)，实现实时处理膈肌肌电信号；其中，k为当前时刻。

2.根据权利要求1所述的基于线性预测的膈肌肌电实时处理方法，其特征在于：包括以下步骤：

S1、参数初始化：设定时间序号k＝0，心电峰值序号j＝0，初始化心电周期t_ECG(0)、心电片段长度t_p(0)和膈肌肌电滤波器系数w_i(0)，i＝0，...，M_conv-1；

其中，M_conv为滤波器阶数；

S2、确定心电峰值x_P及其搜索方式，设定置零阈值初始值Th(0)：

以采样频率f_s采集被心电干扰的膈肌肌电信号x(k)，搜索连续N秒内x(k)的最大值和最小值，其中最大值和最大值时刻分别记为x_max和k_max，最小值和最小值时刻分别记为x_min和k_min；若x_max＞|x_min|，则令心电峰值x_P＝x_max，且设置心电峰值搜索方式标志flag＝1；否则令心电峰值x_P＝|x_min|，且设置心电峰值搜索方式标志flag＝0；N的取值应保证期间至少包括3个心电周期；

设定置零阈值的初始值为Th(0)∈[0.3x_P，0.7X_P]；

S3、以采样频率f_s采集被心电干扰的膈肌肌电信号x(k)；

S4、若k＜M_conv-1，则k＝k+1，并跳转S3；否则，采用膈肌肌电滤波器w_i(j)，i＝0，...，M_conv-1对信号x(k)作卷积运算，卷积后的信号记为计算公式为：

其中K_conv为卷积放大倍数；膈肌肌电滤波器阶数M_conv≥2，卷积放大倍数K_conv∈[1，30]；

S5、阈值处理：对进行超阈值置零处理，所得信号记为y(k)，并自适应更新置零阈值；

对进行超阈值置零处理，所得信号记为y(k)，计算公式为：

自适应更新置零阈值的计算公式为：

Th(j+1)＝α_Th·Th(j)+(1-α_Th)·y(k)

其中，α_Th为阈值调节因子；阈值调节因子α_Th∈[0.9，0.98]；

S6、定位心电峰值时刻：

若flag＝1，则判断当前时刻k是否满足条件x(k-1)≥x(k-2)且x(k-1)≥x(k)且|x(k-1)-x_P|＜δ；

若flag＝0，则判断当前时刻k是否满足条件x(k-1)≤x(k-2)且x(k-1)≤x(k)且|x(k-1)+x_P|＜δ；

若满足则：

记录并保存k_p(j)＝k-1；

若j＞0，则更新心电周期t_ECG(j)与心电片段长度t_p(j)，计算公式为：

t_p2p＝k_p(j)-k_p(j-1)

t_ECG(j)＝α_ECG·t_ECG(j-1)+(1-α_ECG)·t_p2p

其中，为四舍五入取整运算，α_ECG为心电周期调节因子，α_p为心电片段在整个心电周期中的时间占比；

更新心电峰值序号j+1→j；

δ∈[0.01x_P，0.2x_P]，α_ECG∈[0.9，0.98]，α_p∈[0.5，0.7]；

S7、判断当前时刻是否满足若满足则：

记含心电的膈肌肌电片段为计算s(j)的线性预测模型系数

若j＞0，则更新膈肌肌电滤波器系数w_i(j)，计算公式为：

其中，α_w为滤波器系数调节因子；

线性预测模型系数的求解采用自相关法实现；

S8、对y(k)依次进行截止频率为f_HPc的二阶高通滤波和截止频率为f_LPc的M_conv阶低通滤波，滤波结果记为z(k)；

S9、对z(k)进行包络计算，得到包络信号z_e(k)：

对z(k)的绝对值进行对数运算，运算结果记为v(k)，计算公式为：

v(k)＝ln(|z(k)|+ξ)

对v(k)作截止频率为f_LPe的M₂阶低通滤波，滤波结果记为o(k)；

对o(k)进行指数运算，运算结果即为包络信号z_e(k)，计算公式为：

z_e(k)＝e^o(k)

其中，ξ为一很小的正数，ξ＝0.001；M₂≥2；

S10、分别输出S8带通滤波后的信号z(k)(即降噪后膈肌肌电信号)及S9的包络信号z_e(k)(即降噪后膈肌肌电信号的包络)；更新时间序号k+1→k，并跳转至S3。