CN110916636B - 一种基于动态二阶差分阈值的bcg信号心率计算方法及系统 - Google Patents

一种基于动态二阶差分阈值的bcg信号心率计算方法及系统 Download PDF

Info

Publication number
CN110916636B
CN110916636B CN201911156128.6A CN201911156128A CN110916636B CN 110916636 B CN110916636 B CN 110916636B CN 201911156128 A CN201911156128 A CN 201911156128A CN 110916636 B CN110916636 B CN 110916636B
Authority
CN
China
Prior art keywords
wave
new
amplitude
value
threshold
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201911156128.6A
Other languages
English (en)
Other versions
CN110916636A (zh
Inventor
李润超
魏春雨
宿天赋
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ennova Health Technology Co ltd
Original Assignee
Ennova Health Technology Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Ennova Health Technology Co ltd filed Critical Ennova Health Technology Co ltd
Priority to CN201911156128.6A priority Critical patent/CN110916636B/zh
Publication of CN110916636A publication Critical patent/CN110916636A/zh
Application granted granted Critical
Publication of CN110916636B publication Critical patent/CN110916636B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/02—Detecting, measuring or recording pulse, heart rate, blood pressure or blood flow; Combined pulse/heart-rate/blood pressure determination; Evaluating a cardiovascular condition not otherwise provided for, e.g. using combinations of techniques provided for in this group with electrocardiography or electroauscultation; Heart catheters for measuring blood pressure
    • A61B5/024—Detecting, measuring or recording pulse rate or heart rate
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/103—Detecting, measuring or recording devices for testing the shape, pattern, colour, size or movement of the body or parts thereof, for diagnostic purposes
    • A61B5/11—Measuring movement of the entire body or parts thereof, e.g. head or hand tremor, mobility of a limb
    • A61B5/1102—Ballistocardiography
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A61B5/725—Details of waveform analysis using specific filters therefor, e.g. Kalman or adaptive filters

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Public Health (AREA)
  • Molecular Biology (AREA)
  • Veterinary Medicine (AREA)
  • General Health & Medical Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Animal Behavior & Ethology (AREA)
  • Biophysics (AREA)
  • Pathology (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Physiology (AREA)
  • Surgery (AREA)
  • Signal Processing (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Artificial Intelligence (AREA)
  • Psychiatry (AREA)
  • Cardiology (AREA)
  • Dentistry (AREA)
  • Oral & Maxillofacial Surgery (AREA)
  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)

Abstract

本发明公开了一种基于动态二阶差分阈值的BCG信号心率计算方法及系统，所述方法包括：采集获得BCG信号；根据预设规则对BCG信号进行滤波预处理，获得滤波后的BCG信号；根据动态二阶差分阈值算法处理滤波后的BCG信号，确定每个信号周期中的J波波峰位置，获得J波波峰序列；根据所述J波波峰序列中的每个J波波峰位置，计算BCG信号主波的平均长度；根据所述平均长度以及采样率计算获得心率；所述方法及系统的算法适应性强，在不同信号上都具有较好的去噪效果；所述方法及系统参数结构简单且计算复杂度低，解决了现有BCG信号提取方法计算时间长、时间延迟大的问题，具有较高的计算精度且适用于实时监测心率。

Description

一种基于动态二阶差分阈值的BCG信号心率计算方法及系统
技术领域
本发明涉及医疗技术领域，更具体地，涉及一种基于动态二阶差分阈值的BCG信号心率计算方法及系统。
背景技术
心冲击信号(Ballistocardiogram,BCG)是一种记录心脏泵血引起身体震动的方法，是心脏监测方法之一，主要由血液循环过程中人体重力变化引起。相对于现有的心血管检测手段，心冲击图信号检测具有无创性、无接触式、检测方便等优势。BCG信号是一种人体生物信号，具有频率低、强度低的特性，而且易受到呼吸、体动、工频噪声的干扰，导致直接测量得到的BCG信号含有太多噪声，不能直接获得准确的生理特征信息，所以需要对采集到的BCG信号做滤波处理，去除噪声干扰。
目前，对于BCG信号的特征分析方法较复杂，主要结合脉搏波信号来提取BCG信号的特征。通过获取BCG信号J波位置检测心率，目前常用的J波检测法有阈值法、伪周期法、自适应模版匹配法等，这些方法使用波形时窗分析，计算时间长、时间延迟大，不利于实时监测。
发明内容
为了解决背景技术存在的现有的对于BCG信号的提取方法计算时间长、时间延迟大、不适于实时监测心率的问题，本发明提供了一种基于动态二阶差分阈值的BCG信号心率计算方法及系统，所述方法及系统通过低通滤波以及经验模态分解算法优化采集信号，通过动态二阶差分阈值算法计算获得精确的J波波峰位置，进而获得实时准确的心率计算结果；所述一种基于动态二阶差分阈值的BCG信号心率计算方法包括：
采集获得BCG信号；
根据预设规则对BCG信号进行滤波预处理，获得滤波后的BCG信号；
根据动态二阶差分阈值算法处理滤波后的BCG信号，确定每个信号周期中的J波波峰位置，获得J波波峰序列；
根据所述J波波峰序列中的每个J波波峰位置，计算BCG信号主波的平均长度；根据所述平均长度以及采样率计算获得心率。
进一步的，所述根据预设规则对BCG信号进行滤波预处理，获得滤波后的BCG信号，包括：
通过截止频率为30Hz的低通滤波器滤除工频干扰及高频噪声；
根据经验模态分解算法，通过预设的多个阈值对原始的BCG信号进行经验模态分解，获得多阶IMF函数；
将含有噪声成分的IMF函数去除，对剩余的含有信号成分的IMF函数进行重构，获得滤波后的BCG信号。
进一步的，在所述通过截止频率为30Hz的低通滤波器滤除工频干扰及高频噪声前，所述方法还包括：
通过抖动标记方法，将由于体动而导致的振幅过大的区域去除；
通过滑动窗口对BCG信号进行平滑处理，去除异常值。
进一步的，根据动态二阶差分阈值算法处理滤波后的BCG信号，包括：
对滤波后的BCG信号进行离散化处理，获得K个信号点；按照预设的周期长度将所述K个信号点等分为N个周期，使每一段中包含一个BCG周期信号；
根据预设规则计算获得初始阈值；所述初始阈值包括初始差分阈值th1、初始振幅阈值上限th2以及初始振幅阈值下限th3；
根据初始阈值、预设的一阶差分条件以及二阶差分条件确定第一个周期的J波峰值点；
根据所述第一个周期的J波峰值点对初始阈值进行更新，获得新的阈值，即新的差分阈值、振幅阈值上限以及振幅阈值下限；并根据所述新的阈值以及预设的一阶差分条件以及二阶差分条件确定下一个周期的J波峰值点，直至获得所有J波峰值点。
进一步的，所述根据预设规则计算获得初始阈值，包括：
对每个周期中所包括的信号点进行差分处理，获得每个周期对应的差分最大值；比较所述N个周期中每个周期对应的差分最大值，去掉其中的最大值及最小值，并对剩余的差分最大值求算术平均值，获得平均差分最大值M0；取初始差分阈值th1＝a*M0，其中0.5≤a≤1；
获得每个周期对应的振幅最大值，比较所述N个周期中每个周期对应的振幅最大值，去掉其中的最大值及最小值，并对剩余的振幅最大值求算术平均值，获得平均振幅最大值H0；取初始振幅阈值下限为th2＝b*H0，初始振幅阈值上限为th3＝c*H0，其中0≤b≤1；1≤c≤2。
进一步的，根据初始阈值、预设的一阶差分条件以及二阶差分条件确定第一个周期的J波峰值点，包括：
在第一个周期内获得主波升支上的一个离散点yi，所述yi通过满足yi+1-yi>th1且yi-yi-1>th1的条件；
判断从yi开始后的相邻4个点为yk，yk+1，yk+2，yk+3，是否满足差分预设条件；若满足所述差分预设条件，则记录yk+2为所述第一个周期的J波峰值点，其幅度记为Hnew；所述差分预设条件为：yk+1-yk>0；yk+2-yk+1>0；yk+3-yk+2<0；且同时满足二阶差分条件：yn＝yk+1-yk；yn+1＝yk+2-yk+1；yn+2＝yk+3-yk+2；yn+1-yn<0且yn+2-yn+1<0；
若不满足所述差分预设条件，则重新获得主波升支上的新的离散点，直至获得满足所述差分预设条件的情形，获得J波峰值点。
进一步的，根据所述第一个周期的J波峰值点对初始阈值进行更新，包括：
获得从yk+1到yk+2这一过程的差分最大值，记为Mnew；使用所述最大差分值Mnew更换计算初始差分阈值中该周期对应的差分最大值，并重新计算新的平均差分最大值M0；根据新的M0计算获得新的差分阈值；
将所述J波峰值点yk+2对应的振幅Hnew更换计算初始振幅阈值中对应该周期的振幅最大值，并重新计算新的平均振幅最大值H0；根据新的H0计算获得新的振幅阈值上限以及振幅阈值下限。
进一步的，在确定每个信号周期中的J波波峰位置后，所述方法还包括：
确定每个J波波峰的前一波谷的幅值是否为该波峰所在周期内的最低值，将前一波谷幅值非最低值对应的J波波峰删除，获得新的J波波峰序列；
取所述新的J波波峰序列中的每一个J波波峰的点的左右各P个点进行局部极值检测；将局部极值检测出的振幅最高值对应的点替换该J波波峰，作为该周期内寻优后的J波波峰；
将新的J波波峰序列中所有寻优后的J波波峰作为最终的J波波峰序列。
进一步的，所述心率的计算公式为：
Pulse＝60/(Fs*L)
其中，Fs为采样率，L为主波的平均长度。
所述一种基于动态二阶差分阈值的BCG信号心率计算系统包括：
BCG信号采集单元，所述BCG信号采集单元用于采集获得BCG信号；
滤波预处理单元，所述滤波预处理单元用于根据预设规则对BCG信号进行滤波预处理，获得滤波后的BCG信号；
J波波峰计算单元，所述J波波峰计算单元根据动态二阶差分阈值算法处理滤波后的BCG信号，确定每个信号周期中的J波波峰位置，获得J波波峰序列；
心率计算单元，所述心率计算单元用于根据所述J波波峰序列中的每个J波波峰位置，计算BCG信号主波的平均长度；所述心率计算单元用于根据所述平均长度以及采样率计算获得心率。
进一步的，所述滤波预处理单元用于通过截止频率为30Hz的低通滤波器滤除工频干扰及高频噪声；
所述滤波预处理单元用于根据经验模态分解算法，通过预设的多个阈值对原始的BCG信号进行经验模态分解，获得多阶IMF函数；
所述滤波预处理单元用于将含有噪声成分的IMF函数去除，对剩余的含有信号成分的IMF函数进行重构，获得滤波后的BCG信号。
进一步的，所述滤波预处理单元用于通过抖动标记方法，将由于体动而导致的振幅过大的区域去除；
所述滤波预处理单元用于通过滑动窗口对BCG信号进行平滑处理，去除异常值。
进一步的，所述J波波峰计算单元用于对滤波后的BCG信号进行离散化处理，获得K个信号点；按照预设的周期长度将所述K个信号点等分为N个周期，使每一段中包含一个BCG周期信号；
所述J波波峰计算单元用于根据预设规则计算获得初始阈值；所述初始阈值包括初始差分阈值th1、初始振幅阈值上限th2以及初始振幅阈值下限th3；
所述J波波峰计算单元用于根据初始阈值、预设的一阶差分条件以及二阶差分条件确定第一个周期的J波峰值点；
所述J波波峰计算单元用于根据所述第一个周期的J波峰值点对初始阈值进行更新，获得新的阈值，即新的差分阈值、振幅阈值上限以及振幅阈值下限；并根据所述新的阈值以及预设的一阶差分条件以及二阶差分条件确定下一个周期的J波峰值点，直至获得所有J波峰值点。
进一步的，所述J波波峰计算单元用于对每个周期中所包括的信号点进行差分处理，获得每个周期对应的差分最大值；比较所述N个周期中每个周期对应的差分最大值，去掉其中的最大值及最小值，并对剩余的差分最大值求算术平均值，获得平均差分最大值M0；取初始差分阈值th1＝a*M0，其中0.5≤a≤1；
所述J波波峰计算单元用于获得每个周期对应的振幅最大值，比较所述N个周期中每个周期对应的振幅最大值，去掉其中的最大值及最小值，并对剩余的振幅最大值求算术平均值，获得平均振幅最大值H0；取初始振幅阈值下限为th2＝b*H0，初始振幅阈值上限为th3＝c*H0，其中0≤b≤1；1≤c≤2。
进一步的，所述J波波峰计算单元用于在第一个周期内获得主波升支上的一个离散点yi，所述yi通过满足yi+1-yi>th1且yi-yi-1>th1的条件；
所述J波波峰计算单元用于判断从yi开始后的相邻4个点为yk，yk+1，yk+2，yk+3，是否满足差分预设条件；若满足所述差分预设条件，则记录yk+2为所述第一个周期的J波峰值点，其幅度记为Hnew；所述差分预设条件为：yk+1-yk>0；yk+2-yk+1>0；yk+3-yk+2<0；且同时满足二阶差分条件：yn＝yk+1-yk；yn+1＝yk+2-yk+1；yn+2＝yk+3-yk+2；yn+1-yn<0且yn+2-yn+1<0；
若不满足所述差分预设条件，则重新获得主波升支上的新的离散点，直至获得满足所述差分预设条件的情形，获得J波峰值点。
进一步的，所述J波波峰计算单元用于获得从yk+1到yk+2这一过程的差分最大值，记为Mnew；使用所述最大差分值Mnew更换计算初始差分阈值中该周期对应的差分最大值，并重新计算新的平均差分最大值M0；根据新的M0计算获得新的差分阈值；
所述J波波峰计算单元用于将所述J波峰值点yk+2对应的振幅Hnew更换计算初始振幅阈值中对应该周期的振幅最大值，并重新计算新的平均振幅最大值H0；根据新的H0计算获得新的振幅阈值上限以及振幅阈值下限。
进一步的，所述系统还包括波峰寻优单元；
所述波峰寻优单元用于确定每个J波波峰的前一波谷的幅值是否为该波峰所在周期内的最低值，将前一波谷幅值非最低值对应的J波波峰删除，获得新的J波波峰序列；
所述波峰寻优单元用于取所述新的J波波峰序列中的每一个J波波峰的点的左右各P个点进行局部极值检测；将局部极值检测出的振幅最高值对应的点替换该J波波峰，作为该周期内寻优后的J波波峰；
所述波峰寻优单元用于将新的J波波峰序列中所有寻优后的J波波峰作为最终的J波波峰序列。
进一步的，所述心率计算单元计算心率的公式为：Pulse＝60/(Fs*L)
其中，Fs为采样率，L为主波的平均长度。
本发明的有益效果为：本发明的技术方案，给出了一种基于动态二阶差分阈值的BCG信号心率计算方法及系统；所述方法及系统通过多种滤波方法的结合后实现BCG信号的有效滤波；通过动态二阶差分阈值方法对J波进行检测，并基于J波的波峰计算心率，达到较高的计算精度。所述方法及系统的算法适应性强，在不同信号上都具有较好的去噪效果；所述方法及系统参数结构简单且计算复杂度低，解决了现有BCG信号提取方法计算时间长、时间延迟大的问题，具有较高的计算精度且适用于实时监测心率。
附图说明
通过参考下面的附图，可以更为完整地理解本发明的示例性实施方式：
图1为本发明具体实施方式的一种基于动态二阶差分阈值的BCG信号心率计算方法的流程图；
图2为本发明具体实施方式的一种基于动态二阶差分阈值的BCG信号心率计算系统的结构图。
具体实施方式
现在参考附图介绍本发明的示例性实施方式，然而，本发明可以用许多不同的形式来实施，并且不局限于此处描述的实施例，提供这些实施例是为了详尽地且完全地公开本发明，并且向所属技术领域的技术人员充分传达本发明的范围。对于表示在附图中的示例性实施方式中的术语并不是对本发明的限定。在附图中，相同的单元/元件使用相同的附图标记。
除非另有说明，此处使用的术语(包括科技术语)对所属技术领域的技术人员具有通常的理解含义。另外，可以理解的是，以通常使用的词典限定的术语，应当被理解为与其相关领域的语境具有一致的含义，而不应该被理解为理想化的或过于正式的意义。
图1为本发明具体实施方式的一种基于动态二阶差分阈值的BCG信号心率计算方法的流程图；如图1所示，所述方法包括：
步骤110，采集获得BCG信号；
所述采集获得BCG信号是根据现有的BCG信号采集方法对BCG信号进行采集，获得BCG的模拟或数字信号；若获得模拟信号，再后面步骤将进行离散化的取样处理，转换为数字信号。
步骤120，根据预设规则对BCG信号进行滤波预处理，获得滤波后的BCG信号；
BCG信号是复杂的人体内部运动产生的低频微弱信号，主要频率范围在0.6～20Hz之间，对其产生干扰的噪声源很多，比如体动干扰、0.2～0.7Hz的呼吸干扰、30～300Hz范围内的肌电干扰，以及50Hz的工频干扰等；故需要多种信号滤波方式结合来获得更清晰平滑的BCG信号；
首先，通过抖动标记方法，将由于体动而导致的振幅过大的区域去除；并通过滑动窗口对BCG信号进行平滑处理，去除异常值。
其次，通过截止频率为30Hz的低通滤波器滤除工频干扰及高频噪声；
最后根据经验模态分解算法，通过预设的多个阈值对原始的BCG信号进行经验模态分解，获得多阶IMF函数；
将含有噪声成分的IMF函数去除，对剩余的含有信号成分的IMF函数进行重构，获得滤波后的BCG信号。
步骤130，根据动态二阶差分阈值算法处理滤波后的BCG信号，确定每个信号周期中的J波波峰位置，获得J波波峰序列；
步骤131，对滤波后的BCG信号进行离散化处理，获得K个信号点；按照预设的周期长度将所述K个信号点等分为N个周期，使每一段中包含一个BCG周期信号；假设有10000个BCG信号，可以以1000个点为时长将10000个BCG信号等分为10段；
步骤132，根据预设规则计算获得初始阈值；所述初始阈值包括初始差分阈值th1、初始振幅阈值上限th2以及初始振幅阈值下限th3；
进一步的，对每个周期中所包括的信号点进行差分处理，获得每个周期对应的差分最大值；比较所述N个周期中每个周期对应的差分最大值，去掉其中的最大值及最小值，并对剩余的差分最大值求算术平均值，获得平均差分最大值M0；取初始差分阈值th1＝a*M0，其中0.5≤a≤1；
获得每个周期对应的振幅最大值，比较所述N个周期中每个周期对应的振幅最大值，去掉其中的最大值及最小值，并对剩余的振幅最大值求算术平均值，获得平均振幅最大值H0；取初始振幅阈值下限为th2＝b*H0，初始振幅阈值上限为th3＝c*H0，其中0≤b≤1；1≤c≤2。
本实施例中，a＝0.75；b＝0.5；c＝1.5；
步骤133，根据初始阈值、预设的一阶差分条件以及二阶差分条件确定第一个周期的J波峰值点；
在第一个周期内获得主波升支上的一个离散点yi，所述yi通过满足以下条件：
yi+1-yi>th1；
yi-yi-1>th1；
判断从yi开始后的相邻4个点为yk，yk+1，yk+2，yk+3，是否满足差分预设条件；若满足所述差分预设条件，则记录yk+2为所述第一个周期的J波峰值点，其幅度记为Hnew；
所述差分预设条件为：
yk+1-yk>0；
yk+2-yk+1>0；
yk+3-yk+2<0；
且同时满足二阶差分条件：
yn＝yk+1-yk；
yn+1＝yk+2-yk+1；
yn+2＝yk+3-yk+2；
yn+1-yn<0且yn+2-yn+1<0；
若不满足所述差分预设条件，则重新获得主波升支上的新的离散点，直至获得满足所述差分预设条件的情形，获得J波峰值点。
步骤134，根据所述第一个周期的J波峰值点对初始阈值进行更新，获得新的阈值，即新的差分阈值、振幅阈值上限以及振幅阈值下限；
具体的，获得从yk+1到yk+2这一过程的差分最大值，记为Mnew；使用所述最大差分值Mnew更换计算初始差分阈值中该周期对应的差分最大值，并重新计算新的平均差分最大值M0；根据新的M0计算获得新的差分阈值；
将所述J波峰值点yk+2对应的振幅Hnew更换计算初始振幅阈值中对应该周期的振幅最大值，并重新计算新的平均振幅最大值H0；根据新的H0计算获得新的振幅阈值上限以及振幅阈值下限。
步骤135，根据所述新的阈值以及预设的一阶差分条件以及二阶差分条件确定下一个周期的J波峰值点，直至获得所有J波峰值点。
进一步的，在确定每个信号周期中的J波波峰位置后，所述方法还包括：
确定每个J波波峰的前一波谷的幅值是否为该波峰所在周期内的最低值，将前一波谷幅值非最低值对应的J波波峰删除，获得新的J波波峰序列；
取所述新的J波波峰序列中的每一个J波波峰的点的左右各P个点进行局部极值检测；将局部极值检测出的振幅最高值对应的点替换该J波波峰，作为该周期内寻优后的J波波峰；
将新的J波波峰序列中所有寻优后的J波波峰作为最终的J波波峰序列。
步骤140，根据所述J波波峰序列中的每个J波波峰位置，计算BCG信号主波的平均长度；根据所述平均长度以及采样率计算获得心率。
所述心率的计算公式为：
Pulse＝60/(Fs*L)
其中，Fs为采样率，L为主波的平均长度。
图2为本发明具体实施方式的一种基于动态二阶差分阈值的BCG信号心率计算系统的结构图。如图2所示，所述系统包括：
BCG信号采集单元210，所述BCG信号采集单元210用于采集获得BCG信号；
滤波预处理单元220，所述滤波预处理单元220用于根据预设规则对BCG信号进行滤波预处理，获得滤波后的BCG信号；
所述滤波预处理单元220用于通过抖动标记方法，将由于体动而导致的振幅过大的区域去除；
所述滤波预处理单元220用于通过滑动窗口对BCG信号进行平滑处理，去除异常值。
所述滤波预处理单元220用于通过截止频率为30Hz的低通滤波器滤除工频干扰及高频噪声；
所述滤波预处理单元220用于根据经验模态分解算法，通过预设的多个阈值对原始的BCG信号进行经验模态分解，获得多阶IMF函数；
所述滤波预处理单元220用于将含有噪声成分的IMF函数去除，对剩余的含有信号成分的IMF函数进行重构，获得滤波后的BCG信号。
J波波峰计算单元230，所述J波波峰计算单元230根据动态二阶差分阈值算法处理滤波后的BCG信号，确定每个信号周期中的J波波峰位置，获得J波波峰序列；
所述J波波峰计算单元230用于对滤波后的BCG信号进行离散化处理，获得K个信号点；按照预设的周期长度将所述K个信号点等分为N个周期，使每一段中包含一个BCG周期信号；
所述J波波峰计算单元230用于根据预设规则计算获得初始阈值；所述初始阈值包括初始差分阈值th1、初始振幅阈值上限th2以及初始振幅阈值下限th3；
所述J波波峰计算单元230用于对每个周期中所包括的信号点进行差分处理，获得每个周期对应的差分最大值；比较所述N个周期中每个周期对应的差分最大值，去掉其中的最大值及最小值，并对剩余的差分最大值求算术平均值，获得平均差分最大值M0；取初始差分阈值th1＝a*M0，其中0.5≤a≤1；
所述J波波峰计算单元230用于获得每个周期对应的振幅最大值，比较所述N个周期中每个周期对应的振幅最大值，去掉其中的最大值及最小值，并对剩余的振幅最大值求算术平均值，获得平均振幅最大值H0；取初始振幅阈值下限为th2＝b*H0，初始振幅阈值上限为th3＝c*H0，其中0≤b≤1；1≤c≤2。
所述J波波峰计算单元230用于根据初始阈值、预设的一阶差分条件以及二阶差分条件确定第一个周期的J波峰值点；
所述J波波峰计算单元230用于在第一个周期内获得主波升支上的一个离散点yi，所述yi通过满足yi+1-yi>th1且yi-yi-1>th1的条件；
所述J波波峰计算单元230用于判断从yi开始后的相邻4个点为yk，yk+1，yk+2，yk+3，是否满足差分预设条件；若满足所述差分预设条件，则记录yk+2为所述第一个周期的J波峰值点，其幅度记为Hnew；所述差分预设条件为：yk+1-yk>0；yk+2-yk+1>0；yk+3-yk+2<0；且同时满足二阶差分条件：yn＝yk+1-yk；yn+1＝yk+2-yk+1；yn+2＝yk+3-yk+2；yn+1-yn<0且yn+2-yn+1<0；
若不满足所述差分预设条件，则重新获得主波升支上的新的离散点，直至获得满足所述差分预设条件的情形，获得J波峰值点。
所述J波波峰计算单元230用于根据所述第一个周期的J波峰值点对初始阈值进行更新，获得新的阈值，即新的差分阈值、振幅阈值上限以及振幅阈值下限；并根据所述新的阈值以及预设的一阶差分条件以及二阶差分条件确定下一个周期的J波峰值点，直至获得所有J波峰值点。
所述J波波峰计算单元230用于获得从yk+1到yk+2这一过程的差分最大值，记为Mnew；使用所述最大差分值Mnew更换计算初始差分阈值中该周期对应的差分最大值，并重新计算新的平均差分最大值M0；根据新的M0计算获得新的差分阈值；
所述J波波峰计算单元230用于将所述J波峰值点yk+2对应的振幅Hnew更换计算初始振幅阈值中对应该周期的振幅最大值，并重新计算新的平均振幅最大值H0；根据新的H0计算获得新的振幅阈值上限以及振幅阈值下限。
心率计算单元240，所述心率计算单元240用于根据所述J波波峰序列中的每个J波波峰位置，计算BCG信号主波的平均长度；所述心率计算单元240用于根据所述平均长度以及采样率计算获得心率。
所述心率计算单元240计算心率的公式为：Pulse＝60/(Fs*L)
其中，Fs为采样率，L为主波的平均长度。
所述系统还包括波峰寻优单元；
所述波峰寻优单元用于确定每个J波波峰的前一波谷的幅值是否为该波峰所在周期内的最低值，将前一波谷幅值非最低值对应的J波波峰删除，获得新的J波波峰序列；
所述波峰寻优单元用于取所述新的J波波峰序列中的每一个J波波峰的点的左右各P个点进行局部极值检测；将局部极值检测出的振幅最高值对应的点替换该J波波峰，作为该周期内寻优后的J波波峰；
所述波峰寻优单元用于将新的J波波峰序列中所有寻优后的J波波峰作为最终的J波波峰序列。
在此处所提供的说明书中，说明了大量具体细节。然而，能够理解，本公开的实施例可以在没有这些具体细节的情况下实践。在一些实例中，并未详细示出公知的方法、结构和技术，以便不模糊对本说明书的理解。
本领域那些技术人员可以理解，可以对实施例中的设备中的模块进行自适应性地改变并且把它们设置在与该实施例不同的一个或多个设备中。可以把实施例中的模块或单元或组件组合成一个模块或单元或组件，以及此外可以把它们分成多个子模块或子单元或子组件。除了这样的特征和/或过程或者单元中的至少一些是相互排斥之外，可以采用任何组合对本说明书(包括伴随的权利要求、摘要和附图)中公开的所有特征以及如此公开的任何方法或者设备的所有过程或单元进行组合。除非另外明确陈述，本说明书(包括伴随的权利要求、摘要和附图)中公开的每个特征可以由提供相同、等同或相似目的的替代特征来代替。本说明书中涉及到的步骤编号仅用于区别各步骤，而并不用于限制各步骤之间的时间或逻辑的关系，除非文中有明确的限定，否则各个步骤之间的关系包括各种可能的情况。
此外，本领域的技术人员能够理解，尽管在此所述的一些实施例包括其它实施例中所包括的某些特征而不是其它特征，但是不同实施例的特征的组合意味着处于本公开的范围之内并且形成不同的实施例。例如，在权利要求书中所要求保护的实施例的任意之一都可以以任意的组合方式来使用。
本公开的各个部件实施例可以以硬件实现，或者以在一个或者多个处理器上运行的软件模块实现，或者以它们的组合实现。本公开还可以实现为用于执行这里所描述的方法的一部分或者全部的设备或者系统程序(例如，计算机程序和计算机程序产品)。这样的实现本公开的程序可以存储在计算机可读介质上，或者可以具有一个或者多个信号的形式。这样的信号可以从因特网网站上下载得到，或者在载体信号上提供，或者以任何其他形式提供。
应该注意的是上述实施例对本公开进行说明而不是对本公开进行限制，并且本领域技术人员在不脱离所附权利要求的范围的情况下可设计出替换实施例。单词“包含”不排除存在未列在权利要求中的元件或步骤。位于元件之前的单词“一”或“一个”不排除存在多个这样的元件。本公开可以借助于包括有若干不同元件的硬件以及借助于适当编程的计算机来实现。在列举了若干系统的单元权利要求中，这些系统中的若干个可以是通过同一个硬件项来具体体现。
以上所述仅是本公开的具体实施方式，应当指出的是，对于本领域的普通技术人员来说，在不脱离本公开精神的前提下，可以作出若干改进、修改、和变形，这些改进、修改、和变形都应视为落在本申请的保护范围内。

Claims (12)

1.一种基于动态二阶差分阈值的BCG信号心率计算方法，所述方法包括：
采集获得BCG信号；
根据预设规则对BCG信号进行滤波预处理，获得滤波后的BCG信号；
根据动态二阶差分阈值算法处理滤波后的BCG信号，确定每个信号周期中的J波波峰位置，获得J波波峰序列；
根据所述J波波峰序列中的每个J波波峰位置，计算BCG信号主波的平均长度；根据所述平均长度以及采样率计算获得心率；
其中，根据动态二阶差分阈值算法处理滤波后的BCG信号，包括：
对滤波后的BCG信号进行离散化处理，获得K个信号点；按照预设的周期长度将所述K个信号点等分为N个周期，使每一段中包含一个BCG周期信号；
根据预设规则计算获得初始阈值；所述初始阈值包括初始差分阈值th1、初始振幅阈值上限th2以及初始振幅阈值下限th3；
根据初始阈值、预设的一阶差分条件以及二阶差分条件确定第一个周期的J波峰值点；
根据所述第一个周期的J波峰值点对初始阈值进行更新，获得新的阈值，即新的差分阈值、振幅阈值上限以及振幅阈值下限；并根据所述新的阈值以及预设的一阶差分条件以及二阶差分条件确定下一个周期的J波峰值点，直至获得所有J波峰值点；
其中，根据初始阈值、预设的一阶差分条件以及二阶差分条件确定第一个周期的J波峰值点，包括：
在第一个周期内获得主波升支上的一个离散点yi，所述yi通过满足yi+1-yi>th1且yi-yi-1>th1的条件；
判断从yi开始后的相邻4个点为yk，yk+1，yk+2，yk+3，是否满足差分预设条件；若满足所述差分预设条件，则记录yk+2为所述第一个周期的J波峰值点，其幅度记为Hnew；所述差分预设条件为：一阶差分条件：yk+1-yk>0；yk+2-yk+1>0；yk+3-yk+2<0；且同时满足二阶差分条件：yn＝yk+1-yk；yn+1＝yk+2-yk+1；yn+2＝yk+3-yk+2；yn+1-yn<0且yn+2-yn+1<0；
若不满足所述差分预设条件，则重新获得主波升支上的新的离散点，直至获得满足所述差分预设条件的情形，获得J波峰值点；
其中，根据所述第一个周期的J波峰值点对初始阈值进行更新，包括：
获得从yk+1到yk+2这一过程的差分最大值，记为Mnew；使用所述差分最大值Mnew更换计算初始差分阈值中该周期对应的差分最大值，并重新计算新的平均差分最大值M0；根据新的M0计算获得新的差分阈值；
将所述J波峰值点yk+2对应的振幅Hnew更换计算初始振幅阈值中对应该周期的振幅最大值，并重新计算新的平均振幅最大值H0；根据新的H0计算获得新的振幅阈值上限以及振幅阈值下限。
2.根据权利要求1所述的方法，其特征在于，所述根据预设规则对BCG信号进行滤波预处理，获得滤波后的BCG信号，包括：
通过截止频率为30Hz的低通滤波器滤除工频干扰及高频噪声；
根据经验模态分解算法，通过预设的多个阈值对原始的BCG信号进行经验模态分解，获得多阶IMF函数；
将含有噪声成分的IMF函数去除，对剩余的含有信号成分的IMF函数进行重构，获得滤波后的BCG信号。
3.根据权利要求2所述的方法，其特征在于，在所述通过截止频率为30Hz的低通滤波器滤除工频干扰及高频噪声前，所述方法还包括：
通过抖动标记方法，将由于体动而导致的振幅过大的区域去除；
通过滑动窗口对BCG信号进行平滑处理，去除异常值。
4.根据权利要求1所述的方法，其特征在于，所述根据预设规则计算获得初始阈值，包括：
对每个周期中所包括的信号点进行差分处理，获得每个周期对应的差分最大值；比较所述N个周期中每个周期对应的差分最大值，去掉其中的最大值及最小值，并对剩余的差分最大值求算术平均值，获得平均差分最大值M0；取初始差分阈值th1＝a*M0，其中0.5≤a≤1；
获得每个周期对应的振幅最大值，比较所述N个周期中每个周期对应的振幅最大值，去掉其中的最大值及最小值，并对剩余的振幅最大值求算术平均值，获得平均振幅最大值H0；取初始振幅阈值下限为th2＝b*H0，初始振幅阈值上限为th3＝c*H0，其中0≤b≤1；1≤c≤2。
5.根据权利要求1所述的方法，其特征在于，在确定每个信号周期中的J波波峰位置后，所述方法还包括：
确定每个J波波峰的前一波谷的幅值是否为该波峰所在周期内的最低值，将前一波谷幅值非最低值对应的J波波峰删除，获得新的J波波峰序列；
取所述新的J波波峰序列中的每一个J波波峰的点的左右各P个点进行局部极值检测；将局部极值检测出的振幅最高值对应的点替换该J波波峰，作为该周期内寻优后的J波波峰；
将新的J波波峰序列中所有寻优后的J波波峰作为最终的J波波峰序列。
6.根据权利要求1所述的方法，其特征在于：所述心率的计算公式为：
Pulse＝60/(Fs*L)
其中，Fs为采样率，L为主波的平均长度。
7.一种基于动态二阶差分阈值的BCG信号心率计算系统，所述系统包括：
BCG信号采集单元，所述BCG信号采集单元用于采集获得BCG信号；
滤波预处理单元，所述滤波预处理单元用于根据预设规则对BCG信号进行滤波预处理，获得滤波后的BCG信号；
J波波峰计算单元，所述J波波峰计算单元根据动态二阶差分阈值算法处理滤波后的BCG信号，确定每个信号周期中的J波波峰位置，获得J波波峰序列；
心率计算单元，所述心率计算单元用于根据所述J波波峰序列中的每个J波波峰位置，计算BCG信号主波的平均长度；所述心率计算单元用于根据所述平均长度以及采样率计算获得心率；
其中，所述J波波峰计算单元用于对滤波后的BCG信号进行离散化处理，获得K个信号点；按照预设的周期长度将所述K个信号点等分为N个周期，使每一段中包含一个BCG周期信号；
所述J波波峰计算单元用于根据预设规则计算获得初始阈值；所述初始阈值包括初始差分阈值th1、初始振幅阈值上限th2以及初始振幅阈值下限th3；
所述J波波峰计算单元用于根据初始阈值、预设的一阶差分条件以及二阶差分条件确定第一个周期的J波峰值点；
所述J波波峰计算单元用于根据所述第一个周期的J波峰值点对初始阈值进行更新，获得新的阈值，即新的差分阈值、振幅阈值上限以及振幅阈值下限；并根据所述新的阈值以及预设的一阶差分条件以及二阶差分条件确定下一个周期的J波峰值点，直至获得所有J波峰值点；
其中，所述J波波峰计算单元用于在第一个周期内获得主波升支上的一个离散点yi，所述yi通过满足yi+1-yi>th1且yi-yi-1>th1的条件；
所述J波波峰计算单元用于判断从yi开始后的相邻4个点为yk，yk+1，yk+2，yk+3，是否满足差分预设条件；若满足所述差分预设条件，则记录yk+2为所述第一个周期的J波峰值点，其幅度记为Hnew；所述差分预设条件为：一阶差分条件：yk+1-yk>0；yk+2-yk+1>0；yk+3-yk+2<0；且同时满足二阶差分条件：yn＝yk+1-yk；yn+1＝yk+2-yk+1；yn+2＝yk+3-yk+2；yn+1-yn<0且yn+2-yn+1<0；
若不满足所述差分预设条件，则重新获得主波升支上的新的离散点，直至获得满足所述差分预设条件的情形，获得J波峰值点；
其中，所述J波波峰计算单元用于获得从yk+1到yk+2这一过程的差分最大值，记为Mnew；使用所述差分最大值Mnew更换计算初始差分阈值中该周期对应的差分最大值，并重新计算新的平均差分最大值M0；根据新的M0计算获得新的差分阈值；
所述J波波峰计算单元用于将所述J波峰值点yk+2对应的振幅Hnew更换计算初始振幅阈值中对应该周期的振幅最大值，并重新计算新的平均振幅最大值H0；根据新的H0计算获得新的振幅阈值上限以及振幅阈值下限。
8.根据权利要求7所述的系统，其特征在于：
所述滤波预处理单元用于通过截止频率为30Hz的低通滤波器滤除工频干扰及高频噪声；
所述滤波预处理单元用于根据经验模态分解算法，通过预设的多个阈值对原始的BCG信号进行经验模态分解，获得多阶IMF函数；
所述滤波预处理单元用于将含有噪声成分的IMF函数去除，对剩余的含有信号成分的IMF函数进行重构，获得滤波后的BCG信号。
9.根据权利要求8所述的系统，其特征在于：
所述滤波预处理单元用于通过抖动标记方法，将由于体动而导致的振幅过大的区域去除；
所述滤波预处理单元用于通过滑动窗口对BCG信号进行平滑处理，去除异常值。
10.根据权利要求7所述的系统，其特征在于：
所述J波波峰计算单元用于对每个周期中所包括的信号点进行差分处理，获得每个周期对应的差分最大值；比较所述N个周期中每个周期对应的差分最大值，去掉其中的最大值及最小值，并对剩余的差分最大值求算术平均值，获得平均差分最大值M0；取初始差分阈值th1＝a*M0，其中0.5≤a≤1；
所述J波波峰计算单元用于获得每个周期对应的振幅最大值，比较所述N个周期中每个周期对应的振幅最大值，去掉其中的最大值及最小值，并对剩余的振幅最大值求算术平均值，获得平均振幅最大值H0；取初始振幅阈值下限为th2＝b*H0，初始振幅阈值上限为th3＝c*H0，其中0≤b≤1；1≤c≤2。
11.根据权利要求7所述的系统，其特征在于：所述系统还包括波峰寻优单元；
所述波峰寻优单元用于确定每个J波波峰的前一波谷的幅值是否为该波峰所在周期内的最低值，将前一波谷幅值非最低值对应的J波波峰删除，获得新的J波波峰序列；
所述波峰寻优单元用于取所述新的J波波峰序列中的每一个J波波峰的点的左右各P个点进行局部极值检测；将局部极值检测出的振幅最高值对应的点替换该J波波峰，作为该周期内寻优后的J波波峰；
所述波峰寻优单元用于将新的J波波峰序列中所有寻优后的J波波峰作为最终的J波波峰序列。
12.根据权利要求7所述的系统，其特征在于：所述心率计算单元计算心率的公式为：Pulse＝60/(Fs*L)
其中，Fs为采样率，L为主波的平均长度。
CN201911156128.6A 2019-11-22 2019-11-22 一种基于动态二阶差分阈值的bcg信号心率计算方法及系统 Active CN110916636B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911156128.6A CN110916636B (zh) 2019-11-22 2019-11-22 一种基于动态二阶差分阈值的bcg信号心率计算方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911156128.6A CN110916636B (zh) 2019-11-22 2019-11-22 一种基于动态二阶差分阈值的bcg信号心率计算方法及系统

Publications (2)

Publication Number Publication Date
CN110916636A CN110916636A (zh) 2020-03-27
CN110916636B true CN110916636B (zh) 2023-05-26

Family

ID=69850692

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911156128.6A Active CN110916636B (zh) 2019-11-22 2019-11-22 一种基于动态二阶差分阈值的bcg信号心率计算方法及系统

Country Status (1)

Country Link
CN (1) CN110916636B (zh)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111887858B (zh) * 2020-08-04 2021-05-04 西安电子科技大学 基于跨模态映射的心冲击图信号心率估计方法
CN112869733B (zh) * 2021-01-08 2021-12-24 广州中科新知科技有限公司 一种心冲击图实时心搏间期测算方法
CN113017613B (zh) * 2021-03-03 2022-05-06 四川大学华西医院 基于人工智能的心冲击波信号处理方法及计算机设备
CN113288072A (zh) * 2021-04-28 2021-08-24 新绎健康科技有限公司 一种脉搏波的主波检测方法及装置
CN113100727B (zh) * 2021-05-12 2023-09-19 深圳市通久电子有限公司 实时分析识别脉搏波峰的方法
CN113499059B (zh) * 2021-06-01 2022-07-05 武汉理工大学 基于光纤传感非接触式的bcg信号处理系统及方法
CN113647925A (zh) * 2021-07-05 2021-11-16 新绎健康科技有限公司 基于心冲击信号的心率确定方法及装置
CN114176548A (zh) * 2021-12-03 2022-03-15 新绎健康科技有限公司 基于模板匹配的心冲击信号心率计算方法及系统
CN114287919A (zh) * 2021-12-14 2022-04-08 深圳数联天下智能科技有限公司 基于心冲击信号的j波定位方法、装置、设备及介质
CN114010186B (zh) * 2022-01-11 2022-03-18 华南师范大学 一种心冲击图信号定位方法和计算机设备

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2010067297A1 (en) * 2008-12-11 2010-06-17 Koninklijke Philips Electronics N.V. Method and apparatus for the analysis of ballistocardiogram signals
CN104545863A (zh) * 2013-10-10 2015-04-29 上海宽带技术及应用工程研究中心 基于模糊模式识别的bcg心率提取方法及系统
CN107970590A (zh) * 2016-10-25 2018-05-01 四川理工学院 一种基于Android平台的跑步健身数据系统及方法
CN109602234A (zh) * 2019-01-02 2019-04-12 南京信息工程大学 一种用于检测睡眠质量的智能枕头及远程自动调整方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20190094725A (ko) * 2018-02-05 2019-08-14 삼성전자주식회사 식이 습관 관리 장치 및 방법

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2010067297A1 (en) * 2008-12-11 2010-06-17 Koninklijke Philips Electronics N.V. Method and apparatus for the analysis of ballistocardiogram signals
CN104545863A (zh) * 2013-10-10 2015-04-29 上海宽带技术及应用工程研究中心 基于模糊模式识别的bcg心率提取方法及系统
CN107970590A (zh) * 2016-10-25 2018-05-01 四川理工学院 一种基于Android平台的跑步健身数据系统及方法
CN109602234A (zh) * 2019-01-02 2019-04-12 南京信息工程大学 一种用于检测睡眠质量的智能枕头及远程自动调整方法

Also Published As

Publication number Publication date
CN110916636A (zh) 2020-03-27

Similar Documents

Publication Publication Date Title
CN110916636B (zh) 一种基于动态二阶差分阈值的bcg信号心率计算方法及系统
CN109907752B (zh) 一种去除运动伪影干扰与心电特征检测的心电诊断与监护系统
Strasser et al. Motion artifact removal in ECG signals using multi-resolution thresholding
CN103405227B (zh) 基于双层形态学滤波的心电信号预处理方法
Han et al. Electrocardiogram signal denoising based on a new improved wavelet thresholding
CN110680308B (zh) 基于改进emd与阈值法融合的心电信号去噪方法
CN107361764B (zh) 一种心电信号特征波形r波的快速提取方法
CN111160090B (zh) 一种bcg信号降噪方法及系统
Abbaspour et al. Evaluation of wavelet based methods in removing motion artifact from ECG signal
US20170086752A1 (en) Devices, systems, and methods for determining heart rate of a subject from noisy electrocardiogram data
Xin et al. ECG baseline wander correction based on mean-median filter and empirical mode decomposition
Phinyomark et al. EMG feature extraction for tolerance of white Gaussian noise
Almalchy et al. Noise removal from ECG signal based on filtering techniques
CN110507317B (zh) 一种心电信号r波的自适应ca-cfar定位方法
Swathi et al. R peak detection and feature extraction for the diagnosis of heart diseases
Janušauskas et al. Ensemble empirical mode decomposition based feature enhancement of cardio signals
Jain et al. An adaptive method for shrinking of wavelet coefficients for phonocardiogram denoising
CN111493821A (zh) 一种基于modwt及中值滤波的ppg信号实时去噪方法
CN110881958A (zh) 一种用于中医脉诊仪的脉搏信号非生理信号去除方法
Awodeyi et al. On the filtering of photoplethysmography signals
Suchetha et al. Denoising and arrhythmia classification using EMD based features and neural network
CN113925482A (zh) 心率计算方法、可穿戴电子设备和存储介质
Franchevska et al. The Method and Algorithm for Detecting the Fetal ECG Signal in the Presence of Interference
Georgieva-Tsaneva A novel photoplethysmographic noise removal method via wavelet transform to effective preprocessing
Hu et al. A robust beat-to-beat artifact detection algorithm for pulse wave

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant