CN117414117A - Heart rate calculation method, device, equipment and medium - Google Patents
Heart rate calculation method, device, equipment and medium Download PDFInfo
- Publication number
- CN117414117A CN117414117A CN202311369772.8A CN202311369772A CN117414117A CN 117414117 A CN117414117 A CN 117414117A CN 202311369772 A CN202311369772 A CN 202311369772A CN 117414117 A CN117414117 A CN 117414117A
- Authority
- CN
- China
- Prior art keywords
- wrist
- heart rate
- data
- angle
- ppg signal
- 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.)
- Pending
Links
- 238000004364 calculation method Methods 0.000 title claims abstract description 95
- 210000000707 wrist Anatomy 0.000 claims abstract description 160
- 230000033001 locomotion Effects 0.000 claims abstract description 104
- 230000001133 acceleration Effects 0.000 claims abstract description 36
- 238000000034 method Methods 0.000 claims abstract description 36
- 238000001914 filtration Methods 0.000 claims abstract description 33
- 238000005259 measurement Methods 0.000 claims abstract description 17
- 238000010183 spectrum analysis Methods 0.000 claims abstract description 16
- 238000010586 diagram Methods 0.000 claims description 43
- 239000011159 matrix material Substances 0.000 claims description 34
- 238000012545 processing Methods 0.000 claims description 20
- 238000010606 normalization Methods 0.000 claims description 10
- 230000004927 fusion Effects 0.000 claims description 8
- 238000004590 computer program Methods 0.000 claims description 7
- 238000013213 extrapolation Methods 0.000 claims description 7
- 238000000354 decomposition reaction Methods 0.000 claims description 5
- 238000001514 detection method Methods 0.000 abstract description 10
- 230000000694 effects Effects 0.000 abstract description 2
- 238000013186 photoplethysmography Methods 0.000 description 111
- 238000004422 calculation algorithm Methods 0.000 description 21
- 239000013598 vector Substances 0.000 description 15
- 230000000875 corresponding effect Effects 0.000 description 10
- 230000008569 process Effects 0.000 description 6
- 230000003044 adaptive effect Effects 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 5
- 238000000605 extraction Methods 0.000 description 5
- 230000006870 function Effects 0.000 description 5
- 230000008859 change Effects 0.000 description 4
- 238000006243 chemical reaction Methods 0.000 description 4
- 230000005484 gravity Effects 0.000 description 4
- 238000012544 monitoring process Methods 0.000 description 4
- 238000012935 Averaging Methods 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 238000005370 electroosmosis Methods 0.000 description 2
- 239000013388 element-organic framework Substances 0.000 description 2
- 230000001121 heart beat frequency Effects 0.000 description 2
- 238000012880 independent component analysis Methods 0.000 description 2
- 230000033764 rhythmic process Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 238000013473 artificial intelligence Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000017531 blood circulation Effects 0.000 description 1
- 210000004204 blood vessel Anatomy 0.000 description 1
- CLOMYZFHNHFSIQ-UHFFFAOYSA-N clonixin Chemical compound CC1=C(Cl)C=CC=C1NC1=NC=CC=C1C(O)=O CLOMYZFHNHFSIQ-UHFFFAOYSA-N 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000000295 complement effect Effects 0.000 description 1
- 230000008602 contraction Effects 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000013178 mathematical model Methods 0.000 description 1
- 230000004089 microcirculation Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 239000013307 optical fiber Substances 0.000 description 1
- 238000000053 physical method Methods 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000012552 review Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 239000002023 wood Substances 0.000 description 1
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
- A61B5/02416—Detecting, measuring or recording pulse rate or heart rate using photoplethysmograph signals, e.g. generated by infrared radiation
-
- 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/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
- A61B5/7207—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal of noise induced by motion artifacts
-
- 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
-
- 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/7253—Details of waveform analysis characterised by using transforms
- A61B5/7257—Details of waveform analysis characterised by using transforms using Fourier transforms
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Signal Processing (AREA)
- Physics & Mathematics (AREA)
- Medical Informatics (AREA)
- General Health & Medical Sciences (AREA)
- Biophysics (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Physiology (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- Pathology (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Psychiatry (AREA)
- Cardiology (AREA)
- Mathematical Physics (AREA)
- Measurement Of The Respiration, Hearing Ability, Form, And Blood Characteristics Of Living Organisms (AREA)
Abstract
The application relates to a heart rate calculation method, a device, equipment and a medium, and relates to the technical field of heart rate detection, wherein the method comprises the steps of fusing acceleration data and first angular velocity data of a wrist of a target user under a local horizontal coordinate system to obtain first angular data of the wrist under an IMU coordinate system; judging whether the pitch angle in the first angle data is smaller than a first threshold value or larger than a second threshold value; the first threshold is less than the second threshold; if yes, determining the baseline heart rate as a heart rate calculation result of the target user; if not, judging whether the pulse signal in the photoelectric measurement pulse wave signal PPG signal and the motion artifact of the wrist are mutually independent; if yes, filtering the PPG signal to obtain a pulse signal, and obtaining a heart rate calculation result of the target user according to the pulse signal; if not, carrying out singular spectrum analysis on the PPG signal to obtain a heart rate calculation result of the target user. The heart rate calculation method and device have the effect of improving heart rate calculation accuracy.
Description
Technical Field
The present disclosure relates to the field of heart rate detection technologies, and in particular, to a heart rate calculation method, device, apparatus, and medium.
Background
The Photoplethysmography (PPG) signal is an optical signal reflecting the heart beat and microcirculation of the human body, and is obtained by placing a light source and a photosensitive detector on the skin and measuring the change in light intensity after the light source is reflected or transmitted by the skin tissue. When the heart contracts, the blood flow of skin tissues also changes due to the expansion and contraction of arterial blood vessels, so that the amplitude and waveform of the PPG signal change, and therefore, the PPG signal can be used for monitoring physiological indexes such as heart rate and heart rhythm, and the PPG signal is widely applied to the fields such as medical treatment, health monitoring and exercise monitoring.
In PPG signals, motion artifacts are a common noise source, which may cause errors in heart rate extraction, and in order to remove the motion artifacts, researchers have proposed various algorithms, such as a method based on peak detection, a method based on adaptive filtering, a method based on independent component analysis, and the like. Some algorithms are only suitable for specific motion patterns and motion intensities, and may not accurately remove motion artifacts for different motion types and motion intensities, resulting in large errors in heart rate calculation.
Disclosure of Invention
In order to solve the problem of large heart rate calculation error caused by motion artifact, the application provides a heart rate calculation method, device, equipment and medium.
In a first aspect, the present application provides a heart rate calculation method, which adopts the following technical scheme:
a heart rate calculation method, comprising:
according to the acceleration data and the first angular velocity data of the wrist of the target user in a local horizontal coordinate system, fusion is carried out, and first angle data of the wrist in an Inertial Measurement Unit (IMU) coordinate system is obtained; the first angle data comprises a pitch angle of rotation of the wrist around a Y axis under an IMU coordinate system;
judging whether the pitch angle in the first angle data is smaller than a first threshold value or larger than a second threshold value; the first threshold is less than the second threshold;
if yes, determining a baseline heart rate as a heart rate calculation result of the target user; if not, judging whether the pulse signal in the photoelectric measurement pulse wave PPG signal and the motion artifact of the wrist are mutually independent;
if yes, filtering the PPG signal to obtain the pulse signal, and obtaining a heart rate calculation result of the target user according to the pulse signal; if not, carrying out singular spectrum analysis on the PPG signal to obtain a heart rate calculation result of the target user.
By adopting the technical scheme, the obtained pitch angle is more accurate by fusing the acceleration data and the first angular velocity data, so that the accuracy of a subsequent judgment result based on the pitch angle is ensured, and the accuracy of heart rate calculation is further improved. When the pitch angle is smaller than the first threshold value or larger than the second threshold value, the baseline heart rate is directly determined as a heart rate calculation result of the target user, and therefore heart rate calculation efficiency can be improved. When the pitch angle is larger than or equal to the first threshold value and smaller than or equal to the second threshold value, judging whether the pulse signal in the PPG signal and the motion artifact of the wrist are mutually independent, if so, calculating the heart rate by adopting a filtering algorithm, and if not, calculating the heart rate by adopting a singular spectrum analysis technology. Different heart rate calculation methods are provided aiming at different characteristics of motion artifacts, so that the accuracy of heart rate calculation can be improved while the efficiency of heart rate calculation is improved.
Optionally, the fusing is performed according to acceleration data and first angular velocity data of the wrist of the target user in a local horizontal coordinate system, to obtain first angle data of the wrist in an IMU coordinate system, including:
acquiring second angle data of the wrist under an IMU coordinate system according to the acceleration data; the second angle data comprises a roll angle of the wrist rotating around an X axis and a pitch angle of the wrist rotating around a Y axis under an IMU coordinate system;
Obtaining second angular velocity data of the wrist under the IMU coordinate system according to the first angular velocity data and the second angular velocity data; the first angular velocity data comprises angular velocities at which the wrist rotates about an X axis, a Y axis and a Z axis respectively in a local horizontal coordinate system; the second angular velocity data comprises angular velocities at which the wrist rotates about an X axis, a Y axis and a Z axis respectively in an IMU coordinate system;
fusing the second angular velocity data and the second angle data to obtain the first angle data; the first angle data also includes a roll angle at which the wrist rotates about an X-axis in an IMU coordinate system.
Optionally, fusing the second angular velocity data and the second angle data to obtain the first angle data, including:
the first angle data is obtained through a Kalman filtering equation, and the calculation formula is as follows:
wherein k represents time, x k,k For the state at time k after filtering, i.e. the first angle data, r k Is the roll angle at time k, i.e. the roll angle p in the first angle data k Is the pitch angle at time k, namely the pitch angle x in the first angle data k，k-1 The state at time k, which is derived from the state at time k-1, is represented by the state extrapolation equation:
r k+1 Roll angle at time k+1, p k+1 Is the pitch angle at time k+1, [ r ] k ,p k ]The initial value is set to 0,0]The method comprises the steps of carrying out a first treatment on the surface of the dr/dt is the angular velocity of rotation around the X axis in the second angular velocity data, dp/dt is the angular velocity of rotation around the Y axis in the second angular velocity data; k (K) k Is a Kalman filtering matrix; z is Z k For the observation equation, the formula is:r acc for the roll angle, p, in the second angle data acc And the pitch angle in the second angle data.
By adopting the technical scheme, the Kalman filter equation is designed, so that the filtered roll angle and pitch angle are more accurate, the subsequent judgment result based on the pitch angle is also more accurate, and the accuracy of the subsequent heart rate calculation result is further improved.
Optionally, determining whether the pulse signal in the PPG signal and the motion artifact of the wrist are independent of each other includes:
performing discrete Fourier transform and normalization processing on the PPG signal to obtain a normalized waveform diagram of the PPG signal;
performing discrete Fourier transform and normalization processing on the waveform diagram of the wrist movement degree to obtain a normalized waveform diagram of the wrist movement degree; the oscillogram of the wrist motion degree is used for indicating the motion amplitude of the wrist at each moment;
And judging whether the pulse signal in the PPG signal and the motion artifact of the wrist are mutually independent or not according to the normalized waveform diagram of the PPG signal and the normalized waveform diagram of the wrist motion degree.
By adopting the technical scheme, whether the pulse signal in the PPG signal is independent of the motion artifact of the wrist can be accurately judged according to the normalized waveform diagram of the PPG signal and the normalized waveform diagram of the wrist motion degree through discrete Fourier transform and normalization processing, so that the accuracy of a subsequent heart rate calculation result is improved.
Optionally, determining whether the pulse signal in the PPG signal and the motion artifact of the wrist are independent of each other includes:
determining a first heart rate value corresponding to a frequency component with a maximum peak value from a normalized waveform diagram of the PPG signal, and determining a second heart rate value corresponding to the frequency component with the maximum peak value from a normalized waveform diagram of the wrist movement degree;
judging whether the absolute value of the difference value between the first heart rate value and the second heart rate value is greater than or equal to a third threshold value;
If yes, determining that a pulse signal in the PPG signal and the motion artifact of the wrist are mutually independent; if not, determining that the pulse signal in the PPG signal and the motion artifact of the wrist are not mutually independent.
By adopting the technical scheme, whether the pulse signal in the PPG signal is independent of the motion artifact of the wrist can be accurately judged according to the first heart rate value corresponding to the frequency component with the largest peak value in the normalized waveform diagram of the PPG signal and the second heart rate value corresponding to the frequency component with the largest peak value in the normalized waveform diagram of the wrist motion degree, so that the accuracy of the subsequent heart rate calculation result is improved.
Optionally, the waveform of the wrist motion degree is obtained through the following steps:
according to the second angular velocity data, calculating the square sum of the angular velocities of the wrist rotating around an X axis, a Y axis and a Z axis respectively under an IMU coordinate system according to the second angular velocity data;
determining the sum of squares as a motion amplitude value of the wrist at each moment;
and fitting to obtain a waveform diagram of the wrist motion degree according to the motion amplitude value of the wrist at each moment.
Optionally, performing singular spectrum analysis on the PPG signal to obtain a heart rate calculation result of the target user, including:
Converting the normalized waveform diagram of the PPG signal into a track matrix, and carrying out singular value decomposition on the track matrix to obtain different mode groups of the track matrix;
performing anti-diagonal equalization processing on the different mode groups to obtain a plurality of mode sequences with the same length;
and according to the mode sequences with the same length, obtaining a reconstructed PPG signal, performing discrete Fourier transform on the reconstructed PPG signal, and determining the main frequency of the transformed PPG signal as a heart rate calculation result of the target user.
By adopting the technical scheme, the singular spectrum analysis technology is adopted to accurately remove the motion artifact components from the PPG signal, and a more accurate heart rate calculation result can be obtained based on the reconstructed PPG signal.
In a second aspect, the present application provides a heart rate calculating device, which adopts the following technical scheme:
a heart rate computing device, comprising:
an obtaining module for: acquiring first angle data of the wrist under an Inertial Measurement Unit (IMU) coordinate system according to acceleration data and first angular velocity data of the wrist of a target user under a local horizontal coordinate system; the first angle data comprises a pitch angle of rotation of the wrist around a Y axis under an IMU coordinate system;
The judging module is used for: judging whether the pitch angle in the first angle data is smaller than a first threshold value or larger than a second threshold value; the first threshold is less than the second threshold;
a heart rate calculation module for: if the pitch angle in the first angle data is smaller than a first threshold value or larger than a second threshold value, determining a baseline heart rate as a heart rate calculation result of the target user;
the judging module is further configured to: if the pitch angle in the first angle data is larger than or equal to a first threshold value and smaller than or equal to a second threshold value, judging whether a pulse signal in a photoelectric measurement pulse wave PPG signal and the motion artifact of the wrist are mutually independent;
the heart rate calculation module is further configured to filter the PPG signal if the pulse signal in the PPG signal and the motion artifact of the wrist are independent of each other, obtain the pulse signal, and obtain a heart rate calculation result of the target user according to the pulse signal; if the pulse signal in the PPG signal and the motion artifact of the wrist are not mutually independent, performing singular spectrum analysis on the PPG signal to obtain a heart rate calculation result of the target user.
In a third aspect, the present application provides a computer device, which adopts the following technical scheme:
a computer device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, the processor implementing the method of the first aspect when executing the program.
In a fourth aspect, the present application provides a computer readable storage medium, which adopts the following technical scheme:
a computer readable storage medium storing a computer program capable of being loaded by a processor and executing any one of the methods of the first aspect.
In summary, the present application includes at least one of the following beneficial technical effects:
different heart rate calculation methods are provided according to different characteristics of pitch angle and motion artifacts, when the pitch angle is too small or too large, a baseline heart rate is determined to be a heart rate calculation result of a target user, when the pitch angle is in a proper range, the characteristics of the motion artifacts are considered, when the motion artifacts are uncorrelated with heart beat frequency, a filtering algorithm is adopted to calculate the heart rate, and when the motion artifacts are correlated with the heart beat frequency, a singular spectrum analysis technology is adopted to calculate the heart rate. Different heart rate calculation methods are adopted according to actual conditions, so that the accuracy of heart rate calculation can be improved while the efficiency of heart rate calculation is improved. And the acceleration data and the first angular velocity data of the wrist of the target user under the local horizontal coordinate system are fused, so that the accuracy of the filtered pitch angle is ensured, and the accuracy of the subsequent heart rate calculation is further improved.
Drawings
Fig. 1 is a flowchart of a heart rate calculation method according to an embodiment of the present application.
Fig. 2 is a block diagram of a heart rate calculating device according to an embodiment of the present application.
Reference numerals illustrate: 201. obtaining a module; 202. a judging module; 203. a heart rate calculation module; 204. a processing module; 205 determine a module.
Detailed Description
In order to make the objects, technical solutions and advantages of the present application more apparent, the present application will be further described in detail with reference to fig. 1-2 and the embodiments. It should be understood that the specific embodiments described herein are for purposes of illustration only and are not intended to limit the present application.
The terms related to the embodiments of the present application will be explained first.
1. MPU6050 is a commonly used small inertial measurement unit (Inertial Measurement Unit, IMU) chip, manufactured by Invensense. The chip integrates a triaxial accelerometer and a triaxial gyroscope, can measure linear acceleration and angular velocity of an object, and is widely applied to the fields of robots, aircrafts, intelligent bracelets and the like. The MPU6050 supports an I2C interface, can communicate with a main controller through an I2C protocol, is internally integrated with a digital signal processor, and can process and filter acquired data, so that the accuracy and stability of the data are improved.
2. Singular spectrum analysis (Singular Spectrum Analysis, SSA), SAS decomposes the signal into some basic modal functions (Empirical Orthogonal Functions, EOFs) and corresponding time series (Principal Component Time Series, tfts) by singular value decomposition. EOFs represent different frequency, amplitude and phase information in the signal, while tfts represent the variation in time of each fundamental mode function.
3. Kalman Filter (KF), KF comprises two main steps: and predicting and updating, wherein the predicting step predicts the state at the next moment by using information such as a state equation, control input, a state transition matrix and the like of the system. And in the updating step, the prediction result is corrected by utilizing information such as observation data, an observation matrix and the like so as to obtain the calculation of the current state. In the updating step, the Kalman filtering determines the trade-off ratio between the predicted result and the observed result by calculating the Kalman gain, thereby obtaining the optimal calculation of the current state.
In PPG signals, motion artifacts are a common noise source, which may cause errors in heart rate extraction, and in order to remove the motion artifacts, researchers have proposed various algorithms, such as a method based on peak detection, a method based on adaptive filtering, a method based on independent component analysis, and the like.
Chengyu Liu, yanning Zhang, and Hua Yang, "An Improved PPG-based Heart Rate Estimation method," In 2016IEEE International Conference on Bioinformatics and Biomedicine (BIBM), pp.607-610,2016. This paper proposes An adaptive threshold-based peak detection method for extracting pulse wave peaks In the PPG signal and calculating the heart rate from the time interval between the peaks.
Wang, r., blackburn, g., desai, m., phelan, d., gillinov, l., houghting, p., gillinov, a (2015), accuracy of wrist-wood heart rate monitors, jama Cardiology,1 (3), 331-334, which evaluates the accuracy of a plurality of wrist-type heart rate monitors, some of which use PPG signals for heart rate calculation.
Li, q., rajagopalan, c., & Clifford, g.d. (2014) & Signal quality indices and data fusion for determining clinical acceptability of electrical diagnostics, physical Measurement,35 (9), 1745 the paper proposes a method based on adaptive filtering and peak detection for extracting the pulse wave peak in the PPG signal and calculating the heart rate.
Elgendi, m. (2012) On the analysis of fingertip photoplethysmogram signs.current Cardiology Reviews,8 (1), 14-25 the paper proposes a method based on adaptive filtering and peak detection for extracting pulse wave peaks in the PPG signal and calculating the heart rate from the time interval between the peaks.
The paper An Improved Method of Heart Rate Extraction Algorithm Based on Photoplethysmography for Sports Bracelet mainly describes a heart rate extraction algorithm based on a photoelectric measurement pulse wave signal (PPG) for heart rate monitoring in a sport wristband. The algorithm utilizes the idea of sliding window to segment PPG signals, and the heart rate value in each segment of signals is extracted by a peak detection and time interval calculation method.
From the foregoing, it can be seen that the methods presented in these papers mainly include the steps of filtering algorithms, peak detection, and heart rate calculation. The filtering algorithm is used for removing high-frequency noise and baseline drift in the PPG signal, the peak detection is used for extracting pulse wave peak values in the PPG signal, and the pulse wave peak values are used for calculating indexes such as heart rate and heart rhythm. These algorithms still suffer from a number of drawbacks, including mainly the following:
(1) Algorithm complexity: some high-precision motion artifact removal algorithms require complex mathematical models and calculation methods, so that the implementation difficulty is high, and the calculation resource requirement is high.
(2) The application range is limited: some algorithms are only suitable for specific motion patterns and motion intensities, and may not accurately remove motion artifacts for different motion types and motion intensities.
(3) The accuracy is limited: some algorithms do not work well for low amplitude artifacts, easily resulting in errors in heart rate extraction.
(4) The real-time performance is poor: some algorithms require long computation time and cannot meet the real-time requirements.
In view of this, the embodiments of the present application disclose a heart rate calculation method, which may be performed by a heart rate calculation device, which may be implemented by a terminal, such as a mobile terminal, a fixed terminal, or a portable terminal, such as a mobile phone, a multimedia computer, a multimedia tablet, a desktop computer, a notebook computer, a tablet computer, a smart bracelet, etc. The server may be an independent physical server, a server cluster or a distributed system formed by a plurality of physical servers, or a cloud server providing cloud services, cloud databases, cloud computing, cloud functions, cloud storage, network services, cloud communication, middleware services, domain name services, security services, CDNs, basic cloud computing services such as big data and artificial intelligence platforms, but is not limited thereto.
In one possible embodiment, the wrist of the target user wears a heart rate computing device with an inertial measurement unit (Inertial Measurement Unit, IMU) chip, such as MPU6050, embedded therein, in which a tri-axial accelerometer may collect acceleration data of the wrist of the target user. The tri-axis gyroscope in the chip may collect first angular velocity data of the wrist of the target user. The heart rate computing device is embedded with a photoelectric sensor, and PPG signal data of a target user can be acquired.
In another possible embodiment, the wrist of the target user wears other data acquisition devices, the other data acquisition devices are embedded with a photoelectric sensor and an IMU chip, the data acquisition devices acquire PPG signal data of the target user through the photoelectric sensor, acquire acceleration data of the wrist through a triaxial accelerometer in the IMU chip, and after acquiring first angular velocity data of the wrist through a triaxial gyroscope in the IMU chip, the acquired data are sent to the heart rate calculation device.
Further, the heart rate computing device performs heart rate computation based on the acceleration data, the first angular velocity data, and the PPG signal data. Fig. 1 is a flowchart of a heart rate calculating method according to an embodiment of the present application. The specific steps of the heart rate calculation method are described below in connection with fig. 1.
S101, fusing acceleration data and first angular velocity data of the wrist of the target user in a local horizontal coordinate system to obtain first angular data of the wrist in an IMU coordinate system.
The acceleration data includes accelerations of the wrist rotating about the X-axis, the Y-axis and the Z-axis, respectively, in the local horizontal coordinate system, the first angular velocity data includes angular velocities of the wrist rotating about the X-axis, the Y-axis and the Z-axis, respectively, in the local horizontal coordinate system, the first angular data includes roll angles of the wrist rotating about the X-axis, pitch angles of the wrist rotating about the Y-axis, in the IMU coordinate system, also called the carrier coordinate system.
S102, judging whether the pitch angle in the first angle data is smaller than a first threshold value or larger than a second threshold value.
After the heart rate computing device obtains the first angle data, it is determined whether the pitch angle in the first angle data is less than a first threshold or greater than a second threshold, the first threshold being less than the second threshold. If yes, S103 is executed, and if no, S104 is executed.
And S103, determining the baseline heart rate as a heart rate calculation result of the target user.
If the pitch angle in the first angle data is smaller than the first threshold value, the wrist movement amplitude is too small, and if the pitch angle in the first angle data is larger than the second threshold value, the data failure is indicated. The baseline heart rate, which is the heart rate measured while the body of the target user is stationary, is determined as the heart rate calculation of the target user.
S104, judging whether the pulse signal in the PPG signal and the motion artifact of the wrist are mutually independent.
If the pitch angle in the first angle data is greater than or equal to the first threshold value and the pitch angle in the first angle data is less than or equal to the second threshold value, judging whether the pulse signal in the PPG signal and the motion artifact of the wrist are mutually independent, if so, executing S105, and if not, executing S106.
S105, filtering the PPG signal to obtain a pulse signal, and obtaining a heart rate calculation result of the target user according to the pulse signal.
Specifically, the heart rate computing device may extract a pulse signal from the PPG signal by wiener filtering, and perform heart rate computation according to the pulse signal, to obtain a heart rate computation result of the target user.
S106, performing singular spectrum analysis on the PPG signal to obtain a heart rate calculation result of the target user.
In one possible embodiment, after the heart rate computing device obtains the data collected by the IMU chip (e.g., MPU 6050), the acceleration data and the first angular velocity data are preprocessed as follows:
s1.1, data calibration.
Because of the offset, scale factor error and non-orthogonality error of the IMU chip, the acceleration data and the first angular velocity data are calibrated separately using a 6-position method.
S1.2, filtering data.
The data acquired by the IMU chip has the problem of high-frequency noise, and the acceleration data and the first angular velocity data are respectively subjected to downsampling processing by using low-pass filtering to remove the high-frequency noise.
S1.3, aligning data.
The IMU chip generally collects data of 6 channels in total of the accelerometer and the gyroscope, different sampling rates and time stamps exist in the channels, and alignment processing is carried out on the acceleration data and the angular velocity data respectively, so that subsequent data fusion is facilitated.
In one implementation of S101, the data fusion steps are as follows:
s2.1, acquiring second angle data of the wrist under an IMU coordinate system according to the acceleration data.
The second angle data comprises a roll angle of the wrist rotating around the X axis and a pitch angle of the wrist rotating around the Y axis under an IMU coordinate system, and the calculation formula is as follows:
wherein r is acc For roll angle, p in the second angle data acc For pitch angle, a, in the second angle data x 、a y And a z Acceleration data measured for an accelerometer, a x Acceleration, a, of the wrist in a local horizontal coordinate system rotating around the X axis y A is acceleration of wrist rotation around Y axis under local horizontal coordinate system z Is the acceleration of the wrist rotating around the Z axis under the local horizontal coordinate system.
S2.2, obtaining second angular velocity data of the wrist under the IMU coordinate system according to the first angular velocity data and the second angular velocity data.
The first angular velocity data comprises the angular velocities of the wrist rotating around the X axis, the Y axis and the Z axis respectively under the local horizontal coordinate system, and the second angular velocity data comprises the angular velocities of the wrist rotating around the X axis, the Y axis and the Z axis respectively under the IMU coordinate system, and the calculation formula is as follows:
wherein [ dr/dt dp/dt dy/dt] T For the second angular velocity data, dr/dt is the angular velocity of the wrist rotating about the X-axis in the IMU coordinate system, dp/dt is the angular velocity of the wrist rotating about the Y-axis in the IMU coordinate system, and dy/dt is the angular velocity of the wrist rotating about the Z-axis in the IMU coordinate system. [ g ] x g y g z ] T G is the first angular velocity data x G is the angular velocity of the wrist rotating around the X axis under the local horizontal coordinate system y G is the angular velocity of the wrist rotating around the Y axis under the local horizontal coordinate system z Is the angular velocity of the wrist rotating about the Z-axis in the local horizontal coordinate system. r is the roll angle in the second angle data and p is the pitch angle in the second angle data.
S2.3, fusing the second angular velocity data and the second angle data to obtain first angle data.
The first angle data includes roll angle of the wrist rotating about the X-axis and pitch angle of the wrist rotating about the Y-axis in the IMU coordinate system. The first angle data is obtained through a Kalman filtering equation, and the calculation formula is as follows:
wherein k represents time, x k,k For the state at time k after filtering, i.e. the first angle data, r k For roll angle at time k, i.e. in the first angle data, p k For pitch angle at time k, i.e. pitch angle in the first angle data, x k，k-1 The state at time k, which is derived from the state at time k-1, is represented by the state extrapolation equation:
r k+1 roll angle at time k+1, p k+1 Is the pitch angle at time k+1, [ r ] k ,p k ]The initial value is set to 0,0]The method comprises the steps of carrying out a first treatment on the surface of the dr/dt is the angular velocity of rotation around the X axis in the second angular velocity data, and dp/dt is the angular velocity of rotation around the Y axis in the second angular velocity data; k (K) k Is a Kalman filtering matrix; z is Z k For the observation equation, the formula is:r acc for roll angle, p in the second angle data acc Is the pitch angle in the second angle data.
The above describes how the data fusion is performed by kalman filtering, which involves the conversion of the local horizontal coordinate system and the IMU coordinate system, and how the first angular velocity data is obtained is described in detail below in connection with the conversion of the two coordinate systems.
S3.1, defining Euler angles.
The euler angles are defined using the order of rotation Z, Y, X and assuming that the IMU coordinate system (carrier coordinate system) coincides with the local horizontal coordinate system at the initial time. The rotation is sequentially performed around a Z, Y, X axis of the IMU chip, an angle of rotation around a Z axis of the IMU chip is defined as a heading angle, an angle of rotation around a Y axis of the IMU chip is defined as a pitch angle, and an angle of rotation around an X axis of the IMU chip is defined as a roll angle.
S3.2, defining a rotation matrix.
The rotation matrix for rotation about the Z axis is:y is the heading angle.
The rotation matrix for rotation about the Y axis is:p is pitch angle.
The rotation matrix for rotation about the X axis is:r is the roll angle.
And S3.3, calculating the attitude Euler angle through an accelerometer.
When the posture of the IMU chip is not aligned with the local horizontal coordinate system, the projection of the gravity acceleration on the IMU chip is correspondingly changed, and the essence of the projection is that the coordinates of the gravity vector (0, g) under the local horizontal coordinate system under the IMU coordinate system are changed, and 3 values measured by the accelerometer are new projection coordinates of the gravity vector (0, g).
According to the definition of the rotation matrix, when the course angle, the pitch angle and the roll angle exist at the same time, describing the projection process of the gravity vector as follows according to a preset gesture rotation sequence:
C x 、C y 、C z rotation matrix, a, defined for S3.2 respectively x 、a y And a z Accelerometer measurements, i.e. accelerations in local horizontal coordinate system rotating around the X-axis, Y-axis and Z-axis, respectively, can be developed as:
in the process of calculating the measured value of the accelerometer, the heading angle y is not used, the heading angle cannot be calculated by using the accelerometer, and the formula is expressed as a form of an equation set:
the derivation process is as follows:
a y 2 ＝g 2 sinr 2 cosp 2
a z 2 ＝g 2 cosr 2 cosp 2
a y 2 +a z 2 ＝g 2 cosp 2
a y 2 +a z 2 ＝g 2 cosp 2
thus, the roll angle and pitch angle are represented by the acceleration values measured by the accelerometer:
wherein r is acc For the roll angle, p acc Is pitch angle, a x 、a y And a z Acceleration data measured for the accelerometer, i.e. acceleration of the wrist in a local horizontal coordinate system about the X-axis, Y-axis and Z-axis.
And S3.4, calculating an attitude Euler angle through a gyroscope.
Let the Euler angles of the IMU chip at the time t be r, p and y, and let the Euler angles of the IMU chip at the time t+Δt be r+Deltar, p+Deltap and y+Deltay, the calculation mode of the change amount of the attitude angle is as follows:
the attitude euler angles r+Δr, p+Δp, y+Δy at time t+Δt are expressed as:
Where dy/dt represents the angular velocity around the Z-axis, vector 0 [ 0dy/dt] T An angular velocity vector representing rotation about the Z axis in the IMU coordinate system at time t, [ 0dy/dt ]] T After two rotations, the coordinate of the Z axis in the time t+Deltat is obtained, dp/dt represents the angular velocity of the Y axis of the new coordinate system after rotating by Deltay angle around the Z axis under the IMU coordinate system at the time t, and the vector [ 0dy/dt ]] T Representing the angular velocity vector, which also requires a rotational transformation to obtain the Y-axis coordinate at time t+Δt, dr/dt representingAfter the IMU coordinate system rotates by delta Y angle around the Z axis at the moment t, the new IMU coordinate system rotates by delta p angle around the Y axis to obtain the angular velocity of the IMU coordinate system around the X axis, and the vector [ dr/dt 0 0] T Representing the angular velocity vector, and obtaining the posture at the moment t+deltat after rotating for delta r again x 、C y 、C z A rotation matrix defined in S3.2;
thus, the angular velocity vector measured using the gyroscope represents the euler angular velocity vector:
wherein [ dr/dt dp/dt dy/dt] T For angular velocity vectors of wrist rotation about X-axis, Y-axis and Z-axis respectively in IMU coordinate system, [ g ] x g y g z ] T The angular velocity vectors measured for the gyroscope, i.e. the angular velocity vectors of the wrist rotating about the X-axis, Y-axis and Z-axis, respectively, in a local horizontal coordinate system. r is the roll angle r deduced from S3.3 acc P is the pitch angle p derived from S3.3 acc 。
S3.5, designing a Kalman filtering equation.
The accelerometer can calculate the roll angle r at the rest moment according to the sensed gravitational acceleration acc And pitch angle p acc And the angle calculation is related only to the current pose. And the angular velocity of the gyroscope in the time interval is integrated to obtain the angle conversion quantity of each time, the angle conversion quantity is added to the last attitude angle to obtain a new attitude angle, and the roll angle, the pitch angle and the yaw angle can be calculated according to the angular velocity data measured by the gyroscope.
The acceleration can obtain more accurate gesture only at the stationary moment, the gyroscope is only sensitive to gesture change during rotation, and if the gyroscope has errors, the errors can be continuously increased after continuous time integration. Thus, a complementary fusion is required in combination with the two calculated poses. Because the accelerometer does not obtain the yaw angle, only the roll angle and the pitch angle can be fused, the yaw angle is insensitive when the PPG signal is acquired, the pitch angle has the largest influence on the acquisition of the PPG signal, and the roll frequency is calculated, so that the following updating mode is adopted:
(1) Designing a state extrapolation equation: x is x k+1,k ＝Fx k,k +Gu k ；
The state extrapolation equation is developed as:
Wherein r is k Roll angle at time k, p k Is the pitch angle at the moment k,the initial value is set to->I.e. assuming that the initial moment is horizontally placed.
(2) Designing an observation equation: z is Z k ＝Hx k The method comprises the steps of carrying out a first treatment on the surface of the H is the observation matrix, here the identity matrix, and the observation equation is developed as:
(3) Designing a state covariance matrix at an initial moment:the above indicates that the roll angle variance at the initial time is 1 radian, the pitch angle variance is 1 radian, and the covariance between the roll angle and the pitch angle is 0.
(4) The design noise covariance matrix Q is:
(5) Designing a covariance extrapolation equation: p (P) k+1,k ＝FP k,k F T +Q；
(6) Calculating a Kalman filter matrix:
(7) The equation formula for the state update is: x is x k,k ＝x k,k-1 +K k (Z k -x k,k-1 )；
(8) The covariance update equation is: p (P) k,k ＝(I-K k H)P k,k-1 ；
Kalman filtering is carried out through the 8 equations to obtain the latest stateNamely roll angle and pitch angle in the first angular velocity data.
In one possible embodiment, after the heart rate computing device obtains the PPG signal, the PPG signal may be preprocessed as follows:
and S4.1, removing frequency components with frequency components larger than a fourth threshold value from the original PPG signal by using a low-pass filter, and removing frequency components smaller than a fifth threshold value by using a high-pass filter to obtain a filtered PPG signal. Wherein the fourth threshold is, for example, 3.5HZ and the fifth threshold is, for example, 0.8HZ.
And S4.2, processing the filtered PPG signal by using a first average filter to obtain an average filter output PPG signal. Wherein the window size of the first averaging filter is for example 6, by filtering the baseline drift in the PPG signal can be removed.
And S4.3, subtracting the average filtered PPG signal from the filtered PPG signal to obtain a normalized output PPG signal.
And S4.4, processing the normalized output PPG signal by using a second average filter to obtain a preprocessed output PPG signal. Wherein the window size of the second averaging filter is for example 3.
In one implementation manner of S102, the specific step of determining, by the heart rate computing device, whether the pulse signal in the PPG signal and the motion artifact of the wrist are independent of each other is as follows:
s5.1, performing discrete Fourier transform and normalization processing on the PPG signal to obtain a normalized waveform diagram of the PPG signal.
Specifically, the heart rate computing device performs discrete fourier transform on the PPG signal output by preprocessing to obtain a waveform diagram of the PPG signal, and performs normalization operation on the waveform diagram of the PPG signal to obtain a PPG normalized waveform diagram.
S5.2, performing discrete Fourier transform and normalization processing on the waveform diagram of the wrist motion degree to obtain a normalized waveform diagram of the wrist motion degree.
The waveform diagram of the degree of wrist movement is used to indicate the magnitude of movement of the wrist at various moments in time. The heart rate calculating device can calculate the square sum of the angular velocities of the wrist rotating around the X axis, the Y axis and the Z axis respectively under the IMU coordinate system according to the second angular velocity data, determine the square sum as a motion amplitude value of the wrist at each moment, and obtain a waveform diagram of the wrist motion degree by fitting according to the motion amplitude value of the wrist at each moment.
The second angular velocity data is obtained by referring to the foregoing discussion, and will not be described herein. The formula for calculating the motion amplitude value is as follows:
mag is the magnitude of the degree of movement of the wrist at time t,for angular velocity of the wrist rotation about the X-axis in the IMU coordinate system, i.e. angular velocity of roll angle, +.>For the angular velocity of the wrist rotation about the Y-axis in IMU coordinate system, i.e. the angular velocity of the pitch angle, +.>Is the angular velocity of the wrist rotating about the Z-axis in the IMU coordinate system, i.e. the angular velocity of the yaw angle.
S5.3, judging whether the pulse signal in the PPG signal and the motion artifact of the wrist are mutually independent according to the normalized waveform diagram of the PPG signal and the normalized waveform diagram of the wrist motion degree.
In one possible embodiment, the heart rate computing device determines a first heart rate value corresponding to the frequency component with the largest peak from the normalized waveform of the PPG signal, and a second heart rate value corresponding to the frequency component with the largest peak from the normalized waveform of the wrist motion level. The heart rate computing device may determine whether an absolute value of a difference between the first heart rate value and the second heart rate value is greater than or equal to a third threshold, if so, determine that a pulse signal in the PPG signal and a motion artifact of the wrist are independent of each other, and if not, determine that the pulse signal in the PPG signal and the motion artifact of the wrist are not independent of each other.
In another possible embodiment, the heart rate computing device may determine whether a ratio between the first heart rate value and the second heart rate value is greater than or equal to a preset threshold, if so, determine that the pulse signal in the PPG signal and the motion artifact of the wrist are independent of each other, and if not, determine that the pulse signal in the PPG signal and the motion artifact of the wrist are not independent of each other.
In one implementation of S106, the step of performing, by the heart rate computing device, a singular spectrum analysis on the PPG signal is as follows: s5.1, converting the normalized waveform diagram of the PPG signal into a track matrix, and carrying out singular value decomposition on the track matrix to obtain different mode groups of the track matrix.
Specifically, the normalized waveform of the PPG signal has a one-dimensional time sequence of [ x ] 1 、x 2 、...、x N ]N is the sequence length, and the formula for calculating the track matrix X is as follows:
where L is a given window.
Calculate X T All characteristic vectors of X are arranged in columns to obtain a unit orthogonal array V, and XX is calculated T Is arranged in a column to obtain a unit orthogonal array U, and X is calculated T X, a plurality of singular values are obtained by squaring the characteristic values, the plurality of singular values form a diagonal matrix P, and a diagonal element lambda of the diagonal matrix P 1 、λ 2 、...、λ L Ranging from large to small.
It can be verified that the following holds: x=upv T ；
At the same time, the following formula can be verified to be established:wherein lambda is i X is the ith diagonal element of the diagonal array P 1 、X 2 、...、X L Is a different set of modes of the track matrix X.
S5.2, performing anti-diagonal line equalization processing on different mode groups to obtain a plurality of mode sequences with the same length.
For different mode groups X 1 、X 2 、...、X L Performing anti-diagonal equalization to obtain L equal-length mode sequences with length N, respectively denoted as t 1 、t 2 、...、t L 。
S5.3, according to a plurality of mode sequences with the same length, obtaining a reconstructed PPG signal, performing discrete Fourier transform on the reconstructed PPG signal, and determining the main frequency of the transformed PPG signal as a heart rate calculation result of the target user.
Specifically, a pattern sequence t is selected 1 、t 2 、...、t L The first Q mode sequences in (a) are used as main modes of PPG signals, PPG signal construction operation is carried out, and a construction formula is as follows:y represents the reconstructed PPG signal, t i Representation ofAn i-th sequence. The reconstructed PPG signal is subjected to a discrete fourier transform and its dominant frequency is calculated, which is determined as the heart rate calculation of the target user.
In summary, to the problem of limited motion range, the application fuses the data of multiple sensors, so as to improve the application range of the motion artifact removal algorithm, and in particular, the application fuses the data of the accelerometer and the gyroscope with the PPG data, so as to improve the accuracy and stability of the motion artifact removal algorithm. Aiming at the problem of limited precision, the method comprehensively uses Fourier transform and singular spectrum analysis to predict the heart rate of the PPG signal at two layers of a frequency domain and a time domain. Aiming at the problem of poor instantaneity, the method and the device downsamples the high-frequency PPG signal and the high-frequency IMU signal through mean value filtering, improves the operation efficiency of the algorithm by adopting the measures such as fast Fourier transform technology, optimizing the algorithm structure and the like, and improves the instantaneity of the algorithm.
The embodiment of the application also discloses a heart rate calculation device which can be arranged in the heart rate calculation equipment. Referring to fig. 2, a heart rate calculation apparatus includes:
an obtaining module 201, configured to: acquiring first angle data of the wrist under an Inertial Measurement Unit (IMU) coordinate system according to acceleration data and first angular velocity data of the wrist of a target user under a local horizontal coordinate system; the first angle data comprises a pitch angle of rotation of the wrist around a Y axis under an IMU coordinate system;
a judging module 202, configured to: judging whether the pitch angle in the first angle data is smaller than a first threshold value or larger than a second threshold value; the first threshold is less than the second threshold;
heart rate calculation module 203 for: if the pitch angle in the first angle data is smaller than a first threshold value or larger than a second threshold value, determining the baseline heart rate as a heart rate calculation result of the target user;
the judging module 202 is further configured to: if the pitch angle in the first angle data is larger than or equal to a first threshold value and smaller than or equal to a second threshold value, judging whether a pulse signal in the photoelectric measurement pulse wave PPG signal and the motion artifact of the wrist are mutually independent;
the heart rate calculation module 203 is further configured to filter the PPG signal to obtain a pulse signal if the pulse signal in the PPG signal and the motion artifact of the wrist are independent of each other, and obtain a heart rate calculation result of the target user according to the pulse signal; if the pulse signal in the PPG signal and the motion artifact of the wrist are not mutually independent, singular spectrum analysis is carried out on the PPG signal, and a heart rate calculation result of the target user is obtained.
Optionally, the obtaining module 201 is specifically configured to: acquiring second angle data of the wrist under an IMU coordinate system according to the acceleration data; the second angle data comprises a roll angle of the wrist rotating around the X axis and a pitch angle of the wrist rotating around the Y axis under an IMU coordinate system; obtaining second angular velocity data of the wrist under an IMU coordinate system according to the first angular velocity data and the second angular velocity data; the first angular velocity data comprises angular velocities at which the wrist rotates about an X axis, a Y axis and a Z axis, respectively, in a local horizontal coordinate system; the second angular velocity data comprises angular velocities at which the wrist rotates about an X axis, a Y axis and a Z axis respectively in an IMU coordinate system; fusing the second angular velocity data and the second angle data to obtain first angle data; the first angle data also includes a roll angle at which the wrist rotates about the X-axis under the IMU coordinate system.
Optionally, the obtaining module 201 is specifically configured to:
the first angle data is obtained through a Kalman filtering equation, and the calculation formula is as follows:
wherein k represents time, x k,k For the state at time k after filtering, i.e. the first angle data, r k For roll angle at time k, i.e. in the first angle data, p k For pitch angle at time k, i.e. pitch angle in the first angle data, x k,k-1 The state at time k, which is derived from the state at time k-1, is represented by the state extrapolation equation:
r k+1 roll angle at time k+1, p k+1 Is the pitch angle at time k+1, [ r ] k ,p k ]The initial value is set to 0,0]The method comprises the steps of carrying out a first treatment on the surface of the dr/dt is the angular velocity of rotation about the X axis in the second angular velocity data, dp/dt is the angular velocity of rotation about the Y axis in the second angular velocity data; k (K) k Is a Kalman filtering matrix; z is Z k For the observation equation, the formula is:r acc for roll angle, p in the second angle data acc Is the pitch angle in the second angle data.
With continued reference to fig. 2, the apparatus further includes a processing module 204, the processing module 204 configured to:
performing discrete Fourier transform and normalization processing on the PPG signal to obtain a normalized waveform diagram of the PPG signal; performing discrete Fourier transform and normalization processing on the waveform diagram of the wrist movement degree to obtain a normalized waveform diagram of the wrist movement degree; the oscillogram of the wrist motion degree is used for indicating the motion amplitude of the wrist at each moment;
the judging module 202 is configured to: and judging whether the pulse signal in the PPG signal and the motion artifact of the wrist are mutually independent according to the normalized waveform diagram of the PPG signal and the normalized waveform diagram of the wrist motion degree.
With continued reference to FIG. 2, the apparatus further includes a determination module 205;
The determining module 205 is configured to: determining a first heart rate value corresponding to a frequency component with a maximum peak value from a normalized waveform diagram of the PPG signal, and determining a second heart rate value corresponding to the frequency component with the maximum peak value from a normalized waveform diagram of the wrist movement degree; the judging module 202 is configured to: judging whether the absolute value of the difference value between the first heart rate value and the second heart rate value is larger than or equal to a third threshold value; if yes, determining that a pulse signal in the PPG signal and motion artifact of the wrist are mutually independent; if not, determining that the pulse signal in the PPG signal and the motion artifact of the wrist are not mutually independent.
Optionally, the determining module 205 is further configured to: according to the second angular velocity data, calculating the square sum of angular velocities of the wrist rotating around the X axis, the Y axis and the Z axis respectively under the IMU coordinate system, and determining the square sum as a motion amplitude value of the wrist at each moment;
the obtaining module 201 is configured to: and fitting to obtain a waveform diagram of the wrist movement degree according to the movement amplitude value of the wrist at each moment.
Optionally, the heart rate calculation module 203 is specifically configured to:
converting the normalized waveform diagram of the PPG signal into a track matrix, and carrying out singular value decomposition on the track matrix to obtain different mode groups of the track matrix; performing anti-diagonal equalization processing on different mode groups to obtain a plurality of mode sequences with the same length; and according to the mode sequences with the same length, obtaining a reconstructed PPG signal, performing discrete Fourier transform on the reconstructed PPG signal, and determining the main frequency of the transformed PPG signal as a heart rate calculation result of the target user.
The heart rate calculating device of the embodiment of the application can realize any one of the above-mentioned heart rate calculating methods, and the specific working process of each module in the heart rate calculating device can refer to the corresponding process in the above-mentioned method embodiment.
In several embodiments provided herein, it should be understood that the provided methods and apparatus may be implemented in other ways. For example, the device embodiments described above are merely illustrative; for example, a division of a module is merely a logical function division, and there may be another division manner in actual implementation, for example, multiple modules may be combined or may be integrated into another device, or some features may be omitted or not performed.
The embodiment of the application also discloses a computer device. Computer device comprising a memory, a processor, and a computer program stored on the memory and executable on the processor, the processor implementing a heart rate calculation method as described above when executing the computer program.
The embodiment of the application also discloses a computer readable storage medium. A computer readable storage medium storing a computer program that can be loaded by a processor and that performs any of the heart rate calculation methods described above. A computer readable storage medium may be any tangible medium that can contain, or store a program for use by or in connection with an instruction execution system, apparatus, or device; program code embodied on a computer readable medium may be transmitted using any appropriate medium, including but not limited to wireless, wireline, optical fiber cable, RF, etc., or any suitable combination of the foregoing.
In the foregoing embodiments, the descriptions of the embodiments are focused on, and for those portions of one embodiment that are not described in detail, reference may be made to the related descriptions of other embodiments.
The foregoing description of the preferred embodiments of the present application is not intended to limit the scope of the application, in which any feature disclosed in this specification (including abstract and drawings) may be replaced by alternative features serving the same, equivalent or similar purpose, unless expressly stated otherwise. That is, each feature is one example only of a generic series of equivalent or similar features, unless expressly stated otherwise.
Claims (10)
1. A heart rate calculation method, the method comprising:
according to the acceleration data and the first angular velocity data of the wrist of the target user in a local horizontal coordinate system, fusion is carried out, and first angle data of the wrist in an Inertial Measurement Unit (IMU) coordinate system is obtained; the first angle data comprises a pitch angle of rotation of the wrist around a Y axis under an IMU coordinate system;
judging whether the pitch angle in the first angle data is smaller than a first threshold value or larger than a second threshold value; the first threshold is less than the second threshold;
if yes, determining a baseline heart rate as a heart rate calculation result of the target user; if not, judging whether the pulse signal in the photoelectric measurement pulse wave PPG signal and the motion artifact of the wrist are mutually independent;
If yes, filtering the PPG signal to obtain the pulse signal, and obtaining a heart rate calculation result of the target user according to the pulse signal; if not, carrying out singular spectrum analysis on the PPG signal to obtain a heart rate calculation result of the target user.
2. The heart rate calculation method according to claim 1, wherein the obtaining the first angle data of the wrist in the IMU coordinate system according to the fusion of the acceleration data and the first angular velocity data of the wrist of the target user in the local horizontal coordinate system includes:
acquiring second angle data of the wrist under an IMU coordinate system according to the acceleration data; the second angle data comprises a roll angle of the wrist rotating around an X axis and a pitch angle of the wrist rotating around a Y axis under an IMU coordinate system;
obtaining second angular velocity data of the wrist under the IMU coordinate system according to the first angular velocity data and the second angular velocity data; the first angular velocity data comprises angular velocities at which the wrist rotates about an X axis, a Y axis and a Z axis respectively in a local horizontal coordinate system; the second angular velocity data comprises angular velocities at which the wrist rotates about an X axis, a Y axis and a Z axis respectively in an IMU coordinate system;
Fusing the second angular velocity data and the second angle data to obtain the first angle data; the first angle data also includes a roll angle at which the wrist rotates about an X-axis in an IMU coordinate system.
3. The heart rate calculation method according to claim 2, wherein the fusing the second angular velocity data and the second angle data to obtain the first angle data includes:
the first angle data is obtained through a Kalman filtering equation, and the calculation formula is as follows:
wherein k represents time, x k,k For the state at time k after filtering, i.e. the first angle data, r k Roll at time kThe angle, i.e. the roll angle, p in the first angle data k Is the pitch angle at time k, namely the pitch angle x in the first angle data k,k-1 The state at time k, which is derived from the state at time k-1, is represented by the state extrapolation equation:
r k+1 roll angle at time k+1, p k+1 Is the pitch angle at time k+1, [ r ] k ,p k ]The initial value is set to 0,0]The method comprises the steps of carrying out a first treatment on the surface of the dr/dt is the angular velocity of rotation around the X axis in the second angular velocity data, dp/dt is the angular velocity of rotation around the Y axis in the second angular velocity data; k (K) k Is a Kalman filtering matrix; z is Z k For the observation equation, the formula is:r acc for the roll angle, p, in the second angle data acc And the pitch angle in the second angle data.
4. A heart rate calculation method according to claim 2, wherein determining whether the pulse signal in the PPG signal and the motion artifact of the wrist are independent of each other comprises:
performing discrete Fourier transform and normalization processing on the PPG signal to obtain a normalized waveform diagram of the PPG signal;
performing discrete Fourier transform and normalization processing on the waveform diagram of the wrist movement degree to obtain a normalized waveform diagram of the wrist movement degree; the oscillogram of the wrist motion degree is used for indicating the motion amplitude of the wrist at each moment;
and judging whether the pulse signal in the PPG signal and the motion artifact of the wrist are mutually independent or not according to the normalized waveform diagram of the PPG signal and the normalized waveform diagram of the wrist motion degree.
5. The method according to claim 4, wherein determining whether the pulse signal in the PPG signal and the motion artifact of the wrist are independent of each other comprises:
Determining a first heart rate value corresponding to a frequency component with a maximum peak value from a normalized waveform diagram of the PPG signal, and determining a second heart rate value corresponding to the frequency component with the maximum peak value from a normalized waveform diagram of the wrist movement degree;
judging whether the absolute value of the difference value between the first heart rate value and the second heart rate value is greater than or equal to a third threshold value;
if yes, determining that a pulse signal in the PPG signal and the motion artifact of the wrist are mutually independent; if not, determining that the pulse signal in the PPG signal and the motion artifact of the wrist are not mutually independent.
6. The heart rate calculation method according to claim 5, wherein the waveform of the degree of wrist movement is obtained by:
calculating the square sum of the angular speeds of the wrist rotating around the X axis, the Y axis and the Z axis respectively under an IMU coordinate system according to the second angular speed data;
determining the sum of squares as a motion amplitude value of the wrist at each moment;
and fitting to obtain a waveform diagram of the wrist motion degree according to the motion amplitude value of the wrist at each moment.
7. A heart rate calculation method according to any one of claims 4-6, wherein performing singular spectrum analysis on the PPG signal to obtain a heart rate calculation result of the target user comprises:
Converting the normalized waveform diagram of the PPG signal into a track matrix, and carrying out singular value decomposition on the track matrix to obtain different mode groups of the track matrix;
performing anti-diagonal equalization processing on the different mode groups to obtain a plurality of mode sequences with the same length;
and according to the mode sequences with the same length, obtaining a reconstructed PPG signal, performing discrete Fourier transform on the reconstructed PPG signal, and determining the main frequency of the transformed PPG signal as a heart rate calculation result of the target user.
8. A heart rate computing device, comprising:
an obtaining module for: acquiring first angle data of the wrist under an Inertial Measurement Unit (IMU) coordinate system according to acceleration data and first angular velocity data of the wrist of a target user under a local horizontal coordinate system; the first angle data comprises a pitch angle of rotation of the wrist around a Y axis under an IMU coordinate system;
the judging module is used for: judging whether the pitch angle in the first angle data is smaller than a first threshold value or larger than a second threshold value; the first threshold is less than the second threshold;
a heart rate calculation module for: if the pitch angle in the first angle data is smaller than a first threshold value or larger than a second threshold value, determining a baseline heart rate as a heart rate calculation result of the target user;
The judging module is further configured to: if the pitch angle in the first angle data is larger than or equal to a first threshold value and smaller than or equal to a second threshold value, judging whether a pulse signal in a photoelectric measurement pulse wave PPG signal and the motion artifact of the wrist are mutually independent;
the heart rate calculation module is further configured to filter the PPG signal if the pulse signal in the PPG signal and the motion artifact of the wrist are independent of each other, obtain the pulse signal, and obtain a heart rate calculation result of the target user according to the pulse signal; if the pulse signal in the PPG signal and the motion artifact of the wrist are not mutually independent, performing singular spectrum analysis on the PPG signal to obtain a heart rate calculation result of the target user.
9. A computer device, characterized by: comprising a memory and a server, said memory having stored thereon a computer program for loading and executing by the server a method according to any of claims 1-7.
10. A computer readable storage medium, characterized in that a computer program is stored which can be loaded by a server and which performs the method according to any of claims 1-7.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311369772.8A CN117414117A (en) | 2023-10-20 | 2023-10-20 | Heart rate calculation method, device, equipment and medium |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311369772.8A CN117414117A (en) | 2023-10-20 | 2023-10-20 | Heart rate calculation method, device, equipment and medium |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117414117A true CN117414117A (en) | 2024-01-19 |
Family
ID=89529560
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311369772.8A Pending CN117414117A (en) | 2023-10-20 | 2023-10-20 | Heart rate calculation method, device, equipment and medium |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117414117A (en) |
-
2023
- 2023-10-20 CN CN202311369772.8A patent/CN117414117A/en active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Roetenberg et al. | Estimating body segment orientation by applying inertial and magnetic sensing near ferromagnetic materials | |
| US8712508B2 (en) | State detection device, electronic apparatus, measurement system and program | |
| Sessa et al. | A methodology for the performance evaluation of inertial measurement units | |
| CN108836346A (en) | A kind of Human Body Gait Analysis method and system based on inertial sensor | |
| CN104757976A (en) | Human gait analyzing method and system based on multi-sensor fusion | |
| Guo et al. | A low-cost body inertial-sensing network for practical gait discrimination of hemiplegia patients | |
| Steins et al. | A smart device inertial-sensing method for gait analysis | |
| CN109846487B (en) | Thigh movement posture measuring method and device based on MIMU/sEMG fusion | |
| Wang et al. | Estimate spatial-temporal parameters of human gait using inertial sensors | |
| CN109893111B (en) | Dynamic blood pressure measurement mode selection method and device | |
| Sabatini et al. | Fourier-based integration of quasi-periodic gait accelerations for drift-free displacement estimation using inertial sensors | |
| CN109284006B (en) | Human motion capturing device and method | |
| US11833392B2 (en) | Methods and systems for swimming performance analysis | |
| CN108836344A (en) | Step-length cadence evaluation method and device and gait detector | |
| US11395633B2 (en) | Systems and methods for determining engagement of a portable device | |
| Segarra et al. | Orientation estimation using filter-based inertial data fusion for posture recognition | |
| Young et al. | Distributed estimation of linear acceleration for improved accuracy in wireless inertial motion capture | |
| Huang et al. | Wearable Sensor for Measurement of Gait Walking and Running Motion. | |
| CN117414117A (en) | Heart rate calculation method, device, equipment and medium | |
| Martin | Real time patient's gait monitoring through wireless accelerometers with the wavelet transform | |
| US20220292694A1 (en) | Information processing apparatus, method, and non-transitory computer-readable storage medium | |
| Alcaraz et al. | Mobile quantification and therapy course tracking for gait rehabilitation | |
| Okkalidis et al. | A review of foot pose and trajectory estimation methods using inertial and auxiliary sensors for kinematic gait analysis | |
| Tinh et al. | Image-based gramian angular field processing for pedestrian stride-length estimation using convolutional neural network | |
| Biswas et al. | On fall detection using smartphone sensors |
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 |