SPECMAR: Fast Heart Rate Estimation from PPG Signal using a Modified Spectral Subtraction Scheme with Composite Motion Artifacts Reference Generation
Mohammad Tariqul Islam, Sk. Tanvir Ahmed, Celia Shahnaz, Shaikh Anowarul Fattah
aa r X i v : . [ ee ss . SP ] N ov SPECMAR: Fast Heart Rate Estimation from PPG Signal using aModified Spectral Subtraction Scheme with Composite MotionArtifacts Reference Generation
Mohammad Tariqul Islam ∗ , Sk. Tanvir Ahmed † , Celia Shahnaz ‡ , and Shaikh AnowarulFattah § BUET, Dhaka, Bangladesh
Abstract
The task of heart rate estimation using photoplethysmographic (PPG) signal is challenging due tothe presence of various motion artifacts in the recorded signals. In this paper, a fast algorithm for heartrate estimation based on modified SPEctral subtraction scheme utilizing Composite Motion ArtifactsReference generation (SPECMAR) is proposed using two-channel PPG and three-axis accelerometersignals. First, the preliminary noise reduction is obtained by filtering unwanted frequency componentsfrom the recorded signals. Next, a composite motion artifacts reference generation method is developedto be employed in the proposed SPECMAR algorithm for motion artifacts reduction. The heart rate isthen computed from the noise and motion artifacts reduced PPG signal. Finally, a heart rate trackingalgorithm is proposed considering neighboring estimates. The performance of the SPECMAR algorithmhas been tested on publicly available PPG database. The average heart rate estimation error is found tobe 2.09 BPM on 23 recordings. The Pearson correlation is 0.9907. Due to low computational complexity,the method is faster than the comparing methods. The low estimation error, smooth and fast heart ratetracking makes SPECMAR an ideal choice to be implemented in wearable devices.
Accepted in Medical and Biological Engineering and Computing. This is a pre-copyedit version.http://dx.doi.org/10.1007/s11517-018-1909-x ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] Introduction
Monitoring life signals using wearable devices have become a ubiquitous technology in recent years. Oneof the most common way to monitor life signals in such devices is by using photoplethysmographic (PPG)signals [8]. The photoplethysmography is an electro-optic technique, where a PPG sensor is placed aboveskin and the light reflected from the tissues is captured by a photo detector [16]. This reflected light varieswith respect to the amount of blood present in the tissue which varies with cardiac cycle. This signal whichvaries as cardiac cycle is called the PPG signal. PPG signals have various usage for wearable electronics,which include measurement of heart rate [1], arterial blood oxygen saturation [10], blood pressure [12] andreducing false arrythmia alarm reduction in ICU [4]. Traditionally, electrocardiogram (ECG) is used formonitoring the heart rate. However, ECG is hardly suitable for day to day usage for its cumbersome sensorplacement in the body [6]. In this context, PPG is in a favorable position for its low cost and simpleacquisition mechanism of the PPG sensor. However, depending on the mounting mechanism of the sensorthe captured PPG signal can be heavily corrupted by motion artifacts which is generally considered as animpediment to estimation of heart rate [20]. In this section, the recent methods have been discussed andthen the scope in the domain is discussed along with our contribution.A number of methods have been developed to reduce motion artifacts in the PPG signal and to computeheart rate. Kim et al. [9] employed independent component analysis to remove motion artifacts in PPGsignals. A Kalman filter based method is proposed in [14] for adaptive reduction of motion artifacts from thePPG signals. Fu et al. [3] employed wavelet method to extract heart rate from PPG signal. In a normalizedleast mean square adaptive scheme a difference of two channel PPG signal was used as reference to denoisethe PPG signals [18]. Zhang et al. [20] introduced a large PPG dataset with ground truth to perform heartrate estimation from PPG signal during exercise, which inspired numerous researchers to develop and testdifferent methods. Briefly, the TROIKA [20] and JOSS [19] are the first two methods to report the resultson part of the database. After that a number of researchers worked to solve the problem efficiently. In thiscontext, a detailed review on the state-of-the-art of PPG based heart rate estimation techniques is presentedin [13, 15]. Additionally, Temko [17] proposed a heart rate estimation method where the wiener filter andphase vocoder are employed to remove motion artifacts from the PPG signals. Instead of performing onlysequential feed-forward operations to remove noise, Islam et al. [7] removed the motion artifacts parallelyand independently in time and frequency domain. Later, network of adaptive filters is explored in timedomain using cascade and parallel combination (CPC) of adaptive filters to achieve faster and more accurateresults [6]. One major concern in time domain filtering based approach is that, these methods requirefiltering in each time step thus causing considerable delay. As a result, a method that avoids such timedomain computation yet provides less computational burden is still in great demand.We present a spectral subtraction based motion artifacts removal scheme from PPG signals utilizingaccelerometer data as motion artifacts reference. In this paper, a modified SPEctral subtraction method isproposed along with a Composite Motion Artifacts Reference generation, referred to as SPECMAR. First,the signals are pre-filtered to make them band limited and to remove unwanted frequency components. Next,the three channel accelerometer data are combined to produce a composite motion artifacts spectrum. Amodified spectrum subtraction scheme is then introduced to obtain motion artifacts reduced PPG signal.The heart rate is computed from the noise reduced spectrum. Finally, a tracking algorithm is proposed tofinalize heart rate estimate from the obtained spectrum. The proposed method has been tested on a publiclyavailable database from [20] and has been compared with some state-of-the-art methods.
The data used in this study was introduced by Zhang et al. in the TROIKA framework [20], where results ofonly 12 subjects were presented. Each recording consists of two channel PPG signals, three-axis accelerationsignal and single channel ECG signal sampled at 125 Hz for approximately 5 minutes. Later, data ofadditional 11 recordings have been published where only the 13th recording contains the ECG data. Thetwo channel PPG signals were recorded from wrist using a pulse oximeter with green LEDs of 515 nmwavelength which was used in reflection PPG setup. The distance from center to center of the PPG sensors2as 2 cm. The three-axis acceleration signals were also and one channel chest ECG signal were recordedsimultaneously. Both the pulse oximeters and the accelerometer were embedded in a wristband which wasasked to be worn comfortably. During data recording subjects walked or ran on a treadmill with a speed of1-2 km/hour for 0.5 minute, followed by a speed of 6-8 km/hour for 1 minute, a speed of 12-15 km/hour for1 minute, a speed of 6-8 km/hour for 1 minutes, a speed of 12-15 km/hour for 1 minute, and a speed of 1-2km/hour for 0.5 minute [20]. The subjects were asked to use the hand with the wristband to perform varioustasks, such as, to pull clothes, wipe sweat on forehead, and push buttons on the treadmill, in addition tofreely swing. The 13th recording is of a subject who was also for running on treadmill using the same exerciseroutine as of the first 12 recording. However, the subject had abnormal heart rhythm. In the recordings13, 14, 15, 20 and 23, the subjects performed various forearm and upper arm exercise (e.g. shake hands,stretch, push), running, jumping, or push-ups. In the rest of the recordings, the subjects performed intensiveforearm and upper arm movements (e.g. boxing). The 23rd recording is of a female subject with abnormalheart rhythm and blood pressure.In addition to the signals each recordings contain ground truth heart rates computed from the ECGsignal. The ground truth heart rate has been computed from a 8s time frame. Each time frame has anoverlap of 6s with the previous time frame.
The proposed heart rate estimation method consists of three major stages. They are pre-processing of thesignals, performing spectral subtraction using composite motion artifacts data to remove motion artifactsand estimation of heart rate. Each of the steps are described in the subsections that follow.
The PPG signal has its biological origin that corresponds to the cardiac activity and thus widely used forheart rate measurements. Human heart rate generally varies between 0 . . . . . . N F F T which is also the number of frequency bins in the spectrum. In the magnitudespectrum of a frame of PPG data upto samples corresponding to 3 . M FFT pointscorresponding to H BPM, are termed as the pre-processed PPG spectrum ( X P P G ) and acceleration spectrum( C X , C Y , and C Z ) are obtained. This action does not hamper the heart rate estimation due to the fact thatno useful information related to heart rate lies beyond this specified range. The value of M is computed as M = H × N F F T f s (1)where f s is the sampling frequency. The greatest advantage of spectral subtraction algorithm lies in its simplicity. While using spectral sub-traction, the underlying assumption is that noise is additive and stationary or slowly time varying signal,3 (a) Spectrum of the PPG Signal (b) Waveform of the PPG Signal (c) Spectrum of the Acceleration Signal (Y-axis) (d) Waveform of the Acceleration Signal (Y-axis) (e) Spectrum after Spectral Subtraction
Heart Rate (BPM)
Time (seconds) -1000100 (f) Waveform of the Signal after Spectral Subtraction
Figure 1: An example of Case I where motion artifacts and heart rate peaks are separable and motion artifactspeaks are also clearly visible in PPG signal. The blue circle and red diamond correspond to location of trueheart rate peak and dominant peak, respectively.whose spectrum does not change significantly between the updating periods. However, the main assumptionof noise being stationary does not hold in PPG and sometimes the noise spectrum changes significantlybetween two successive epochs. In case of PPG, improvement of signal intelligibility is the main concernalong with quality. Thus special cares have to be taken to implement spectral subtraction in PPG signalprocessing. In what follows, challenges and observations made while using the traditional spectral subtrac-tion for PPG signals are presented followed by the considerations for generating composite motion artifactsspectrum, and finally the proposed modified spectral subtraction method is described.
Generalized Spectral Subtraction:
By denoting the spectrum of noisy observed PPG signal as Y ( k )and the noise spectrum as N ( k ), a generalized direct spectral subtraction algorithm performs the followingoperation | X ( k ) | p = ( | Y ( k ) | p − | N ( k ) | p , if | Y ( k ) | p > | N ( k ) | p , otherwise (2)where X is the spectrum of noise free signal, p is the power exponent and k = 0 , , . . . , M − p = 1 , (a) Spectrum of the PPG Signal (b) Waveform of the PPG Signal (c) Spectrum of the Acceleration Signal (Y-axis) (d) Waveform of the Acceleration Signal (Y-axis) (e) Spectrum after Spectral Subtraction Heart Rate (BPM)
Time (seconds) -50050 (f) Waveform of the Signal after Spectral Subtraction
Figure 2: An example of Case II, where heart rate peaks and motion artifacts peaks are overlapping, thusnon-separable by naked eye. However, after spectral subtraction a peak near to the true heart rate isdominant. The blue circle and red diamond correspond to location of true heart rate peak and dominantpeak, respectively.spectral subtraction. In order to construct the signal after spectral subtraction the phase information ofthe PPG signal is used. In time frames following Case I, where the motion artifacts peaks’ locations in thePPG spectrum exactly match with the peaks’ location of acceleration spectrum and they locate at distantneighborhood of heart rate peak, calculating difference spectrum by direct spectral subtraction of (2) willgenerate a cleansed PPG spectrum. An example of such case is shown in fig. 1. In this case, an accelerationsignal ( C Y ) is chosen as motion artifacts reference signal. It is obvious that, the motion artifacts dominantpeaks present in the acceleration spectrum, are also dominant in the pre-processed PPG spectrum and thetrue heart rate peak is dominant and at a separable distance. Thus after the basic spectral subtraction, themotion artifacts peaks are significantly reduced and only dominant peak corresponds to true heart rate. Itis seen from the waveforms that the high frequency oscillations corresponding to the motion artifacts aresignificantly reduced in the PPG signal after spectral subtraction.For cases II and III, the spectral subtraction of one or all of the acceleration signals from the pre-processedPPG signal cannot always recover the true heart rate peak. The examples of these cases are shown in Figs.2 and 3. In Case II shown in Fig. 2, though the dominant peak in the PPG represent the true heart rate,the motion artifacts peak in acceleration spectrum are at very close vicinity of the true heart rate peak(the movement rate of hand closely matches the heart rate). Thus performing basic spectral subtraction,the magnitude of the true heart rate peak is reduced and the dominant peak in the spectrum no longercorresponds to true heart rate. Since, heart rate and motion artifacts peaks overlap, the waveforms of thePPG signal before and after spectral subtraction remain mostly unchanged. In Fig. 3, the acceleration signaldoes not represent any motion artifacts. By subtracting it from the PPG spectrum, true heart rate peaklocation cannot be recovered. It is observed from the waveforms of the signal that the morphology of theobtained signal after spectral subtraction is similar to that of the original PPG signal.However from these figures, it is observable that the true heart rate peak can be extracted by alteringthe basic spectral subtraction algorithm in (2) in such a way that motion artifacts removal is facilitated. Inorder to improve the spectral subtraction for these above mentioned cases, beside modification to the actualalgorithm, a careful choice of motion artifacts spectrum is also required. Motion artifacts spectrum can beconstructed as one of the acceleration signals or a composite spectrum constructed from the three channelacceleration spectra. Such propositions are discussed in the next subsection. Motion Artifact Spectrum Generation:
Though motion at a certain time can be differently capturedby any of the acceleration signals, a careful observation shows that, in most time frames, the motion artifactspeaks are available simultaneously in all acceleration signals. This is largely due to the periodic movement5 (a) Spectrum of the PPG Signal (b) Waveform of the PPG Signal (c) Spectrum of the Acceleration Signal (Y-axis) (d) Waveform of the Acceleration Signal (Y-axis) (e) Spectrum after Spectral Subtraction
Heart Rate (BPM)
Time (seconds) -2000200 (f) Waveform of the Signal after Spectral Subtraction
Figure 3: An example of Case III, where the acceleration signal does not bear much information aboutmotion artifacts in the PPG signal. The blue circle and red diamond correspond to location of true heartrate peak and dominant peak, respectively.of hands during the physical activity and such movement may not necessarily be aligned in any particulardirection of the accelerometer. For example, during running the hand movement appears as back and forthif the subject is seen from a side, appears as up and down if the subject is observed from behind and circularif seen from top. As a result, all of the acceleration signals capture similar frequencies of motion. Hence,by using only one of the acceleration signals in spectral subtraction may result in partial removal of motionartifacts from the PPG signal. On the other hand, using all three acceleration signals successively may resultin diminishing the PPG signal. Hence, a carefully designed composite of three axis acceleration signal mayprovide satisfactory or desired performance.It has already been discussed that, the channels of accelerometer experience movements with same periodduring a physical exercise. Since, all three acceleration signals are expected to have same frequency com-ponents of motion artifacts of different amplitude, the question thus arise is, for each frequency componentwhich one to choose among the three components. Naively, the largest component for a frequency bin ofthe three spectra may be considered. However, sometimes the large component in one of the channels isnot due to motion artifacts but due to noise in accelerometer. Since, motion artifacts is present in all threechannels, choosing the minimum magnitude of the three spectra for a frequency bin captures the motionartifacts most efficiently. Thus, for the k th frequency bin of composite motion artifact reference (CMAR)spectrum N CMAR is formed by, N CMAR ( k ) = min( C X ( k ) , C Y ( k ) , C Z ( k )) . (3)The composite spectrum is normalized in range [0 , C X captures motion artifact with approximately equal dominance in twofrequencies, one fundamental and another at its double frequency. C Y and C Z capture these frequenciesas well, however, only one of them is dominant. It is also seen from the figure that, the dominant motionartifact in the PPG signal includes only the higher frequency of the accelerometer signal. The compositemotion artifact reference spectrum constructed by following (3) is dominated by the higher frequency fromthe motion artifacts whereas the lower frequency motion artifact component is not dominant. Thus, in thespectrum after spectral subtraction, none of the motion artifacts frequencies are dominant and the location ofthe true heart rate and dominant peak coincide. In a similar way, Fig. 5 shows a case where the accelerationsignals are noisy and does not correspond to motion artifact. In addition to previously considered referencesan alternative to (3) is considered where instead of min( · ) function, the max( · ) function is considered for6 (a) Spectrum of the PPG Signal (b) Spectrum of the PPG Signal (c) Spectrum of X-ch. Acceleration (d) Spectral Subtraction Employing X-ch. Acceleration (e) Spectrum of Y-ch. Acceleration (f) Spectral Subtraction Employing Y-ch. Acceleration (g) Spectrum of Z-ch. Acceleration (h) Spectral Subtraction Employing Z-ch. Acceleration (i) Spectrum of CMAR Heart Rate (BPM) (j) Spectral Subtraction Employing CMAR
Heart Rate (BPM)
Figure 4: Comparison of the spectra generated after spectral subtraction using different motion artifactsreferences. (a),(b) Spectrum of the PPG signal. Spectrum of the motion artifacts reference using (c) X-axisacceleration signal C X , (e) Y-axis acceleration signal C Y , (g) Z-axis acceleration signal C Z , and (i) usingthe proposed motion artifacts reference generation method N CMAR . (d), (f), (h), (j) The spectra generatedafter employing spectrum subtraction on the references (c), (e), (g), and (i), respectively. The blue circleand red diamond correspond to location of true heart rate peak and dominant peak, respectively.motion artifacts reference. In this case, the large components in the spectra of the acceleration signals aredue to noise of the accelerometer and do not correspond to the motion artifacts. It is seen from the figurethat after spectral subtraction, the max( · ) function causes degradation in the peak near the heart rate. Onthe other hand, using the proposed CMAR in the spectral subtraction, the dominant peak and the locationof the heart rate coincide. Thus, the proposed CMAR can be an efficient motion artifacts reference forspectral subtraction denoising mechanism for PPG signals. Modified Spectral Subtraction (MSS) Algorithm:
In order for the spectral subtraction to workbetter with PPG, a modification of the basic spectral subtraction is required. Direct subtraction using (2)sometimes produces unwanted results, as seen in Figs. 2 and 3. In formulation of (2), | Y ( k ) | is assumedto be the sum of | X ( k ) | and | N ( k ) | . However, true | N ( k ) | is not available and a scaled version of | N ( k ) | isestimated as the motion artifacts spectrum. Thus we can replace | N ( k ) | with [ α N | N ( k ) | ], where α N is thescaling factor of the noise. Moreover, the | Y ( k ) | which is found from the PPG sensor is also an estimatedor the scaled version of the internal sum of | X ( k ) | and α N | N ( k ) | . Thus the expression of Y ( k ) is given by | Y ( k ) | = α Y ( | X ( k ) | + α N | N ( k ) | ) , (4)where α Y is the scaling factor of the observed output. Hence, we can write the estimation of | X ( k ) | as | X ( k ) | = 1 α Y | Y ( k ) | − α N α Y | N ( k ) | . (5)7 (a) Spectrum of the PPG Signal (b) Spectrum of the PPG Signal (c) Spectrum of X-ch. Acceleration (d) Spectral Subtraction Employing X-ch. Acceleration (e) Spectrum of Y-ch. Acceleration (f) Spectral Subtraction Employing Y-ch. Acceleration (g) Spectrum of Z-ch. Acceleration (h) Spectral Subtraction Employing Z-ch. Acceleration (i) Spectrum using Max Function (j) Spectral Subtraction Employing Max Function (k) Spectrum of CMAR Heart Rate (BPM) (l) Spectral Subtraction Employing CMAR
Heart Rate (BPM)
Figure 5: Comparison of the spectra generated after spectral subtraction using different motion artifactsreferences. (a),(b) Spectrum of the PPG signal. Spectrum of the motion artifacts reference using (c) X-axis acceleration signal C X , (e) Y-axis acceleration signal C Y , (g) Z-axis acceleration signal C Z , (i) usingthe max( C X [ k ] , C Y [ k ] , C Z [ k ]) fucntion and (k) using the proposed motion artifacts reference generationmethod N CMAR . (d), (f), (h), (j) and (l) The spectra generated after employing spectrum subtraction onthe references (c), (e), (g), (i) and (k), respectively. The blue circle and red diamond correspond to locationof true heart rate peak and dominant peak, respectively.From the discussion above, the modified spectral subtraction (MSS) is proposed as | X ( k ) | p = ( α | Y ( k ) | p − α | N ( k ) | p , if α | Y ( k ) | p > α | N ( k ) | p , (6)where, α and α essentially works as the weight parameters. As both PPG signal spectrum ( X = X P P G )and noise signal spectrum ( N = N CMAR ) values are in the range of [0 , α and α , canalso be chosen in the range of [0 , α = α = 1 results in α Y = 1 and α N = 1, which is the case considered in (2) with N = N CMAR , is termed as SPECMAR without scaling (SPECMARWS).8 (a) Spectrum of the PPG Signal (b) Spectrum of CMAR (c) Spectral Subtraction Employing CMAR (d) Modified Spectral Subtraction Employing CMAR
Heart Rate (BPM)
Figure 6: Comparison of the spectra generated from spectral subtraction and proposed modified spectralsubtraction using CMAR. (a) Spectrum of the PPG signal. (b) Spectrum of motion artifacts referenceemploying composite motion artifacts reference generation method. Spectrum obtained after (c) emloyingspectral subtration without scaling and (d) employing modified spectral subtraction with α = 0 .
88 and α = 0 .
70. The blue circle and red diamond correspond to location of true heart rate peak and dominantpeak, respectively.
Once the result of an individual frame is available, it is more logical to consider the results obtained innearby frames and by incorporating an efficient tracking algorithm for estimating heart rate. The trackingalgorithm proposed here consists of two main parts: initial estimation of heart rate at the beginning of themethod and estimation and tracking of current heart rate.
Initial Estimate:
The method requires an initial estimate for the first time frame. The initial estimationis performed using the output of the modified spectral subtraction ( X MSS ). The highest peak of the spectrumis selected as the initial estimate. This estimated heart rate is used as the heart rate of the previous timeframe ( N ) for the next frame. Estimation and Tracking of Heart Rate Using Best Candidate Approach:
Even though thenoise suppression using modified spectral subtraction is highly effective, there are some frames where itcannot resolve the heart rate peak. Most of them are due to motion artifacts components being in closevicinity of the heart rate component. In such cases, an error in heart rate estimation may occur. A properestimation technique has to be employed to mitigate such error. One key point in estimating the heart ratesuccessfully is that, the change of heart rate is generally not abrupt, rather gradual. If the time overlappingbetween two successive time frames is small, the true heart rate for latest frame will be close to that of theprevious frame. Considering this, both pre-processed PPG spectrum ( X P P G ) and processed PPG spectrum( X MSS ) are employed to track BPM correctly. The candidate peaks are selected from the X P P G signalconsidering the search region, P = [ N − ∆ s , ..., N + ∆ s ]. That is, the neighborhood, ± ∆ s bins, of thefrequency bin of the previous time frame ( N ) are selected to be the most likely location of the heart ratecomponent. Then peaks and their locations inside this search region are extracted. The location of thepeaks that have magnitude greater than 25% of the magnitude of maximum peak of X P P G are put in theset S = [ N , N , ... ]. Now selection of location of heart rate component can be considered for three cases,Case I: If only one peak is found in S , that is cardinality of the set is 1, this is the dominant peak, D ofthe frame. Then in the close vicinity, ± δ frequency bins, of the location of that peak, a search is performed9n X MSS to know whether that peak is present or a different peak has surfaced due to the motion artifactremoval. The peak with highest magnitude in this search region is selected as the location of the heart ratepeak ( N cur ).Case II: If more than one peak is found, it is suspected that motion artifacts or some other kind of noiseis present in this spectrum and so emphasis is given on the dominant peak of the X P P G . If the dominantpeak is within the range, ± ∆ t frequency bins, of N , a search is performed in X MSS within the ± δ rangeof N to find a satisfactory peak, that is the magnitude of the peak is at least 10% of the dominant peak.The extracted location is set to N cur .Case III: If a satisfactory peak is not found in the two cases above, then all the peaks and their locationin the search region P is extracted from X MSS . The set containing the normalized magnitude of the peaksis A and the set containing the location of the peaks is S . Then a closeness to the heart rate location ofprevious frame is stored in set C as C = (cid:18) − | S − N | κN (cid:19) , (7)where κ is a normalizing parameter chosen such that the values inside C is in the range [0 , A and the closeness set C provides theestimated location of heart rate i = argmax( β A + β C ) (8)and the i th member of S is set to N cur . Here, β and β are positive weight parameters satisfying β + β = 1.After that, heart rate in BPM is calculated to beˆ B est = N cur − N F F T × × f s . (9)Next in order to obtain a smooth heart rate tracking, a three point moving average is employed wherethe estimate of BPM in current frame is computed as B ′ est = γ ˆ B est + γ B − + γ B − , (10)where B − and B − represent estimated BPM of previous two successive time frames and sum of the weightparameters satisfies γ + γ + γ = 1 and these values are empirically chosen.Finally in order to prevent extremely high or low estimated values of BPM in comparison to previousBPM estimate moving averaged value B ′ est is modified as B est = B ′ est + 5 , if B ′ est − B − ≥ λ inc B ′ est − , if B ′ est − B − ≤ λ dec B ′ est , otherwise , (11)where constants λ inc and λ dec are set empirically to small integer values. The proposed method requires a number of values to be initialized. Most of these parameters are chosenempirically with reference to previous literature. The number of points to generate spectrum, N F F T ischosen to be 4096. In order to produce truncated spectrum, H is set to 240. The weight parameters α and α of modified spectral subtraction are set to 0 .
88 and 0 .
70, respectively and the exponent p is set to 1. Forlarge searching neighborhood, the parameters ∆ s and ∆ t are both set to 30, whereas small search regions δ and δ are both set to 3. The weight paramers β and β of (8) are empirically set to 0 . .
3. Theweights γ , γ and γ , of the moving average operation are set to 0 .
9, 0 .
05 and 0 .
05, respectively and theparameters λ inc and λ dec are set to 5 and −
3, respectively as in [7] considering the rate of change of heartrate. 10able 1: Performance comparison in terms of AAE for 23 Data Recordings.Data TROIKA [20] MISPT-125 [11] WFPV [17] TFD [7] CPC [6] SPECMARWS SPECMAR1 2 .
29 1 .
58 1 . . .
56 1 .
22 1 .
222 2 .
19 1 .
80 1 . . .
25 1 .
82 1 .
513 2 . . .
72 0 .
79 0 .
66 0 .
78 0 .
754 2 .
15 0 .
99 0 .
98 0 . . .
26 1 .
265 2 .
01 0 .
74 0 .
75 0 . . .
83 0 .
756 2 .
76 0 . . .
14 0 .
96 1 .
88 1 .
877 1 .
67 0 . . .
71 0 .
80 1 .
06 0 .
808 1 . . .
91 0 .
73 0 .
56 1 .
09 1 .
079 1 . . .
54 0 .
64 0 .
56 0 .
71 0 . .
70 3 .
60 2 .
61 3 .
09 2 .
62 3 . .
11 1 . . .
94 1 .
34 1 .
27 1 .
46 1 . . . .
98 1 .
54 0 .
79 1 .
46 1 . Mean 12 2 .
34 1 .
11 1 .
05 1 .
16 1 .
12 1 .
33 1 . .
47 2 .
33 1 .
33 1 .
74 2 .
30 1 .
77 1 .
13 - - 3 .
58 4 . . .
98 3 . .
66 6 .
69 14 .
58 6 . .
15 - - 2 .
31 2 .
29 1 .
86 2 . .
16 - - 4 .
93 16 .
24 10 .
19 3 . .
17 - - 3 .
07 3 .
02 7 .
62 2 . .
18 - - 2 .
67 2 . . .
18 2 . . .
27 8 .
39 5 .
53 4 . . .
87 2 .
53 2 .
56 2 . . .
99 10 .
45 4 .
06 3 . .
35 1 . . .
83 1 . .
75 0 . . .
82 0 . Mean 11 - - .
61 4 .
22 5 .
63 3 .
35 3 . - - −− .
91 12 .
98 5 .
10 4 . - - .
28 2 .
73 3 .
26 2 .
30 2 . - - .
33 4 .
01 4 .
25 3 .
81 3 . .4 Performance Measurement The evaluation of performance is carried out in terms of few different performance indices, i.e. averageabsolute error (AAE), pearson correlation ( r ) and Bland-Altman plot [2].For each data recording the AAE is defined as AAE = 1 N N X i =1 | B est ( i ) − B true ( i ) | , (12)where B est ( i ) and B true ( i ), respectively, are estimated and ground truth value of heart rate in BPM of the i th frame and N is the total number of time frames. Pearson correlation is a measure of degree of similaritybetween true and estimated values of heart rate. Higher the value of r better the estimates. The Bland-Altman plot measures the agreement between true and estimates of heart rate. Here limit of agreement(LOA) is computed using the average difference ( µ ) and the standard deviation ( σ ), which is defined as[ µ − . σ, µ + 1 . σ ]. LOA is a measurement which shows that 95% data exist within 1 . σ from mean.Development of the method and the evaluation are performed in MATLAB 2012b with a 2.2 GHz IntelProcessor and 4 GB of RAM. Few recently reported heart rate estimation methods have been taken into consideration for the purposeof performance comparison with the proposed method. The methods are TROIKA [20], MISPT [11],WFPV [17], TFD [6], and CPC [7]. Moreover, in order to show the effectiveness of the proposed modifiedspectral subtraction (referred to as SPECMAR), the algorithm’s performance with direct spectral subtrac-tion (referred to as SPECMARWS) is also compared. In both cases, proposed SPECMARWS (using (2))and proposed SPECMAR (using (6)) estimated noise spectrum is obtained using proposed motion artifactsspectrum generation in (3). Table 1 shows the average absolute error on 23 recordings when implementingdifferent methods. The summary of the first 12 recordings are denoted as
Mean 12 and
Std 12 for mean ofAAE and standard deviation of absolute errors, respectively. For the rest of the 11 recordings the summaryof metrics are denoted as
Mean 11 and
Std 11 . And finally, for all 23 recordings of the database, the meanof AAE and standard deviation of absolute errors are denoted as
Mean All and
Std All , respectively. Itis to be noted that in case of TROIKA method, for performance evaluation, some initial time frames areexcluded. However, in the proposed method all the time frames are included and none are excluded. It canbe observed that, the proposed method provides a better AAE in 5 of the 23 recordings and overall meanAAE is lowest among the compared method.The Pearson correlation between ground truth and estimated heart rate is computed and shown in Fig. 7.The value of the Pearson coefficient is found 0 . . − . , . − . , . α is optimal, the value of α has been varied by ± .
88 while setting the value of α to 0 . α has been varied from its nominal value of 0 . α has been set to 0 .
88. It is seen from Fig. 10 that as the parameters are varied, the AAE firstdecreases and then increases again. The values α = 0 .
88 and α = 0 .
70 provides the lowest AEE whichhas been marked using a black circle. Additionally, a similar experiment has been conducted to investigatethe effect of the number of frequency bins in the spectra N F F T . In this experiment, the N F F T has been12
Ground Truth of Heart Rate (BPM) E s ti m a t e d H ea r t R a t e ( B P M ) Pearson Correlation Between Ground Truth and Estimated Heart RatePearson Correlation Coefficient is 0.9907
Figure 7: The Pearson correlation of the estimation results on the 23 data recordings is 0 . . . . . . . . .
96, 2 .
65, 2 .
09, 2 . .
41, respectively. The fact that, the value 4096 provides the lowest AAE is consistent with the resultsobtained in [20]. Note that, the number of frequency bins in the search region has been scaled according tothe value of N F F T in the experiments.Finally, the computation time required by the proposed method is compared with those required bysome state of the art methods. In Table 2, mean runtime for TROIKA, MISPT, WFPV, TFD, CPC and theproposed SPECMAR methods are provided. It can be observed from the table that, the proposed methodhas lowest runtime compared to any of the methods in the tables.
In this paper, an efficient method, for motion artifacts removal from PPG signal based on modified spectralsubtraction and composite motion artifacts reference generation is proposed. It has been shown in Figs. 1,2 and 3 that only using a single accelerometer channel may not be beneficial depending on the informationthe channel provides. Thus, construction of a composite motion artifacts reference has been explored.It has been observed that the motion artifacts frequencies are generally present in all three axes of theaccelerometer albeit having different amplitudes. As a result, motion artifacts can be captured efficientlyfrom the accelerometer signals. It is found that the minimum magnitude of the three accelerometer spectrafor a frequency bin efficiently captures the motion artifacts. It is shown in Figs. 4 and 5 that the proposedcomposite motion artifacts reference generation method can efficiently capture the motion artifacts. It hasalso been stated that the captured PPG signal contains a scaled version of both the noise-free PPG signaland motion artifacts. With this view, the spectral subtraction algorithm has been modified to be a scaledspectral subtraction. It is shown in Fig. 6 that the proposed modified spectral subtraction can improve the13
Average of Ground Truth and Estimate (BPM) -25-20-15-10-50510152025 D i ff e r e n ce o f t h e T w o M ea s u r e m e n t s Bland Altman Plot -1.96+1.96
Figure 8: The Bland-Altman plot of the estimation on the 23 recordings. The LOA is [ − . , . − . , .
86] BPM.spectrum quality and facilitate the extraction of true heart rate.In the proposed SPECMAR method, the motion artifacts generation and modified spectral subtractionhave been combined. It can be observed from Table 1 that SPECMAR has lower AAE than SPECMARWSwhich uses spectral subtraction without scaling. The high value of Pearson correlation of 0 . − . , .
18] shows that the method can track the heart rate satisfactorily. Moreover,in order to reduce computational complexity, only the relevant frequency bins (corresponding to 0 to 3.5Hz)have been used in the computation. The overall computation pipeline includes only few filtering operations,spectral subtraction, and a decision mechanism for estimating and tracking the heart rate. As a result, theproposed SPECMAR is faster than many recently published methods as shown in Table 2.With respect to previous methods the proposed SPECMAR provides several significant improvements.Similar to the proposed method, the methods, namely, TROIKA, WFPV and TFD employ frequency domainnoise removal. TROIKA and TFD use the singular spectrum analysis [5] to decompose the PPG signal.Then each of the decomposed signals are processed in spectral domain where the spectral peaks of thedecomposed signals are matched with the spectral peaks of the acceleration signal in order to discard someof the decomposed signals. Such a process requires large number of computations for each of the frames.Moreover, the TROIKA method employs computationally expensive MFOCUSS algorithm to construct thespectra of the signals which is the primary cause of the long runtime of the algorithm. The WFPV methodperforms the simplest Wiener adaptive filtering in the frequency domain to remove motion artifacts from thePPG signals and the phase vocoder part requires extensive time-frequency analysis. The methods, namely,MISPT, TFD and CPC employ time domain noise removal using adaptive filters. The time domain adaptivefilter generally suffer from slower convergence rate compared to its frequency domain counterparts and theoutput has to be computed for each time steps which causes the methods to require extensive computation.In addition, MISPT algorithm requires the spectra of all the previous time steps for computing trajectoryand thus causes large memory requirement. On the other hand, TFD method requires cascaded adaptivefiltering for noise removal where three adaptive filters are processed serially and the CPC method requirestwo of such cascaded adaptive filters. As a result, the number of computation is high in these algorithms. Onthe other hand, the proposed SPECMAR is simplest using only modified spectrum subtraction as the noiseremoval part and overall employs fewest computations. The computations are further reduced by discardingthe frequency bins which are outside the region of interest (greater than 3 .
50 100 150
Time Windows B ea t s P e r M i nu t e ( B P M ) (a) Estimated Heart Rate for Recording 9 Ground TruthProposed SPECMAR
Time Windows B ea t s P e r M i nu t e ( B P M ) (b) Estimated Heart Rate for Recording 15 Ground TruthProposed SPECMAR
Figure 9: Tracking performance of the proposed SPECMAR for recording (a) 9 (AAE: 0 .
64) and (b) 15(AAE: 1 . -10 -8 -6 -4 -2 0 2 4 6 8 10 Variation of the Parameter (%) AA E Effect of and varied, =0.70 varied, =0.88 Figure 10: Effet of values of α and α on AEE. One of the values has been varied keeping the other onefixed. The lowest AAE has been marked using a black circle. In this paper, an efficient method, for motion artifacts removal from PPG signal based on modified spectralsubtraction and composite motion artifacts reference generation is proposed. The motion artifact is capturedby taking minimum amplitude of each of the frequency bins of the three channel accelerometer. By usingthis composite motion artifacts in the modified spectral subtraction algorithm, significant noise reductionis achieved which helps in obtaining accurate heart rate estimation. Finally, an efficient tracking scheme isprovided to track the heart rate computed from the noise removed spectrum. From extensive experimentson recorded PPG data, it is shown that the proposed SPECMAR method is capable of tracking the groundtruth with high estimation accuracy and comparable to recently published methods. Also the method isquantitatively faster than these methods. As a result, it can be used in real-time application withoutsacrificing estimation accuracy. The proposed spectral subtraction method can be employed in other areassuch as noise removal in speech processing. An improvement of the modified spectral subtraction can beobtained by studying different techniques for setting the parameters α and α automatically. Additionally,in future alternative tracking algorithm such as the adaptive notch filter can be explored for tracking the15eart rates in the algorithm. Acknowledgement
Authors would like to thank Andriy Temko, Research Fellow, Irish Centre for Fetal and Neonatal Transla-tional Research, Department of EEE, University College Cork for his correspondence regarding computationtime of his method. The authors would also like to thank the anonymous reviewers for their valuablecomments that were useful to improve the quality of the paper.
References [1] Allen, J.: Photoplethysmography and its application in clinical physiological measurement. Physiologicalmeasurement (3), R1 (2007)[2] Bland, J.M., Altman, D.G.: Comparing methods of measurement: why plotting difference againststandard method is misleading. The lancet (8982), 1085–1087 (1995)[3] Fu, T.H., Liu, S.H., Tang, K.T.: Heart rate extraction from photoplethysmogram waveform usingwavelet multi-resolution analysis. J. Medical and Biological Engineering (4), 229–232 (2008)[4] Gambarotta, N., Aletti, F., Baselli, G., Ferrario, M.: A review of methods for the signal quality assess-ment to improve reliability of heart rate and blood pressures derived parameters. Medical & biologicalengineering & computing (7), 1025–1035 (2016)[5] Golyandina, N., Nekrutkin, V., Zhigljavsky, A.A.: Analysis of time series structure: SSA and relatedtechniques. Chapman and Hall/CRC (2001)[6] Islam, M.T., Ahmed, S.T., Zabir, I., Shahnaz, C., Fattah, S.A.: Cascade and parallel combination(CPC) of adaptive filters for estimating heart rate during intensive physical exercise from photoplethys-mographic signal. Healthcare Technology Letters (1), 18–24 (2018)[7] Islam, M.T., Zabir, I., Ahamed, S.T., Yasar, M.T., Shahnaz, C., Fattah, S.A.: A time-frequency do-main approach of heart rate estimation from photoplethysmographic (PPG) signal. Biomedical SignalProcessing and Control , 146–154 (2017)[8] Jain, P.K., Tiwari, A.K.: Heart monitoring systems a review. Computers in Biology and Medicine ,1–13 (2014)[9] Kim, B.S., Yoo, S.K.: Motion artifact reduction in photoplethysmography using independent componentanalysis. IEEE Trans. Biomedical Engineering (3), 566–568 (2006)[10] Kyriacou, P., Powell, S., Langford, R., Jones, D.: Investigation of oesophageal photoplethysmographicsignals and blood oxygen saturation measurements in cardiothoracic surgery patients. Physiologicalmeasurement (3), 533 (2002)[11] Lakshminarasimha Murthy, N.K., Madhusudana, P.C., Suresha, P., Periyasamy, V., Ghosh, P.K.: Multi-ple spectral peak tracking for heart rate monitoring from photoplethysmography signal during intensivephysical exercise. IEEE Signal Processing Letters (12), 2391–2395 (2015)[12] McCombie, D.B., Reisner, A.T., Asada, H.H.: Adaptive blood pressure estimation from wearable PPGsensors using peripheral artery pulse wave velocity measurements and multi-channel blind identificationof local arterial dynamics. In: in Proc. 28th Annual Int. Conf. IEEE Engineering in Medicine andBiology Society, pp. 3521–3524. IEEE, New York City, NY (2006)[13] Periyasamy, V., Pramanik, M., Ghosh, P.K.: Review on heart-rate estimation from photoplethysmog-raphy and accelerometer signals during physical exercise. J. Indian Institute of Science pp. 1–12 (2017)1614] Seyedtabaii, S., Seyedtabaii, L.: Kalman filter based adaptive reduction of motion artifact from photo-plethysmographic signal. World Academy of Science, Engineering and Technology , 173–176 (2008)[15] Tamura, T., Maeda, Y.: Photoplethysmogram. In: Seamless Healthcare Monitoring, pp. 159–192.Springer (2018)[16] Tamura, T., Maeda, Y., Sekine, M., Yoshida, M.: Wearable photoplethysmographic sensorspast andpresent. Electronics (2), 282–302 (2014)[17] Temko, A.: Estimation of heart rate from photoplethysmography during physical exercise using wienerfiltering and the phase vocoder. In: 37th Annual International Conference of the IEEE Engineering inMedicine and Biology Society, pp. 1500–1503. IEEE, Milano, Italy (2015)[18] Yousefi, R., Nourani, M., Ostadabbas, S., Panahi, I.: A motion-tolerant adaptive algorithm for wearablephotoplethysmographic biosensors. IEEE J. Biomedical and Health Informatics (2), 670–681 (2014)[19] Zhang, Z.: Photoplethysmography-based heart rate monitoring in physical activities via joint sparsespectrum reconstruction. IEEE Trans. Biomed. Eng. (8), 1902–1910 (2015)[20] Zhang, Z., Pi, Z., Liu, B.: TROIKA: A general framework for heart rate monitoring using wrist-type photoplethysmographic signals during intensive physical exercise. Biomedical Engineering, IEEETransactions on62