技术领域Technical Field
本发明属于呼吸频率检测领域，具体涉及一种基于移动开窗法的自适应卡尔曼滤波的呼吸频率计算方法。The invention belongs to the field of respiratory frequency detection, and in particular relates to a respiratory frequency calculation method based on an adaptive Kalman filter of a moving windowing method.
背景技术Background technique
呼吸是人类维持生命的必需生理活动之一。呼吸频率(Respiratory Rate)作为重要的生理参数，直观反映人的即时身体状况。通过监测病患的呼吸频率，可以有效推算心率、细胞乳酸及二氧化碳的含量等生理指标，为后续医疗操作提供理论指导。Breathing is one of the essential physiological activities for humans to maintain life. Respiratory rate is an important physiological parameter that directly reflects a person's immediate physical condition. By monitoring the patient's respiratory rate, physiological indicators such as heart rate, cellular lactate and carbon dioxide content can be effectively calculated, providing theoretical guidance for subsequent medical operations.
目前呼吸检测的方法主要包括针对气体成分、光学、声学、肌电、压力及胸部阻抗进行检测等，在野外救援、特种训练、战场作战等复杂环境下，现有的常用于呼吸频率检测方法存在如下缺点：1)基于呼吸道传感器的方法，需要将测量装置内置到上呼吸道，严重干扰复杂环境下的救援工作；2)二氧化碳检测法或容积描记术需要专用设备支持，在复杂环境下较难部署；3)基于声学、光学或电磁学的识别方法，容易受到环境干扰，稳定性和可靠性较低；4)通过心电信号的分析获取呼吸波形的方法，需要心电图的支持，而在复杂环境下心电图的实时跟踪和动态测量难以实现；5)使用压力传感器，必须保持对胸部的压力；6)使用阻抗测量的方法需要针织、镀银等复杂工艺，不便于患者行动或大规模量产推广。At present, the methods of respiratory detection mainly include detection of gas composition, optics, acoustics, myoelectricity, pressure and chest impedance. In complex environments such as field rescue, special training, battlefield operations, etc., the existing commonly used respiratory rate detection methods have the following disadvantages: 1) Methods based on respiratory sensors require the measurement device to be built into the upper respiratory tract, which seriously interferes with rescue work in complex environments; 2) Carbon dioxide detection or volume volumography requires special equipment support and is difficult to deploy in complex environments; 3) Recognition methods based on acoustics, optics or electromagnetics are easily affected by environmental interference and have low stability and reliability; 4) Methods for obtaining respiratory waveforms through analysis of ECG signals require the support of electrocardiograms, but real-time tracking and dynamic measurement of electrocardiograms are difficult to achieve in complex environments; 5) When using pressure sensors, pressure on the chest must be maintained; 6) Impedance measurement methods require complex processes such as knitting and silver plating, which are not convenient for patient movement or large-scale mass production and promotion.
随着传感技术的发展和惯性测量单元(Inertial Measurement Unit，IMU)分辨率的提高，使用其作为运动传感器测量呼吸运动的数据进行呼吸频率检测逐渐变得可行。惯性测量单元集成了三轴陀螺仪、三轴加速度计以及三轴磁力计，各测量元器件使用微机电方案：加速度计集成弹簧质量块和压电或压阻式器件；角速度计集成科里奥利力的测量器件和机械结构；磁力计集成霍尔效应传感器。由于各测量计分别将测量部件集成到一个芯片上，因此测量过程具有频率高、体积小、质量轻、成本低的特点。利用IMU测量在静止或运动状态下的物体惯性物理量，如姿态角速度、加速度等，可以用于解算胸廓运动情况，实现直观反应呼吸运动的波形曲线以及对呼吸频率的计算。With the development of sensor technology and the improvement of the resolution of inertial measurement units (IMUs), it has become increasingly feasible to use them as motion sensors to measure respiratory movement data for respiratory rate detection. The inertial measurement unit integrates a three-axis gyroscope, a three-axis accelerometer, and a three-axis magnetometer. Each measuring component uses a micro-electromechanical solution: the accelerometer integrates a spring mass block and a piezoelectric or piezoresistive device; the angular velocity meter integrates the Coriolis force measurement device and mechanical structure; the magnetometer integrates a Hall effect sensor. Since each measuring meter integrates the measuring components onto a chip, the measurement process has the characteristics of high frequency, small size, light weight, and low cost. Using IMU to measure the inertial physical quantities of objects in a static or moving state, such as attitude angular velocity, acceleration, etc., can be used to resolve the movement of the chest, realize a waveform curve that intuitively reflects the respiratory movement, and calculate the respiratory rate.
现有技术中针对IMU解算呼吸频率的方法比较少见，主要技术壁垒是对呼吸信号的捕捉。呼吸信号频率范围较小，而人作为传感器载体在运动时姿态范围显著，运动姿态变化可能淹没呼吸信号。另外，针对IMU的姿态解算算法以传统卡尔曼滤波为主，传统卡尔曼滤波的噪声协方差阵均需要在滤波前预设并保持固定，因此其对初始位置敏感，不能针对各种环境的实时变化进行自动适应。现有针对呼吸曲线计算呼吸频率的方法中，主要采用峰值检测的方式，寻找峰值点的方法一般使用的是简单的阈值法或差分法，这两种方法都无法针对运动场景的变化进行自动适应，也容易受到噪声的干扰。In the prior art, methods for calculating respiratory frequency from IMU are relatively rare, and the main technical barrier is the capture of respiratory signals. The frequency range of respiratory signals is relatively small, and people as sensor carriers have a significant range of postures when in motion, and changes in motion postures may drown out respiratory signals. In addition, the posture calculation algorithm for IMU is mainly based on traditional Kalman filtering. The noise covariance matrix of traditional Kalman filtering needs to be preset and kept fixed before filtering, so it is sensitive to the initial position and cannot automatically adapt to real-time changes in various environments. Among the existing methods for calculating respiratory frequency from respiratory curves, peak detection is mainly used. The method for finding peak points generally uses a simple threshold method or a differential method. Both methods cannot automatically adapt to changes in motion scenes and are easily affected by noise.
发明内容Summary of the invention
本发明针对现有技术中传统呼吸频率检测方案对检测设备要求高，在变化的运动场景或复杂环境下稳定性和可靠性低等缺陷，提出一种基于移动开窗法的自适应卡尔曼滤波呼吸频率计算方法，以达到在特种训练等复杂环境下测量呼吸频率的目的。In view of the defects of traditional respiratory rate detection schemes in the prior art, such as high requirements on detection equipment and low stability and reliability in changing motion scenes or complex environments, the present invention proposes an adaptive Kalman filter respiratory rate calculation method based on the moving window method, so as to achieve the purpose of measuring respiratory rate in complex environments such as special training.
本发明是采用以下的技术方案实现的：一种基于移动开窗法的自适应卡尔曼滤波的呼吸频率方法，包括以下步骤：The present invention is implemented by adopting the following technical solution: a respiratory frequency method based on adaptive Kalman filtering of moving window method, comprising the following steps:
步骤A、运动姿态数据获取：Step A: Acquisition of motion posture data:
采用差分测量方案，使用胸带在前胸中心和后背中心固定两个IMU(前IMU和后IMU)，两个IMU同时测量载体的运动参数，将地理坐标系作为导航坐标系(n系)，记前IMU所在的载体坐标系为b1，后IMU所在的载体坐标系为b2，并取后IMU所在的载体坐标系为标准载体坐标系，在胸廓运动时将前加速度计的值向后加速度计对准。记在载体坐标系(b系)中前加速度计采集的值向后加速度计对准所需的差分值为将加速度的差分值映射到后IMU的坐标系，记从载体坐标系(b系)向导航坐标系(n系)转换时姿态转换矩阵为和则加速度的差分过程表示为：A differential measurement scheme is adopted. Two IMUs (front IMU and rear IMU) are fixed at the center of the front chest and the center of the back using a chest strap. The two IMUs measure the motion parameters of the carrier at the same time. The geographic coordinate system is used as the navigation coordinate system (n system). The carrier coordinate system where the front IMU is located is b1 , and the carrier coordinate system where the rear IMU is located is b2 . The carrier coordinate system where the rear IMU is located is taken as the standard carrier coordinate system. When the chest moves, the value of the front accelerometer is aligned with the rear accelerometer. The differential value required to align the value collected by the front accelerometer with the rear accelerometer in the carrier coordinate system (b system) is Map the differential value of acceleration to the coordinate system of the rear IMU, and record the attitude transformation matrix when converting from the carrier coordinate system (b system) to the navigation coordinate system (n system) as: and The differential process of acceleration is expressed as:
其中，为前IMU在其所在的载体坐标系下测量的加速度；为后IMU在其所在的载体坐标系下测量的加速度，和是过程量，分别指代前后IMU在导航坐标系下的加速度值，也是过程量，代表在导航坐标系下前IMU测量的加速度向后IMU测量的加速度对准所需的差分值。in, It is the acceleration measured by the front IMU in the carrier coordinate system; is the acceleration measured by the rear IMU in the carrier coordinate system, and are process quantities, referring to the acceleration values of the front and rear IMUs in the navigation coordinate system. It is also a process quantity, representing the differential value required to align the acceleration measured by the front IMU with the acceleration measured by the rear IMU in the navigation coordinate system.
步骤B、运动姿态数据解算：Step B: Motion posture data calculation:
根据IMU获得的传感器差分数据，基于自适应噪声及其协方差阵的卡尔曼滤波算法进行姿态解算，并结合移动开窗法和自适应补偿因子算法对协方差矩阵进行迭代最优估计，获得呼吸波形曲线；According to the sensor differential data obtained by IMU, the attitude is solved based on the Kalman filter algorithm of adaptive noise and its covariance matrix, and the covariance matrix is iteratively estimated in combination with the moving window method and adaptive compensation factor algorithm to obtain the breathing waveform curve;
在进行姿态解算时，具体采用以下方式。When performing attitude solving, the following method is specifically adopted.
(1)以自适应卡尔曼滤波为基准，对于系统噪声协方差阵，在系统开始运行时初始化系统噪声为环境值，然后固定系统噪声，每隔一定的时间(如1个呼吸窗口)更新一次，而非实时更新；对于测量噪声协方差阵，在近似静止以及呼吸频率过高时不进行实时补偿；(1) Based on the adaptive Kalman filter, the system noise covariance matrix is initialized to the ambient value when the system starts running, and then the system noise is fixed and updated once every certain period of time (such as one breathing window) instead of real-time update; for the measurement noise covariance matrix, no real-time compensation is performed when the system is nearly stationary or the breathing rate is too high;
(2)为预防量测产生的异常或对异常进行纠正，对IMU测量的原始信号进行低通滤波，以降低高频噪声分量的影响；滤波后利用移动开窗法对协方差矩阵进行迭代最优估计，提高对单次量测异常的健壮性；使用自适应补偿因子对协方差矩阵进行修正；在残差可能失配时调整自适应补偿因子纠正协方差阵，具体的姿态解算过程如下：(2) To prevent or correct measurement anomalies, the original signal measured by the IMU is low-pass filtered to reduce the impact of high-frequency noise components; after filtering, the covariance matrix is iteratively estimated using the moving window method to improve the robustness to single measurement anomalies; the covariance matrix is corrected using an adaptive compensation factor; when the residual may be mismatched, the adaptive compensation factor is adjusted to correct the covariance matrix. The specific attitude solution process is as follows:
1)首先对系统噪声和观测噪声协方差阵进行估计：1) First, estimate the covariance matrix of system noise and observation noise:
系统噪声协方差矩阵的改变量如下式所示：The change in the system noise covariance matrix is shown below:
系统噪声协方差矩阵的更新公式如下所示：The update formula of the system noise covariance matrix is as follows:
观测噪声协方差矩阵的改变量如下式所示：The change in the observation noise covariance matrix is shown as follows:
观测噪声协方差矩阵的更新公式如下所示：The update formula of the observation noise covariance matrix is as follows:
其中，dk是更新参量，计算式为(0.95≤b≤0.995)，b是遗忘因子，在滤波前进行设置；W(k)和V(k)表示的是第k个历元的系统的过程噪声和系统的观测噪声，记它们的统计特性关系为E(W(k))＝qk，E(V(k))＝rk；K(k)表示的是第k个历元的卡尔曼增益系数，A是对系统状态的参数矩阵，在卡尔曼滤波中它和对系统控制的参数矩阵B满足如下的关系X(k)＝AX(k-1)+BU(k-1)+W(k)；H是测量系统的参数矩阵，它和在k时刻的系统测量值Z(k)满足Z(k)＝HX(k)+V(k)；是每一轮的预测协方差阵，是每一轮的更新协方差阵。Among them, dk is the update parameter, and the calculation formula is (0.95≤b≤0.995), b is the forgetting factor, which is set before filtering; W(k) and V(k) represent the process noise and observation noise of the system at the kth epoch, and their statistical characteristics are E(W(k))＝qk , E(V(k))＝rk ; K(k) represents the Kalman gain coefficient of the kth epoch, A is the parameter matrix of the system state, and in the Kalman filter, it and the parameter matrix B of the system control satisfy the following relationship X(k)＝AX(k-1)+BU(k-1)+W(k); H is the parameter matrix of the measurement system, and it and the system measurement value Z(k) at time k satisfy Z(k)＝HX(k)+V(k); is the prediction covariance matrix for each round, is the updated covariance matrix for each round.
2)引入自适应补偿因子αk，对预测协方差矩阵进行补偿，具体补偿过程如下：2) Introduce the adaptive compensation factor αk to compensate the prediction covariance matrix. The specific compensation process is as follows:
其中，Q'是系统噪声协方差阵，此时已固定。αk的表示式如下：Among them, Q' is the system noise covariance matrix, which is fixed at this time. The expression of αk is as follows:
对于使用自适应的指数加权法，其求解公式如下：for Using the adaptive exponential weighting method, the solution formula is as follows:
其中，ek表示残差，满足关系式是残差理论协方差的估计值，tr()表示对矩阵进行求迹运算，E[]表示对矩阵求期望。Among them, ek represents the residual, satisfying the relationship is the estimate of the residual theoretical covariance, tr() represents the trace operation of the matrix, and E[] represents the expectation of the matrix.
步骤C、呼吸频率求解：Step C, respiratory rate solution:
利用快速离散傅里叶变换处理后的频谱图得到信号频率的先验估计值，并据此调整阈值与采样间隔等自适应参数；根据调整好的自适应参数，利用基于步长的自适应差分法计算呼吸频率，实现更加准确的进行呼吸频率的求解，实现不同情景下的呼吸频率计算。The frequency spectrum after fast discrete Fourier transform processing is used to obtain a priori estimated value of the signal frequency, and adaptive parameters such as the threshold and sampling interval are adjusted accordingly; according to the adjusted adaptive parameters, the respiratory frequency is calculated using the step-based adaptive difference method to achieve a more accurate solution to the respiratory frequency and realize the calculation of the respiratory frequency in different scenarios.
具体的步骤如下：The specific steps are as follows:
步骤C1、使用自适应差分算法求呼吸波形曲线的下一个峰谷(峰谷即波峰或波谷)；Step C1, using an adaptive differential algorithm to find the next peak or valley (peak or valley is a peak or valley) of the respiratory waveform curve;
步骤C2、利用快速离散傅里叶变换得到频谱图，并由此计算信号频率的先验估计值，根据该先验估计值调整有效数据阈值与差分间隔(以下简称自适应参数)，以给下一步自适应差分法提供相对应的自适应参数；Step C2, using fast discrete Fourier transform to obtain a spectrum diagram, and thereby calculating a priori estimated value of the signal frequency, and adjusting the effective data threshold and the differential interval (hereinafter referred to as adaptive parameters) according to the prior estimated value, so as to provide corresponding adaptive parameters for the next step of adaptive differential method;
1)计算姿态解算得到的呼吸波形曲线上各点数据的移动开窗平均值；1) Calculate the moving window average of the data at each point on the respiratory waveform curve obtained by posture solution;
其中，nL,nR为待计算的点两边的点数，cn为窗口大小，pi为姿态解算得到的数据，gi为经由移动开窗法计算得到的结果，此处计算窗长与自适应卡尔曼滤波处的窗长一致；WherenL ,nR are the number of points on both sides of the point to be calculated,cn is the window size,pi is the data obtained by attitude solution, and giis the result calculated by the moving window method. The calculation window length here is consistent with the window length at the adaptive Kalman filter;
2)对gi进行1024点快速离散傅里叶变换，得到频谱数组，用下述公式计算最大幅值索引对应的信号频率作为呼吸频率的先验估计值；2) Perform a 1024-point fast discrete Fourier transform on gi to obtain a spectrum array, and use the following formula to calculate the signal frequency corresponding to the maximum amplitude index as a priori estimate of the respiratory frequency;
其中，n是最大幅值的索引值，F是信号的实际采样频率，N是采集点数；Among them, n is the index value of the maximum amplitude, F is the actual sampling frequency of the signal, and N is the number of acquisition points;
步骤C3、采用基于步长的自适应差分法对数据进行差分、计数和筛选，最终计算呼吸频率，其具体步骤包括：Step C3: using the step-based adaptive difference method to differentiate, count and filter the data, and finally calculate the respiratory rate. The specific steps include:
(1)根据上一步设置的采样间隔和自适应的差分阈值等自适应参数，对数据进行差分，以判断呼吸波形数据的每个离散点是否到达波峰和波谷；(1) Differentiate the data according to the adaptive parameters such as the sampling interval and the adaptive differential threshold set in the previous step to determine whether each discrete point of the respiratory waveform data reaches the peak and the trough;
(2)采用基于步长的呼吸峰谷判断法遍历计数呼吸窗口内呼吸曲线的呼吸波峰谷数，并分别计算当前呼吸步长、预测呼吸步长；根据实际呼吸步长和预测呼吸步长的偏差范围为呼吸峰谷赋值当前置信度和累积置信度，根据置信度的变化情况判断呼吸波峰谷的有效性；(2) Using the step-based breathing peak and valley judgment method to traverse and count the number of breathing peaks and valleys of the breathing curve in the breathing window, and respectively calculate the current breathing step length and the predicted breathing step length; assigning the current confidence and cumulative confidence to the breathing peak and valley according to the deviation range between the actual breathing step length and the predicted breathing step length, and judging the validity of the breathing peak and valley according to the change of the confidence;
(3)筛选得到有效的呼吸波峰谷数量，并据此计算呼吸频率。(3) Screen and obtain the effective number of respiratory peaks and valleys, and calculate the respiratory rate based on this.
与现有技术相比，本发明的优点和积极效果在于：Compared with the prior art, the advantages and positive effects of the present invention are:
本方案采用差分测量方案，将两个IMU传感器分别固定在前胸中心与后背中心进行同时测量，取加速度相对值(差分值)作为胸廓运动变化的依据，有效克服了运动姿态变化可能淹没呼吸信号的问题。另外，以自适应噪声及其协方差阵的卡尔曼滤波算法作为基准进行姿态解算，克服传统卡尔曼滤波对初始位置敏感和固定噪声协方差阵的限制，提高滤波解算对环境的适应性；This solution adopts a differential measurement solution, fixing two IMU sensors at the center of the front chest and the center of the back for simultaneous measurement, taking the relative value of acceleration (differential value) as the basis for the change of chest movement, effectively overcoming the problem that the change of movement posture may drown out the breathing signal. In addition, the Kalman filter algorithm with adaptive noise and its covariance matrix is used as a benchmark for posture solution, overcoming the limitations of traditional Kalman filter on initial position sensitivity and fixed noise covariance matrix, and improving the adaptability of filter solution to the environment;
采用结合移动开窗法、自适应指数加权和自适应补偿因子的算法，一方面突出呼吸过程中的趋势变化，另一方面在检测到新息失配时调整补偿因子以抑制滤波发散，提高滤波算法在复杂环境中的稳定性，有效防止了自适应噪声滤波过程中可能出现的滤波发散问题；利用基于步长的自适应差分法和傅里叶变换方法实现根据呼吸曲线有效地计算呼吸频率，一方面根据步长预测和自适应差分减小了呼吸波形曲线中的毛刺影响和小幅震荡的干扰，有效提高了测量的灵活性和准确性；另一方面根据傅里叶变换对信号频率进行先验估计、阈值调整等方法，进一步增强计算结果的准确性。An algorithm combining moving windowing, adaptive exponential weighting and adaptive compensation factor is adopted. On the one hand, it highlights the trend changes in the breathing process. On the other hand, it adjusts the compensation factor to suppress the filter divergence when the new information mismatch is detected, thereby improving the stability of the filtering algorithm in complex environments and effectively preventing the filter divergence problem that may occur in the adaptive noise filtering process. The step-based adaptive difference method and Fourier transform method are used to effectively calculate the respiratory frequency according to the breathing curve. On the one hand, the step-based prediction and adaptive difference reduce the influence of glitches and interference of small oscillations in the respiratory waveform curve, effectively improving the flexibility and accuracy of the measurement. On the other hand, the Fourier transform is used to perform a priori estimation of the signal frequency, threshold adjustment and other methods to further enhance the accuracy of the calculation results.
本方案可以在室外运动训练、抢险救灾、战场作战等复杂环境场景下使用，适应复杂环境或特种训练等复杂或变化场景，测量方法具有较好的适应性和鲁棒性。This solution can be used in complex environment scenarios such as outdoor sports training, disaster relief, battlefield operations, etc. It can adapt to complex or changing scenarios such as complex environments or special training, and the measurement method has good adaptability and robustness.
附图说明BRIEF DESCRIPTION OF THE DRAWINGS
图1为本发明实施例所述呼吸频率检测方法流程示意图；FIG1 is a schematic flow chart of a respiratory rate detection method according to an embodiment of the present invention;
图2为本发明实施例IMU传感器差分测量示意图；FIG2 is a schematic diagram of differential measurement of an IMU sensor according to an embodiment of the present invention;
图3为本发明实施例呼吸段各组成部分示意图；FIG3 is a schematic diagram of the components of the breathing segment according to an embodiment of the present invention;
图4为本发明实施例静坐时KF、AKF与MAKF三种方法解算呼吸曲线比较图；FIG4 is a comparison diagram of breathing curves calculated by three methods, namely, KF, AKF and MAKF, during meditation according to an embodiment of the present invention;
图5为本发明实施例走路时KF、AKF与MAKF三种方法解算呼吸曲线比较图；FIG5 is a comparison diagram of breathing curves calculated by three methods, namely, KF, AKF and MAKF, when walking according to an embodiment of the present invention;
图6为本发明实施例跑步时KF与MAKF两种方法解算呼吸曲线比较图。FIG6 is a comparison diagram of the breathing curves calculated by the KF and MAKF methods when running according to an embodiment of the present invention.
具体实施方式Detailed ways
为了能够更加清楚地理解本发明的上述目的、特征和优点，下面结合附图及实施例对本发明做进一步说明。在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用不同于在此描述的其他方式来实施，因此，本发明并不限于下面公开的具体实施例。In order to more clearly understand the above-mentioned purpose, features and advantages of the present invention, the present invention is further described below in conjunction with the accompanying drawings and embodiments. In the following description, many specific details are set forth to facilitate a full understanding of the present invention, but the present invention can also be implemented in other ways different from those described herein, and therefore, the present invention is not limited to the specific embodiments disclosed below.
本实施例提出一种基于移动开窗法的自适应卡尔曼滤波的呼吸频率计算方法，其流程如图1所示，包括以下步骤：This embodiment proposes a method for calculating respiratory frequency based on an adaptive Kalman filter using a moving windowing method, the process of which is shown in FIG1 and includes the following steps:
(1)首先通过差分测量方案获得测量参数，基于自适应噪声及其协方差阵的卡尔曼滤波算法进行姿态解算获得呼吸波形曲线，克服原始卡尔曼滤波对初始位置敏感、噪声协方差阵固定的缺点；(1) First, the measurement parameters are obtained through a differential measurement scheme, and the posture is solved based on the Kalman filter algorithm of adaptive noise and its covariance matrix to obtain the breathing waveform curve, which overcomes the shortcomings of the original Kalman filter that is sensitive to the initial position and the noise covariance matrix is fixed;
(2)以自适应噪声及其协方差阵的卡尔曼滤波算法作为基准进行姿态解算，针对滤波发散的问题，使用结合移动开窗法、自适应指数加权和自适应补偿因子的算法对协方差矩阵进行迭代最优估计，在检测到新息失配时调整补偿因子以抑制发散，提高滤波算法对复杂环境的鲁棒性和稳定性；(2) The Kalman filter algorithm with adaptive noise and its covariance matrix is used as the benchmark for attitude solution. To address the problem of filter divergence, an algorithm combining moving windowing, adaptive exponential weighting and adaptive compensation factor is used to iteratively optimally estimate the covariance matrix. When a new information mismatch is detected, the compensation factor is adjusted to suppress divergence, thereby improving the robustness and stability of the filtering algorithm in complex environments.
(3)根据离散呼吸波形曲线，使用基于步长的自适应差分法计算呼吸频率，通过估计呼吸步长和自适应差分的方法有效减小呼吸波形曲线中的毛刺影响及小幅震荡的干扰；使用傅里叶变换求信号频率的估计值调整自适应阈值和自适应差分间隔，以适应不同情景下的呼吸频率计算。(3) According to the discrete respiratory waveform curve, the respiratory frequency is calculated using the step-based adaptive difference method. The influence of glitches and interference from small oscillations in the respiratory waveform curve is effectively reduced by estimating the respiratory step length and the adaptive difference method. The adaptive threshold and adaptive difference interval are adjusted by using the Fourier transform to estimate the signal frequency to adapt to the respiratory frequency calculation in different scenarios.
具体的实施步骤与原理如下：The specific implementation steps and principles are as follows:
1、数据采集1. Data Collection
1.1差分采集方法1.1 Differential Acquisition Method
物体在静止时三轴加速度矢量和约等于重力加速度，三轴角速度矢量和为地球自转角速度，约为7.292×10-5rad/s，数量级显著小于便携运动传感器的分辨率，因此可以安全地直接使用运动传感器对运动物体的六轴数据进行测量。常规的测量方法是将运动传感器通过胸带绑定到胸部前廓，在人与IMU组成的载体右手坐标系(b系)中，初始情况下，运动传感器Z轴负向对准胸部，正向朝载体正前方，X轴和Y轴正向则分别对应载体正右方和载体正下方。When an object is stationary, the sum of the three-axis acceleration vector is approximately equal to the acceleration of gravity, and the sum of the three-axis angular velocity vector is the angular velocity of the earth's rotation, which is approximately 7.292×10-5 rad/s, which is significantly smaller than the resolution of the portable motion sensor. Therefore, the six-axis data of the moving object can be safely measured directly using the motion sensor. The conventional measurement method is to bind the motion sensor to the front contour of the chest through a chest strap. In the right-handed coordinate system (b system) of the carrier composed of the person and the IMU, initially, the negative direction of the Z axis of the motion sensor is aimed at the chest, the positive direction is toward the front of the carrier, and the positive directions of the X axis and Y axis correspond to the right and bottom of the carrier, respectively.
然而，由于呼吸信号频率范围较小，变化较为微弱，而运动时姿态范围更加显著，变化较为强烈，因此运动时人本身的姿态变化可能淹没呼吸信号。为达到在特种训练等复杂环境下测量呼吸频率的目的，本实施例采用了一种差分测量方案，即将运动传感器(IMU)A、B分别固定在前胸中心与后背中心进行同时测量，根据A和B运动传感器的陀螺仪和加速度计的差值解算胸部运动情况，以有效减小外部运动对呼吸时胸廓的舒张与收缩的干扰，测量过程如图2所示，实线代表呼吸前，虚线代表呼吸后，呼吸过程中AB的相对位置会发生变化，变成A’和B’，可以根据相对位置变化计算差分值。However, since the frequency range of the respiratory signal is small and the changes are relatively weak, while the posture range during exercise is more significant and the changes are relatively strong, the posture changes of the person during exercise may drown out the respiratory signal. In order to achieve the purpose of measuring the respiratory frequency in complex environments such as special training, this embodiment adopts a differential measurement scheme, that is, the motion sensors (IMU) A and B are fixed at the center of the front chest and the center of the back for simultaneous measurement, and the chest movement is calculated based on the difference between the gyroscope and accelerometer of the A and B motion sensors to effectively reduce the interference of external movement on the expansion and contraction of the chest during breathing. The measurement process is shown in Figure 2, the solid line represents before breathing, and the dotted line represents after breathing. The relative position of AB will change during breathing and become A' and B', and the differential value can be calculated based on the relative position change.
1.2建立姿态解算模型1.2 Establishing the attitude solution model
欧拉角是表示三维旋转的描述方法，针对三维空间里的一个参考系，任何坐标系的取向，都可以用三个欧拉角来表现。通过旋转矩阵可以表示坐标系间转换关系，并计算欧拉角。设惯性测量单元IMU在第n个时刻的姿态角度分别为r、p、y，它们分别代表IMU从b系初始位置，绕Z轴旋转yaw角度，绕Y轴旋转pitch角度，绕X轴旋转roll角度得到的最终姿态。Euler angles are a method of describing three-dimensional rotation. For a reference system in three-dimensional space, the orientation of any coordinate system can be expressed by three Euler angles. The transformation relationship between coordinate systems can be represented by the rotation matrix, and the Euler angles can be calculated. Assume that the attitude angles of the inertial measurement unit IMU at the nth moment are r, p, and y, which represent the final attitude of the IMU from the initial position of the b system, rotating around the Z axis by the yaw angle, rotating around the Y axis by the pitch angle, and rotating around the X axis by the roll angle.
姿态旋转矩阵为：The attitude rotation matrix is:
由加速度计测量3轴方向上感受到的加速度，当物体静止时，其本身仅受重力加速度，根据相对运动理论，其读取到的加速度数据为[0,0,g]，当加速度计旋转时，重力加速度会在加速度计的3轴上产生分量，加速度计读取到的值即为初始静止时的向量在b系下的新坐标。故加速度计与姿态角对应关系如式2所示：The acceleration felt in the three-axis direction is measured by the accelerometer. When the object is stationary, it is only subject to the acceleration of gravity. According to the theory of relative motion, the acceleration data read is [0,0,g]. When the accelerometer rotates, the acceleration of gravity will produce components on the three axes of the accelerometer. The value read by the accelerometer is the new coordinate of the vector in the b system at the initial stationary state. Therefore, the corresponding relationship between the accelerometer and the attitude angle is shown in Formula 2:
由式2可得roll和pitch角的更新公式：From formula 2, we can get the update formula of roll and pitch angle:
由陀螺仪测量绕3个轴转动的角速度积分可解算得到角度。设n+1时刻的姿态角为r+△r,p+△p,y+△y，则陀螺仪的姿态更新公式如下所示：The angle can be calculated by integrating the angular velocity around the three axes measured by the gyroscope. Assuming the attitude angle at time n+1 is r+△r,p+△p,y+△y, the attitude update formula of the gyroscope is as follows:
1.3建立加速度标定及差分模型1.3 Establish acceleration calibration and differential model
系统启动时记录前后IMU的加速度计差值(下面将A称为前IMU，B称为后IMU)，此时认为加速度计处于同一坐标系中，记前加速度计采集的值向后加速度计对准所需的差分值表示为：When the system starts, the accelerometer difference between the front and rear IMUs is recorded (hereinafter referred to as A front IMU, B rear IMU). At this time, the accelerometers are considered to be in the same coordinate system. The difference value required to align the value collected by the front accelerometer with the rear accelerometer is expressed as:
前后IMU测量出人体前后的三轴加速度和三轴角速度后，需要对前后加速度求取差分值。由于呼吸时后背的变化较缓，而前胸的变化较为明显，因此将加速度的差分值映射到后IMU的坐标系。设前IMU的坐标系为b1，后IMU的坐标系为b2，姿态转换矩阵为和同时使用式5进行加速度的零偏补偿，则加速度的差分过程表示为：After the front and rear IMUs measure the three-axis acceleration and three-axis angular velocity of the human body, it is necessary to calculate the difference between the front and rear accelerations. Since the back changes slowly during breathing, while the chest changes more obviously, the difference in acceleration is mapped to the coordinate system of the rear IMU. Assume that the coordinate system of the front IMU is b1 , the coordinate system of the rear IMU is b2 , and the attitude transformation matrix is and At the same time, using equation 5 to compensate for the acceleration bias, the acceleration differential process is expressed as:
其中，为前IMU在其所在的载体坐标系下测量的加速度；为后IMU在其所在的载体坐标系下测量的加速度，和是过程量，分别指代前后IMU在地理坐标系下的加速度值，也是过程量，代表在导航坐标系(n系)下前IMU测量的加速度向后IMU测量的加速度对准所需的差分值。求解出加速度的差分值之后，根据后IMU的姿态角和加速度的差分值进行自适应卡尔曼滤波即可得到呼吸离散曲线。in, It is the acceleration measured by the front IMU in the carrier coordinate system; is the acceleration measured by the rear IMU in the carrier coordinate system, and is a process quantity, which refers to the acceleration values of the front and rear IMUs in the geographic coordinate system. It is also a process quantity, representing the differential value required to align the acceleration measured by the front IMU with the acceleration measured by the rear IMU in the navigation coordinate system (n system). After solving the differential value of the acceleration, an adaptive Kalman filter is performed based on the differential value of the attitude angle and acceleration of the rear IMU to obtain the breathing discrete curve.
2、改进自适应卡尔曼滤波进行姿态解算2. Improved adaptive Kalman filter for attitude calculation
卡尔曼滤波器最优估计过程预设了系统噪声和测量噪声，通常情况下难以根据经验预设精确匹配噪声参数的值。自适应卡尔曼滤波通过引入自适应性的噪声及协方差阵，对噪声进行在线调整、估计和补偿，降低了传统卡尔曼滤波器的误差，较好地解决了传统卡尔曼滤波忽略的噪声时变问题。然而自适应卡尔曼滤波在高频情况下可能过补偿导致残差失配、滤波发散。因此，本实施例以Sage-Husa自适应卡尔曼滤波为基准，对其进行改进，以得到适配复杂环境下且更加稳定更加可靠的滤波算法。具体地：The Kalman filter optimal estimation process presets the system noise and measurement noise. It is usually difficult to preset the value of the noise parameter based on experience to accurately match it. The adaptive Kalman filter reduces the error of the traditional Kalman filter by introducing adaptive noise and covariance matrix to adjust, estimate and compensate for the noise online, and better solves the problem of noise time variation ignored by the traditional Kalman filter. However, the adaptive Kalman filter may over-compensate in high-frequency conditions, resulting in residual mismatch and filter divergence. Therefore, this embodiment takes the Sage-Husa adaptive Kalman filter as a benchmark and improves it to obtain a filtering algorithm that is more stable and reliable and adaptable to complex environments. Specifically:
2.1本实施例以Sage-Husa自适应卡尔曼滤波为基准，对其进行简化。2.1 This embodiment is based on the Sage-Husa adaptive Kalman filter and simplifies it.
1)固定系统噪声：系统噪声一般变化不大，实时补偿意义不大，反而会加重计算负担，并可能导致残差更新失配。因此对于系统噪声协方差阵，在系统开始运行时初始化一次后固定，仅在每隔一定的时间(如1个呼吸窗口)更新一次，以适应环境变化或运动场景变化；1) Fixed system noise: System noise generally does not change much, so real-time compensation is not very meaningful. Instead, it will increase the computational burden and may cause residual update mismatch. Therefore, the system noise covariance matrix is initialized once when the system starts running and then fixed. It is only updated once at regular intervals (such as a breathing window) to adapt to environmental changes or motion scene changes;
2)自适应固定测量噪声：对于测量噪声协方差阵，由于在近似静止时IMU存在零漂，在高频运动时IMU的误差会逐渐增大，这两种情况都可能导致滤波发散。因此在系统近似静止时或在呼吸频率过高时(显著大于60次每分)不进行实时补偿，保持测量噪声及其协方差稳定不变，以保持姿态解算结果稳定。2) Adaptive fixed measurement noise: For the measurement noise covariance matrix, since the IMU has zero drift when it is approximately stationary, the error of the IMU will gradually increase during high-frequency movement. Both of these situations may cause filter divergence. Therefore, when the system is approximately stationary or when the breathing rate is too high (significantly greater than 60 times per minute), no real-time compensation is performed to keep the measurement noise and its covariance stable and unchanged, so as to keep the attitude solution result stable.
2.2本实施例对滤波进行高频噪声的量测异常进行预防或纠正2.2 This embodiment prevents or corrects abnormal measurement of high-frequency noise in filtering
1)对IMU测量的原始信号进行低通滤波，以尽量去除实时估计测量噪声时受到的高频噪声分量影响；1) Perform low-pass filtering on the original signal measured by the IMU to remove the influence of high-frequency noise components when estimating measurement noise in real time;
2)替代Sage-Husa自适应卡尔曼滤波中普通迭代最优估计，而是利用移动开窗法对协方差矩阵进行迭代最优估计，原公式与替代公式如下：2) Instead of the ordinary iterative optimal estimation in the Sage-Husa adaptive Kalman filter, the moving window method is used to iteratively optimally estimate the covariance matrix. The original formula and the alternative formula are as follows:
原系统噪声协方差矩阵的更新公式如式7所示：The update formula of the original system noise covariance matrix is shown in Equation 7:
采用移动开窗法后系统噪声协方差矩阵的改变量如式8所示：The change in the system noise covariance matrix after using the moving window method is shown in Equation 8:
采用移动开窗法后系统噪声协方差矩阵的更新公式如式9所示：The updating formula of the system noise covariance matrix after using the moving window method is shown in Equation 9:
原系统观测噪声协方差矩阵更新公式如式10所示：The update formula of the original system observation noise covariance matrix is shown in formula 10:
采用移动开窗法后观测噪声协方差的改变量如式11所示。The change in the observation noise covariance after using the moving window method is shown in Equation 11.
采用移动开窗法后观测噪声协方差矩阵的更新公式如式12所示。The update formula of the observation noise covariance matrix after using the moving window method is shown in Equation 12.
其中，dk是更新参量，计算式为(0.95≤b≤0.995)，b是遗忘因子，在滤波前进行设置；W(k)和V(k)表示的是第k个历元的系统的过程噪声和系统的观测噪声，记它们的统计特性关系为E(W(k))＝qk，E(V(k))＝rk；K(k)表示的是第k个历元的卡尔曼增益系数，A是对系统状态的参数矩阵；AT表示对矩阵A进行转置，在卡尔曼滤波中它和对系统控制的参数矩阵B满足如下的关系X(k)＝AX(k-1)+BU(k-1)+W(k)；H是测量系统的参数矩阵，它和在k时刻的系统测量值Z(k)满足Z(k)＝HX(k)+V(k)；是每一轮的预测协方差阵，是每一轮的更新协方差阵。Among them, dk is the update parameter, and the calculation formula is (0.95≤b≤0.995), b is the forgetting factor, which is set before filtering; W(k) and V(k) represent the process noise and observation noise of the system at the kth epoch, and their statistical characteristics are E(W(k))＝qk , E(V(k))＝rk ; K(k) represents the Kalman gain coefficient of the kth epoch, and A is the parameter matrix of the system state;AT represents the transposition of the matrix A, and in the Kalman filter, it and the parameter matrix B for system control satisfy the following relationship X(k)＝AX(k-1)+BU(k-1)+W(k); H is the parameter matrix of the measurement system, and it and the system measurement value Z(k) at time k satisfy Z(k)＝HX(k)+V(k); is the prediction covariance matrix for each round, is the updated covariance matrix for each round.
3)引入自适应补偿因子αk，对预测协方差矩阵进行补偿，方法如下：3) Introduce the adaptive compensation factor αk to compensate the prediction covariance matrix. The method is as follows:
其中，αk表示式如下：Among them, αk is expressed as follows:
其中，ek表示残差，满足关系式是残差理论协方差的估计值，tr()表示求迹运算，E[]表示求期望。Among them, ek represents the residual, satisfying the relationship is the estimate of the residual theoretical covariance, tr() represents the trace operation, and E[] represents the expectation.
对于使用自适应的指数加权法，其求解公式如下：for Using the adaptive exponential weighting method, the solution formula is as follows:
3、频率计算3. Frequency calculation
3.1呼吸步长预测偏差与置信度关系3.1 Relationship between breathing step length prediction deviation and confidence
由于呼吸运动具有规律性，正常情况下每个呼吸步长应落在预测呼吸步长的0.9～1.1倍范围内。然而实际测量过程中可能存在运动场景变化或量测过程异常等情况致使呼吸频率产生渐变趋势，因此实际呼吸步长与预测呼吸步长的理想范围存在差距。Due to the regularity of respiratory movement, under normal circumstances, each breathing step should fall within the range of 0.9 to 1.1 times the predicted breathing step. However, in the actual measurement process, there may be changes in the motion scene or abnormalities in the measurement process, which may cause the respiratory frequency to change gradually. Therefore, there is a gap between the actual breathing step and the ideal range of the predicted breathing step.
本实施例根据呼吸频率渐变的统计规律设置置信度，捕捉潜在的呼吸频率渐变趋势，提高了对毛刺、振荡、局部极值等情况的容忍程度，减小了呼吸计数偏差。呼吸步长预测偏差范围的值表示的是实际呼吸步长与预测呼吸步长的倍数，其相对应的置信度和最大出现次数如表1所示。This embodiment sets the confidence level according to the statistical law of the gradual change of the respiratory frequency, captures the potential gradual change trend of the respiratory frequency, improves the tolerance for glitches, oscillations, local extreme values, etc., and reduces the respiratory count deviation. The value of the respiratory step prediction deviation range represents the multiple of the actual respiratory step length and the predicted respiratory step length, and its corresponding confidence level and maximum number of occurrences are shown in Table 1.
表1呼吸步长预测偏差与置信度关系Table 1 Relationship between breathing step length prediction deviation and confidence
3.2基于步长的自适应差分法3.2 Adaptive difference method based on step size
传统计算呼吸频率的方法是差分法或阈值法，差分法是通过斜率变化来判断波峰和波谷，该方式特别容易受噪声影响，计数结果严重偏大；阈值法是通过计数大于或小于阈值的点作为峰值点，该方式在运动场景变化时难以确定阈值大小，计数结果一般偏小。The traditional methods for calculating respiratory rate are differential method or threshold method. The differential method determines the peaks and troughs by the change in slope. This method is particularly susceptible to noise and the counting results are seriously biased. The threshold method uses the points that are greater than or less than the threshold as the peak points. This method is difficult to determine the threshold size when the motion scene changes, and the counting results are generally biased small.
为有效计算呼吸频率，本实施例提出一种基于步长的自适应差分法。In order to effectively calculate the respiratory frequency, this embodiment proposes an adaptive difference method based on step length.
首先，将传统的阈值设置方式改为自适应阈值，即根据呼吸频率的先验估计值设置阈值和差分间隔，呼吸频率的先验估计值是通过快速离散傅里叶变换求解得到；其次，算法在运行时通过上一步设置的差分间隔，对数据进行差分以判断当前是否到达波峰或波谷，但是是否计数峰值点取决于算法通过移动开窗法预测的“呼吸自适应步长”以及此时的累积置信度。通过计算步长可以预测下一次呼吸的动作完成时刻落点，结合步长进行差分运算能有效排除波形曲线上的小振动、毛刺、局部伪峰谷，在计算窗口到达终点时也可以为下一个计算窗口的起始提供较好的拟合。First, the traditional threshold setting method is changed to an adaptive threshold, that is, the threshold and differential interval are set according to the prior estimate of the respiratory frequency. The prior estimate of the respiratory frequency is obtained by fast discrete Fourier transform. Secondly, the algorithm uses the differential interval set in the previous step to differentiate the data during operation to determine whether the current peak or trough has been reached. However, whether to count the peak point depends on the "respiratory adaptive step" predicted by the algorithm through the moving window method and the cumulative confidence at this time. By calculating the step length, the landing point of the next breathing action completion time can be predicted. The differential operation combined with the step length can effectively eliminate small vibrations, burrs, and local pseudo-peaks and valleys on the waveform curve. When the calculation window reaches the end, it can also provide a better fit for the start of the next calculation window.
本实施例将呼吸过程中IMU解算得到的姿态信号波形作为数据源，认为信号波形具有局部最大值和局部最小值。对于呼吸波形曲线来说，本实施例定义最小有效数据阈值Ymin和最大有效数据阈值Ymax用于区分姿态变化幅度，姿态变化幅度小于Ymin时说明胸廓此时完全吸入气体或完全排出气体，运动平缓，人体呼吸动作处于完成状态；姿态变化幅度介于Ymin和Ymax之间时说明胸廓此时正在吸入气体或正在排出气体，运动稍显剧烈，人体呼吸动作处于切换阶段；姿态变化幅度超过Ymax时认为是异常状态，可能是传感器固定松动、前后传感器未对齐或受试者呼吸严重过速、严重过载等情况。This embodiment uses the posture signal waveform obtained by IMU solution during breathing as the data source, and considers that the signal waveform has a local maximum value and a local minimum value. For the breathing waveform curve, this embodiment defines the minimum valid data threshold value Ymin and the maximum valid data threshold value Ymax to distinguish the amplitude of posture change. When the amplitude of posture change is less than Ymin , it means that the chest is completely inhaling or expelling gas at this time, the movement is gentle, and the human body's breathing action is in a completed state; when the amplitude of posture change is between Ymin and Ymax , it means that the chest is inhaling or expelling gas at this time, the movement is slightly violent, and the human body's breathing action is in a switching stage; when the amplitude of posture change exceeds Ymax , it is considered to be an abnormal state, which may be due to loose sensor fixation, misalignment of the front and rear sensors, or severe over-breathing and severe overload of the subject.
基于步长的自适应差分法分为两个算法，具体算法内容如下所示：The step-size-based adaptive difference method is divided into two algorithms. The specific algorithm contents are as follows:
算法1：求呼吸波形曲线的下一个峰谷(峰谷即波峰或波谷)。Algorithm 1: Find the next peak or valley (peak or valley is the peak or valley) of the respiratory waveform curve.
输入：呼吸波形曲线Input: Respiration waveform curve
输出：下一个呼吸峰谷Output: Next breathing peak
1)设起始历元i＝1,设当前呼吸状态为ti，起始呼吸状态为t0，(呼吸状态包括动作切换状态和动作完成状态)。1) Let the starting epoch i=1, let the current breathing state beti , and the starting breathing state bet0 (the breathing state includes the action switching state and the action completion state).
2)若遍历结束，算法向上返回遍历结束状态，算法结束；2) If the traversal is completed, the algorithm returns to the traversal end state and the algorithm ends;
3)求第i个历元数据与第i+M个历元数据的差分值△x；/*历元差△x为正则处于波谷到波峰段，历元差为负则处于波峰到波谷段；M表示差分间隔，此处算法实际使用的是自适应差分间隔，见3.2节*/3) Calculate the difference value △x between the i-th epoch data and the i+M-th epoch data; /*If the epoch difference △x is positive, it is in the trough to peak section, and if the epoch difference is negative, it is in the peak to trough section; M represents the difference interval. The algorithm actually uses the adaptive difference interval here, see Section 3.2*/
4)如果△x超出Ymax，呼吸频率测量出现问题，抛出异常，算法结束；4) If △x exceeds Ymax , there is a problem with the respiratory rate measurement, an exception is thrown, and the algorithm ends;
5)如果△x满足Ymin<|△x|<Ymax，指示器记录当前呼吸处于峰谷变化段，ti为动作切换状态；转到步骤3；5) If △x satisfies Ymin <|△x|<Ymax , the indicator records that the current breathing is in the peak-to-valley change segment, and tiis the action switching state; go to step 3;
6)如果△x满足0<|△x|<Ymin，指示器记录当前呼吸接近峰谷，ti为动作完成状态；6) If △x satisfies 0<|△x|<Ymin , the indicator records that the current breathing is close to the peak and trough, and tiis the action completion state;
6.1)如果ti-1是动作完成状态，则重新向后寻找下一个峰谷，i＝i+1，转到2)；6.1) If ti-1 is the action completion state, then search for the next peak and valley again, i=i+1, and go to 2);
6.2)如果ti-1是动作切换状态，则认为当前位置是峰谷部分，输出呼吸峰谷位置，算法结束；6.2) If ti-1 is in the action switching state, the current position is considered to be the peak-valley part, the respiratory peak-valley position is output, and the algorithm ends;
算法2：使用基于步长的呼吸峰谷判断法计数呼吸窗口内呼吸曲线的呼吸波峰谷数。呼吸曲线分为起始不完整呼吸段(若存在)、中间的完整呼吸段、结尾不完整呼吸段(若存在)，记起始不完整呼吸的呼吸进度为n0(0<n0<1)，结尾不完整呼吸的呼吸进度为n1(0<n1<1)，如图3所示；则上一个呼吸窗口的结尾未完成呼吸的呼吸进度为1-n0。Algorithm 2: Use the step-based breathing peak-valley judgment method to count the number of breathing peaks and valleys of the breathing curve within the breathing window. The breathing curve is divided into the initial incomplete breathing segment (if any), the middle complete breathing segment, and the final incomplete breathing segment (if any). The breathing progress of the initial incomplete breathing is n0 (0<n0 <1), and the breathing progress of the final incomplete breathing is n1 (0<n1 <1), as shown in Figure 3; the breathing progress of the unfinished breathing at the end of the previous breathing window is 1-n0 .
输入：呼吸波形曲线Input: Respiration waveform curve
输出：呼吸波的峰谷数NOutput: Number of peaks and valleys of the respiratory wave N
1)拟合本呼吸窗口的起始呼吸峰谷数N＝N+n0；1) Fitting the number of initial respiratory peaks and valleys of this respiratory window N=N+n0 ;
2)调用算法1向后寻找下一个呼吸峰谷；2) Call Algorithm 1 to search for the next respiratory peak and valley;
2.1)如果算法1遍历完成整个呼吸窗口，转到步骤5)2.1) If Algorithm 1 traverses the entire breathing window, go to step 5)
2.2)如果算法1抛出异常，则本算法将该异常向上抛出，算法结束并退出；2.2) If Algorithm 1 throws an exception, this algorithm will throw the exception upward, and the algorithm ends and exits;
2.3)如果该峰谷为本窗口的第一个峰谷，则重新调用算法1向后寻找下一个峰谷；/*第一个峰谷以前的呼吸进度已拟合，本窗其他峰谷数从第一个呼吸峰谷开始向后计算*/2.3) If the peak valley is the first peak valley of this window, then re-call Algorithm 1 to search for the next peak valley; /*The breathing progress before the first peak valley has been fitted, and the number of other peak valleys in this window is calculated backward from the first breathing peak valley*/
2.4)如果算法1返回峰谷位置，由当前呼吸步长si与预测呼吸步长psi差分得到预测偏差△si(△si＝si-psi)和预测偏差范围2.4) If Algorithm 1 returns the peak and valley position, the prediction deviation △s i( △si = si - pi ) and the prediction deviation range are obtained by the difference between the current breathing step si and the predicted breathing step ps i
3)根据预测偏差范围Si和表1得到当前呼吸峰谷的置信度confi；3) Obtain the confidence level confi of the current respiratory peak and valley according to the prediction deviation range Si and Table 1;
4)将上一次呼吸峰谷的累积置信度∏confi-1与当前呼吸峰谷的置信度confi相乘得到当前呼吸峰谷的累积置信度∏confi；/*该步骤假设当前为呼吸峰谷，并根据呼吸步长偏差验证正确性*/；4) Multiply the cumulative confidence ∏confi-1 of the previous respiratory peak and valley by the confidence confi of the current respiratory peak and valley to obtain the cumulative confidence ∏confi of the current respiratory peak and valley; /*This step assumes that the current is a respiratory peak and valley, and verifies the correctness based on the respiratory step deviation*/;
4.1)如果∏confi>50％，认为当前确为呼吸峰谷段，N＝N+1，记录si，对本窗口的所有历史呼吸步长进行自适应指数加权得到下一次呼吸的预测呼吸步长psi，计算公式4.1) If ∏confi >50%, it is considered that the current breathing peak and valley segment is indeed, N = N + 1, recordsi , and perform adaptive exponential weighting on all historical breathing steps in this window to obtain the predicted breathing step lengthpsi of the next breath, and the calculation formula is
为for
计算完成后转到步骤2)；After the calculation is completed, go to step 2);
4.2)如果∏confi<50％，认为当前不在呼吸峰谷段，忽略当前呼吸峰谷，回退4.2) If ∏confi <50%, it is considered that the current breathing peak and valley are not in progress, the current breathing peak and valley are ignored, and the process is reversed.
∏confi为∏confi-1；∏confi is ∏confi-1 ;
4.2.1)若预测偏差△si小于0，转到步骤2；/*偏差小于0时，偏差范围从0.5往1变化，此时该次呼吸峰谷的置信度会增大，因此存在下一个峰谷使它的置信度大于50％，算法继续*/4.2.1) If the prediction deviation △si is less than 0, go to step 2; /* When the deviation is less than 0, the deviation range changes from 0.5 to 1. At this time, the confidence of the peak and valley of this breath will increase. Therefore, there is a next peak and valley that makes its confidence greater than 50%, and the algorithm continues*/
4.2.2)若预测偏差△si大于0，该次呼吸出现异常，蜂鸣提示用户呼吸存在异常，可能是呼吸暂停、低通气或其他异常，算法结束并退出；/*偏差大于0时，偏差范围从1往1.5变化，此时该次呼吸峰谷的置信度只会减小，由于累积置信度已经低于最低置信度阈值且不可能再上升，算法结束并退出，并报告异常原因*/4.2.2) If the prediction deviation △siis greater than 0, the breathing is abnormal, and the buzzer reminds the user that the breathing is abnormal, which may be apnea, hypopnea or other abnormalities. The algorithm ends and exits; /* When the deviation is greater than 0, the deviation range changes from 1 to 1.5. At this time, the confidence of the peak and valley of the breathing will only decrease. Since the cumulative confidence is already lower than the minimum confidence threshold and cannot rise again, the algorithm ends and exits, and reports the cause of the abnormality*/
5)拟合本呼吸窗口的结尾呼吸峰谷数N＝N+n15) Fit the number of end breathing peaks and valleys of this breathing window N = N + n1
6)输出呼吸峰谷数N,算法结束并退出。6) Output the number of respiratory peaks and valleys N, and the algorithm ends and exits.
算法1中，由于呼吸动作的切换和完成状态的姿态变化幅度取决于呼吸频率的大小，因此Ymin与Ymax使用的是自适应阈值来计算，见3.2节。算法2中呼吸进度等于不完全呼吸步长(即si)除以呼吸窗口的自适应指数加权平均呼吸步长(即当前预测步长psi)；In Algorithm 1, since the switching of breathing action and the change amplitude of posture in the completed state depend on the size of breathing frequency, Ymin and Ymax are calculated using adaptive thresholds, see Section 3.2. In Algorithm 2, the breathing progress is equal to the incomplete breathing step length (i.e., si ) divided by the adaptive exponential weighted average breathing step length of the breathing window (i.e., the current predicted step length pi );
算法1与算法2中，为了降低数据在区间波动时引起的计算冗余，在算法中所有求取差分值的地方均采取自适应差分间隔，即取其与前第M个数据的差值。In Algorithm 1 and Algorithm 2, in order to reduce the computational redundancy caused by data fluctuations in the interval, an adaptive differential interval is adopted in all places where differential values are obtained in the algorithm, that is, the difference between the differential value and the previous M-th data is taken.
3.3自适应阈值与自适应差分间隔3.3 Adaptive Threshold and Adaptive Differential Interval
本实施例从经过快速离散傅里叶变换得到的频谱图计算信号频率的先验估计值，并据此调整阈值与差分间隔，以提高基于步长的自适应差分法计算呼吸频率的准确度。This embodiment calculates a priori estimated value of the signal frequency from the spectrum diagram obtained by fast discrete Fourier transform, and adjusts the threshold and the difference interval accordingly to improve the accuracy of calculating the respiratory frequency by the step-based adaptive difference method.
本实施例采用基于移动开窗法的1024点快速离散傅里叶变换，首先计算经过姿态解算得到的各点数据的移动开窗平均值，计算公式如17所示。This embodiment adopts a 1024-point fast discrete Fourier transform based on the moving windowing method. First, the moving windowing average value of each point data obtained after posture solution is calculated. The calculation formula is shown in 17.
其中nL,nR为待计算的点两边的点数，cn为窗口大小，pi为姿态解算得到的数据，gi为经由移动开窗法计算得到的结果。此处计算窗长与自适应卡尔曼滤波处的窗长保持一致。WherenL ,nR are the number of points on both sides of the point to be calculated,cn is the window size,pi is the data obtained by attitude solution, andgi is the result calculated by the moving window method. The calculation window length here is consistent with the window length at the adaptive Kalman filter.
接下来对计算结果进行快速离散傅里叶变换，由于变换结果具有对称性，因此只需要保存512个频率值到频谱序列。其中，数据波形可能包含有低频或高频噪声，因此需要忽略频谱中的起始低频部分(第0-9频率值)和结尾高频部分(第502-511频率值)，而后取最大频值，并根据其索引经过式18计算得到信号频率估计值。Next, the calculation results are subjected to fast discrete Fourier transform. Since the transformation results are symmetrical, only 512 frequency values need to be saved in the spectrum sequence. Among them, the data waveform may contain low-frequency or high-frequency noise, so it is necessary to ignore the starting low-frequency part (0-9 frequency values) and the ending high-frequency part (502-511 frequency values) in the spectrum, and then take the maximum frequency value, and calculate the signal frequency estimate value according to its index through formula 18.
其中n是最大幅值的索引值，F是信号的实际采样频率，N是采集点数。Where n is the index value of the maximum amplitude, F is the actual sampling frequency of the signal, and N is the number of acquisition points.
根据香农定理，采样频率大于信号最高频率的两倍时可以完全还原原始信号，正常人的呼吸频率应在12-20次每分钟，通过控制呼吸深度，呼吸频率可以降低到6次每分钟；在剧烈运动后，呼吸频率可以上升到60次每分钟(或以上)，但理论最大呼吸频率一般不会超过120次/分钟。因此本实施例设置数据采集频率为5Hz，采集点数为1024点。为了简化计算，增强系统的实时测量能力，自适应间隔、自适应阈值与信号频率的对应关系使用历史呼吸数据经验进行设置，具体将呼吸分为如下四类：According to Shannon's theorem, the original signal can be completely restored when the sampling frequency is greater than twice the highest frequency of the signal. The normal breathing rate should be 12-20 times per minute. By controlling the depth of breathing, the breathing rate can be reduced to 6 times per minute; after strenuous exercise, the breathing rate can rise to 60 times per minute (or more), but the theoretical maximum breathing rate generally does not exceed 120 times per minute. Therefore, this embodiment sets the data acquisition frequency to 5Hz and the number of acquisition points to 1024 points. In order to simplify the calculation and enhance the real-time measurement capability of the system, the corresponding relationship between the adaptive interval, the adaptive threshold and the signal frequency is set using historical breathing data experience, and breathing is specifically divided into the following four categories:
1)正常呼吸，呼吸频率fn≤201) Normal breathing, respiratory rate fn ≤ 20
2)慢速运动呼吸，呼吸频率21≤fn≤302) Slow exercise breathing, respiratory rate 21 ≤f n ≤ 30
3)快速运动呼吸，呼吸频率31≤fn≤603) Rapid exercise breathing, respiratory rate 31≤fn ≤60
4)过速呼吸，呼吸频率fn>604) Rapid breathing, respiratory rate fn >60
其中，最小有效阈值、最大有效阈值分别对应的设置为各类呼吸场景中平均姿态数据变化的1％、10％；差分间隔设置为各类呼吸场景中呼吸接近峰谷时历元差大于最小有效阈值但小于最大有效阈值的平均历元间隔。Among them, the minimum effective threshold and the maximum effective threshold are set to 1% and 10% of the average posture data change in various breathing scenarios respectively; the differential interval is set to the average epoch interval when the epoch difference is greater than the minimum effective threshold but less than the maximum effective threshold when the breathing is close to the peak and valley in various breathing scenarios.
根据以上思想，本实施例中信号频率估计值与自适应参数的对应关系如表2所示。表2信号频率估计值与自适应参数对应关系According to the above ideas, the corresponding relationship between the signal frequency estimation value and the adaptive parameter in this embodiment is shown in Table 2. Table 2 Corresponding relationship between the signal frequency estimation value and the adaptive parameter
其中，呼吸频率超过60次时，呼吸严重过速，平均呼吸深度水平会显著下降，但是剧烈运动中测量误差会增大，因此自适应阈值不宜作过多回调。Among them, when the respiratory rate exceeds 60 times, the breathing is seriously too fast and the average breathing depth level will drop significantly. However, the measurement error will increase during strenuous exercise, so the adaptive threshold should not be adjusted too much.
4、仿真实验与结果分析4. Simulation experiment and result analysis
4.1实验设备4.1 Experimental equipment
仿真实验选择GY-9250型IMU作为呼吸运动检测的数据源，其集成了三轴陀螺仪、三轴加速度计以及三轴磁力计。GY-9250重量仅0.01kg，数据传输速率较高，适合进行仿真实验。The simulation experiment selected the GY-9250 IMU as the data source for respiratory motion detection, which integrates a three-axis gyroscope, a three-axis accelerometer, and a three-axis magnetometer. The GY-9250 weighs only 0.01kg and has a high data transmission rate, making it suitable for simulation experiments.
4.1实验方法4.1 Experimental methods
仿真实验首先采集三轴加速度值和陀螺仪角速度值，而后将其传递给Cortex-M7处理器核心进行计算。仿真实验使用卡尔曼滤波(下称KF)、Sage-Husa自适应卡尔曼滤波(下称AKF)以及本实施例提出的改进的自适应卡尔曼滤波(下称MAKF)三种方法同时进行解算，然后分别对静坐、正常走路、跑步等三个场景的呼吸运动曲线进行呼吸频率计算。相关参数设定如下：GY-9250工作频率为50Hz、GY-9250低通滤波频率为25Hz、数据记录频率(实际采样频率)为5Hz、采样点数为1024点。The simulation experiment first collects the three-axis acceleration value and the gyroscope angular velocity value, and then passes them to the Cortex-M7 processor core for calculation. The simulation experiment uses three methods, Kalman filter (hereinafter referred to as KF), Sage-Husa adaptive Kalman filter (hereinafter referred to as AKF) and the improved adaptive Kalman filter (hereinafter referred to as MAKF) proposed in this embodiment to solve at the same time, and then calculates the respiratory frequency of the respiratory motion curves of three scenes such as sitting still, normal walking and running. The relevant parameters are set as follows: GY-9250 operating frequency is 50Hz, GY-9250 low-pass filter frequency is 25Hz, data recording frequency (actual sampling frequency) is 5Hz, and the number of sampling points is 1024 points.
4.3数据分析4.3 Data Analysis
本实施例就静坐、正常走路和跑步后三种呼吸场景进行了实验。This embodiment conducts experiments on three breathing scenarios: sitting quietly, walking normally, and after running.
1)静坐深呼吸场景：静坐深呼吸时，KF、AKF、MAKF解算的呼吸曲线如图3所示。1) Sitting quietly and taking deep breath: When sitting quietly and taking deep breath, the breathing curves solved by KF, AKF and MAKF are shown in Figure 3.
总体上看，深呼吸时呼吸频率较为缓慢，此时三种解算方法得出的呼吸曲线差异较为明显：在低频率呼吸时，由于KF的噪声阵与环境不匹配，因此解算产生的漂移较大；AKF虽然能够绘制大致的波形，但受噪声过补偿影响，有一部分呼吸曲线不明显；而MAKF则能够绘制较为规律的呼吸曲线。In general, the breathing frequency is relatively slow during deep breathing, and the breathing curves obtained by the three solution methods are quite different: during low-frequency breathing, the drift caused by the solution is relatively large because the noise array of KF does not match the environment; although AKF can draw a rough waveform, some breathing curves are not obvious due to the influence of noise over-compensation; while MAKF can draw a more regular breathing curve.
利用本实施例提出的基于步长的自适应差分法分别对三个呼吸曲线进行计算，求出的呼吸频率与误差其结果如表3所示。The three respiratory curves are calculated respectively using the step-size-based adaptive difference method proposed in this embodiment, and the respiratory frequency and error obtained are shown in Table 3.
表3静坐时深呼吸频率与误差Table 3 Deep breathing frequency and error during meditation
2)走路场景：正常走路时，KF、AKF、MAKF解算的呼吸曲线如图4所示。2) Walking scene: When walking normally, the breathing curves solved by KF, AKF, and MAKF are shown in Figure 4.
走路时呼吸频率在正常区间，KF和AKF解算得到的曲线存在一些毛刺或“异常”(见图中红圈部分)，而本实施例提出的MAKF得到的结果较为平滑。计算所得的呼吸频率与误差结果如表4所示。When walking, the breathing frequency is in the normal range. The curves obtained by KF and AKF have some burrs or "abnormalities" (see the red circle in the figure), while the results obtained by MAKF proposed in this embodiment are relatively smooth. The calculated breathing frequency and error results are shown in Table 4.
表4走路时呼吸频率与误差Table 4 Respiratory frequency and error during walking
3)运动场景：跑步时，KF和MAKF解算的呼吸曲线如图5所示。3) Sports scene: When running, the breathing curves solved by KF and MAKF are shown in Figure 5.
根据呼吸频率的相关研究结果，呼吸频率在运动时会随时间近似线性增加，并在运动结束时接近最大值。因此本实施例取跑步接近结束时的数据，以测试高频呼吸下算法的效果。According to the relevant research results of respiratory frequency, the respiratory frequency will increase approximately linearly with time during exercise and approach the maximum value at the end of exercise. Therefore, this embodiment takes the data near the end of running to test the effect of the algorithm under high-frequency breathing.
由于跑步时呼吸频率较快，AKF会因为残差失配中止计算，因此无法绘制呼吸曲线。而在此情景下，本实施例提出的简化思路和自适应补偿因子保证了MAKF不会产生残差失配现象。利用本实施例提出的基于步长的自适应差分法分别对两个呼吸曲线进行计算，求出的呼吸频率与误差其结果如表5所示。Since the breathing rate is faster during running, AKF will terminate the calculation due to residual mismatch, so it is impossible to draw the breathing curve. In this scenario, the simplified idea and adaptive compensation factor proposed in this embodiment ensure that MAKF will not produce residual mismatch. The two breathing curves are calculated respectively using the step-based adaptive difference method proposed in this embodiment, and the breathing frequency and error are shown in Table 5.
表5跑步时呼吸频率与误差Table 5 Respiratory frequency and error during running
由上述三部分内容可以看出，相比KF和AKF来说，本实施例提出的MAKF与基于步长的自适应差分法相结合计算呼吸频率的方法更加有效，误差始终保持在3％及以下，能够较为准确地实现静态、动态等场景下的呼吸状态检测。It can be seen from the above three parts that compared with KF and AKF, the method of calculating the respiratory frequency by combining MAKF with the step-based adaptive difference method proposed in this embodiment is more effective, and the error is always kept at 3% or less, which can more accurately realize the respiratory state detection in static and dynamic scenarios.
对不同运动场景分别测试3次呼吸频率，并用均方根误差法评估结果，如表6所示。准确性测试中，每次测量的呼吸窗长为200秒，呼吸次数即在该呼吸窗长中的呼吸周期数；小数为呼吸窗口前后两端未能构成一个完整呼吸次数的拟合值；准确性测试按呼吸频率递增的顺序排列，其中，1，2，3为静坐场景；4，5，6为走路场景；7，8，9为跑步场景。The breathing frequency was tested three times for different sports scenes, and the results were evaluated by the root mean square error method, as shown in Table 6. In the accuracy test, the breathing window length for each measurement was 200 seconds, and the number of breaths was the number of breathing cycles in the breathing window length; the decimal was the fitting value that failed to form a complete breathing number at the front and back ends of the breathing window; the accuracy test was arranged in the order of increasing breathing frequency, where 1, 2, and 3 were sitting scenes; 4, 5, and 6 were walking scenes; and 7, 8, and 9 were running scenes.
表6呼吸频率准确性测试Table 6 Respiratory rate accuracy test
由准确性测试可知，KF在场景变换情况下均方根误差较大，AKF在呼吸频率大于16后会因残差失配中止计算，而MAKF对呼吸频率的检测具有较高准确度，在各种呼吸场景下保持了较好的稳定性，其呼吸频率均方根误差小于1，证明了MAKF和基于步长的自适应差分算法的有效性和准确性。From the accuracy test, it can be seen that the root mean square error of KF is large when the scene changes, and AKF will terminate the calculation due to residual mismatch when the breathing frequency is greater than 16. MAKF has a high accuracy in detecting the breathing frequency and maintains good stability in various breathing scenarios. Its root mean square error of breathing frequency is less than 1, which proves the effectiveness and accuracy of MAKF and the step-based adaptive difference algorithm.
以上所述，仅是本发明的较佳实施例而已，并非是对本发明作其它形式的限制，任何熟悉本专业的技术人员可能利用上述揭示的技术内容加以变更或改型为等同变化的等效实施例应用于其它领域，但是凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与改型，仍属于本发明技术方案的保护范围。The above description is only a preferred embodiment of the present invention and does not limit the present invention in other forms. Any technician familiar with the profession may use the technical content disclosed above to change or modify it into an equivalent embodiment with equivalent changes and apply it to other fields. However, any simple modification, equivalent change and modification made to the above embodiment based on the technical essence of the present invention without departing from the content of the technical solution of the present invention still falls within the protection scope of the technical solution of the present invention.
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311244889.3ACN117958792B (en) | 2023-09-26 | 2023-09-26 | Self-adaptive Kalman filtering respiratory frequency calculation method based on mobile windowing method |
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311244889.3ACN117958792B (en) | 2023-09-26 | 2023-09-26 | Self-adaptive Kalman filtering respiratory frequency calculation method based on mobile windowing method |
| Publication Number | Publication Date |
|---|---|
| CN117958792A CN117958792A (en) | 2024-05-03 |
| CN117958792Btrue CN117958792B (en) | 2024-07-23 |
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311244889.3AActiveCN117958792B (en) | 2023-09-26 | 2023-09-26 | Self-adaptive Kalman filtering respiratory frequency calculation method based on mobile windowing method |
| Country | Link |
|---|---|
| CN (1) | CN117958792B (en) |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN119469741B (en)* | 2025-01-16 | 2025-04-18 | 青岛理工大学 | Concealed conduit detection method based on track reconstruction and storage medium |
| CN120477940B (en)* | 2025-07-21 | 2025-09-23 | 西安交通大学医学院第一附属医院 | Real-time image registration and path planning system in bronchoscope operation navigation |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108784703A (en)* | 2018-07-05 | 2018-11-13 | 西南石油大学 | A kind of wearable monitoring of respiration method of the middle-aged and the old |
| CN110132271A (en)* | 2019-01-02 | 2019-08-16 | 中国船舶重工集团公司第七0七研究所 | A kind of adaptive Kalman filter Attitude estimation algorithm |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN100444785C (en)* | 2006-11-10 | 2008-12-24 | 南京航空航天大学 | Air cushion type automatic monitoring method for human respiration, heartbeat and turning over |
| JP2015217143A (en)* | 2014-05-19 | 2015-12-07 | パナソニックＩｐマネジメント株式会社 | Heart rate measurement device |
| JP2016013239A (en)* | 2014-07-01 | 2016-01-28 | 国立大学法人 筑波大学 | Cardiopulmonary function monitoring device |
| CN114545355A (en)* | 2015-04-20 | 2022-05-27 | 瑞思迈传感器技术有限公司 | Detection and identification of humans from characteristic signals |
| US11338158B2 (en)* | 2018-03-15 | 2022-05-24 | Safran Aerotechnics Sas | System and a method for delivering breathing gas to passengers on-board an aircraft |
| CN110327036B (en)* | 2019-07-24 | 2021-11-30 | 东南大学 | Method for extracting respiratory signal and respiratory frequency from wearable electrocardiogram |
| CN110811572B (en)* | 2019-10-18 | 2020-10-27 | 西安交通大学 | Simulation synthesis method and device for photoelectric volume wave signal |
| CN115429251A (en)* | 2021-06-03 | 2022-12-06 | 安徽华米健康科技有限公司 | Wearable device and its monitoring method and monitoring device |
| US11954862B2 (en)* | 2021-09-20 | 2024-04-09 | Nvidia Corporation | Joint estimation of heart rate and respiratory rate using neural networks |
| CN114983373B (en)* | 2022-06-02 | 2023-03-28 | 谢俊 | Method for detecting human heart rate |
| CN116269298B (en)* | 2023-02-21 | 2023-11-10 | 郑州大学 | Non-contact sleep breathing monitoring method and system based on millimeter wave radar |
| CN116350206A (en)* | 2023-04-12 | 2023-06-30 | 南京大学 | Method for deriving respiration waves from cardiac-related signals and for calculating respiration rate from the respiration waves |
| CN116651534B (en)* | 2023-07-31 | 2023-10-20 | 延边大学 | Separated noise-reducing and vibration-absorbing optical platform for patch clamp experiment |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108784703A (en)* | 2018-07-05 | 2018-11-13 | 西南石油大学 | A kind of wearable monitoring of respiration method of the middle-aged and the old |
| CN110132271A (en)* | 2019-01-02 | 2019-08-16 | 中国船舶重工集团公司第七0七研究所 | A kind of adaptive Kalman filter Attitude estimation algorithm |
| Publication number | Publication date |
|---|---|
| CN117958792A (en) | 2024-05-03 |
| Publication | Publication Date | Title |
|---|---|---|
| CN117958792B (en) | Self-adaptive Kalman filtering respiratory frequency calculation method based on mobile windowing method | |
| JP6555692B2 (en) | Method for measuring respiration rate and system for measuring respiration rate | |
| US20210219925A1 (en) | Apparatus and method for detection of physiological events | |
| US10758164B2 (en) | Processing apparatus and processing method for determining a respiratory signal of a subject | |
| CN108478206B (en) | Heart rate monitoring method based on pulse wave in exercise state | |
| Aly et al. | Zephyr: Ubiquitous accurate multi-sensor fusion-based respiratory rate estimation using smartphones | |
| US7634379B2 (en) | Newtonian physical activity monitor | |
| RU2515404C2 (en) | Respiration monitors and methods of respiration monitoring | |
| US8614630B2 (en) | Fall detection using sensor fusion | |
| CN100362963C (en) | Portable health-care monitoring device capable of performing motion compensation and compensation method thereof | |
| US20130144130A1 (en) | System method and device for monitoring a person's vital signs | |
| CN105662345B (en) | heartbeat signal processing method, device and system | |
| US9572521B2 (en) | Monitoring biometric characteristics of a user of a user monitoring apparatus | |
| EP3651643A1 (en) | A wearable device for the continuous monitoring of the respiratory rate | |
| EP2563218A1 (en) | Method, apparatus, computer program and system for measuring oscillatory motion | |
| CN115429251A (en) | Wearable device and its monitoring method and monitoring device | |
| CN108992053A (en) | A method of real-time chainless detection heart rate and eartbeat interval | |
| EP4216812A1 (en) | A processor and method for determining a respiratory signal | |
| Erfianto et al. | IMU‐Based respiratory signal processing using cascade complementary filter method | |
| CN105852826B (en) | The method that terminal and terminal determine physiologic information | |
| CN113551687B (en) | Step counting method, device, step counting device, computer storage medium and chip | |
| CN115227225A (en) | Heart rate monitoring method, device and system based on BCG signal and intelligent mattress | |
| Li et al. | A Respiratory Rate Detection Scheme Based on Improved Adaptive Kalman Filtering and Gait Matching | |
| László et al. | Extracting Physiological Signals From Smartphone Sensors. | |
| CN116965800A (en) | Respiratory state evaluation method based on electrocardiographic data |
| 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 | ||
| CB03 | Change of inventor or designer information | Inventor after:Jin Hua Inventor after:Li Qinghan Inventor before:Li Qinghan Inventor before:Jin Hua | |
| CB03 | Change of inventor or designer information |