ECG-Derived Respiration Using a Real-Time QRS

Detector Based on Empirical Mode Decomposition

Christina Kozia Randa Herzallah David Lowe
Mathematics Department Mathematics Department Mathematics Department
Aston University Aston University Aston University
Birmingham, UK Birmingham, UK Birmingham, UK
koziac@aston.ac.uk r.herzallah@aston.ac.uk d.lowe@aston.ac.uk

Abstract—Respiration Rate (RR) is an important physiological been identified as high risk up to 24 hours before the event
indicator and plays a major role in health deterioration monitor- with a specificity of over 95%. The respiration signal can be
ing. Despite that, it has been neglected in hospital wards due to recorded using the following methods: spirometry, pneumogra-
inadequate nursing skills and insufficient equipment. ECG signal,
which is always monitored in a clinical setting, is modulated phy, plethysmography, or capnography. The main disadvantage
by respiration which renders it a highly enticing mean for the of these methods is that they use expensive equipment and
automatic RR estimation. In addition, accurate QRS detection is make the patients’ hospitalisation uncomfortable.
pivotal to RR estimation from the ECG signal. The investigation
of QRS complexes is a continuing concern in ECG analysis A. Respiratory-Induced Modulation of ECG
because current methods are still inaccurate and miss heart
beats. This paper presents a frequency domain RR estimation Research into extraction of respiration signal using the ECG
method which uses a novel real-time QRS detector based on has a long history. In [3] the use of the ECG in respiration
Empirical Mode Decomposition (EMD). Another novelty of the
proposed work stems from the RR estimation in the frequency acquisition was suggested by pointing out that a normal respi-
domain as opposed to some of the current methods which rely ration cycle affects the heart rate, which causes sinus arrhyth-
on a time domain analysis. As will be shown later, the RR mia. More recently, it was observed that the ECG signal is
extraction in the frequency domain provides more accurate affected by events occurring during the breathing process [4].
results compared to the time domain methods. Moreover, our A change of heart-to-electrode distances is observed during
novel QRS detector uses an adaptive threshold over a sliding
window and differentiates large Q- from R-peaks, facilitating a the thoracic expansion which modulates the QRS morphology.
more accurate RR estimation. The performance of our methods Furthermore, a variation of amplitude, frequency, and phase is
was tested on real data from Capnobase dataset. An average observed because the ECG is affected by changes in thoracic
mean absolute error of less than 0.5 breath per minute was impedance as air fills spaces in lungs.
achieved using our frequency domain method, compared to
6 breaths per minute of the time domain analysis. Moreover,
our modified QRS detector shows comparable results to other B. ECG-Derived Respiration Methods
published methods, achieving a detection rate over 99.80%. The correlation between the RR and the ECG signal leads
Index Terms—ECG-derived respiration, Frequency domain
analysis, R-peak detection, Empirical Mode Decomposition
us to RR estimation through ECG. ECG signal is monitored
(EMD), Local Signal Energy routinely in a hospital setting, it is non-invasive, inexpensive,
and safe. In recent years there has been an increasing amount
I. I NTRODUCTION of literature on ECG-derived (EDR) respiration methods [4]
Many of the literature since the mid-1990s emphasises the [5]. The alternation of QRS morphology facilitates the respi-
importance of the respiration signal and its derivative, the ration signal extraction because it is related to the breathing
Respiration Rate (RR). In [1] it was indicated that an RR cycle.
higher than 27 breaths per minute is the most important In [5] an EDR method which is based on peak-to-trough
predictor in failure of the heart to contract effectively in QRS amplitude was suggested. A single lead ECG signal is
hospitals. In [2] the necessity of RR was investigated. It was given as input to the algorithm. The R-peaks of the signal are
claimed that 21% of hospitalised patients with an RR of 25- detected and the peak-to-trough amplitude is measured. As
29 breaths per minute assessed by a critical care outreach soon as all R-peak amplitudes are measured, a cubic spline
service died in hospital. An increase of the mortality rate has interpolation is attempted followed by a filtering stage in order
been reported for patients with higher RR [2]. It has been to extract the respiration signal. In [4] the so-called Respira-
demonstrated that just over the half of unhealthy subjects tory Sinus Arrhythmia (RSA) derived respiration method was
suffering a serious event in the hospital wards has a RR greater proposed. The latter uses the R-peak time locations in order
than 24 breaths per minute and these subjects could have to compute the R-R intervals. The instantaneous heart rate
(IHR) values are then computed. IHR is in fact the inverse of
This work is supported by Isasnys Lifecare. R-R intervals. After successful computation of the IHR values,

978-1-5386-5602-0/18/$31.00 ©2018 IEEE

the latter are interpolated using cubic splines, followed by a needs to decide if the i-th component, hi (t), extracted from the
filtering stage in order to derive the respiration signal. data, is an IMF based on the conditions mentioned above. In
After successful extraction of the respiration signal, the RR order to achieve this, the EMD method uses a systematic way
estimation is attempted, using a time domain analysis method. which is called the sifting process and is described as follows:
In [6] it is suggested to use a basic peak detector in order to for a given signal x(t), the extrema points are first identified,
extract the peaks in the respiration signal, along with their time followed by approximation of the upper, rb(t), and lower, r(t),
locations. The instantaneous respiration rate (IRR) is evaluated envelopes of the signal through a cubic spline interpolation.
from the time intervals of consecutive peaks in the respiration The mean is then obtained, and the i-th component, hi (t),
signal. As soon as the IRR values are computed, they are is computed as the difference between the signal and the
converted to breaths per minute to get the RR. mean. The sifting process has to be repeated as many times
as required to reduce the extracted signal to an IMF. For
C. QRS Detection our implementation in order to terminate the EMD algorithm,
The QRS complex is the most distinguishable feature of the the number of zero-crossings and the number of extrema are
ECG signal and comprises valuable information which facili- checked for equality or whether they differ at most by one. If
tates the estimation of the RR [3] and Heart Rate Variability the final residue, rN (t), is obtained as a monotonic function,
[7]. R-peak detection complexity results from signal contam- the sifting process is stopped, ci = hi , and the signal, x(t),
ination due to noise and the QRS morphological variability can be written as follows:
due to respiration. Furthermore, the complexity of the QRS N
identification stems from the difficulty of differentiating R-
X
x(t) = ci (t) + rN (t) , (1)
peaks from large P or T peaks [8]. i=1
A great deal of previous research has focused on the QRS
detection using a variety of approaches from derivative based where N is the total number of the extracted IMFs.
methods [8] [9] to neural networks [10] [11] methods. The
E. Overview of the Proposed EDR method
majority of the algorithms first pre-process the signal, in order
to enhance the QRS complex and reduce noise, and then apply Our EDR method first locates the QRS complexes. For each
a set of thresholds for the R-peak identification. Signal differ- QRS complex the peak-to-baseline amplitude is computed
entiation [8] [9] and Hilbert transform [12] [13] have been and outliers due to noise and artifact are discarded. As soon
proposed as pre-processing methods for QRS amplification. as the peak-to-baseline values are evaluated, the sequence is
Most recently, Empirical Mode Decomposition (EMD) [14] upsampled using a cubic spline, followed by a filtering stage
has been proven to be an effective QRS pre-processor [15] in order to derive the respiration signal. The latter is then
[16] [17]. The R-peaks are finally located by a set of thresholds divided into one minute windows in order to estimate the
such as: estimators of signal and noise level like the median [8] RR. For each window the Fourier transform of the respiration
or the Root Mean Square (rms) [12], the heart refractory period signal is computed. The next step requires the identification
(200 miliseconds) [18], or thresholds based on the maximum of the frequency that corresponds to the most dominant
amplitude of the signal [16] [17]. A major weakness of the peak in the frequency spectrum, which is assumed to be the
detectors being proposed in [15], [16] and [17] is that the respiration rate in the current window. The novelty of our
threshold is derived from the full length ECG. Ectopic R- EDR method lies on the fact that the RR estimation takes
peaks can make the threshold high which results in the failure place in the frequency domain, whereas time domain analysis
to detect lower R-peaks. Our QRS detector provides a solution is sensitive to dicrotic notches, which makes RR estimation
to the latter by dividing the signal into segments. An additional inaccurate. Moreover, our frequency analysis approach can be
drawback of [15], [16] and [17], which emerges from the use implemented in real-time as the estimation requires windows
of the full length ECG for the threshold computation, is that of one minute duration.
these methods cannot be implemented on-line. In contrast, we Our proposed QRS detector overcomes the aforementioned
propose a real-time QRS detector which only makes use of problems in the current state-of-the-art EMD methods by intro-
the most recent history of the ECG signal. ducing an adaptive threshold which is calculated from the local
energy of the reconstructed ECG signal from the EMD. The
D. Empirical Mode Decomposition signal is first pre-processed by a bandpass filter and EMD in
EMD decomposes the signal, x(t), into a series of narrow- order to eliminate noise, baseline wander and enhance the QRS
band signals, ci (t), which are called IMFs, and fulfill special complexes. Our R-peak detector has a number of attractive
conditions. An oscillatory mode of the signal is an IMF ex- features, compared with the existing results on the topic.
clusively under the conditions that: first, in the whole dataset, Firstly, our adaptive threshold, derived from the local signal
the number of zero-crossings and the number of extrema are energy of a prespecified number of most recent segments,
either equal or differ at most by one; and second, at any prevents the failed detection of small R-peaks. Secondly, the
point, the mean value of the maximum and the minimum present approach explores the differentiation of R-peaks from
envelope is zero. The key advantage of EMD is that it is a large Q peaks in the absolute value of the signal, by using
data-driven analysis method. In each iteration the algorithm the first derivative of the ECG signal. Thirdly, our method

978-1-5386-5602-0/18/$31.00 ©2018 IEEE

proposes an on-line algorithm with a detection delay equal to Plot 1
PBA value
the duration of the segment and a low memory usage. 1

ecg
0
R peaks

Amplitude (mV)
II. M ETHODS -1 EDR spline
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
A. Proposed Frequency Domain ECG-Derived Respiration Plot 2
10 -5
1

The proposed RR estimation method is based on the fact that 0.5

the most dominant peak of the respiration frequency spectrum 0

-0.5
corresponds to the RR. In order to derive the respiration signal, -1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
the R-peaks are located using our proposed QRS detector Time (minutes)

(Section II-C). In contrast to existing EDR methods [4] [5]

Fig. 1. Plot 1 shows our proposed EDR method for a small part (0.7 minutes)
the peak-to-baseline amplitude (PBA) for each QRS complex of an ECG. R-peaks are represented by black dots. The output of the EDR
is computed because the ECG signal is first filtered in order to interpolation is represented in a red curve. Plot 2 shows the derived respiration
remove baseline wander and noise during the R-peak detection signal for a one minute window corresponding to the ECG of Plot 1.
(Fig. 1, Plot 1). In order to increase the accuracy of the RR
estimation, outliers are discarded to minimise the presence of Frequency spectra of one minute respiration signal
0.4

false positives due to R-peak detection. Then the PBA values 0.35
RR = 0.3 * 60 = 18 BPM
are upsampled in order to derive the EDR waveform. The latter
0.3
is then filtered within reasonable respiration frequencies in
order to extract the final respiration signal (Fig. 1, Plot 2). The 0.25

Magnitude
next step is the frequency domain analysis of the respiration 0.2

signal where a Fourier transform is applied to one minute 0.15

windows of the respiration in order to estimate the RR (Fig.

0.1
2). To summarise, the proposed EDR method is as follows:
0.05
1. Extract the QRS complexes from a single-lead ECG
signal, 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

2. For each QRS complex compute the peak-to-baseline Frequency (Hz)

amplitude (PBA),
Fig. 2. The frequency spectrum of the one minute window respiration signal
3. Discard outliers by restricting PBA values to remain of Fig. 1, Plot 2. The dominant peak is located (black dashed line) and it is
within the range of ± 2 SDs from the mean value, converted to breaths per minute (bpm).
4. Upsample the remaining PBA values at 8 Hertz using
10 -5
cubic spline interpolation, 1

5. Filter the output of step 4 within reasonable respiration

Amplitude (mV)

0.5

frequencies (0.0666-0.5 Hertz) to get the respiration sig- 0

nal,
-0.5
6. Divide the respiration signal into one minute windows, respiration
peaks
7. For each window compute the Fourier transform of the -1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

respiration signal, Time (minutes)

8. Locate the frequency that corresponds to the most dom-

Fig. 3. Time domain analysis for one minute window. The peaks used in RR
inant peak in the frequency spectrum, estimation are represented by black triangles.
9. Convert frequency to breaths per minute to get the
respiration rate.
3. Compute the IRR values,
B. Time Domain ECG-Derived Respiration 1
IRRi = 60 · , (2)
One popular approach for the RR estimation is based on ti+1 − ti
the assumption that the time interval between two peaks in where i = 1, . . . , n − 1,
the respiration signal, is the interval between two breaths [6]. 4. For every one minute window, take the average of the
In other words, the time domain analysis of the respiration IRR values within that window in order to get the RR in
attempts to find the sampling interval between the peaks of the breaths per minute.
respiration signal and convert it to breaths per minute (Fig. 3).
To summarise, the time domain analysis for the RR estimation C. Proposed R-peak Detection
is as follows: The main assumption of the proposed R-peak detector is
1. Derive the respiration signal as proposed in Section II-A, that the QRS complex can be amplified by reconstructing the
2. Extract the time locations of the peaks of the respiration signal from the first three IMFs of the EMD. This assumption,
signal, that is t1 , t2 , . . . , tn , where n is the total number which will be discussed in section III, was verified on all
of peaks in the respiration signal, the tested recordings. The signal is first pre-processed by a

978-1-5386-5602-0/18/$31.00 ©2018 IEEE

bandpass filter and EMD to reduce baseline wander, noise and 8. The threshold of the k-th segment is set to be the mean
enhance the QRS complexes. Then the reconstructed signal of the most recent eight RM Sk values,
is divided into a number of segments. The envelope of the k
maxima of each segment is approximated and followed by the 1 X
Tk = RM Sj , (5)
computation of the local signal energy and an averaging step 8
j=k−7
for the evaluation of the threshold. The final R-peak decision
requires a gradient-based and a refractory period checks to 9. The peaks, which exceed threshold Tk in the absolute
prevent false detection. The chosen segment duration provides sequence ak (t), are classified as candidate peaks.
an adequate number of QRS complexes and depends on the 10. In order to segregate large Q peaks from R peaks,
sampling frequency. Moreover, for the averaging step the eight we compute the first derivative of xr (t). Peaks with
most recent segments were used. The number of segments to a negative derivative will be investigated further at the
be averaged is a prespecified parameter which can be decided refractory period check given next,
based on the clinical condition of the patient and whether it 11. Apply a refractory period check. When the R-R interval
is expected that their ECG signal is going to be less or highly of two adjacent peaks is less than 200 milliseconds, keep
variant. However, it is recommended that just the most recent the peak with the maximum amplitude.
history of the ECG vital signs are kept. To summarise, the III. R ESULTS AND D ISCUSSION
proposed QRS complex detection algorithm is as follows:
A. ECG and RR Reference Data
Pre-processing Stage For the evaluation of our proposed methods, 20 children
1. The raw signal, x(t), is first filtered with a band-pass recordings (median age: 9.1, range: 1-16 years) and 5 adult
filter, whose coefficients were designed using a Kaiser- recordings (median age: 46.2, range: 37-64 years) were used
Bessel window [19]. The band-stop frequencies were set from the Capnobase dataset [21]. The database contains ECG,
at 8 and 20 Hertz [20] in order to amplify the QRS pulse oximetry and capnography recordings of 8 minute dura-
complex, eliminate noise and reduce the number of IMFs. tion acquired during elective surgery and routine anaesthesia,
The output of the filter is denoted as xf (t), sampled at 300 Hertz. Moreover, all recordings under study
2. The EMD method is applied on xf (t) to extract the provide annotated R-peaks and reference RR.
IMFs, c1 (t) . . . cN (t), where N is the total number of
B. ECG-Derived Respiration
the extracted IMFs,
3. The signal is reconstructed by adding the first three IMFs, For all the recordings the results are shown in Table I.
In order to evaluate the performance of our proposed EDR
3
X method, the Mean Absolute Error (MAE) in breaths per minute
xr (t) = ci (t) , (3)
(bpm) was calculated as follows:
i=1
N
where the number of IMFs is experimentally selected and 1 X
M AE = |RRrefi − RResti | , (6)
it will be discussed later, N i=1
4. Then, the absolute value of the reconstructed signal is
computed, that is a(t) = |xr (t)|. This makes all data where RRref is the reference RR, RRest is the estimated RR,
points positive and implements a linear amplification of and N is the total number of one minute windows.
the reconstructed signal emphasising the higher frequen- As can be seen from Table I our EDR method shows a
cies. promising performance for most of the recordings, achieving
an average MAE of 0.4150 breaths per minute. Fig. 4 shows
Decision Stage the respiration signal obtained from our method for recording
5. In order to increase the efficiency of the algorithm, we 0009. Each plot corresponds to one minute window. As soon
divide a(t) into k segments of 3 seconds duration, that is as the respiration signal is ready, our method analyses the
k = (total number of samples)/(3∗fs). The starting point latter in the frequency domain in order to estimate the RR.
of the k-th segment should match the last R peak located Fig. 5 shows the frequency analysis of the respiration signal
in the k − 1 segment in order to increase the accuracy, of recording 0009. As we have already mentioned, the ECG
6. Compute the envelope of the maxima, âk (t), of ak (t) for signals duration, from Capnobase dataset, is 8 minutes, thus
each segment k through a cubic spline interpolation of we expect to see the frequency spectra of 8 one minute
the local maxima, windows. Plots 1 to 8 of Fig. 5 show the frequency domain
7. Compute the local signal energy for each segment as, analysis for the full length respiration signal, along with the
v reference RR. It can be seen that each window contains a
u
u 1 X M
dominant peak which corresponds to the RR.
RM Sk = t [âk (t)]2 , (4)
M t=1 The largest MAE, which was obtained for recording 0032, is
about 2.6240 breaths per minute. As can be seen in Fig. 6, the
where k is the current segment and M is the number of respiration signal shows smaller oscillations which can overlap
samples in the segment, that is M = 3 ∗ fs, the actual respiration, and it does not have a smooth shape as

978-1-5386-5602-0/18/$31.00 ©2018 IEEE

 TABLE I
ECG-D ERIVED R ESPIRATION E VALUATION P ERFORMANCE Plot 1 - Window 1 Plot 2 - Window 2
0.4 0.2
18 bpm 17 bpm
Frequency Domain Analysis Time Domain Analysis 0.2 Ref: 17.79 bpm 0.1 Ref: 17.60 bpm
Children Recordings MAE (bpm) MAE (bpm) 0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0 0.1 0.2 0.3 0.4 0.5 0.6
0009 0.3716 0.7478 Plot 3 - Window 3 Plot 4 - Window 4
0015 0.2086 0.0710 0.1 0.04
0016 0.0048 11.2462 18 bpm 22 bpm
0.05 0.02

Magnitude
Ref: 18.36 bpm Ref: 21.46 bpm
0018 1.2627 9.0078 0 0
0023 0.2624 0.5929 0 0.1 0.2 0.3 0.4 0.5 0.6 0 0.1 0.2 0.3 0.4 0.5 0.6
0032 2.6240 13.1101 Plot 5 - Window 5 Plot 6 - Window 6
0.04 0.01
0035 2.0175 1.5300 20 bpm 19 bpm
0.02 0.005
0038 0.6886 0.5697 Ref: 19.66 bpm Ref: 19.30 bpm
0103 0.0079 0.0873 0
0 0.1 0.2 0.3 0.4 0.5 0.6
0
0 0.1 0.2 0.3 0.4 0.5 0.6
0104 0.0049 3.9799 Plot 7 - Window 7 Plot 8 - Window 8
10 -3 10 -3
0105 0.0476 6.4162 4 2
19 bpm 20 bpm
0121 0.0069 11.2100 2 Ref: 19.46 bpm 1
Ref: 19.66 bpm
0122 0.0070 10.6559 0 0
0125 0.3194 3.0763 0 0.1 0.2 0.3 0.4 0.5 0.6 0 0.1 0.2 0.3 0.4 0.5 0.6

0127 0.3776 0.8477 Frequency (Hz) Frequency (Hz)

0128 0.4403 13.8068
0134 0.0062 17.0908 Fig. 5. Frequency domain analysis of the respiration signal obtained from
0142 0.5118 4.2626 our EDR method for recording 0009. Plots 1 to 8 show the frequency spectra
0147 0.0043 13.2784 of each one minute window, along with the reference RR.
0148 0.0051 11.4484
Adult Recordings
0311 0.0263 8.1826
0312 0.6346 13.7058 Plot 1 - Window 1 Plot 2 - Window 2
10 -5 10 -5
0313 0.0110 7.5415 2 2
0322 0.5099 0.5884 0 0
0325 0.0150 14.1616 -2
0 0.2 0.4 0.6 0.8 1
-2
1 1.2 1.4 1.6 1.8 2
Average 0.4150 6.6987 10 -5
Plot 3 - Window 3 10 -5
Plot 4 - Window 4
Amplitude (mV)

2 2
0 0
-2 -2
2 2.2 2.4 2.6 2.8 3 3 3.2 3.4 3.6 3.8 4

10 -5
Plot 5 - Window 5 10 -5
Plot 6 - Window 6
10 -5
Plot 1 - Window 1 10 -5
Plot 2 - Window 2 2 2
1 2
0 0
0 0 -2 -2
4 4.2 4.4 4.6 4.8 5 5 5.2 5.4 5.6 5.8 6
-1 -2
0 0.2 0.4 0.6 0.8 1 1 1.2 1.4 1.6 1.8 2
10 -5
Plot 7 - Window 7 10 -5
Plot 8 - Window 8
Plot 3 - Window 3 Plot 4 - Window 4 2 1
10 -5 10 -5
1 1 0 0
0 0 -2 -1
6 6.2 6.4 6.6 6.8 7 7 7.2 7.4 7.6 7.8 8
Amplitude (mV)

-1 -1
2 2.2 2.4 2.6 2.8 3 3 3.2 3.4 3.6 3.8 4 Time (minutes) Time (minutes)
10 -5
Plot 5 - Window 5 10 -5
Plot 6 -Window 6
2 1

0 0 Fig. 6. The respiration signal obtained from our EDR method for recording
-2 -1
0032. Plots 1 to 8 correspond to the 8 one minute windows obtained from
4 4.2 4.4 4.6 4.8 5 5 5.2 5.4 5.6 5.8 6 our algorithm.
10 -5
Plot 7 - Window 7 10 -5
Plot 8 - Window 8
1 1

0 0

-1 -1
6 6.2 6.4 6.6 6.8 7 7 7.2 7.4 7.6 7.8 8 described in Section II-B [6]. The results obtained from this
Time (minutes) Time (minutes)
approach are shown in Table I. It can be observed that the
MAE obtained from the time domain analysis is high about 6
Fig. 4. The respiration signal obtained from our EDR method for recording
0009. Plots 1 to 8 correspond to the 8 one minute windows obtained from breaths per minute. The main drawback of this method is that
our algorithm. dicrotic notches are also detected as peaks in the respiration
signal. Fig. 8 shows the time domain analysis for the recording
0134. As can be seen a large number of dicrotic notches are
the respiration of recording 0009 (Fig. 4). Our observation detected as peaks, making the RR estimation inaccurate.
becomes more clear from the frequency domain analysis of A major advantage of our proposed EDR method is that
the respiration signal. Fig. 7 shows the frequency spectrum of it shows promising results for both real children and adult
each one minute window for the respiration signal of recording data, without changing any parameters during the experiments.
0032. What stands out in Fig. 7 is that there are several peaks Moreover, as can be seen from Table I our proposed method
close to the dominant, which means that there are overlapping outperforms the proposed time domain analysis, achieving
frequencies, thus we cannot obtain the correct RR for all the an average MAE of 0.4150 breaths per minute compared
windows. to 6.6987 breaths per minute [6]. Furthermore, the practical
In order to compare our RR estimation results, we have advantage of this method is that it can be implemented on-line
implemented the time domain analysis of the respiration signal with an estimation delay equal to one minute, because the

978-1-5386-5602-0/18/$31.00 ©2018 IEEE

 the spectra of the first three IMFs correspond to the frequency
Plot 1 - Window 1 Plot 2 - Window 2 band of the QRS complex. The dominant frequencies in Plots 7
0.2 0.2

0.1 0.1
to 9 (Fig. 9) are about 10-20 Hertz, whereas the dominant fre-
0 0 quencies in Plots 10 and 11 are about 2-5 Hertz, which shows
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
Plot 3 - Window 3 Plot 4 - Window 4 that the last two IMFs correspond to P and T waves, hence
0.1 0.04 they should be discarded before signal reconstruction. Fig. 10
0.05 0.02
shows that the filtered signal, xf (t), can be approximated by
Magnitude

0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 the reconstructed signal, xr (t), because the difference of the
Plot 5 - Window 5 10 -3
Plot 6 - Window 6 two signals (yellow dotted line) is small and the shape of the
0.02 5

0.01
QRS complex is preserved. Therefore, the first three IMFs are
0 0 sufficient to designate the QRS complex. Our assumption was
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
-3 Plot 7 - Window 7 Plot 8 - Window 8
validated for all the recordings of Capnobase dataset, and the
10 10 -3
2 1 first three IMFs were found to be sufficient for reconstructing
1 0.5 the signal, amplifying the QRS complex and reducing low
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 frequency interference.
Frequency (Hz) Frequency (Hz) For all the recordings from Capnobase dataset the results
obtained from our proposed QRS detector are shown in Table
Fig. 7. Frequency domain analysis of the respiration signal obtained from II. Table III shows a comparison of our method’s performance
our EDR method for recording 0032. Plots 1 to 8 show the frequency spectra
of each one minute window.
with other existing detectors. Fig. 11 shows the sequential
steps of the R-peak detection method. The detected R-peaks
are marked by an asterisk ’*’ on the filtered signal, xf (t)
2
10 -5
Plot 1 - Window 1
2
10 -5
Plot 2 - Window 2 (Plot 4). A false negative (FN) occurs when the algorithm fails
0 0 to detect an actual R-peak. A false positive (FP) represents
-2
0 0.2 0.4 0.6 0.8 1
-2
1 1.2 1.4 1.6 1.8 2
a false peak detection. Sensitivity (Se), Positive Predictivity
2
10 -5
Plot 3 - Window 3
2
10 -5
Plot 4 - Window 4 (+P) and Detection Error Rate (DER) were calculated for
0 0 each recordings using the following formulas respectively:
Amplitude (mV)

-2 -2
2 2.2 2.4 2.6 2.8 3 3 3.2 3.4 3.6 3.8 4 TP
10 -5
Plot 5 - Window 5 10 -5
Plot 6 - Window 6 Se(%) = %, (7)
2 2
TP + FN
0 0

-2 -2 TP
4 4.2 4.4 4.6
Plot 7 - Window 7
4.8 5 5 5.2 5.4 5.6
Plot 8 - Window 8
5.8 6
+P(%) = %, (8)
2
10 -5
2
10 -5 TP + FP
0 0
FP + FN
-2
6 6.2 6.4 6.6 6.8 7
-2
7 7.2 7.4 7.6 7.8 8
DER(%) = %, (9)
Time (minutes) Time (minutes)
total number of R peaks
where TP (true positives) is the total number of R-peaks
Fig. 8. Time domain analysis of the respiration signal obtained from our correctly identified.
EDR method for recording 0134. Plots 1 to 8 correspond to the 8 one minute
windows. The detected peaks, which are used in time domain analysis, are As can be seen from Table II and Table III our QRS de-
represented by black triangles. tector shows a better performance for the Capnobase records,
achieving a Se of 100%, a higher +P of 99.80% compared to
99.78% in [7] and 99.70% in [11], as well as a lower DER of
frequency domain analysis of the respiration signal requires 0.24% compared to 0.25% in [7] and 0.31% in [11].
windows of one minute. The reduction of the duration of the During our experiments the following observations were
estimation delay will be part of our future work. found, compared to existing methods. Firstly, an important
observation, which yields high detection error, is that the
C. R-peak Detection absolute amplitude of a Q-peak is larger than the R-peak. This
The proposed QRS detector is based on the assumption that was found to identify the Q-peak as a real R-peak, because
the range of the frequencies of the first three IMFs of the the threshold is applied to the absolute of the reconstructed
EMD correspond to the QRS complex which includes high signal. To overcome this issue the first derivative of the ECG
frequencies in the range 3-40 Hertz [18]. Moreover, P and T signal is computed. The derivative after an R-peak is negative
wave frequencies are about 0.7-10 Hertz [18]. Therefore in because the signal decreases in time. The derivative after a
order to amplify QRS complexes, the IMFs that correspond Q-peak is positive as the signal increases in time. Hence,
to P and T waves should be omitted. The validity of our peaks with a negative derivative were investigated further in
assumption is shown below. Fig. 9 shows a filtered ECG signal the decision stage by applying the refractory period check
and its first five IMFs, obtained after the EMD algorithm. It of 200 milliseconds. Secondly, another significant advantage
is evident from Plots 7 to 11 that as the order of the IMFs of our proposed method is that it can be implemented in
increases, the frequency content decreases. It is shown that real-time with a detection delay equal to the duration of the

978-1-5386-5602-0/18/$31.00 ©2018 IEEE

 Plot 1 - x f(t) Plot 1 - x f(t)
1 1 cand peaks
0 a(t)
-1 0
0 0.5 1 1.5 2 2.5 thresh
Plot 2 - c 1 (t) Plot 7 - F{c 1 (t)} max env
-1
1 4000 0.2 0.4 0.6 0.8 1 1.2 1.4
0 2000
-1 0 Plot 2 - x r (t)
0 0.5 1 1.5 2 2.5 0 5 10 15 20 25 1
Plot 3 - c 2 (t) Plot 8 - F{c 2 (t)}
0.2 2000 0
Amplitude (mV)

0 1000

Amplitude (mV)
-0.2 0
0 0.5 1 1.5 2 2.5 0 5 10 15 20 25 -1
0.2 0.4 0.6 0.8 1 1.2 1.4
Plot 4 - c 3 (t) Plot 9 - F{c 3 (t)}

Magnitude
0.5 2000 Plot 3 - Decision stage
0 1000 1
-0.5 0
0 0.5 1 1.5 2 2.5 0 5 10 15 20 25
0.5
Plot 5 - c 4 (t) Plot 10 - F{c 4 (t)}
0.2 1000
0 500 0
0.2 0.4 0.6 0.8 1 1.2 1.4
-0.2 0
0 0.5 1 1.5 2 2.5 0 5 10 15 20 25 Plot 4 - Identified R peaks
Plot 6 - c 5 (t) Plot 11 - F{c 5 (t)} 1

0.05 1000
0 500 0
-0.05 0
0 0.5 1 1.5 2 2.5 0 5 10 15 20 25
-1
Frequency (Hz) 0.2 0.4 0.6 0.8 1 1.2 1.4
Time (seconds)
Time (seconds)

Fig. 9. The result on the EMD and the spectrum of each IMF. Plot 1
corresponds to the filtered ECG, xf (t). Plots 2 to 6 correspond to the first Fig. 11. Steps of the QRS detector for a part of the record 100 from the MIT-
five IMFs. Plot 7 to 11 correspond to the Fourier transform of each IMF. BIH database. Plot 1, corresponds to the filtered ECG signal, xf (t). Plot 2,
corresponds to the reconstructed signal, xr (t). Plot 3, shows the absolute
sequence, ak (t), (blue line) and its maximum envelope, âk (t), (dotted black
1
line) along with the threshold (dashed black horizontal line) and candidate
x f(t) peaks marked with a red asterisk ‘*’. Plot 4, shows the identified R peaks on
0.8 xf (t) as red asterisk ‘*’.
x r (t)
0.6
x f(t)-xr (t)
0.4
Amplitude (mV)

0.2 estimate respiration rate. Our EDR method shows promising

results for both real children and adult recordings from Cap-
0
nobase dataset, achieving an average mean absolute error of
-0.2
0.4150 breaths per minute, compared to 6.6987 of the time
-0.4 domain analysis [6]. The latter is important because no change
-0.6 of parameters was attempted during our experiments between
-0.8
children and adults. The novelty of our method results from the
assumption that the dominant peak of the frequency spectrum
-1
0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 of the respiration signal can give us an accurate estimation of
Time (seconds) the respiration rate. Additionally, a significant advantage of the
proposed EDR method is that it can be executed in real-time
Fig. 10. Reconstruction of the filtered signal, xf (t), by the summation of
the first three IMFs, xr (t), and their difference, xf (t) − xr (t). with a detection delay of one minute.
Moreover, the present research suggested a new QRS detec-
tor based on the Empirical Mode Decomposition (EMD) using
segment. On-line implementation requires a small alteration an adaptive threshold which relies on the local signal energy.
of our method. The segmentation can be executed at the very This method provides a solution for the detection of small R-
beginning of the pre-processing stage and the sequential steps peaks by establishing a threshold derived from the mean of the
of our algorithm can be implemented for each segment. Hence, Root Mean Square over eight segments. Although our detector
the time between the task request and the availability of the is based on EMD, it uses a gradient-based and refractory
detected R-peaks requires only that we await the part of the period checks to differentiate large Q-peaks and omit false
ECG signal of one segment time duration. For Capnobase R-peaks. Using the Capnobase dataset, the method performed
recordings, for example, a segment of 3 seconds duration effectively with accurate QRS complex detection of 99.80%.
contains an adequate number of QRS complexes. Furthermore, the proposed QRS detector shows comparable
results to other existing detectors using derivative-based [8]
IV. C ONCLUSION and Hilbert methods [12] on real data. In addition, the present
To conclude, a new ECG-Derived (EDR) respiration method method provides a real-time algorithm with a computational
was presented based on the accurate R-peak detection and delay equal to the duration of the segment, whereas proposed
the frequency domain analysis of the respiration in order to methods require the full length ECG [16] [17].

978-1-5386-5602-0/18/$31.00 ©2018 IEEE

 TABLE II [10] N. Maglaveras, N. Stamkopoulos, T. Diamantaras, K. Pappas, and M.
QRS DETECTION PERFORMANCE USING THE C APNOBASE DATASET Strintzis, “ECG pattern recognition and classification using non-linear
transformations and neural networks: A review,” in International journal
Rec Annotated peaks DER (%) Se (%) +P (%) of Medical Informatics, vol. 52(1), pp.191-208, 1998.
9 815 0.00 100 100 [11] B. U. Kohler, C. Hennig, and R. Orglmeister, “The principles of software
15 960 0.00 100 100 QRS detection,” IEEE Engineering in Medicine and Biology Magazine,
16 1012 0.00 100 100 vol. 21(1), pp. 42-57, 2002.
18 1131 0.00 100 100 [12] D. S. Benitez, P. A. Gaydecki, A. Zaidi, and A. P. Fitzpatrick, “A new
23 818 0.00 100 100 QRS detection algorithm based on the Hilbert transform,” in Computers
32 685 0.00 100 100 in Cardiology, vol. 27, pp. 379-382, 2000.
35 900 0.18 100 99.89 [13] D. S. Benitez, P. A. Gaydecki, A. Zaidi, and A. P. Fitzpatrick, “The use
38 956 0.00 100 100 of Hilbert transform in ECG signal analysis,” Computers in Cardiology
103 826 0.00 100 100 and Medicine, vol. 31(5), pp. 399-406, 2001.
104 912 0.00 100 100 [14] N. E. Huang et al., “The Empirical Mode Decomposition and the
105 530 0.37 100 99.62 Hilbert spectrum for nonlinear and non-stationary time series analysis,”
121 579 0.00 100 100 Proceedings of the Royal Society of London A: Mathematical, Physical
122 588 0.00 100 100 and Engineering Sciences, vol. 454, pp. 903-995, 1998.
125 627 0.00 100 100 [15] X. L. Yang, and J. T. Tang, “Hilbert-Huang transform and Wavelet
127 615 0.00 100 100 transform for ECG detection,” 4th International Conference on Wireless
128 541 0.18 100 99.82 Communications, Networking and Mobile Computing, WiCOM’08, pp.
134 578 1.53 100 98.50 1-4, 2008.
142 739 0.00 100 100 [16] M. A. Arafat, and M. K. Hasan, “Automatic detection of ECG wave
147 538 3.58 100 97.52 boundaries using Empirical Mode Decomposition,” in ICASSP, IEEE
148 624 0.00 100 100 International Conference on Acoustics, Speech and Signal Processing,
311 555 0.17 100 99.82 pp. 461-466, 2009.
312 432 0.00 100 100 [17] S. Taouli, and F. Bereksi-Reguig, “Detection of QRS complexes in ECG
313 588 0.00 100 100 signals based on Empirical Mode Decomposition,” in Global Journal of
322 589 0.16 100 99.83 Computer Science and Technology, vol. 11(20), 2011.
325 584 0.00 100 100 [18] J. Malmivuo, and R. Plonsey, Bioelectromagnetism: Principles and
Average 17716 0.24 100 99.80 Applications of Bioelectric and Biomagnetic fields, Oxford University
Press, 1995.
[19] J. G. Proakis, and D. G. Manolakis, “Digital Signal Processing: princi-
ples, algorithms, and applications (3rd edition)”, Prentice hall, 1996.
TABLE III [20] M. Elgendi, M. Jonkman, and F. DeBoer, “Frequency bands effects on
C OMPARISON OF QRS DETECTOR PERFORMANCE WITH OTHER METHODS QRS detection,” in 3rd International Joint Conference on Biomedical
FOR C APNOBASE DATASET Engineering Systems and Technologies (BIOSIGNALS), vol. 5, pp. 428-
431, 2010.
Method DER (%) Se (%) +P (%) [21] W. Karlen, S. Raman, J. M. Ansermino, and G. A. Dumont, “Multi-
Derivative based [7] 0.25 100 99.78 parameter respiratory rate estimation from photoplethysmogram,” IEEE
Hilbert transform [11] 0.31 100 99.70 Transactions on Bio-medical Engineering, vol. 60(7), pp. 1946-1953,
Our method 0.24 100 99.80 2013.

R EFERENCES
[1] J. F Fieselmann, M. S. Hendryx, C. M Helms, and D. S. Wakefield,
“Respiratory rate predicts cardiopulmonary arrest for internal medicine
inpatients,” Journal of General Internal medicine, vol. 8, pp. 354–360,
1993.
[2] D. R. Goldhill, A. F. McNarry, G. Mandersloot, and A. McGinley,
“A physiologically-based early warning score for ward patients: the
association between score and outcome,” Anaesthesia, vol. 6, pp. 547–
553, 2005.
[3] G. B. Moody, R. G. Mark, A. Zoccola, and S. Mantero, “Derivation
of respiratory signals from multi-leas ECGs,” Computers in Cardiology,
vol. 12, pp. 113-116, 1985.
[4] E. Helfenbein, R. Firrozabadi, S. Chien, E. Carlson, and S. Babaeizadeh,
“Development of three methods for extracting respiration from the
surface ECG: a review,” Journal of Electrocardiography, vol. 47, pp.
819–825, 2014.
[5] S. Babaizadeh, S. Zhou, S. D. Stephen, and D. P. White,
“Electrocardiogram-derived respiration in screening of sleep-disordered
breathing,” Journal of Electrocardiology, vol. 44, pp. 700–706, 2011.
[6] S. A. Shah, “Vital sign monitoring and data fusion for paediatric triage,”
PhD Thesis, University of Oxford, 2012.
[7] G. B. Moody, “Spectral analysis of Heart Rate without resampling,” in
Computers in Cardiology, Proceedings, pp. 715-718, 1993.
[8] J. Pan, and W. J. Tompkins, “A real-time QRS detection algorithm,”
IEEE transactions on Biomedical Engineering, vol. 3, pp. 230-236, 1985.
[9] P. S. Hamilton, and W. J. Tompkins, “Quantitative investigation of
QRS detection rules using the MIT/BIH Arrhythmia database,” IEEE
transactions on Biomedical Engineering, vol 12, pp. 1157-1165, 1986.

978-1-5386-5602-0/18/$31.00 ©2018 IEEE