Pulse wave digital processing and beat detection based on wavelet transform and adaptive thresholding

. This paper considers the adaptive algorithm for pulse wave processing in the presence of various physiological interferences, such as motion artefacts and baseline wander. The proposed processing technique consists of two main steps: comprehensive filtering of pulse wave based on wavelet transforms and noise-resistant pulse beats detection based on the set of pass-band filtering, nonlinear transforms and adaptive thresholding. To eliminate baseline wander from pulse wave signal we introduced a new algorithm based on the principles of adaptive noise cancellation, while the reference signal for adaptive filtering was generated by using multiresolution wavelet transform. For reducing motion artefacts and others high frequency distortions of pulse wave we suggested an approach based on soft thresholding of detailed coefficients from wavelet decomposition. The efficiency of the developed algorithm for pulse wave signal processing was assessed in comparison with the existing approaches by using mathematical modeling of pulse wave signal and occurred contaminations. The performance of the pulse beats detector was further verified for different recordings of clinical pulse wave signals from the Physionet database. © 2015 Samara State Aerospace University (SSAU).


Introduction
The recording and processing of the pulse wave signal finds wide application in the instrumental systems of cardiologic diagnostics, aimed at monitoring the heart rate and blood arterial pressure, assessing the saturation of arterial blood haemoglobin with oxygen, and studying haemodynamic processes in the human arterial bed [1][2][3].
In contrast to the ECG signal, for which the processing technique is a subject of substantial studies carried out by the researchers of different scientific profile during the last 40 years, the approaches to processing the pulse wave signals, in spite of their obvious diagnostic value, are much less known and presented in scientific literature.This fact is mainly due to the history and the field of application of the methods and instruments for pulse wave signal recording and processing in clinical practice, first of all, in clinical monitors and instruments for blood arterial pressure measurements, as well as in pulse oximeters.The software in these devices is proprietary, and this is just the reason why the algorithms used in them did not become available to the scientific community and, therefore, were not developed as intensely as the ECG signal processing techniques.The absence or insufficient development of open databases with real clinical records of this biological signal also caused reduced attention of researchers and clinical physicians to the problems of elaborating and comparative analysis of various methods for pulse wave signal automated processing.
However, there are a few publications and researches that describe algorithms for pulse wave filtering and morphological features detection [4][5][6][7][8][9][10].The most successful approaches for pulse wave filtering from baseline wander and motion artefacts are based on different applications of multiresolution wavelet transforms [4][5][6].For wearable devices and long-term monitoring using accelerometer for motions recording and principles of adaptive filtering for eliminating motion artefacts demonstrates very good results [7].
At the same time there is a need to create the novel approach for comprehensive filtering of pulse wave from baseline wander and motion artefacts, as well as for noise-resistant pulse wave beat detection.The present paper is devoted to the development and thorough analysis of high-efficiency techniques for high quality processing of pulse wave signals contaminated by various types of distortions.

Theory
Pulse wave recording by means of plethysmography or sphygmography sensors is accompanied by disturbance and noise of different nature.The strongest signal distortions are due to the artefacts having the physiological origin, i.e., caused by the movements and breathing of the examined patient.Breathing trends, present in the arterial blood pulsation signal, distort the baseline and the shape of the biological signal.Motion artefacts are random and lead to the maximal distortions of the pulse wave morphology [1,2].
The processing of the pulse wave signal against the background of physiologic artefacts faces a number of algorithmic difficulties.They are mainly due to the random character of the noise and the fact that their frequency components overlap with the fundamental band of the biosignal, which makes it difficult to use classical methods of linear frequency filtering and causes the necessity for creating nontrivial processing methods.
At present, the filtering based on the multiresolution wavelet transforms is a promising line of research in the field of processing the biological signals, perturbed by broadband stochastic noises [4-6, 11, 12].The signal processing procedure based on using discrete wavelet transforms includes the signal decomposition and signal reconstruction stage.
The signal decomposition is implemented in accordance with the Mallat algorithm and consists in the decomposition of the initial signal into a sequence of approximation and detailed coefficients.The signal filtering to get rid of the noise can be implemented by direct modification of the detailed coefficients of the wavelet decomposition [12].After the signal decomposition stage and the removal of noise components, the reconstruction of the initial signal is carried out.

Materials and methods
In the present paper, to make the pulse wave signal free of physiologic artefacts, we propose the wavelet filtering technique based on the discrete decomposition in terms of orthogonal wavelets that includes the following processing stages: 1) calculation of the direct discrete wavelet transform of the pulse wave signal; 2) modification of the detailed coefficients of the wavelet decomposition in correspondence with the threshold processing parameters; 3) reconstruction of the pulse wave signal basing on the initial approximation coefficients and modified detailed coefficients by means of the inverse wavelet transform.The efficiency of the proposed pulse wave signal filtering technique based on the multi-resolution wavelet transforms depends upon the form of the mother wavelet function, the level of de-composition, the type of the threshold function for processing the specifying coefficients, and the algorithm for threshold determination.
The adaptive determination of the threshold level T for each level of the wavelet decomposition was performed using the Donoho formula [12]: where σ is the root-mean-square deviation of noise, evaluated for the sequence of detailed coefficients at the given level of the wavelet decomposition; M is the total number of the signal samples.
where x is the input value of the detailed coefficients, y is their output value, sign(x) is the function determining the sign of the number x, and T is the value of the adaptive threshold.
To study the efficiency of the pulse wave signal wavelet filtering technique we used the mathematical modelling of present distortions.The model of the pulse wave signal, distorted by the baseline wander and motion artefacts, was assumed to be additive.To get the model dependences of the pulse wave, we used the simulation model proposed by P.E.McSharry [13] that allowed the formation of biosignal fragments with the required morphology and with the given values of amplitude-time parameters.
The analysis of the structure and factors, affecting the appearance of the pulse wave baseline wander, has shown that this type of noise is a low-frequency signal of mainly stochastic nature, which can be described in the form of additive combination of deterministic and random components: where W max is the amplitude of the baseline wander model signal; ψ(k) is the random component, obtained by filtering the white Gaussian noise with zero mean value and unit variance with a low-frequency filter with the cut-off frequency 1 Hz; f i is a set of harmonic signal frequencies, representing the deterministic component of the baseline wander.In the process of simulation, the following values were used: f1…f4 =0.1 Hz, 0.2 Hz, 0.4 Hz, and 0.6 Hz.
The experiments with the accelerometer sensor, fixed on the arm of the investigated human, followed by the spectral analysis of the recorded motion distortion signals, have shown that the frequency region of motion artefacts, accompanying the pulse wave recording at different motion activity (walking, gesticulation, running) ranges from 0 to 5 Hz [7].
As a mathematical model of motion artefacts, we used the Markov process model, or the first-order autoregressive model, described by the following recursion expression: where z max is the amplitude of the motion artefact model signal; ε(k) is a sequence of random numbers, distributed in accordance with the normal law with the zero mean value and the unit variance; a and c are constants; z(k) is the current sample of the Markov process; z(k-1) is the previous sample of the Markov process.
Since the spectral band of motion artefact frequencies observed under the real conditions is limited, the model signal z(k) was processed with a lowfrequency filter (LFF) having the cut-off frequency 5 Hz.
As one of the criteria for estimating the efficiency of filtering the pulse wave signal, we proposed the distortion coefficient δ of the biosignal after passing the relevant processing stages: where i is the signal sample's number; M is the numbers of samples in the considered signal fragments; Y f (i) is the sample of the distortion-free model pulse wave signal.
The optimal parameters for processing the pulse wave signals were chosen using the criterion of minimizing the value of the signal distortion coefficient within the range of variation of the signal-to-noise ratio, estimated as follows: where S max is the total spectral power of the initial model pulse wave signal; X max is the total spectral power of the additive noise.
One of the basic problems of the pulse wave signal digital processing is to detect the fiducial points, mainly the systolic maxima, for further determination of heart rhythm parameters and measurement of the biosignal amplitudes for calculating the diagnostic indices of the cardiovascular system.Increasing the reliability of the performed diagnostics is possible at the expense of reducing the measurement error for the amplitude-time parameters of the biosignal, depending on the precision and efficiency of the implementation of the algorithms that detect the pulse wave fiducial points.
At present, the issues of minimizing the errors, caused by the inaccuracy of detecting the pulse wave fiducial points, remain to be insufficiently studied.The major requirement to the means of detecting the fiducial points is the possibility of efficient operation under the conditions of high-intensity noises and high variability of the biosignal shape.As fiducial points for the pulse wave signal, they usually choose the most recognizable points against the background of the noise, namely, the systolic maximum of the signal, the minimum of the pulse wave signal, and the maximum of the first derivative signal [8][9][10]14].
To provide high-efficiency detection of the pulse wave beat detection, we propose here a detector with the adaptive threshold algorithm.At the stage of preliminary processing the operations of frequency band-filtering, differentiation, and nonlinear transformation are applied; the maximum of the biosignal first derivative is chosen to be a reference point for further localization of the pulse wave beats or systolic maximums.
To improve the accuracy of pulse wave signal processing and pulse beats detection we preliminary applied digital interpolation to the recorded biosignals by using cubic splines, thus the effective sampling frequency after interpolation was selected at 2000 Hz.
To implement the frequency filtering of the pulse wave signal we used the Butterworth digital filter of the 8-th order with the pass-band from 0.5 to 10 Hz.To linearize the phase characteristic, the signal output from the Butterworth filter was again passed through the filter, but in the reversed direction [15].
Then the pulse wave signal was differentiated using the operator: where n is the number of the signal sample; Δ d is the sampling interval; Ppg(n) is the pulse wave signal after the band-pass frequency filtering.
The signal obtained after the differentiation is raised to the third power, which additionally enhances the large differences of the signal P(n).Then the samples with negative amplitude values are removed from the obtained signal.
Figure 1 presents the time dependences of the pulse wave signal at different stages of processing: curve 1 is the model signal of the pulse wave Ppg(n), curve 2 is the signal after the differentiation P(n), curve 3 is the signal raised to the third power P 2 (n), and curve 4 is the signal after selecting the positive samples P 3 (n).
The output signal P 3 (n) after passing all stages of the preliminary processing arrives at the input of the circuit executing the threshold pulse beat detection.The detection unit forms a moving window with the duration of 2 seconds; for each window the threshold unit determines the value of the threshold (Lev).Only those signal samples are input into the maximum detector, for which the condition P 3 (n) > Lev is fulfilled.
The running window duration is chosen such that at least one pulse beat should be present in the signal during this time.Since the periodicity of the pulse beats is determined by the heart rate, the minimal value of which is not less than 30 beats per minute, the appropriate detection interval should not exceed 2 seconds.
The detector of the maximum determines the position of the maximum of the pulse wave signal first derivative using the three-point detection method:

P n Max if P n Lev
P n P n P n P n If this logical condition is fulfilled, the number of the sample n determines the temporal position of the pulse wave signal derivative maximum.
The efficiency of the proposed pulse beat detection scheme is largely dependent on the choice of the threshold value Lev.The threshold value is independently determined for each moving window basing on the following threshold function: where Ω(i) is the value of the root-mean-square deviation of the signal samples amplitude value within the limits of the i-th moving window; Max(i) is the maximal value of the amplitude of the signal samples within the limits of the i-the moving window; Max(i-1) is the similar maximum within the pre-ceding (i-1)-th moving window.
The numerical values of the threshold function parameters were chosen empirically as a result of the preliminary processing of the noisy model signal using the criterion of maximizing the correct detections of pulse beats and minimizing the false detections and misses.
The described algorithm provides the detection of fiducial point, i.e. the maximum of the pulse wave's first derivative in our case.At the second stage detector uses the information provided by the previous step.A defined width window subset (±100ms form the location of the peak found in the pulse wave's first derivative) is selected in the original filtered pulse wave to locate the real systolic peak.Once again a simple maximum peak locator in the values of this subset sequence is used.
The detection of pulse beats is implemented with a certain error, caused by the imperfection of the detection scheme, as well as by the distortions of the biosignal.This error may produce a random error in measuring the pulse wave amplitude and the duration of the beat-to-beat intervals.As one of the criteria for estimating the efficiency of the scheme for detecting the pulse beats the deviation of measuring beat-to-beat intervals from the true value with the confidence level P=0.9 could be chosen: , 6 , 1 where Δ и is the absolute error of measuring beat-to-beat intervals, σ u is the root-mean-square deviation of the beat-to-beat intervals from the true value, determined as where PP(i) is the true value of the i-th beat-to-beat interval, PP'(i) is the measured value of the i-th beat-tobeat interval, N is the total number of beat-to-beat intervals in the analyzed recording of the biosignal.
One more criterion for evaluating the efficiency of the pulse beat detection scheme is the random deviation of the measured amplitude from the true value, also determined at the confidence level P=0.9: where δ A is the relative error of the pulse wave amplitude measurement; Am is the true value of the pulse wave amplitude; σ is the root-mean-square deviation of the pulse wave amplitude from the true value, determined as where A m (j) is the true value of the j-th pulse wave amplitude; A m '(j) is the measured value of the j-th pulse wave amplitude; N is the total number of pulse wave fragments in the considered recording.

Results
The results of the studies were obtained using the standard functions of Signal Processing Toolbox MATLAB R2013a.
Table 1 presents the values of the pulse wave distortion coefficient δ for different levels of wavelet decomposition using the hard and soft threshold function.The signal-to-noise ratio is K а =8 dB; the mother wavelet function is the Daubechies wavelet of the 6-th order (db6); the heart rate is 70 beats per minute, the pulse wave amplitude is 1; f s =2000 Hz; c=0; a=0.5.
The analysis of the obtained results shows that the number of the wavelet decomposition levels, optimal from the point of view of minimizing the pulse wave distortions, should be equal to 6, when the soft threshold function is used for processing the detailed coefficients.
Figure 2 shows a fragment of the synthesized model pulse wave signal, distorted by strongly expressed motion artefacts, and the result of applying the technique of the noisy biosignal wavelet filtering.The time dependences were obtained for the following values of the model parameters: the heart rate is 70 beats per minute, the amplitude of pulse wave is 1; f s =2000 Hz; c=0; a=0.5;W max =0; z max =0.3.
Table 1 The values of the pulse wave distortion coefficient δ for different levels of wavelet decomposition using soft and hard threshold function.The obtained results demonstrate high quality of processing of the noisy model pulse wave signals using the multi-resolution wavelet transforms.
Figure 3 presents a fragment of clinical record of the finger pulse wave signal, distorted by the baseline wander and motion artefacts, caused by the patient's hand tremor.Figure 4 shows a fragment of the pulse wave signal processed using the proposed wavelet filtering technique.
The analysis of the pulse wave signal fragment shown in Fig. 4 demonstrates that the developed technique of wavelet filtering efficiently removes the motion artefacts, but has practically no effect on the baseline wander present in the signal.Additional studies have shown that the sequence of approximation coefficients, obtained at high levels of wavelet decomposition, has the frequency spectrum, similar to that of the pulse wave baseline wander.
The optimal parameters of wavelet decomposition for removing the pulse wave baseline wander were determined basing on the analysis of the coefficient of correlation between the model signal of baseline wander and the sequence of coefficients in the wavelet decomposition of the model pulse wave signal with the additive baseline wander present.The numerical values of the correlation coefficient are presented in Table 2.
The analysis of the results presented in Table 2 shows that in order to correct the baseline wander of the pulse wave it is optimal to use the sequence of approximation coefficients of the 6-th decomposition level.Table 2 The values of the coefficient of correlation between the model baseline wander of the pulse wave and the sequence of wavelet decomposition coefficients.Figure 5 illustrates the correction of the pulse wave baseline wander basing on the subtraction of the interpolated sequence of approximation coefficients of the 6-th level of wavelet decomposition from the biosignal, distorted by the baseline wander (Fig. 4).The interpolation of the sequence of approximation coefficients is necessary due to the specific feature of the Mallat algorithm implementation, according to which at each stage of the wavelet decomposition the decimation with the factor 2 is carried out 2 [4].The efficiency of the pulse beat detection was studied for different types of detectors, namely, the detector based on the moving average filtering (MAF) [12]; the detector based on band-pass filtering (BPF) and quantile detection thresholds [13]; the detector based on wavelet transforms (WT) [8]; the adaptive detector (AD) proposed in the present paper.
Figures 6 and 7 present the error of measuring the beat-to-beat intervals (Δ и ) and the relative error of measuring the pulse wave amplitude (δ A ), respectively, versus the signal-to-noise ratio K а for different detectors of the pulse beats: 1 -BPF, 2 -MAF, 3 -WT, and 4 -AD.The dependences are obtained at the following parameters of the model: the heart rate 80 beats per minute; f p =50 Hz; f s =2000 Hz; a=0.5; c=0; the noise amplitudes W max and z max were specified in correspondence with the value of the signal-to-noise ratio K а .Fig. 6 The error of beat-to-beat intervals measurement versus the signal-to-noise ratio K а (1 -BF, 2 -RAF, 3 -WT, and 4 -AD).where N T is the number of correctly detected pulse beats (systolic maxima); N F is the number of false detected pulse beats; Nm is the number of missed pulse beats.
The computer processing of the signals was carried out for multichannel records of biological polysomnography signals from the database MIT-BIH Polysomnographic Database (http://physionet.org),containing 18 fragments of pulse wave signals, each 2 hours long.Separately, one recording with low noise (Slp45) and one recording with very intensive noise (Slp67x) were selected for the studies.
Table 3 presents the results of the quantitative evaluation of the efficiency of different detectors of systolic maxima in application to the processing of clinical pulse wave signals.Table 3 Efficiency evaluation for the detection of pulse wave signal systolic maxima.

Conclusions
As a result of the performed studies, the method for complex digital processing of the pulse wave signal was developed and studied.The method is based on multiresolution wavelet transforms and includes the filtering of baseline wander and motion artefacts, as well as the noise resistant fiducial point detection.
The studies have shown that the adaptive detector of systolic maxima, proposed here, provides higher accuracy as compared to the existing approaches, namely, the probability of correct detection amounts to not less than 99%, the probability of erroneous detection does not exceed 0.5%, the random error for the measurement of the pulse wave amplitude is less than 5%, and the random error for the measurements of beatto-beat intervals is less than 5 ms.
The high quality of the pulse wave signal filtering and perfect accuracy of detecting the pulse beats can serve as a reliable base for developing high-efficiency algorithms and hardware-software cardiology diagnostic complexes.

Fig. 1
Fig. 1 Time dependences of the pulse wave signal at different stages of the pulse beat detection; amplitude values are expressed in relative units.

Fig. 2 A
Fig. 2 A fragment of the pulse wave model signal with the motion artefacts (top); the result of wavelet filtering (bottom); amplitude values are expressed in relative units.

Fig. 3 A
Fig. 3 A fragment of pulse wave signal, distorted by the baseline wander and motion artefacts; amplitude values are expressed in relative units.

Fig. 4
Fig.4The fragment of pulse wave signal after the application of wavelet filtering procedures; amplitude values are expressed in relative units.

Fig. 5 A
Fig. 5 A fragment of the pulse wave signal after correcting the baseline wander (A); the estimate of the baseline wander obtained basing on the wavelet decomposition (B); amplitude values are expressed in relative units.

Fig. 7
Fig. 7 The amplitude measurement error versus the signal-to-noise ratio K а (1 -BPF, 2 -MAF, 3 -WT, and 4 -AD).As follows from the obtained dependences, the reduction of the signal-to-noise ratio leads to the exponential growth of errors for the beat-to-beat intervals and pulse wave signal amplitude measurements.The minimal error of measuring the amplitude-time parameters of the pulse wave under the conditions of baseline wander and motion artefacts of different intensity is demonstrated by the pro-posed adaptive detector (AD), with preliminary processing based on the multi-resolution wavelet transforms.To analyze the efficiency of the pulse beat detection in the processing of clinically recorded biosignals we used the Physionet open database of the Massachusetts Institute of Technology.The following statistical indicators were used for quantitative evaluation of the detection efficiency: 1) the probability of correct detection of pulse beats (PT): %; 100 ⋅ = N N P T T 2) the probability of erroneous detection of pulse beats (PF): %; 100 ⋅ = N N P F F 3) the index of detection error level Per: %, 100 ⋅ + = N N N P F m er

Figure 8
Figure 8 presents the block diagram of the considered method of the pulse wave signal digital processing method.The notation is as follows: 1 -initial pulse wave signal; 2 -output pulse wave signal after the removal of physiological artefacts; 3 -array of detected fiducial points of the pulse wave; WDU -wavelet decomposition unit; АC -approximation coefficients; DC -detailed coefficients; TP -threshold processing; IWT -inverse wavelet transform; AC 6 -approximation coefficients of the 6-th decomposition level; A -adder; BF -bandpass filtering; DU -differentiation unit; NO -nonlinear operation unit; MWF -moving window formation; ATF -adaptive threshold formation unit; TU -threshold unit; MD -maximum detector.

Fig. 8
Fig. 8 Block-diagram of pulse wave signal processing