Removal of Motion Artifacts in Photoplethysmograph Sensors during Intensive Exercise for Accurate Heart Rate Calculation Based on Frequency Estimation and Notch Filtering

With photoplethysmograph (PPG) sensors showing increasing potential in wearable health monitoring, the challenging problem of motion artifact (MA) removal during intensive exercise has become a popular research topic. In this study, a novel method that combines heart rate frequency (HRF) estimation and notch filtering is proposed. The proposed method applies a cascaded adaptive noise cancellation (ANC) based on the least mean squares (LMS)-Newton algorithm for preliminary motion artifacts reduction, and further adopts special heart rate frequency tracking and correction schemes for accurate HRF estimation. Finally, notch filters are employed to restore the PPG signal with estimated HRF based on its quasi-periodicity. On an open source data set that features intensive running exercise, the proposed method achieves a competitive mean average absolute error (AAE) result of 0.92 bpm for HR estimation. The practical experiments are carried out with the PPG evaluation platform developed by ourselves. Under three different intensive motion patterns, a 0.89 bpm average AAE result is achieved with the average correlation coefficient between recovered PPG signal and reference PPG signal reaching 0.86. The experimental results demonstrate the effectiveness of the proposed method for accurate HR estimation and robust MA removal in PPG during intensive exercise.


Introduction
Photoplethysmography (PPG) has proven effective in monitoring cardiovascular-related physiological signs, especially heart rate (HR), oxygen saturation (SpO 2 ) and blood pressure (BP) [1]. Due to the advantages of low cost and convenience, PPG sensors are widely applied in wearable healthcare. Though having a great potential for wearable healthcare, the accuracy of PPG sensors during motion such as exercise by the user is still unsatisfactory due to motion artifacts (MA) [2]. The MA is typically caused by the change of blood flow velocity induced by the motion [3] or the relative movement between PPG sensors and human skin [4]. The wide frequency range of MA with time-varying nature makes it difficult to use traditional filtering techniques for the removal of motion artifacts [5]. Thus, to improve the reliability and accuracy of health monitors based on PPG sensors, the removal of MA continues to be a technical challenge that needs to be tackled.

Proposed Method
Shown in Figure 1 is the overall signal flow of the proposed method. The method requires PPG signal(s) and three-axis acceleration signal(s) as its inputs. The 4th order 0.4 Hz-4.0 Hz Butterworth band-pass filters (BPF) are firstly employed to remove the out-of-band noise of the input PPG and acceleration signals. Then a cascaded ANC stage is adopted to adaptively cancel the motion artifacts in the PPG signals using respectively x, y and z-axis acceleration signals as reference noise. After that, a simple HRFT scheme is applied to give a preliminary estimation for the HRF from the spectrum of the ANC-de-noised PPG signal. Then a novel HRFC stage is used to correct possibly wrong HRF estimation and give the final estimated HRF value. The HRF and its second harmonics frequency are used as notch frequencies to construct two notch filters. The notch filters remove the PPG components from the ANC-de-noised PPG signal and extract the motion artifacts. And the MA-removed PPG signal is finally obtained through subtracting the motion artifacts from the ANC-de-noised PPG signal.
de-noised PPG signal and extract the motion artifacts. And the MA-removed PPG signal is finally obtained through subtracting the motion artifacts from the ANC-de-noised PPG signal. A widely-utilized open source data set provided for the 2015 IEEE Signal Processing Cup [12] is used for preliminary evaluation of different stages in the proposed method. The data set contains signals collected from 12 subjects (11 males and one female). For each subject, two channels of PPG signals and three-axis acceleration signals are recorded from one wrist, with one channel electrocardiograph (ECG) signal recorded simultaneously for about 5 min with a sampling rate of 125 Hz. As the ECG signal is barely influenced by MA, it is used for calculating the reference HR. During the 5 min recording, the subjects exercise on a treadmill in the following pattern: 1-2 km/h walking for 30 s, 6-8 km/h running for 60 s, 12-15 km/h running for 60 s, 6-8 km/h running for 60 s, 12-15 km/h running for 60 s and 1-2 km/h walking for 30 s. For the first 30 s, subjects are in a relatively static initialization stage. Average Absolute Error (AAE) between reference HR and estimated HR from the proposed method [12] is used for evaluating performance of the proposed method, which denoted as HRAAE. The HRAAE for each subject is calculated using Equation (1). HREST(i) and HRTRUE(i) are respectively the i-th HR estimation results and reference HR value. N is the number of total estimations for a subject.

Cascaded Adaptive Noise Cancellation
Adaptive noise cancellation is an effective technique to de-noise a noise-corrupted signal given a reference signal that is highly correlated with the noise [17]. Figure 2 shows the signal flow of typical ANC. In the proposed method, the least mean squares-Newton (LMS-Newton) algorithm [18] is applied for adaptive filtering, which features a faster convergence speed than conventional LMSbased algorithms. The LMS-Newton algorithm is given in Equation (2) and Equation (3). The parameter x(k) is the reference signal of MA, which is the acceleration signal. W(k) is the weight vector of the finite impulse response (FIR) filter, and e(k) is the output for de-noised PPG signal. The parameter α and μ in Equation (3) determine the speed of convergence. A widely-utilized open source data set provided for the 2015 IEEE Signal Processing Cup [12] is used for preliminary evaluation of different stages in the proposed method. The data set contains signals collected from 12 subjects (11 males and one female). For each subject, two channels of PPG signals and three-axis acceleration signals are recorded from one wrist, with one channel electrocardiograph (ECG) signal recorded simultaneously for about 5 min with a sampling rate of 125 Hz. As the ECG signal is barely influenced by MA, it is used for calculating the reference HR. During the 5 min recording, the subjects exercise on a treadmill in the following pattern: 1-2 km/h walking for 30 s, 6-8 km/h running for 60 s, 12-15 km/h running for 60 s, 6-8 km/h running for 60 s, 12-15 km/h running for 60 s and 1-2 km/h walking for 30 s. For the first 30 s, subjects are in a relatively static initialization stage. Average Absolute Error (AAE) between reference HR and estimated HR from the proposed method [12] is used for evaluating performance of the proposed method, which denoted as HR AAE . The HR AAE for each subject is calculated using Equation (1). HR EST (i) and HR TRUE (i) are respectively the i-th HR estimation results and reference HR value. N is the number of total estimations for a subject.

Cascaded Adaptive Noise Cancellation
Adaptive noise cancellation is an effective technique to de-noise a noise-corrupted signal given a reference signal that is highly correlated with the noise [17]. Figure 2 shows the signal flow of typical ANC. In the proposed method, the least mean squares-Newton (LMS-Newton) algorithm [18] is applied for adaptive filtering, which features a faster convergence speed than conventional LMS-based algorithms. The LMS-Newton algorithm is given in Equations (2) and (3). The parameter x(k) is the reference signal of MA, which is the acceleration signal. W(k) is the weight vector of the finite impulse response (FIR) filter, and e(k) is the output for de-noised PPG signal. The parameter α and µ in Equation (3) determine the speed of convergence.
Because MA can be decomposed into motions in different directions, a cascaded topology of ANC is used here to reduce the motion artifacts better, as is illustrated in Figure 1. The cascaded adaptive noise cancellers use respectively x, y and z-axis acceleration signals as their MA reference inputs. Band-pass filtered PPG and acceleration signals are streamed into this step at a specified time interval ltv. In order to avoid an unmatched scale between the PPG signal and the acceleration signal, Z-score standardization is adopted to normalize the inputs of each adaptive noise canceller. Especially, if there are multiple PPG signals at the input end of the cascaded ANC, the PPG signals will be averaged into a single PPG signal before the first adaptive noise canceller.
Shown in Table 1 is the experimental results of HRAAE before and after applying the cascaded ANC stage of the proposed method. When only BPF is applied to the MA-corrupted signals, the estimation errors for HR values are very large, reaching an average HRAAE of 11.47 bpm over the 12 subjects, which indicates that the MA frequencies are strongly deviated from the HRFs and the spectral magnitudes of MA are larger. After applying the cascade ANC stage for MA reduction, the HR estimation errors significantly decreased.

Heart Rate Frequency Tracking
After the cascaded ANC stage, the de-noised PPG signal SPPG is used for a preliminary estimation of current HRF. FFT is first applied to acquire the spectrum of SPPG, where only the magnitude information is considered. A relatively large number of points for FFT LFFT is used to increase frequency resolution. Then a simple HRFT scheme from [15] is used for preliminary HRF estimation.
Flowchart of the HRFT scheme is shown in Figure 3. The variable i corresponds to the i-th PPG signal sequence being processed by the proposed algorithm since it starts. fHRlow and fHRhigh are the low and high boundaries of typical human HRF. The partameter fPHR denotes the value of preliminary HRF estimation result and fHR is the final estimated HRF value. For the first signal sequence, fPHR is directly selected to be the frequency that corresponds to the largest spectral magnitude of SPPG within the human HRF range. Starting from the second signal sequence, fPHR is tracked from a range of ±Δf Because MA can be decomposed into motions in different directions, a cascaded topology of ANC is used here to reduce the motion artifacts better, as is illustrated in Figure 1. The cascaded adaptive noise cancellers use respectively x, y and z-axis acceleration signals as their MA reference inputs. Band-pass filtered PPG and acceleration signals are streamed into this step at a specified time interval ltv. In order to avoid an unmatched scale between the PPG signal and the acceleration signal, Z-score standardization is adopted to normalize the inputs of each adaptive noise canceller. Especially, if there are multiple PPG signals at the input end of the cascaded ANC, the PPG signals will be averaged into a single PPG signal before the first adaptive noise canceller.
Shown in Table 1 is the experimental results of HR AAE before and after applying the cascaded ANC stage of the proposed method. When only BPF is applied to the MA-corrupted signals, the estimation errors for HR values are very large, reaching an average HR AAE of 11.47 bpm over the 12 subjects, which indicates that the MA frequencies are strongly deviated from the HRFs and the spectral magnitudes of MA are larger. After applying the cascade ANC stage for MA reduction, the HR estimation errors significantly decreased.

Heart Rate Frequency Tracking
After the cascaded ANC stage, the de-noised PPG signal S PPG is used for a preliminary estimation of current HRF. FFT is first applied to acquire the spectrum of S PPG , where only the magnitude information is considered. A relatively large number of points for FFT L FFT is used to increase frequency resolution. Then a simple HRFT scheme from [15] is used for preliminary HRF estimation.
Flowchart of the HRFT scheme is shown in Figure 3. The variable i corresponds to the i-th PPG signal sequence being processed by the proposed algorithm since it starts. f HRlow and f HRhigh are the low and high boundaries of typical human HRF. The partameter f PHR denotes the value of preliminary HRF estimation result and f HR is the final estimated HRF value. For the first signal sequence, f PHR is directly selected to be the frequency that corresponds to the largest spectral magnitude of S PPG within the human HRF range. Starting from the second signal sequence, f PHR is tracked from a range of ±∆f around f HR of the previous sequence. When i is smaller or equivalent to N init , ∆f is set to be an initial constant ∆f0. N init corresponds to the initialization stage of the proposed algorithm that typically takes 30 s to 1 min, where the users are required to stay motionless. If the initialization stage has finished, ∆f is updated according to be the sum of the largest absolute difference of previous consecutive HRF estimation results and a bias b. Then, f PHR is selected to be corresponding to the largest spectral magnitude within a range of ±∆f around previous f HR and also within human HRF range. The calculated HR AAE results after applying the HRFT stage have been improved, as shown in Table 2.
magnitude within a range of ±Δf around previous fHR and also within human HRF range. The calculated HRAAE results after applying the HRFT stage have been improved, as shown in Table 2.

Heart Rate Frequency Correction
After the preliminary HRF estimation, a heuristic HRF correction scheme is further designed to correct possibly wrong fPHR results and give the final HRF estimation result fHR. The flow of this scheme is illustrated in Figure 4. The value of fHR is initialized as fPHR. For a relatively small time interval, the change in human HRF is not expected to be very large. So if the difference between current fHR and HRF estimation result of the previous PPG signal sequence fHR,pre becomes larger than a threshold Th0, it is highly possible that either fHR or fHR;pre is a wrong estimation result. Under that circumstance, the HRFC scheme tries to judge whether fHR is estimated to be wrong and will correct it if so (As the proposed method targets real-time online application, fHR,pre should not be corrected).
To correct a possibly wrong fHR result, the HRFC stage first tries to find fN, which corresponds to the largest current PPG spectral peak between fHR and fPHR. If fN is found, it is highly possible that fN is the correct HRF estimation result for fHR. However, it need to further check its fidelity to make the final decision. It first needs to check whether fACC is within a small range of ±Δ around fHR, where fACC is the frequency of the largest spectral magnitude of SACC, and SACC is the magnitude spectrum of acceleration signal calculated from the average of the three-axis acceleration sequences. If that is true, fHR is highly possible to be a wrong result which actually corresponds to motion artifacts, and HRFC further checks whether the spectral magnitude of fN is large enough when compared to the one of fHR so that fN can be the correct estimation result of fHR.
If fACC is not found to be around fHR, HRFC cannot decide whether fHR is a wrong result, and it further checks whether fN can be a faked HRF peak which is actually caused by motion artifacts. If fACC is found to be around fN and the spectral magnitude of fN is not large enough, HRFC keeps the value of fHR unchanged. If HRFC is not able to make a decision after all those checks, it simply compares the spectral magnitude of fN and fHR and set fHR to be fN if the magnitude of fN is large enough. For the spectral magnitude comparisons between fN and fHR, the range of thresholds Th1, Th2 and Th3

Heart Rate Frequency Correction
After the preliminary HRF estimation, a heuristic HRF correction scheme is further designed to correct possibly wrong f PHR results and give the final HRF estimation result f HR . The flow of this scheme is illustrated in Figure 4. The value of f HR is initialized as f PHR . For a relatively small time interval, the change in human HRF is not expected to be very large. So if the difference between current f HR and HRF estimation result of the previous PPG signal sequence f HR,pre becomes larger than a threshold Th0, it is highly possible that either f HR or f HR;pre is a wrong estimation result. Under that circumstance, the HRFC scheme tries to judge whether f HR is estimated to be wrong and will correct it if so (As the proposed method targets real-time online application, f HR,pre should not be corrected).
To correct a possibly wrong f HR result, the HRFC stage first tries to find f N , which corresponds to the largest current PPG spectral peak between f HR and f PHR . If f N is found, it is highly possible that f N is the correct HRF estimation result for f HR . However, it need to further check its fidelity to make the final decision. It first needs to check whether f ACC is within a small range of ±∆ around f HR , where f ACC is the frequency of the largest spectral magnitude of S ACC , and S ACC is the magnitude spectrum of acceleration signal calculated from the average of the three-axis acceleration sequences. If that is true, f HR is highly possible to be a wrong result which actually corresponds to motion artifacts, and HRFC further checks whether the spectral magnitude of f N is large enough when compared to the one of f HR so that f N can be the correct estimation result of f HR .
If f ACC is not found to be around f HR , HRFC cannot decide whether f HR is a wrong result, and it further checks whether f N can be a faked HRF peak which is actually caused by motion artifacts. If f ACC is found to be around f N and the spectral magnitude of f N is not large enough, HRFC keeps the value of f HR unchanged. If HRFC is not able to make a decision after all those checks, it simply compares the spectral magnitude of f N and f HR and set f HR to be f N if the magnitude of f N is large enough. For the spectral magnitude comparisons between f N and f HR , the range of thresholds Th1, Th2 and Th3 is (0, 1). After final estimation of HRF, the HR value (60 f HR bpm) for current PPG sequence is also calculated as an output of the proposed method. is (0, 1). After final estimation of HRF, the HR value (60 fHR bpm) for current PPG sequence is also calculated as an output of the proposed method. Figure 4. Flowchart of the heart rate frequency correction.
As shown in Table 3, the final average HRAAE over the 12 subjects is 0.92 bpm for the complete HRF estimation algorithm, with an average standard deviation (SD) of 1.50 bpm (SDs are calculated within the estimation results of every subjects and then averaged over the 12 subjects).

Notch Filtering
Cascaded notch filters are used to restore the PPG signal based on the estimated HRF. As is discussed above, the PPG signal can be treated as a quasi-periodic signal with spectral components mainly distributed as peaks around its HRF and the harmonic frequencies of HRF. Thus notch filters or comb filters can be used to eliminate unwanted components and keep only the PPG-related ones. Shown in Equation (5) is the typical transfer function of a digital notch filter [19]. In Equation (4), fn is the notch frequency and Fs is signal sampling rate. Shown in Figure 5 is the frequency domain amplitude response of the notch filter, where the parameter r controls its bandwidth. It can be observed from Figure 5 that notch filters actually suppress the signal around notch frequencies while keep the components of other frequencies close to their original amplitudes.  As shown in Table 3, the final average HR AAE over the 12 subjects is 0.92 bpm for the complete HRF estimation algorithm, with an average standard deviation (SD) of 1.50 bpm (SDs are calculated within the estimation results of every subjects and then averaged over the 12 subjects).

Notch Filtering
Cascaded notch filters are used to restore the PPG signal based on the estimated HRF. As is discussed above, the PPG signal can be treated as a quasi-periodic signal with spectral components mainly distributed as peaks around its HRF and the harmonic frequencies of HRF. Thus notch filters or comb filters can be used to eliminate unwanted components and keep only the PPG-related ones. Shown in Equation (5) is the typical transfer function of a digital notch filter [19]. In Equation (4), f n is the notch frequency and F s is signal sampling rate. Shown in Figure 5 is the frequency domain amplitude response of the notch filter, where the parameter r controls its bandwidth. It can be observed from Figure 5 that notch filters actually suppress the signal around notch frequencies while keep the components of other frequencies close to their original amplitudes. For every processing iteration (every Itv s), two notch filters are constructed using respectively current HRF estimation result fHR and the second harmonic frequency of fHR, because the main features of time domain PPG signal in one period are its main pulse and dicrotic pulse. For the first notch filter, notch frequency is set to be fHR directly, while for the second one, its notch frequency is set to be the one that corresponds to the largest spectral peak around 2 fHR due to the quasi-periodicity of PPG signals. The output of the cascaded ANC stage is filtered using the two notch filters, where the PPG components are removed while the MA is kept as it is. Through finally subtracting the notch filtering output from the ANC-de-noised PPG signal, the restored PPG signal are achieved. Figure 6 shows the performance of the proposed method to restore the PPG signal from MAcorrupted signal. The PPG signal was seriously distorted during intensive exercise, and the amplitude of the PPG component in the frequency domain was not obvious. After the proposed method was applied, the waveform of the PPG signal was restored and only the fundamental and second harmonic frequency of the PPG signal were found in the frequency domain. Because there was no reference PPG signal in the open source data set, a preliminary judgment was made through ECG that the peak number and interval of the recovered PPG signal are consistent with that of ECG.  For every processing iteration (every Itv s), two notch filters are constructed using respectively current HRF estimation result f HR and the second harmonic frequency of f HR , because the main features of time domain PPG signal in one period are its main pulse and dicrotic pulse. For the first notch filter, notch frequency is set to be f HR directly, while for the second one, its notch frequency is set to be the one that corresponds to the largest spectral peak around 2 f HR due to the quasi-periodicity of PPG signals. The output of the cascaded ANC stage is filtered using the two notch filters, where the PPG components are removed while the MA is kept as it is. Through finally subtracting the notch filtering output from the ANC-de-noised PPG signal, the restored PPG signal are achieved. Figure 6 shows the performance of the proposed method to restore the PPG signal from MA-corrupted signal. The PPG signal was seriously distorted during intensive exercise, and the amplitude of the PPG component in the frequency domain was not obvious. After the proposed method was applied, the waveform of the PPG signal was restored and only the fundamental and second harmonic frequency of the PPG signal were found in the frequency domain. Because there was no reference PPG signal in the open source data set, a preliminary judgment was made through ECG that the peak number and interval of the recovered PPG signal are consistent with that of ECG. For every processing iteration (every Itv s), two notch filters are constructed using respectively current HRF estimation result fHR and the second harmonic frequency of fHR, because the main features of time domain PPG signal in one period are its main pulse and dicrotic pulse. For the first notch filter, notch frequency is set to be fHR directly, while for the second one, its notch frequency is set to be the one that corresponds to the largest spectral peak around 2 fHR due to the quasi-periodicity of PPG signals. The output of the cascaded ANC stage is filtered using the two notch filters, where the PPG components are removed while the MA is kept as it is. Through finally subtracting the notch filtering output from the ANC-de-noised PPG signal, the restored PPG signal are achieved. Figure 6 shows the performance of the proposed method to restore the PPG signal from MAcorrupted signal. The PPG signal was seriously distorted during intensive exercise, and the amplitude of the PPG component in the frequency domain was not obvious. After the proposed method was applied, the waveform of the PPG signal was restored and only the fundamental and second harmonic frequency of the PPG signal were found in the frequency domain. Because there was no reference PPG signal in the open source data set, a preliminary judgment was made through ECG that the peak number and interval of the recovered PPG signal are consistent with that of ECG.
Shown in Table 4 is the HR AAE experimental results of the proposed method, with comparison to other published state-of-art methods from 2016 to 2018 that use the same data set for evaluation. And the proposed method achieves the lowest error, which shows the strength of the proposed method for accurate HR estimation during intensive exercise.  [20] IEEE TMBE'16 [14] MUARD [21] WFPV [15] HSUM [22] This Work  Figure 7a shows the Pearson correlation plot for all the HR estimation results. The Pearson correlation coefficient reached 0.997, showing that the estimated HR values were highly correlated with the reference HR values. Figure 7b shows the Bland-Altman plot [23] of all the HR estimation results for the proposed method, which illustrates the relationship between HR estimation errors and the mean of estimated HR and reference HR. It can be shown that though the largest HR estimation error reaches around 25 bpm, 95% of the errors are within the range of [-4 bpm, 4 bpm], which indicated a high agreement between the estimated HR values and reference HR values.

Evaluation Under Practical Experiment
Due to the lack of reference PPG signals in open source data set, we cannot evaluate the results of signal recovery well. To further verify the performance of the proposed method, a PPG evaluation platform is developed to set up the self-built data set including reference PPG signals, and then

Evaluation Under Practical Experiment
Due to the lack of reference PPG signals in open source data set, we cannot evaluate the results of signal recovery well. To further verify the performance of the proposed method, a PPG evaluation platform is developed to set up the self-built data set including reference PPG signals, and then demonstrate the effectiveness of the proposed algorithm under practical environment.
3.2.1. Implementation of Evaluation Hardware and Software Figure 8 shows the developed PPG evaluation platform, which consists of two signal acquisition finger-bands and a mobile software application. The finger-bands can simultaneously record one channel PPG signal and three-axis acceleration signals with 25 Hz sampling rate, in which the PPG sensor (SEN0203, DFRobot, China) was used to obtain PPG signal with green LEDs and the main control panel (CurieNano, Intel, USA) integrates the functions of 32-bit MCU, 12-bit ADC, six-axis motion sensor and Bluetooth Low Energy (BLE). The Figure 8a shows the signal acquisition finger-band. The software application developed with JAVA was installed on the android smartphone equipped with Octa-core (4 × 2.4 GHz Cortex-A73 and 4 × 1.8 GHz Cortex-A53), which was used to receive data from the two finger-bands and run the proposed algorithm. The restored PPG signals and estimated HR values were displayed on the mobile phone. After testing, the average running time of the proposed algorithm in the developed platform is 83 ms (the average result of 60 runs of the algorithm), which shows a great real-time performance. The parameters of the algorithm realized in the software application are set as the following:

Experimental Process
During experiment, the subject will wear two finger-bands with each of them on one thumb, as shown in Figure 8b. The left hand and arm are kept motionless, where only one channel PPG signal is recorded as the reference PPG signal with no motion artifacts. The right hand and arm performed certain movements, where one channel MA-corrupted PPG signal and corresponding three-axis acceleration signals were recorded. As left and right hands and arms were symmetric, the pure PPG signals collected from them should be highly similar. Thus the PPG signal collected from motionless left hand can be used as a reference PPG signal, which can indicate the effectiveness of our proposed method.
All the raw data mentioned above and processed results of estimated HR values and restored PPG signals are displayed on the mobile phone in real time and saved as self-built data set. A total of 20 subsets of data from 13 male subjects aging from 20 to 23 are collected, with their duration time and motion patterns shown in Table 5. The first subset is a control group where the subject stays motionless for 120 s. The second motion pattern featured relatively regular large-joint movements, while the third pattern was mainly composed of random small-joint movements.

Experimental Process
During experiment, the subject will wear two finger-bands with each of them on one thumb, as shown in Figure 8b. The left hand and arm are kept motionless, where only one channel PPG signal is recorded as the reference PPG signal with no motion artifacts. The right hand and arm performed certain movements, where one channel MA-corrupted PPG signal and corresponding three-axis acceleration signals were recorded. As left and right hands and arms were symmetric, the pure PPG signals collected from them should be highly similar. Thus the PPG signal collected from motionless left hand can be used as a reference PPG signal, which can indicate the effectiveness of our proposed method.
All the raw data mentioned above and processed results of estimated HR values and restored PPG signals are displayed on the mobile phone in real time and saved as self-built data set. A total of 20 subsets of data from 13 male subjects aging from 20 to 23 are collected, with their duration time and motion patterns shown in Table 5. The first subset is a control group where the subject stays motionless for 120 s. The second motion pattern featured relatively regular large-joint movements, while the third pattern was mainly composed of random small-joint movements. Table 5. Detailed description for the motion patterns of the practical experiment.

Subset
Number

Experiment Results
The correlation coefficient between MA-removed PPG and reference PPG signal [24] was adopted for evaluating the performance of the proposed method to restore the PPG signal. The correlation between two PPG signal sequences S 1 and S 2 can be calculated using Equation (6), where L is the length of the output PPG sequence. µ 1 and µ 2 are the averages of S 1 and S 2 . σ 1 and σ 2 are the standard deviations. The typical range of correlation coefficient was [-1, 1]. A higher value of the absolute correlation means a higher similarity in two signal sequences. In this study, we denote the correlation between MA-removed PPG and reference PPG as CORR aMAR . Denoted as CORR bMAR , the correlation coefficient between MA-corrupted PPG signal and reference PPG signal is also used for showing the improvement in PPG signal purity. Table 6 shows the experimental results for the 20 data subsets of self-built data set. The reference PPG signal or MA-corrupted PPG signal are first filtered using the BPF in Figure 1 to avoid the influence of out-of-band noise to the experimental results. In the respect of PPG signal recovery, the result of control group (subset 1) shows the correlation between the PPG signals from left hand and right hand is quite good, which means the left-hand signal can be treated as a reference for the right-hand indeed. The average CORR bMAR over the 20 subsets was 0.57, which is significantly lower than unity, showing that right-hand PPG signals were strongly corrupted by motion artifacts and the waveforms are distorted. After applying the proposed MA-removal method, the average correlation is increased to 0.86, which is close to the value of control group, indicating that PPG signal waveforms have been recovered well. In respect of HR estimation, the reference HR values were directly calculated from the magnitude spectrums of reference PPG signals. With the complete proposed method, the overall average HR AAE is 0.89 bpm, which shows the robustness of the proposed method for accurate HR estimation during different types of intensive exercise.
To illustrate the performance of the proposed method more clearly, the subset 19 (the best case) and the subset 8 (the worst case) were taken for examples to further analysis, as shown in Figures 9 and 10. For the best case, Figure 9a shows the estimated HR stays close to the reference HR and is able to follow the instantaneous changes in reference HR. Figure 9b shows the two correlation values against time for subset 19. For the first 20 s where both hands/arms were motionless, the correlation values were both high and close to unity. After motion began (walking on the treadmill for subset 19), CORR bMAR dropped significantly to the range between 0.5 and 0.6. While CORR aMAR also dropped at the transition, it recovered quickly and always stayed around 0.9 during the intensive exercise. For the worst case, the estimated HR still stayed close to the reference most of the time, as shown in the Figure 10a. The wide differences at some time are caused by the sudden change of HR. Because of the excessive movement of right arm during the data acquisition, the left arm cannot be motionless completely, which results in the corrupted PPG reference signal. and is able to follow the instantaneous changes in reference HR. Figure 9b shows the two correlation values against time for subset 19. For the first 20 s where both hands/arms were motionless, the correlation values were both high and close to unity. After motion began (walking on the treadmill for subset 19), CORRbMAR dropped significantly to the range between 0.5 and 0.6. While CORRaMAR also dropped at the transition, it recovered quickly and always stayed around 0.9 during the intensive exercise. For the worst case, the estimated HR still stayed close to the reference most of the time, as shown in the Figure 10a. The wide differences at some time are caused by the sudden change of HR. Because of the excessive movement of right arm during the data acquisition, the left arm cannot be motionless completely, which results in the corrupted PPG reference signal.

Discussions and Conclusions
PPG sensors have been widely used in wearable health tracking, from the early HR monitoring to the present BP and SpO2 monitoring. However, the accuracy of PPG sensors during exercise is unsatisfactory due to the motion artifact. At present, many studies have proposed methods to improve the accuracy of HR monitoring, which ignore the morphological recovery of the PPG signal. It limits

Discussions and Conclusions
PPG sensors have been widely used in wearable health tracking, from the early HR monitoring to the present BP and SpO 2 monitoring. However, the accuracy of PPG sensors during exercise is unsatisfactory due to the motion artifact. At present, many studies have proposed methods to improve the accuracy of HR monitoring, which ignore the morphological recovery of the PPG signal. It limits the application of PPG sensors. In this paper, a novel method is proposed for PPG MA removal during intensive exercise. On the one hand, an LMS-Newton-based cascaded ANC and unique HRF estimation scheme are applied to increase the accuracy of HR calculation. The mean AAE of HR during intensive exercise is improved from 11.47 bpm to 0.92 bpm, which verified preliminary on a widely-used open source data set of 12 subjects. On the other hand, the notch filters are employed here to recover the PPG components with the calculated HRF, which maintaining the morphological characteristics of PPG signals during intensive exercise. In the practical experiment, the proposed method achieves mean AAE result of 0.89 bpm, which shows the stability and strength of the proposed method for accurate HR calculation during intensive exercise. At the same time, the average correlation coefficient between recovered PPG signal and reference PPG signal reaching 0.86. The main features of PPG signal are restored, which is helpful to the future study of PPG sensors, such as wearable BP monitoring based on PPG. In addition, the proposed method has the advantage of real-time performance. It only costs an average of 0.83 ms to process the PPG signal with 25 Hz sampling rate. The proposed method provides a good basis for the improvement of the wearable PPG sensors.