EP2954840A1

EP2954840A1 - Method for the estimation of the heart-rate and corresponding system

Info

Publication number: EP2954840A1
Application number: EP15170706.4A
Authority: EP
Inventors: Stefano Cervini
Original assignee: STMicroelectronics SRL
Current assignee: STMicroelectronics SRL
Priority date: 2014-06-09
Filing date: 2015-06-04
Publication date: 2015-12-16
Anticipated expiration: 2035-06-04
Also published as: US9936886B2; US20150351646A1; CN109381175A; EP2954840B1; CN109381175B; CN205083467U; CN105125198A; CN105125198B

Abstract

Method for the estimation of the heart-rate using photoplethysmography on a body organ, in particular a wrist of a user, comprising acquiring optically (13a, 13b, 17) from said body organ a heart beat signal (o), acquiring (21) an acceleration signal (ax, ay, az) representative of the acceleration of said body organ, selecting data blocks (Bo, Bx, By, Bz) of said acquired heart beat signal (o) and acceleration signal (ax, ay, az), compensating said heart beat signal (o) by the acceleration signal (ax, ay, az), calculating the heart rate value (r) on the basis of said compensated heart beat signal (O') .

The method provides
obtaining (130) from said selected data blocks (Bo, Bx, By, Bz) corresponding frequency domain data blocks for the heart beat signal (O) and for the acceleration signal (AX, AY, AZ),
said compensating operation includes performing a motion compensation (210) in the frequency domain, compensating the frequency domain data blocks for the heart beat (O) with the corresponding frequency domain data blocks for the acceleration signal (AX, AY, AZ).

Description

Technical field

The present description relates to techniques for the estimation of the heart-rate using photoplethysmography on a body organ and an information on acceleration of said body organ.
Various embodiments may apply e.g. in wearable, in particular wrist-wearable devices, for continuous monitoring of the heart-rate in fitness/wellness applications.

Description of the prior art

Among the methods for the estimation of the heart-rate in a subject, one of the most employed makes use of optical means to detect heart beat to evaluate the heart-rate, based on photoplethysmography (PPG). Photoplethysmography involves obtaining optically a volumetric measurement of an organ (plethysmogram). A photoplethysmogram is often obtained by using a pulse oximeter which illuminates the skin and measures changes in light absorption. A conventional pulse oximeter monitors the perfusion of blood to the dermis and subcutaneous tissue of the skin.
PPG was historically first employed with finger clips in medical applications. Lately PPG has been employed also on wrist, arm, forearm, to make it suitable for fitness applications.
The PPG techniques typically allow an easy estimation of the heart rate at rest. In motion conditions the PPG signal is affected however by a very low SNR (Signal to Noise Ratio) and SIR (Signal to Interference Ratio. In this specific context SNR means the ratio of the heart rate, i.e. the signal S, to all the other signals, i.e. the noise N. It can be measured by having a subject wearing a cardio frequency meter and a PPG device. The power of the frequency component (peak) measured by the PPG device closest to cardio frequency measured by the cardio frequency meter is the signal S. The total power T is T=(S+N). Thus, for the SNR it is S/N=S/(T-S). The SIR value is obtained as ratio of the signal S over I, the power of the strongest non-cardiac frequency component.
Poor performances are sometimes obtained by processing a PPG signal with known frequency-domain techniques or time-domain techniques, including adaptive filtering techniques.
Thus problems affecting the PPG are strong motion artifacts leading to low SNR, low SIR due to motion artifacts, spikes in the detected signal due to motion.
In the state of the art it is known, in order to compensate the above discussed effects of motion on SNR and SIR and other aspect of the signal detection, to acquire optically from the body organ a heart rate signal, acquire an acceleration signal representative of the acceleration of such body organ, selecting data blocks of said acquired heart rate signal and acceleration signal, compensating in the time domain the heart rate signal by the acceleration signal.
The following publications described similar prior art techniques for the reduction of the artifacts:

Yuta Kuboyama, "Motion Artifact Cancellation for Wearable Photoplethismography Sensor", MIT, 2009;
K. Ashoka eddy, V. Jagadeesh Kumar, "Motion Artifact Reduction in Photoplethysmographic Signals using Singular Value Decomposition", Instrumentation and Measurement Technology Conference, Poland, May 1-3, 2007;
H. Han, M. Kim, J. Kim, "Development of real-time motion artifact reduction algorithm for a wearable photoplethismography", Proceedings of the 29 th Annual International Conference of the IEEE EMBS Cite Internationale, Lyon, France, August 23-26, 2007;
P. Wei, R. Guo, J. Zhang, Y.T. Zhang, "A New Wristband Sensor Using Adaptive Reduction Filter to Reduce Motion Artifact", Proceedings of the 5th International Conference on Information Technology and Application in Biomedicine, Shenzen, China, May 30-31, 2008.

Object and summary

An object of one or more embodiments is to provide a communication system that solves the drawbacks of the prior art and in particular allows:

to overcome noise by means of high processing gain;
to overcome interference due to motion in the frequency domain (motion compensation);
to build a secondary estimator to be used whenever the primary estimator (optical) fails.

According to one or more embodiments, that object is achieved thanks to method having the characteristics specified in Claim 1. One or more embodiments may refer to a corresponding system, to a corresponding measuring device and as well as to a computer program product that can be loaded into the memory of at least one computer and comprises parts of software code that are able to execute the steps of the method when the product is run on at least one computer. As used herein, reference to such a computer program product is understood as being equivalent to reference to a computer-readable means containing instructions for controlling the processing system in order to coordinate implementation of the method according to the embodiments. Reference to "at least one computer" is evidently intended to highlight the possibility of the present embodiments being implemented in modular and/or distributed form.
The claims form an integral part of the technical teaching provided herein in relation to the various embodiments.
According to the solution described herein, the method comprises that said compensation operation includes obtaining from said selected data blocks corresponding frequency domain data blocks for the heart rate signal and acceleration signal, performing a motion compensation of the frequency domain data blocks for the heart rate using the frequency domain data blocks for the acceleration signal.
In various embodiments, the method comprises performing the motion compensation by subtracting from the frequency domain data blocks for the heart the frequency domain data blocks in the frequency by a scalar weight to obtain said compensated heart rate signal.
In various embodiments, the method comprises performing a detection operation on said compensated heart rate signal to obtain an estimate of the heart rate.
In various embodiments, the method comprises computing a non-linear predictor value on the basis of said frequency domain data blocks, performing a correction operation of the detected estimate of the heart rate including performing a decision between the detected estimate of the heart rate and a predicted estimate to select a heart rate value, said predicted estimate being obtained as a linear function of said predictor.
In various embodiments, the method comprises performing a decision which includes smoothing operations to obtain a final heart rate estimate.
In various embodiments, the system implementing the method comprises sensors configured for acquiring optically from said body organ a signal representative of the heartbeat and acquiring an acceleration signal representative of the acceleration of said body organ and a processor module configured for selecting data blocks of said acquired heart beat signal and acceleration signal, compensating said heart beat signal by the acceleration signal, calculating the heart rate value on the basis of said compensated heart beat signal.
In various embodiments, the system includes such sensors and processing module are comprised in a same photoplethysmographic heart rate measuring device.
In various embodiments, the PPG heart rate measuring device is comprised in an apparatus which is wearable on said body organ, in particular the apparatus is wearable on the wrist or the arm, in particular by means of a bracelet.
In various embodiments, it is provided a photoplethysmographic heart rate measuring device capable of operating in the system here described.
In various embodiments such device is associated to a remote processing device connected to it by a communication link.
In various embodiments such device includes a display.

Brief description of the drawings

The invention will now be described purely by way of a non-limiting example with reference to the annexed drawings, in which:

Figures 1A and 1B show a system and a device operating according to an embodiment of the method for estimation of the heart rate;
Figure 2 represents a block diagram of the operations performed by said method;
Figure 3 represents a diagram detailing an operation of the method described in figure 2;
Figures 4A and 4B detail a filtering operation of the method described in figure 2;
Figure 5 shows a diagram representing a model applied by the method described in figure 2;
Figures 6, 7 and 8 represents flow diagrams detailing operations of the method of Figure 2.

Detailed description of embodiments

The ensuing description illustrates various specific details aimed at an in-depth understanding of the embodiments. The embodiments may be implemented without one or more of the specific details, or with other methods, components, materials, etc. In other cases, known structures, materials, or operations are not illustrated or described in detail so that various aspects of the embodiments will not be obscured.
Reference to "an embodiment" or "one embodiment" in the framework of the present description is meant to indicate that a particular configuration, structure, or characteristic described in relation to the embodiment is comprised in at least one embodiment. Likewise, phrases such as "in an embodiment" or "in one embodiment", that may be present in various points of the present description, do not necessarily refer to the one and the same embodiment. Furthermore, particular conformations, structures, or characteristics can be combined appropriately in one or more embodiments.
The references used herein are intended merely for convenience and hence do not define the sphere of protection or the scope of the embodiments.
In figure 1A it is shown a block diagram representing a photoplethysmographic (PPG) device 10 for measuring the heart rate. Such PPG device 10 includes two lighting emitting diodes, i.e. LEDs, 13a and 13b, in particular green LEDs, which emit light in the direction of the user's skin surface. Such LEDs 13a and 13b are driven by a LED driver 15 which is in its turn controlled by a digital to analog converter 16 to control the current value. The skin surface is the skin surface of a body organ which in the preferred example is the wrist. Other suitable body organs can be for instance the arm, the forearm, the finger, the forehead or the ear. Light is absorbed, reflected, scattered, and transmitted as it travels through the user's tissues and encounters one or more blood vessels. Heart beats cause the amount of light reflected to drop during systoles and to increase during diastoles. The reflected light is sensed by a photodiode 17 which converts the amount of light into a current. The two LEDs 13a and 13b are placed respectively immediately above and below the photodiode 17. The current at the output of photodiode 17 is converted into a voltage by means of a trans-impedance amplifier 18. Finally the voltage is converted into a digital word by means of an analog-to-digital converter 18, for instance a 12-bit analog-to-digital converter 18 operating at a sampling rate f_s of 100Hz and acquired by a microcontroller 11, in particular a 32-bit microcontroller STM32L1 (or STM32F), as a optically acquired digital signal representative of the heart beat, or digital optical heart beat signal o. The microcontroller 11 also controls by means of one of its digital outputs connected to digital-to-analog converter 16 the current of the LEDs 13a and 13b, i.e. their light intensity.
The PPG device 10 includes also a three-axis accelerometer 21, which is in particular equipped with an internal analog-to-digital converter operating at the same sampling rate f_s, 100Hz, which provides three accelerometric signals to the microcontroller, namely ax, ay, az, according to axis x, y and z respectively.
An internal memory, in particular a SRAM, of the microcontroller 11, in particular the 48KB SRAM of the STM32L1 microcontroller, is used to implement the method for the estimation of the heart rate here described, indicated as a whole with reference 1000 in figure 2. Such method for the estimation of the heart rate receives as input the optically acquired heart beat signal o and the accelerometric signals ax, ay, az and supplies as output a series of heart rate estimate values r. If more memory is available, for example by using the alternative STM32F4 microcontroller, an operation of decimation during a filtering stage 110 shown in figure 2 is not necessary.
The PPG heart rate measuring device 10 in a preferred embodiment, shown in figure 1B, is arranged in a way which is known per se in a wrist-wearable device 40, including a wrist bracelet 41, with the LEDs 13a and 13b and the photodiode PD oriented in direction of the wrist skin surface. In figure 1B for simplicity's sake only LED 13a is shown, in dashed line.
The PPG device 10 is equipped with an own local display 12 to visualize the series of heart rate estimate values r, in particular, as better detailed in the following, the time series r(n) of the heart rate estimates r. The microcontroller 11 itself can further use the data of the heart estimate r for statistics and other type of data manipulation, for instance in order to supply different types of data presentations to the user, such as statistical analyses, the microcontroller 11 itself can further perform processing of the heart rate estimate r, acting as a local application processor if it has enough computing power. Also, the PPG device 10 comprises a transmitter 22, in the example a Bluetooth transmitter, to send wirelessly the series of heart rate estimate values r to a remote device 30, also shown in figure 1, which receives the series of heart rate estimate values r by means of a corresponding receiver 31, which sends the data to an application processor 32, which in its turn shows the results on a display 33. The remote device 30 can be for instance a PC, as shown in figure 1B, or a smart phone or another device capable of processing and displaying data.
Such option of having the time series of the heart rate estimate values r either visualized on a local display 12, used by the microcontroller 11 itself or sent wirelessly to the remote device 30 depends on the product/application that the system is integrated onto.
In figure 2 it is shown a block diagram representing an embodiment of a method for the estimation of the heart rate, indicated as a whole with the reference 1000, which is performed by the microcontroller 11.
The digital optical heart beat signal o acquired by the PPG device 10 and the digital accelerometric signals ax, ay, az pertaining the three axes x, y z are sent in parallel to a filtering stage 110, then to a selective block generation stage 120, then to a frequency analysis stage 130, which outputs a frequency domain heart beat signal O and frequency domain accelerometric signals AX, AY, AZ pertaining the three axes x, y, z respectively.
The purpose of the filtering stage 110 is to remove the unwanted frequency bands, that is, all the bands not included in the band of interest of the heart rate. Since the band of interest of the heart rate is, for example, the interval of frequencies between 40 bpm (beats per minute) and 200 bpm, the unwanted bands include a low-frequency band (between 0 bpm and 40 bpm) and a high-frequency band (between 200 bpm and half the sampling frequency). Both the low-frequency and the high-frequency unwanted bands need to be removed to clear the heart signal from noise. In addition, the low-frequency unwanted band needs to be removed because it contains a strong sub-band (between 0 bpm and a few bpm) that negatively interferes with the frequency-analysis stage. More specifically, since the frequency-analysis stage is not limited to orthogonal frequency points, but it calculates also non-orthogonal frequency points, the presence of a strong sub-band in the very low end of the spectrum would cause this sub-band to interfere with the calculation of the non-orthogonal frequency points. In addition, the high-frequency unwanted band needs to be removed in order to prevent aliasing, in case decimation is performed.
Each of such stages 110, 120, 130 include four banks in parallel to perform substantially the same operations respectively on the digital optical heart beat signal o and the three digital accelerometric signals ax, ay, az. Of course adjustements can be possible to take in account the different nature of the optically acquired heart signal and of the accelerometric signal. For instance the four banks of the filtering stages 110 can have the same structure, i.e. the one described with refererence to figures 4A and $B. This allows to have signals at the output of the banks of the filtering stage 110 with basically the same delay and provides simplicity of implementation or design. Nevertheless, in various embodiments the filtering banks can have differences to take in account differences among the input signals, provided they perform the task of removing the frequency bands as indicated above.
Downstream the frequency analysis stage 130 the frequency domain heart beat signal O and the frequency domain accelerometers signals, AX, AY, AZ are fed to an estimation module 200 which includes a motion compensation block 210 and a detection block 220. The motion compensation block 210 compensates the frequency domain heart beat signal O using the values of the frequency domain accelerometers signals, AX, AY, AZ. The compensated output, O', of such motion compensation block 210, namely a compensated frequency domain heart beat signal O' is fed to the detection block 220 which computes a detected heart rate estimate d, i.e. identifies the main frequency of the compensated frequency domain heart beat signal O', which is then outputted as detected heart rate estimate d.
Such detected heart rate estimate d is fed to a correction module 400, performing a model estimation, which also receives the frequency domain accelerometers signals, AX, AY, AZ. In the correction module 400, the detected heart rate estimate d is sent as input to a decision block 410 which decides which value, between the detected heart rate estimate d and a predicted estimate m, send as decided heart rate r. The succession of decided heart rate values r form the series of heart rate estimate values r, i.e. the final output of the method described by the block diagram of figure 2.
Also, a non linear prediction calculation block 300 receives the frequency domain accelerometers signals, AX, AY, AZ and calculates a predictor p. Such predictor p is fed to the correction module 400, specifically to a model learning block 430 which also receives from the decision block 410 the detected heart rate estimate d. The model learning block 430, on the basis of the predictor p and of the frequency domain accelerometers signals, AX, AY, AZ calculates linear parameters α and β, respectively the constant term and the first degree coefficient of a linear function, which are fed to a linear function calculation block 420, together with the predictor p from the non linear prediction calculation block 300. The linear function calculation block 420 outputs the predicted estimate m on the basis of the linear function, as better detailed in the following. Thus, summarizing, a parametric model in the form of a linear function calculates a predicted estimate m of the heart rate from the predictor p and the linear parameters α and β.
As mentioned, the decision block 410 chooses which value between the detected heart rate estimate d and the predicted estimate m represents the decided heart rate r. The decision block 410 in various embodiment can simply select to send either the detected heart rate estimate d or the predicted estimate m as output. However, in a preferred embodiment, better detailed in the following the decision block 410 performs a decision which operate in a smoothed manner, to filter out abrupt changes due to the transitions between the two estimators, the detected heart rate estimate d, which represent a first optical estimator, and the secondary predicted estimator, m, supplying as an output the decided heart rate estimate r and, ultimately, the series of heart rate estimate values r.
When the decision block 410 operates its choice, if the optical estimate d is selected, such optical estimate d and the predictor p evaluated at the same time are stored, in particular in respective FIFO memories c_x , c_y, the pair of values identifying an observation point (p,d) on a predictor-estimate diagram, as better shown with reference to figure 5, which is used to refine the model in block 430. Figure 5 shows a scatter plot of such observation points, (p,d), i.e. the detected heart rate d and the predictor p values. The line shown in figure 5 corresponds to predicted estimate m=α + β * p.
The method attempts at learning the value of linear regression parameters α and β on the basis of points (p, d). For efficiency reasons, the number of points on the basis of which the linear regression parameters α and β are maintained limited in number, since the number of the observation points (p, d) grows rapidly, one point being stored each time the decisor 410 chooses to select the estimate d. This is obtained by a procedure which, as better detailed in the following, includes collapsing the observation points (p, d) which are sufficiently similar or near one to the other in a single point, and includes substituting with the most recent observation points the oldest observation points, using the FIFO data structures. The size of such FIFO structures is equal to the maximum number of points on the basis of which the regression line is traced. The observation points are no longer of the type of pairs (p, d), but collapsed unique elements contained in the FIFO (cx, cy).
The output of the method for the estimation of the heart-rate 1000, directed to the local display 12, the microcontroller 11 or the remote device 30, is the numerical time series of estimates r, r(ns), ns being the index of the estimate produced. For instance a decided estimate r is produced every three seconds, so that the index ns is increased every three seconds. The method for the estimation of the heart-rate 1000 makes an attempt to update such series r, with a tentative update period of the estimate T_e, i.e. every T_e seconds from the last estimate, for example T_e=5s or T_e=1s. Occasionally the updating attempt fails, due to poor signal conditions; in this case the subsequent update will be provided after a period longer than the update period T_e .
With reference to the operations shown in figure 2, in various embodiments such method not necessarily includes all the operations and modules shown in such figure and describe above, or, alternatively, some of the operations and module can be different.
By way of example, in various embodiments, the method for the estimation of the heart-rate using photoplethysmography on a body organ, in particular on a wrist of a user, comprises acquiring optically from sensor 13a, 13b, 17 from said body organ a heart beat signal o, acquiring from sensor 21 an acceleration signal ax, ay, az representative of the acceleration of said body organ, selecting by the stage 120 data blocks B_o, B_x, B_y , B_z of said acquired heart beat signal o and acceleration signal ax, ay, az, compensating in block 210 said heart beat signal o by the acceleration signal ax, ay, az, calculating in block 210, 220 a heart-rate value, the heart rate estimate d on the basis of the resulting compensated heart beat signal O', the compensation operation 210, 220 including obtaining in stage 130 from said selected data blocks Bo, Bx, By, Bz corresponding frequency domain data blocks for the heart rate signal O and for the acceleration signal AX, AY, AZ, performing a motion compensation 210 of the frequency domain data blocks for the heart beat O using the frequency domain data blocks for the acceleration signal AX, AY, AZ.
In various embodiments, the correction module 400 operates on the basis of the heart rate estimate d and of the predictor p originated by the non linear prediction calculation block 300 to produce the decided heart rate value r.
In various embodiments the decided heart rate value r is calculated by a smoothed decision operation (as detailed in figure 8).
Now the operation of the specific blocks and stages shown in figure 3 will be detailed.
As regards the banks of the filtering stage 110, each comprises a low-pass FIR (Finite Impulse Response) filter 111, as shown in figure 4A, followed by a high-pass filter 113. Between the low-pass filter 111 and the high pass filter 113 a factor-2 decimation stage 112 is optionally applied, as in the example shown, in order to reduce the memory footprint of the algorithm. The filters are independently applied to each of the four input signals: the optical sequence of values corresponding to the optical signal o, o(n), n being as mentioned the numeric index of the samples acquired, and analogously the three accelerometric sequences ax(n), ay (n) and az(n), corresponding to signals ax, ay and az. The filters of filtering stage 110 produce respectively four filtered time domain signals, one, O_HP(n), pertaining the optical sensor 17, and the other three pertaining each axis of the accelerometric sensor 21, ax_HP(n), ay_HP(n), az_HP(n) respectively.
Since the required transfer function of the high-pass filter 113 is very stringent, preferably the high pass filter 113 is implemented by means of a chain of five IIR (Infinite Impulse Response) filters 113a, as illustrated in figure 4B.
Each of the five IIR filters 113a have the following transfer function: $H (z) = \frac{b_{0} + b_{1} \cdot z^{- 1} + b_{2} \cdot z^{- 2}}{1 + a_{1} \cdot z^{- 1} + a_{2} \cdot z^{- 2}}$

where

a ₁= -1.905850154697027
a ₂=0.911145313595591
b ₀=0.979650817189758
b ₁=-1.907660496005129
b ₂=0.929684155087730

The transfer function resulting from the cascade of the five IIR filters accomplishes the removal of the unwanted low-frequency band as explained before, especially the very first portion of the unwanted low-frequency band.
For what regards the selection block generation stage 120, this stage comprises four parallel banks receiving as input each a respective sequence corresponding to one of the four filtered time domain signals coming from the sensors, O_HP(n), ax_HP(n), ay_HP(n), az_HP(n), and generating four time domain data blocks by operating on such respective filtered optical and accelerometric sequences. The four generated time-domain data blocks, B _o from the data in the optical sequence O_HP(n), and B_x, B_y , B_z from the data in the accelerometric sequences ax_HP(n), ay_HP(n), az_HP(n), are time synchronized, which means that the first samples in each bank of the stage 120 are time aligned to the same temporal reference, which marks the beginning of the time domain data blocks. The selection of the time domain data block is done by means of a procedure, i.e. the selective block-generation better detailed below, that repeatedly tests the data blocks tentatively constructed on the time domain optical sequence o(n) until a given condition is met. The procedure consists in provisionally generating a data block from the optical sequence o(n), testing the condition, and, in case the condition is met, exiting the procedure and forwarding that block B_o to the next stage of the method (along with the three time-aligned accelerometric data blocks B_x, B_y , B_z ), while, in case it is not met, generating a new provisional block that is time shifted relative to the previous block. In this way, the new provisional block contains a sub-block with more recent samples relative to the previous block.
Specifically, the selective block generation procedure operates as follows and as detailed in figure 6, where no indicates an index of the most recent optical sample of the optical sequence o(n) stored in memory out of the filtering stage 110 and B indicates a data block formed from the N most recent optical samples, from no backwards:

a de-trending step 610 is first operated on the block B. A block is said to be de-trended if its samples are replaced with the result of the subtraction of the mean value of its original samples from each original sample. In this way, the mean of the samples of the de-trended block equals zero;
then in a step 620 it is calculated a ratio S on the de-trended block B as the maximum value of B over the mean of the de-trended data block B. N is the number of samples in the block. $S = \frac{\underset{n = 0 \dots N - 1}{\max |B (n)|}}{mean |B (n)|}$
if in a step 630 test condition is evaluated that such ratio S is lower than a determined threshold value S_max:
- ∘ in a step 640 the stage 120 gives as output the optical data block B₀, set equal to block B, and the three time-aligned and de-trended accelerometric blocks B_x , B_y , B_z , which are forwarded to the frequency-analysis stage 130;
- o it is then waited in a step 650 until the filtering stage 110 outputs the sample n₀+N_e and repeated the above procedure starting from the de-trending step 610 to generate the next blocks;
else, if the condition 630 is not met, in a step 660 it is waited until the filtering stage 110 outputs the sample n₀ +N_w (the new most recent optical sample) and then it is repeated the procedure from the de-trending step 610, thus generating a new provisional block that is time shifted relative to the previous block.

The parameter N indicating the number of most recent optical samples taken into account is set to 1024, but other values are also possible (although with some effects on performance and memory requirements).
The throughput parameter N_e is set according to the required throughput of heart-rate estimates, i.e. update period, T_e : $N_{e} = T_{e} \cdot \frac{f_{s}}{2}$
The factor 2 in the denominator of the previous formula refers to the case of the sampling-rate reduction through decimation.
For example, in case a throughput of 0.2Hz is required (one estimate every 5s), and the sample frequency f_s is 100Hz, the throughput parameter N_e is set to 250 samples. The value of the parameter N_w can be set to any amount smaller than throughput parameter N_e, for example to 50 samples.
For what regards the frequency analysis stage 130, this stage operates independently on each of the four blocks B_o, B_x, B_y , B_z at the output of the selective block generation stage 120. The four time domain data blocks B_o, B_x, B_y , B_z pass through the same type of processing. Each of the four time domain data blocks B_o, B_x , B_y , B_z is transformed into the frequency domain in the form of a number K of frequency lines, or frequency bins. The K bins are obtained by interleaving a number M of N-point FFTs (Fast Fourier Transform) as shown in reference to figure 3. In this way a higher (better) frequency resolution can be obtained, as compared to a single FFT. Not all of the N points are necessary for each FFT, but only a subset of K/M points, as described next. K it is assumed to be an integer multiple of M. Each of the M FFTs are calculated respectively on blocks B^(m), with m=0...M-1 obtained by multiplying in a multiplier block 132 the samples of the block U(n) by a complex-valued sequence: $B^{(m)} (n) = U (n) e^{\frac{j 2 π n Δ m}{f_{s}}}$

where m=0...M-1, n=0...N-1 and U(n) is obtained by multiplying in a multiplier block 132 the samples of the generic block B(n) by an apodization function (or window function), for example a Taylor function. The use of apodization functions is well known in the prior art and has the purpose to attenuate the artifacts due to the secondary lobes introduced by the operation of block selection.
The complex argument of the exponential includes the product of index m identifying the m-th FFT, index n identifying the sample and a frequency resolution of the frequency analysis Δ, over the sampling frequency f_s.
The m-th FFT, calculated in block 133 on the block B^(m) is denoted with F_m (k), while S_m (k), indicates the squared modulus of the m-th FFT: S_m (k)=|F_m (k)|² , also calculated in block 133.
The output of the frequency analysis 130 is a block S of frequency domain values, formed with such squared modulus of the FFT S_m (k) by a multiplexer 134 as follows: $S = \{S_{0} (k_{0}), \dots, S_{M - 1} (k_{0}), S_{0} (k_{0} + 1), \dots, S_{M - 1} (k_{0} + 1), \dots, S_{0} (k_{0} + \frac{K}{M} - 1), \dots, S_{M} (k_{0} + \frac{K}{M} - 1)\}$

where k ₀ is the index corresponding to the lower bound of the frequency band of interest, as described later.
The K elements of S are numbered from 0 to K-1 so that: $S (v) = S_{m} (k_{0} + q)$

where v=qM+m, with q and m univocally determined by the conditions: m<M, q integer, and m integer.
Each element S(v) corresponds therefore to the frequency (k ₀+q)·M·Δ+m·Δ, where the quantity Δ is the frequency resolution: Δ=1/(M*D), and D is the duration in seconds of the analysis block. By way of example, if D=20.48s, and M=3, then Δ =0.0162
In order to determine the value of K, that is the number of elements in the set S, it is possible to proceed in the following way. First, a band of frequencies of interest is specified, for example the band 40 bpm - 200 bpm. Secondly, the index k₀ , is determined by imposing $S (0) < \frac{40}{60},$
which leads to $k_{0} \cdot M \cdot Δ < \frac{40}{60},$
and solving with respect to k ₀, it is obtained k ₀=13.
Finally, it is imposed $S (K - 1) > \frac{200}{60},$
and solving it is obtained K=168.
Finally the block S is normalized. Within the context of the present description, the normalization of a block X containing K elements produces the block X^u defined as follows: $X^{u} (k) = \frac{X (k)}{s}, k = 0, \dots, K - 1$

where $S = \sum_{k = 0}^{K - 1} X (k)$
The normalized block S is thus indicated by S^u.
A normalized block S^u is constructed on each of the four time domain data blocks B_o, B_x, B_y , B_z, and these normalized blocks, which are frequency domain data values, are denoted respectively with O^u , A_x ^u, A_y ^u, A_z ^u .
For what regards the motion compensation 210, this operation consists in calculating a motion compensated block O^MC defined as follows: $O^{M C} (k) = O^{u} (k) - [v_{X} \cdot A_{X}^{u} (k) + v_{Y} \cdot A_{Y}^{u} (k) + v_{Z} \cdot A_{Z}^{u} (k)], k = 0, \dots, K - 1$

Where v_X, v_Y, v_Y are three scalar weight values, for example set to 1; O^u(k), A_x ^u(k), A_y ^u(k), A_z ^u(k) are the normalized frequency domain data block outputted by the frequency analysis stage 130 respectively for the optical and the three accelerometric sequences.
The derivation of the motion-compensation formula can be made proceeding in the following way.
The signal O^u (k) can be decomposed into five contributions: a heart-beat signal O_h (k) (the wanted signal), three motion-induced artifact signals, O_x (k), O_y (k), and O_z (k), and a non motion-induced noise O_n (k). So it is possible to write: $O^{u} (k) = O_{h} (k) + O_{x} (k) + O_{y} (k) + O_{z} (k) + O_{n} (k)$

The artifact components, being caused by motion, can be written in terms of transfer function of the accelerometric input: $O_{X} (k) = V_{X} (k) \cdot A_{X}^{u} (k)$
$O_{Y} (k) = V_{Y} (k) \cdot A_{Y}^{u} (k)$
$O_{Z} (k) = V_{Z} (k) \cdot A_{Z}^{u} (k)$

where V_X(k), V_Y (k), and V_Z (k) are three transfer functions defined as: $V_{X} (k) = \frac{O_{X} (k)}{A_{X}^{u} (k)}$
$V_{Y} (k) = \frac{O_{Y} (k)}{A_{Y}^{u} (k)}$
$V_{Z} (k) = \frac{O_{Z} (k)}{A_{Z}^{u} (k)}$
Perfect motion compensation would require the knowledge of the three transfer functions. Unfortunately, the transfer functions are unknown for the following reasons:

they are time variant;
they are motion dependent;
a model for their estimation is not available and it is not easy to derive;
they are subject dependent.

For the reasons listed above, the transfer functions have been approximated to a constant value both in the time and in the frequency axis, denoted respectively with v_X, v_Y, and v_Z . Despite this approximation, the performance of the motion compensation stage is still acceptable in that the purpose of motion compensation is not the complete (over the whole frequency band) estimation of O_h (k), but, rather, is the estimation of only the position of the dominant frequency component of O_h (k), as it will be described next. The optimal values for the three constant (approximated) transfer functions can be obtained, for example, through an optimization process that maximizes the performance of the final heart-rate estimation with respect to v_X' v_Y, and v_Z. The values v_X = 1,v_Y = 1, and v_Z = 1 offer a first guess of the optimal solution.
Another important aspect of the motion compensation stage is the operation of normalization applied to the four signals (one heart signal and three accelerometric signals). The result of this operation is to equalize the energy of the four signals to the unitary value. This operation is necessary in order to bring the magnitude of the spectra of the four signals to the same range of values. Failure to do this would cause the subtraction operation to be not effective, in that it would be performed on terms of non-comparable magnitudes. Through the operation of normalization, a strong frequency component present in the spectra of any of the accelerometric signals can effectively cancel the corresponding artefact on the optical signal.
Thus, the motion compensation substantially includes subtracting from the frequency domain optical signal, O^u , i.e. the heart beat signal in the frequency domain, a linear combination of the three frequency domain acceleration signals A_x ^u(k), Ay^u(k), A_z ^u(k) multiplied by a respective scalar weight v_X , v_Y, v_Y.
For what regards the detection operation 220, the detection consists in finding a argument k_max of the maximum value of the motion compensated block O^MC(k): $k_{\max} = \arg \max_{k} (O^{M C} (k))$
The frequency bin k_det of the detected heart-rate d can then be obtained from such argument k_max as the sum of the index of the first useful bin k₀ and argument k_max : $k_{d s t} = k_{0} + k_{\max}$
Optionally, the estimate of the frequency bin k_det can be refined by a 3-point interpolation: $k_{d s t} = k_{d s t} + f_{int} (O^{MC} (k_{\max} - 1), O^{MC} (k_{\max}), O^{MC} (k_{\max} + 1))$
where f_int (y ₁,y ₂,y ₃) is any interpolation function such that: $f_{int} (y) (_{1} y_{2} y_{3}) = {\begin{matrix} \in (- 1, 0) & i f y_{1} > y_{3} \\ 0 & i f y_{1} = y_{3} \\ \in (0 1) & i f y_{1} < y_{3} \end{matrix}$
Finally, the detected heart rate estimate d is obtained by performing a rescaling the frequency bin k_det of the detected heart-rate d: $d = 60 \cdot Δ \cdot k_{d s t}$
For what regards the predictor calculation in block 300, the frequency analysis stage 130 produces a set of four frequency-domain data blocks at output times t_l =t _l-1+τ _l , where τ _l is the interval between two consecutive calculations and 1 is the index of the output times. This interval τ _l is usually equal to the update period T_e , but it can be longer in case the selective block-generation stage 120 enters in the step 660 of the selective block generation procedure in which it is waited until the filtering stage 110 outputs the sample n₀ +N_w . The l-th set of four frequency domain blocks, calculated at time t_l, is denoted as O^u (n;l), $A_{X}^{u} (n l),$
$A_{Y}^{u} (n l),$
and $A_{Z}^{u} (n l) .$
The calculation of the predictor in block 300 consists of the following procedure comprising the following steps:

a step 710 of summing the normalized frequency domain blocks of the three accelerometric signals to obtain a sum of accelerations A(k;l) $A (k l) = A_{X}^{u} (k l) + A_{Y}^{u} (k l) + A_{Z}^{u} (k l)$
a step 720 of finding the maximum A_max(l) of the sum of accelerations A(k;l) and the corresponding argument n_max: $A_{\max} (l) = \max_{k} A (k l)$
$n_{\max} (l) = \arg \max_{k} A (k l)$
a step 730 of calculating a time sequence q(l) as function of the maximum A_max(l) of the sum of accelerations A(k;l) and of the corresponding argument n_max as follows: $q (l) = 0, 1 \cdot (60 \cdot Δ \cdot (k_{0} + n_{\max} (l)) + 10 \cdot Log (A_{\max} (l)))$
a step 740 of filtering the time sequence q(l) by means of two first-order IIR filters as follows: $p_{f} (l) = α_{f} (l) \cdot p_{f} (l) + (1 - α_{f} (l)) \cdot q (l)$
$p_{s} (l) = α_{s} (l) \cdot p_{s} (l) + (1 - α_{s} (l)) \cdot q (l)$

where: $α_{f} (l) = {0.99}^{5 \cdot τ_{l}}$
$α_{s} (l) = {0.998}^{5 \cdot τ_{l}}$
a step 750 of selecting as predictor at time t_l, p(l), the maximum between the first filtered sequence p_f (l) and the second filtered sequence p_s(l): $p (l) = \max (p_{f} (l), p_{s} (l))$

The previous formula represents the fact that the predictor p(l) is constructed as the maximum between two intermediate functions, p_f(l) and p_s(l). These two intermediate functions are two filtered versions of the time series of q(l). If to the time sequence q(l) is given the interpretation of the level of physical exertion of the user, the two filtered output of the time sequence q(l) represent two different responses to physical exertion. In particular, the series p_f(l) represents a fast response to the time sequence q(l), while p_s(l) represents a slow response to the time sequence q(l). Through the operation of maximum search max, the series p_s(l) becomes relevant during periods of physical recovery (that is when the level of physical exertion is decreasing), while the series p_f(l) becomes relevant during periods of increasing physical activity. This is consistent with the fact that cardiac output (and therefore cardiac frequency) responds differently to increases and decreases of physical activity. Typically, heart rate responds quickly to an increase of physical activity, and responds much more slowly to a decrease of physical activity.
For what concerns the correction stage 400, this stage uses the values of the predictor p(l) and of the detected heart rate estimate d(l) to arrive at a heart-rate estimate r(l).
The correction procedure 400 includes the following steps:

a step 805 in which is verified if an initialization of values is required. This occurs usually only at the power-on of the system. In the affirmative (Y), a step 810 of initialization of the values is performed. Values are by way of example initialized as follows:
- l=1
- α= α₀
- β= β₀

In the following for simplicity of representation index 1 will be used instead of time t_l.
Three FIFO memories, c_X, c_Y, which store, as mentioned, the observation points (p,d) each time the estimate d is selected, and a weight computation FIFO c_W which is used in the aggregation of the points. FIFO memories c_X, c_Y, c_W as described in the following, are initialized as empty.
In particular, the predictor FIFO memory c_x sequentially stores the values of the predictor p at each time instant l at which the estimate d is selected.
The estimate FIFO memory c_y sequentially stores the values of the estimate d at each time instant l.
The weight FIFO memory c_w increments the values of a weight factor at each time instant l.
α₀ and β₀ are initial guesses set for the first order term and for the constant term of the linear regression, for example, α₀=80 and β₀=0.8. Other values are suitable, although they have an effect on the accuracy of the final heart rate estimate;

after step 810 or if the condition verified at step 805 is negative (N), i.e. initialization has been already performed a step of evaluation 820 of a condition c₁(l) at time l is performed: $c_{1} (l) = (d (l) < d_{\min} (l)) o r (d (l) > d_{m a x} (l))$

where: $d_{\min} (l) = m (l) - G$
$d_{\max} (l) = m (l) + G$
$m (l) = α + β \cdot p (l)$

G indicates a safety margin set for instance to 40, therefore d_min and d_max represents two lines parallel to m(l) which identify the lower and upper limit of a safety range. If the value estimate d is outside such range is considered not a good value.
If the condition c₁(l) is true then the estimate r(l) is, in block 830: $r (l) = \frac{r (l - 1) + m (l)}{2}$
The correction is finished, and the procedure wait in a step 840 for the next iteration, updating the time index l, thus it is l=l+1 and going back to operation 820.
Else if condition c₁(l) is false, the condition c₂(l) is evaluated at time l in a step 850: $c_{2} (l) = (d (l) - d (l - 1) < g_{1}) o r (d (l) - d (l - 1) > g_{2}) and not c_{1} (l - 1)$
g₁ and g₂ are parameters set for instance to g₁=-10 and g ₂= 20 which set a confidence range [g1, g2] for the variation between time adjacent estimates. If the estimates d(l) are within said range, c₂(l) is false, meaning that the estimate d is reliable and can be selected as valid output (see step 870 below).
Other values are possible, although they may affect the accuracy of the heart rate estimates.
If condition c₂(l) is true, the final estimate r(l) is evaluated in a step 860 as the median of the latest three detected estimates d(l-2),d(l-1),d(l): $r (l) = median (d (l - 2), d (l - 1), d (l))$
It is noted that if c₂(l) is true the optical estimate d is considered in any case at least partly reliable, since it is used to obtain the final estimate r(l), although through the median function. If c₂(l) is false, the reliability of the optical estimate d(l) is greater, thus there is no need of taking in account the two previous estimates d(l-1), d(l-2) and the estimate d(l) is sent directly as output r(l). Thus operation 860, within the decision operation performed by block 410, represents a smoothing operation applied taking in account a plurality of estimates d at different times, such as l-2,...,l and averaging over them.
Also in this case the correction procedure is finished and the procedure waits in a step 840 for the next iteration, updating the time index l, thus it is l=l+1 and going back to operation 820.
If condition c₂(l) is false, the final estimate is set in a step 870 to the detected estimate d(l): $r (l) = d (l)$
In figure 8 steps 810-870 are indicated by a dashed polygonal as operated specifically by the decision block 410. As mentioned, the decision block 410 in its simpler form decides if the final estimate r(l) is equal to the optical estimate d(l) or to the predicted estimate m(l) on the basis of a given condition, such as c₁(l) or c₂(l)), which evaluates if the optical estimate value is reliable. In the embodiment of figure 8, however, a smoothed decision is implemented, i.e. the final estimate r(l) is at least in some cases the result of an average of the current estimate with previous estimates, optical estimate d and or predicted estimate m. Specifically, while step 870 output simply the optical estimate d(l) as final estimate r(l), step 860 outputs an average over previous optical estimates, step 830 outputs an average with the previous final estimate r(l).It is clear that in various embodiments other choices of reliability conditions such as c₁(l) or c₂(l)) and other type of averaging are possible to the person skilled in the art.
The results of either steps 830, 860 or 870 are collected in a block 896 which supplies the final estimate r(l) outputted by the decision block 410.
Preferably, the update of the linear block 430 includes performing an aggregation operation, described here below, to keep low the number of observation points in the FIFO. Thus, in figure 8 it is indicated, after step 870, a step 881 of evaluating if aggregation is required. In the negative, after the step 870, since the decision block 410 has chosen the detected estimate d(l), the linear model in block 430 is updated, this meaning that observation point (p, d) are stored in a step 875 in the FIFO c_X, c_Y . Then step 896 is performed, outputting the final estimate r(l).
If step 881 indicates that an aggregation procedure must be performed, it is executed a step 882 of finding an element a in the predictor FIFO c_x, storing a time sequence of predictors p, closest to the predictor p(l): $a = \min_{k} |p (l) - c_{x} (k)|$
and evaluating if a<a _max, i.e. the closest element a in the predictor FIFO c_x is closer than a threshold value to the predictor p at time l. This corresponds to evaluate if the new observation point is near to an existing point stored in the FIFO.
In the affirmative in a step 883 it is set: $j_{0} = \arg \min_{j} |p (l) - c_{x} (j)|$
$w = \max (0.2, 1 - 0.1 \cdot c_{w} (j_{0}))$
$c_{x} (j) (_{0}) = c_{x} (j) (_{0}) + w \cdot (p (l) - c_{x} (j_{0}))$
$c_{y} (j) (_{0}) = c_{y} (j) (_{0}) + w \cdot (d (l) - c_{y} (j_{0}))$
$c_{w} (j) (_{0}) = c_{w} (j) (_{0}) + 1$
j₀ is a parameter indicating the value of index j of the value c_x(j) in the predictor FIFO c_x which minimizes the distance with the predictor p(l) at time index l, which identifies an aggregating point (c_x (j₀ ),c_y (j₀ )).
c_w is in general a weight factor FIFO. c_w (j₀ ) represents a weight factor of the aggregating point (c_x (j₀ ),c_y (j₀ )), i.e. a FIFO counter which is incremented by one at each aggregation of a new observation point (p,d) to such aggregating point(c_x (j ₀),c_y (j ₀)).
The operation of aggregation or collapse in a single point is a weighted average performed on both the values in the predictor FIFO c_x and in the estimate FIFO c_y . The aggregating points (c_x (j₀ ),c_y (j₀ )) with the greater weight are those which are the results of a greater number of aggregations in the past. The weight of the aggregating point is used to establish a weight of the newly aggregated point (p,d), which is represented by w. Thus, the relations above implements the operations of assigning to the new point or aggregated (p, d) a relative weight which is inversely proportional to the weight of the aggregating point ((c_x (j ₀),c_y (j ₀)) to which is aggregated.
In this way, if the aggregating point (c_x (j ₀),c_y (j ₀)) is reliable, the aggregated point (p,d) has not to change significantly its value.
Else if condition 882 is not met, it is performed a push operation on the content of the FIFO, before passing to step 875. $c_{x} = push (p (l), c_{x})$
$c_{y} = push (p (l), c_{y})$
$c_{w} = push (1 c_{w})$
The function push(val, list) first shifts one position rightwards all the elements contained in the list (in this way discarding the rightmost element) and then inserts the number val at the leftmost position of the list.
Then in a step 890, following step 883, are calculated the linear regression parameters α₁ , β₁, R₁ of the set of points {(c_x (j),c_y (j))}, for j=0,...,F(l)-1, F(l) being the number of elements in the FIFO memory at time index l.
R indicates the coefficient of determination computed in the classical mode from the linear regression (it is R² properly, R for short). The coefficient of determination R measures the amount of linearity of the aggregated poins, which are used as reliability indication. R1 indicates therefore the coefficient of determination computed over the current aggregated points, i.e. the points (c_x, c_y ) in the FIFO structures, while Rmin indicates a minimum threshold level above which the linear model is considered as reliable.
In a step 892 it is evaluated if the number of elements in the FIFO F(l) is at least a minimum value F _min and if the coefficient of determination R is greater equal than a minimum value R_min.
In the affirmative, in a step 893 the linear parameters are updated, setting them to the linear regression parameters α₁, β₁ calculated in the step 890: $α \leftarrow α_{1}$
$β \leftarrow β_{1}$
Else, in a step 894 the linear parameters are set to the initial values of initialization step 810: $α \leftarrow α_{0}$
$β \leftarrow β_{0}$
The procedure goes to step 840 and then to step 820 to wait for the next iteration.
The parameters F_min and R_min are set, for example to F_min =5 and R_min =0.5. Other values are possible but performance may vary.
Steps 890-894 are indicated by a dashed square as operated specifically by the model learning block 430 to update the linear regression parameters on the basis of the observation points storing step 875. As indicated steps 880-884 represents an accessory aggregation procedure. In various embodiments it is also possible to proceed from the observation points storing step 875 directly to linear parameters calculation step 890.
The solution according to the various embodiments here described allows to obtain the following advantages.
The method according to the various embodiments here described is advantageous since it allows to overcome noise by means of high processing gain (for instance, with 1024 points FFT is obtained a 30dB gain).
The method according to the various embodiments here described is also advantageous since it allows to overcome interference due to motion of the body organ in the frequency domain (motion compensation).
The method according to the various embodiments here described is also advantageous since it to build a secondary estimator, i.e. a prediction, to be used whenever the primary estimator derived by the optical signal fails.
Of course, without prejudice to the principle of the embodiments, the details of construction and the embodiments may vary widely with respect to what has been described and illustrated herein purely by way of example, without thereby departing from the scope of the present embodiments, as defined the ensuing claims.
The method here described can be used in system and device for the continuous monitoring of the heart-rate for fitness/wellness applications. These systems and devices can be for instance wrist-wearable devices for runners monitoring their workout session. These devices can be wristwear or also earwear, smart-watches, smart-headphones.

Claims

Method for the estimation of the heart-rate using photoplethysmography on a body organ, in particular on a wrist of a user, comprising acquiring optically (13a, 13b, 17) from said body organ a signal representative of the heart beat (o), acquiring (21) an acceleration signal (ax, ay, az) representative of the acceleration of said body organ, selecting (120) data blocks (B_o, B_x, B_y , B_z) of said acquired heart beat signal (o) and acceleration signal (ax, ay, az), compensating (210) said heart beat signal (o) by the acceleration signal (ax, ay, az), calculating (220; 410; 500) a heart-rate value (d; r) on the basis of said compensated heart beat signal (O'),
characterized in that said method includes obtaining (130) from said selected data blocks (Bo, Bx, By, Bz) corresponding frequency domain data blocks for the heart beeat signal (O) and for the acceleration signal (AX, AY, AZ),
said compensating operation includes performing a motion compensation (210) in the frequency domain, compensating the frequency domain data blocks for the heart beat (O) with the corresponding frequency domain data blocks for the acceleration signal (AX, AY, AZ).
The method according to claim 1, characterized in that said motion compensation (210) includes subtracting from the frequency domain data blocks for the heart beat (O) the frequency domain data blocks (AX, AY, AZ) multiplied by a respective scalar weight (v_X, v_Y, v_Z ) to obtain said compensated heart beat signal (O').
The method according to claim 1 or claim 2, characterized by performing a detection operation (220) on said compensated heart beat signal (O') to obtain an estimate of the heart rate (d).
The method according to any one of the preceding claims, characterized in that includes computing a non-linear predictor value (p) on the basis of said frequency domain data blocks (AX, AY, AZ),
performing a correction (400) operation of the detected estimate of the heart rate (d) including
performing a decision (410) between the detected estimate of the heart rate (d) and a predicted estimate (m) to select a decided heart rate value (r),
said predicted estimate (m) being obtained (420) as a linear function (α, β) of said predictor (p).
The method according to claim 4, characterized in that said performing a decision (410) between the detected estimate of the heart rate (d) and a predicted estimate (m) to select a decided heart rate value (r) includes a smoothing operation (850, 860) applied taking in account a plurality of estimates (d, m, r) at different times and averaging over them.
The method according to any of the preceding claims, characterized by updating (430) said linear function (α, β) by storing the current detected estimate (d) and predictor value (p) each time the operation of performing a decision selects the optical estimate (d), obtaining a sequence of observation points (p,d), in particular stored in a FIFO (cx, cy), calculating the linear function (α, β) as regression over said sequence of observation points (p,d).
The method according to claim 6, characterized by performing an aggregation operation (880, 882, 884) over said sequence of observation points (p,d) to keep low their number.
The method according to any of the preceding claims, characterized by performing a filtering (110) before said selecting (120) data blocks (B_o, B_x, By, B_z ) of said acquired heart beat signal (o) and acceleration signal (ax, ay, az), in particular said filtering (111) comprising a low-pass FIR (Finite Impulse Response) filter (111) and a high-pass filter (113), in particular implemented by a chain of IIR (Infinite Impulse Response) filters (113a).
The method according to any of the preceding claims, characterized in that said selecting (120) data blocks (B_o, B_x, B_y , B_z) of said acquired heart beat signal (o) and acceleration signal (ax, ay, az) includes provisionally generating a data block (B_o) from the acquired heart beat signal (o), testing (630) a condition (S_max ) on said block (B_o ) and, in case the condition is met (640) forwarding said block (B_o) and the data blocks (B_x, B_y , B_z) of said acquired acceleration signal (ax, ay, az), time-aligned with said data block (B_o) from the acquired heart beat signal (o) to the operation (130) for obtaining corresponding frequency domain data blocks for the heart beat signal (O) and for the acceleration signal (AX, AY, AZ), while, in case said condition is not met (660), generating a new provisional block that is time shifted (N_w) relative to the previous block.
The method according to any of the preceding claims, characterized that said motion compensation operation (210) is followed by a detection operation (220) to identify a main frequency of the compensated frequency domain heart beat signal (O') which is outputted as detected heart rate estimate (d).
The method according to any of the preceding claims, characterized that said operation of obtaining (130) from said selected data blocks (Bo, Bx, By, Bz) corresponding frequency domain data blocks for the heart beat signal (O) and for the acceleration signal (AX, AY, AZ) includes obtaining a plurality (K) of bins by interleaving a plurality (M) of Fast Fourier Transform operations.
A system for the estimation of the heart-rate using photoplethysmography on a body organ, in particular a wrist of a user, comprising sensors (17, 21) configured for acquiring optically (13a, 13b, 17) from said body organ a heart beat signal (o) and acquiring (21) an acceleration signal (ax, ay, az) representative of the acceleration of said body organ and a processor module (11) configured for selecting data blocks (Bo, Bx, By, Bz) of said acquired heart beat signal (o) and acceleration signal (a_x , a_y , a_z ), compensating said heart beat signal (o) by the acceleration signal (ax, ay, az), calculating the heart rate value (r) on the basis of said compensated heart beat signal (O'),characterized in that said module (11) is configured to perform the method according to any of claim 1 to 9.
The system according to Claim 12, wherein said sensors (17, 21) and said processing module (11) are comprised in a same photoplethysmographic heart rate measuring device (11).
The system according to Claim 12, wherein said device (11) is comprised in a apparatus (40) which is wearable on said body organ.
The system according to Claim 14, wherein said wearable apparatus (40) is wearable on the wrist or the arm, in particular by means of a bracelet (41).
The system according to any of Claims 12 to 15, wherein said system comprises a remote processing device (30) connected by a communication link (22, 31) with said photoplethysmographic heart rate measuring device (11).
The system according to Claim 15 or 16, wherein said photoplethysmographic heart rate measuring device (11) includes a display (12).
A photoplethysmographic heart rate measuring device (11) according to any of Claims 12 to 17.
A computer program product that can be loaded into the memory of at least one computer and comprises parts of software code that are able to execute the steps of the method of any of Claims 1 to 11 when the product is run on at least one computer.