Biomedical Signal Processing and Control 7 (2012) 151–172
Contents lists available at ScienceDirect
Biomedical Signal Processing and Control
journal homepage: www.elsevier.com/locate/bspc
A high-speed C++/MEX solution for long-duration arterial blood pressure
characteristic locations detection
M.R. Homaeinezhad a,b,∗ , A. Ghaffari a,b , M. Aghaee a,b , H.N. Toosi b,c , R. Rahmani d
a
Department of Mechanical Engineering, K.N. Toosi University of Technology, Tehran, Iran
CardioVascular Research Group (CVRG), K.N. Toosi University of Technology, Tehran, Iran
Department of Mechatronic Engineering, K.N. Toosi University of Technology, Tehran, Iran
d
Tehran University of Medical Science (TUMS), Cardiovascular Division (Catheter Laboratory) of Imam Khomeini Hospital, Tehran, Iran
b
c
a r t i c l e
i n f o
Article history:
Received 13 January 2011
Received in revised form 17 May 2011
Accepted 24 May 2011
Available online 30 July 2011
Keywords:
ABP events detection–delineation
Discrete wavelet transform
Pulse pressure
Dicrotic pressure
Systolic blood pressure
Diastolic blood pressure
False alarm probability
Mean
Variance
Skewness
Kurtosis
Neyman–Pearson hypothesis test
a b s t r a c t
The major concentration of this study is to describe the structure of a C++/MEX solution for robust detection and delineation of arterial blood pressure (ABP) signal events. Toward this objective, the original
ABP signal was pre-processed by application of à trous discrete wavelet transform (DWT) to extract several dyadic scales. Then, a sliding window with fixed length was moved on the appropriately selected
scale. In each slid, mean, variance, Skewness and Kurtosis values of the excerpted segment were superimposed to generate a newly defined multiple higher order moments (MHOM) metric to be used as
the detection decision statistic (DS). Then, after application of an adaptive-nonlinear transformation for
making the DS baseline static, the histogram parameters of the enhanced DS were used for regulation of
the ␣-level Neyman–Pearson classifier aimed for false alarm probability (FAP)-bounded delineation of
the ABP events. The proposed method was applied to all 18 subjects of the MIT-BIH Polysomnographic
Database (359,000 beats). The end-systolic and end-diastolic locations of the ABP signal as well as the
dicrotic notch pressure were extracted and values of sensitivity Se = 99.86% and positive predictivity
P+ = 99.95% were obtained for the detection of all ABP events. This paper proves the proposed MHOMbased ABP events detection–delineation algorithm as an improvement because of its merits such as: high
robustness against measurement noises, acceptable detection–delineation accuracy of the ABP events in
the presence of severe heart valvular, arrhythmic dysfunctions within a tolerable computational burden
(processing time) and having no parameters dependency on the acquisition sampling frequency.
© 2011 Elsevier Ltd. All rights reserved.
1. Introduction
The constitutive cells (myocytes) of the heart possess two
important characteristics called nervous excitability and force
feedbacked mechanical tension. If under any circumstances, the
Abbreviations: MHOM, multiple higher order moment; EDES, end-diastolic endsystolic; MAP, mean arterial pressure; SBP, systolic blood pressure; DBP, diastolic
blood pressure; SF, sampling frequency; ssn, smoothing stage number; ASF, adaptive
smoothing filter; DS, decision statistic; FAP, false alarm probability; pdf, probability
density function; ABP, arterial blood pressure; DWT, discrete wavelet transform; FP,
false positive; FN, false negative; PSDB, MIT-BIH Polysomnographic Database; SMF,
smoothing function; FIR, finite-duration impulse response; IIR, infinite-duration
impulse response; LE, location error; CHECK#0, procedure of evaluating obtained
results with the QRS annotations of the MIT-BIH Polysomnographic Database;
CHECK#1, procedure of evaluating obtained results consulting with a control cardiologist; CHECK#2, procedure of evaluating obtained results consulting with a control
cardiologist and also at least with 3 cardiology residents.
∗ Corresponding author at: Department of Mechanical Engineering, K.N. Toosi
University of Technology, Tehran, Iran.
E-mail addresses: mrezahomaei@yahoo.com, mrhomaeinezhad@kntu.ac.ir
(M.R. Homaeinezhad).
1746-8094/$ – see front matter © 2011 Elsevier Ltd. All rights reserved.
doi:10.1016/j.bspc.2011.05.011
electro-mechanical function of a region of myocytes fails, the corresponding abnormal effects are reflected in the heart measurable
quantities such as electrocardiogram (ECG), arterial blood pressure
(ABP) and phonocardiogram (PCG) [1]. ABP depends on a variety
of factors including the pumping action of the heart, blood volume, resistance to flow, and the viscosity (thickness) of the blood.
If any of these factors changes, the corresponding effects will be
reflected in the ABP trend behavior. Physiologically, after each electrical trigger originated from the sinoatrial node, due to sequential
operations of heart main muscle and valves (tricuspid–bicuspid,
pulmonary-aortic), the chambers blood is dynamically excited
causing a substantial pressure gradient in the heart chambers and
the arteries. By inserting catheter tip sensors into atrium, ventricle, or into arteries, the measurement of the continuous ABP
waveforms becomes possible which this approach is categorized
as an invasive method associated with some certain risks. Normally, the characteristic locations of the ABP signal are the systolic
blood pressure (SBP), diastolic blood pressure (DBP), and dicrotic
pressures which by their proper detection, the hemodynamic performance of the heart and vascular system can be analyzed. In Fig. 1,
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
152
150
Systole
Diastole
Pressure<mmHg>
140
Systolic
Pressure
130
120
110
100
Dicrotic Notch
Pulse
Pressure
90
Dicrotic
Pressure
80
Diastolic
70 Pressure
60
596.6
596.8
597
597.2
597.4
Time<sec>
Fig. 1. A parameterized ABP waveform including systole and diastole regions, systolic, diastolic, pulse pressures and dicrotic notch.
a detected–delineated ABP trend including the characteristic markers for diastolic, systolic, dicrotic and pulse pressures is depicted.
Although in the normal subjects, detection and delineation of the
ABP events are rather simple and straightforward; however, the
detection and delineation of the non-impulsive events such as P
or T-waves in the ECG or dicrotic notch pressure in the ABP signal
in the presence of motion artifacts, measurement noises, severe
heart arrhythmia(s) and dysfunctions, turn in a difficult and challengeable task. In the area of heart blood pressure analysis, a major
problem is the robust detection of the incidence (detection) and
onset–offset (delineation or segmentation) locations of each genuine event existing in the cardiac cycles in the presence of heart
sever valvular disorders, diseases, arrhythmia, high level measurement noises as well as circumferential disturbances. In such cases,
the detection and delineation of low-power low-frequency events
might be accompanied by unacceptable false-negative (FN) errors,
while due to existence of noises and disturbing potentials, large
false-positive (FP) errors may appear. Assessments and studies
indicate that in the ABP signal, pressure increment from enddiastolic to end-systolic locations, shows an impulsive behavior
while pressure drop from end-systolic to end-diastolic incidences,
presents a smooth behavior and between this interval, the dicrotic
notch might happen which is due to the closure of the aortic valve
[2] (see Fig. 1). Up to now, several methods have been presented
by authors to detect impulsive regions of the cardiac measured
signals such as ECG or ABP quantities. To meet this end, several
methods associated with pre-processing and analysis of the cardiac system quantities have been introduced by immense number
of researchers. For instance, in the area of ECG signal, Chawla et al.
[3,4], Milanesi et al. [5] and He et al. [6], proposed statisticalbased methodologies for the ECG pre-processing (including noise
and motion artifact removal) and detection of this signal major
events. In some other innovative methods, Chouhan et al. [7]
and Mehta et al. [8–13], elaborated efficient algorithms for the
aim of ECG signal impulsive (QRS complex) and non-impulsive
(P and T-waves) events detection as well as their robust delineation (segmentation) based on support vector machine used as
the discrimination method and a signal information-based measure as the feature of the detection–delineation process. Also,
mathematical models [14], Hilbert transform and the first derivative [15–18], multiple higher order moments [19], second order
derivative [20], wavelet transform and the filter banks [21–24],
combination of signal derivatives and multi-resolution digital
filters with non-parametric detection algorithms [25–28], soft
computing (neuro-fuzzy, genetic algorithm) [29], Hidden Markov
Models (HMM) application [30] were used by authors for fulfill-
ment of genuine events detection–delineation process. Analogous
to aforementioned remarks, the major events of the cardiac system
measured quantities such as ECG, ABP and PCG can be divided into
two impulsive and non-impulsive natures. According to the following reasons, detection–delineation of the characteristic events
are accompanied by complexities, ambiguities and challenges making this problem an important critical research field in the area
of cardiovascular engineering; (a) Measurement noises: measurement noises (high-frequency disturbances) are usually generated
by the electronic probes and medical transducers which are important sources for demolition of the delineation algorithm accuracy.
(b) Holter monitoring: because of the fact that a holter system
records patient’s signal during motion and dynamical movements,
strong noises and motion artifacts can be seen in the holter data
causing the accuracy decrement of a detection–delineation algorithm. (c) Heart disease and arrhythmia: sometimes, occurrence
of severe arrhythmias, mechanical disorders and acute dysfunctions, changes drastically the morphology of events which their
detection and delineation might be associated with certain false
detections. In this study, first, the original ABP signal was preprocessed by application of à trous DWT and then, a fixed sample
size sliding window was moved on the appropriately selected scale.
In each slid, a newly defined MHOM metric was generated to be
used as the detection DS. Then, after application of an appropriate
transformation for making the DS baseline static, the histogram
parameters of the enhanced DS were used for regulation of the ␣level Neyman–Pearson tester aimed for FAP-bounded delineation
of the ABP events. The presented algorithm was applied to more
than 359,000 beats of the PSDB and the end-diastolic end-systolic
(EDES) locations of the ABP signal as well as the dicrotic notch
pressure were extracted. The obtained values of Se = 99.86% and
P+ = 99.95% were obtained for the detection of all ABP events.
2. Materials and methods
2.1. Sampling rate conversion of the original signal into a target
sampling frequency: construction of a unified framework
Almost all parameters of any proposed detection delineation
algorithm are highly dependent on the sampling frequency of
the electronic board system. For example, sampling frequencies
128 Hz, 250 Hz, 360 Hz, 500 Hz, 750 Hz, 1 kHz and 10 kHz can be
seen among ECG and ABP databases [31]. In order to introduce a
unified ABP individual events detection framework which is applicable for all sampling frequencies, the original signal in the core
sampling frequency is mapped to a new trend with target sampling frequency 1 kHz. By this operation, once the parameters of
the algorithm are properly regulated for the target sampling frequency, the algorithm can be implemented to the ABP data sampled
at any rate. It should be noticed that this sampling frequency is
the most profitable discretization rate among other rates. In this
sampling frequency, the EDES regions can be detected by application of dyadic scales 21 –27 . Also, for appropriate delineation
of dicrotic notches, scales 23 –25 can be implemented. As a summary, in Table 1, approximate empirical inference for adoption
of appropriate scales helping successful detection–delineation of
EDES-dicrotic events are presented. For resampling a hypothetical signal x(t) from core SF into a target SF, the signal is treated
as a time-series and the new samples of the resampled trend are
created by application of the truncated sinc function interpolation.
More details left to be seen in [32].
2.2. Discrete wavelet transform using à trous method
Generally, it can be stated that the wavelet transform is a quasiconvolution of the hypothetical signal x(t) and the wavelet function
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
153
Table 1
Appropriate dyadic scales versus sampling rate for successful detection–delineation of ABP signal events. The results are based on empirical assessments.
Database sampling
frequency (Hz)
Suitable scales for EDES
detection
Suitable scales for
events delineation
Sensitivity to window length
and other parameters
128
250
360
500
≥1000
21 –22
21 –22 –23
21 –22 –23
21 –22 –23 –24
21 –· · ·–27
21 –22
21 –22
22 –23 –24
23 –24
23 –24 –25
High-undesirable
High-undesirable
Medium-undesirable
Medium-undesirable
Robust-desirable
(t) with the dilation parameter a and translation parameter b, as
follows:
1
Wax (b) = √
a
+∞
t − b
x(t)
a
−∞
dt,
a>0
(1)
The parameter a can be used to adjust the wideness of the basis
function and therefore the transform can be adjusted in temporal
resolutions. Suppose that the function Yax (b) is obtained based on
a quasi-convolution of signal x(t) and function (t), as follows:
Yax (b) =
+∞
x(t)
−∞
t − b
a
(2)
dt
If the derivative of Yax (b) is calculated relative to b, then
1
∂Yax (b)
=−
a
∂b
+∞
x(t) ′
−∞
t − b
a
k,l (t)
On the other hand, if (t) is the derivative of a smoothing function
(t), i.e. (t) = ′ (t), then
sin(˝/4)
˝/4
H(ejω ) = ejω/2 cos
=
4jejω/2
and therefore,
4
h[n] =
1
8
ω 3
sin
2
ω
{ı[n + 2] + 3 ı[n + 1] + 3ı[n] + ı[n − 1]}
(8)
It should be noted that for the main frequency contents, up to
100 Hz, à trous DWT can be used in different sampling frequencies. Therefore, one of the most prominent advantages of the à trous
algorithm is the approximate independency of its results from sampling frequency. This is because of the main frequency contents of
Scale: 21
G(z2)
H(z)H(z2)H(z4)H(z8)H(z16)
(7)
g[n] = 2{ı[n + 1] − ı[n]}
Scale: 22
SMF 21
H(z)H(z2)
G(z4)
Scale: 23
SMF 22
2
8
4
Scale: 24
H(z)H(z ) H(z )
G(z )
SMF 23
2
16
4
8
H(z)H(z ) H(z ) H(z )
G(z )
Scale: 25
SMF 24
H(z)
(6)
2
x[n]
G(z)
(5)
the transfer functions H(z) and G(z) can be obtained from the following equation:
(4)
Accordingly, it can be concluded that wavelet transform at the scale
a is proportional to the quasi-convolution derivative of the signal x(t) and the smoothing function (t). Therefore, if the wavelet
transform of the signal crosses zero, it will be an indicative of
local extremum(s) existence in the smoothed signal and the absolute maximum value of the wavelet transform in different scales
represents a maximum slope in the filtered signal. Thus, useful
information can be obtained using wavelet transform in different
scales. If the scale factor a and the translation parameter b are con-
k, l ∈ Z +
To implement the à trous DWT algorithm, filters H(z) and G(z)
should be used according to the block diagram represented in Fig. 2
[22,23]. According to this block diagram, each smoothing function
(SMF) is obtained by sequential low-pass filtering (convolving with
H(z2k ) filters), while after high-pass filtering of a SMF (convolving
with G(z2k ) filters), the corresponding DWT at appropriate scale is
generated.
For a prototype wavelet (t) with the following quadratic spline
Fourier transform,
G(ejω )
√ ∂Yax (b)
Wax (b) = − a
∂b
= 2−k/2 (2−k t − l);
(˝) = j ˝
(3)
dt
sidered as a = 2k and b = 2k l, the dyadic wavelet with the following
basis function will be generated [33]:
G(z32)
H(z)H(z2) H(z4) H(z8) H(z16) H(z32)
H(z)H(z2)H(z4)H(z8)H(z16)H(z32) H(z64)
Usable only for sampling
frequency more than
1 kHz
Scale: 26
SMF 25
G(z64)
G(z128)
Scale: 27
SMF 26
Scale: 28
SMF 27
Fig. 2. FIR filter-bank implementation to generate discrete wavelet transform based on à trous algorithm.
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
the ABP signal concentrate on the range less than 100 Hz [23]. After
examination of various databases with different sampling frequencies (range between 128 and 1000 Hz), it has been concluded that in
low sampling frequencies (less than 1 kHz), scales 2 ( = 1, 2, . . . , 5)
are feasible while for sampling frequencies more than 1 kHz, scales
2 ( = 1, 2, . . . , 10) contain profitable information that can be used
for the purpose of wave detection, delineation and classification.
In Fig. 3, several dyadic scales obtained via application of à
trous DWT to an arbitrary ABP record belonging to the MIT-BIH
Polysomnographic Database (PSDB) are illustrated. As it can be
Normalized Signals
a
2.3. Multiple higher order moments (MHOM) metric
The general structure of the MHOM measure consists of firstorder (mean), second-order (variance), third-order (Skewness) and
forth-order (Kurtosis) moments of the excerpted signal in the
analysis window. Each of the aforementioned components has its
b
BP
4
seen, by increment of the dyadic dilation parameter 2 , the highfrequency fluctuations of the transform induced by measurement
noises are diminished.
1
DWT:2
2
0
Normalized Signals
154
-2
500
BP
4
DWT:22
2
0
-2
1000
1500
2000
2500
3000
3500
500
1000
Sample Number
3
DWT:2
2
0
-2
500
1000
1500
2000
2500
3000
3500
BP
DWT:24
2
0
500
1000
1500
2000
2500
f
5
DWT:2
2
0
Normalized Signals
Normalized Signals
3500
3000
3500
Sample Number
BP
4
BP
4
DWT:26
2
0
-2
-2
1000
1500
2000
2500
3000
3500
500
1000
2000
2500
3000
3500
h
BP
4
DWT:27
2
0
-2
Normalized Signals
g
1500
Sample Number
Sample Number
Normalized Signals
3000
-2
e
500
2500
4
Sample Number
500
2000
Sample Number
d
BP
4
Normalized Signals
Normalized Signals
c
1500
BP
4
DWT:28
2
0
-2
1000
1500
2000
2500
Sample Number
3000
3500
500
1000
1500
2000
2500
3000
3500
Sample Number
Fig. 3. Illustration of several DWT dyadic scales 2 obtained from application of à trous filter bank to an arbitrary record of PSDB for (a) = 1, (b) = 2, (c) = 3, (d) = 4, (e)
= 5, (f) = 6, (g) = 7 and (h) = 8.
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
specific characteristics so as by implementing them, the identification of ABP events characteristic parameters such as edge and
incidence location becomes possible. In order to show the procedure of the MHOM measure construction, first, presenting adequate
definitions for two simple concepts namely as mth order moment
m and nth order central moment n is necessary. If X is a random
variable with probability density function PX (x), then the mth order
moment m is obtained via application of the mathematical expectation transform as m = E {Xm }. Also, the nth order central moment
n is obtained by taking expectation from the centralized random variable X about its first order moment, i.e., n = E {(X − 1 )n }.
The structural elements of the MHOM metric are described
below.
2.3.1. Mean value
This quantity is first order moment of random variable X showing the value that uttermost samples are distributed around it.
(9)
Cmean = 1 = E{X}
The mean value shows the mean elevation of the signal from its
baseline in the analysis window. The elevation is a factor by which
detection of the local extremums of the original signal in the analysis window can be achieved.
2.3.2. Variance
This measure is the second order central moment of random
variable X demonstrating dispersion of samples in distribution.
From another point of view, the variance can be interpreted as an
indicator of signal X power.
Cvrc = 2 = E{(X − 1 )2 }
(10)
The variance indicates difference between absolute maximum and
minimum values of an excerpted segment. This difference might
not be detected via mean value because it is possible that the mean
of a segment is a small value while the difference between its maximum and minimum values is large.
2.3.3. Skewness
The Skewness statistic Cskw is obtained from following equation:
Cskw =
3
3/2
2
E{(X − 1 ) }
(E{(X − 1 )2 })
(11)
3
This quantity shows the extent of signal distribution asymmetry and in a signal with symmetric pdf such as Normal, Laplace,
and Cauchy, this quantity is equal to zero. A large value of Skewness indicates a sharp transition from low to high or from high to
low value in the excerpted segment. Consequently, this measure
detects the probable edges of the signal in the analysis window
making the MHOM metric sensitive to the behavior of the signal in the edges. In general, the Skewness statistic Cskw shows
the distance of data distribution from normality and also can be
assumed as a measure indicating the Skewness (asymmetry) of data
[34,35].
2.3.4. Kurtosis
The Kurtosis statistic Ckts is obtained from following ratio:
Ckts =
4
22
=
window, either ascending or descending. Generally, the Kurtosis statistic Ckts indicates the extent of flatness (smoothness or
being impulsive) of samples in the analysis window allowing the
detection of sharp ascending/descending regimes occurred in the
excerpted segment [34,35].
2.3.5. Generation of the MHOM Metric
In order to exploit all discriminant characteristics of each abovementioned statistic, the MHOM metric in the kth slid of the analysis
window is defined by superposition of the measures obtained from
Eqs. (9)–(12) as follows:
MHOM(k) = Cmean (k) + Cvrc (k) + Cskw (k) + Ckts (k)
E{(X − 1 )4 }
(E{(X − 1 )2 })
2
(12)
The Kurtosis quantity acts exactly opposite to the Skewness quantity. A large Kurtosis value of a signal points out a sharp ascending or
descending regime. Consequently, this quantity makes the MHOM
sensitive to the high-slope parts of the signal in the analysis
(13)
As a summary of all aforementioned subsections, the MHOM
metric reacts in edges, maximum (minimum) slope locations and
maximum (minimum) elevation locations making this decision
statistic suitable for detection of edges and extremums of an
excerpted segment in the analysis window. Empirical assessments
show that in the MHOM measure, several geometric parameters such as maximum value to minimum value ratio, area,
extent of smoothness or being impulsive and distribution Skewness degree (asymmetry) can be found. The capacities of the
MHOM-based ABP events detection–delineation algorithm will be
discussed in the next sections. For grasping more details associated with the generation of MHOM metric, the corresponding
C++/MEX framework can be seen in the appendix of the present
manuscript.
2.4. Nonlinear (exponential) amplification and normalization of
the decision statistic (DS) to built the static baseline
In order to prepare the MHOM trend for peak and edge detection, a window with length of (20–25 ms) samples is slid sample by
sample on the MHOM trend. In each slid, a nonlinear amplification
is implemented as follows:
Px = MHOM
3
=
155
n − L n + L
:
2
× exp mean
2
dPx
dt
(MHOM)A [n] = mean(Px )
n = 1, 2, . . . , N
(14)
In Eq. (14), three operations namely as mean(·), exp(·), and
abs(diff(·)) can be seen. If signal Px is stationary, therefore, its abs(diff(·)) will be zero, so that the value of
mean(·) × exp(abs(diff(·))) is constant. In this case, static baseline
of the signal can be detected and chosen as measurement reference
of the signal peaks and edges (this property will be discussed in
Section 2.6.2). On the other hand, if by sliding of the analysis
window, the value of the abs(diff(·)) varies with a rate, its exp(·)
will exponentially vary, so that in the edges of MHOMA signal,
very sharp ascending or descending parts can be seen. The aim of
nonlinear (exponential) amplification of the MHOM signal is to
smoothen misleading fluctuations so that minimum slope-value
locations (edges of events) can be discriminated with greater
accuracy. After exponentially amplification of the MHOM signal, a
local search is conducted to determine positions with minimum
value-slope and correspondingly edges of ABP events are localized.
As an illustrating example, suppose that MHOM signal is a function
with the following definition:
3
MHOM =
Ci exp
j=1
−1
2j2
(x − j )2
(15)
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
156
the ASF problem into a mathematical formulation, assume that the
hypothetical signal x(t) is of length N (t = 1, 2, . . . , N). Then, the kth
step ASF is defined as follows:
8
Normalized Signal
Amplified Signal
Original
6
xk (t) = mean[xk−1 (max(1, t − c) : min(N, t + c))] k = 1, . . . , ssn,
4
t = 0, 1, . . . , N,
2
0
-2
0
0.2
0.4
0.6
0.8
1
Time<sec>
Fig. 4. Performance illustration of the nonlinear amplification of a simulated.
excerpted segment indicated by Eq. (15) to improve detection of the edges.
where in Eq. (15), the parameters set is as
C1 = 1, 1 = 0.20, 1 = 3.5
C2 = 4, 2 = 0.10, 2 = 5.0
C3 = 1, 3 = 0.25, 3 = 6.5
(16)
The artificially obtained trend MHOMA is shown in Fig. 4. According to this figure, it can be seen that the ratio between peak to peak
values of the MHOM are C2 /C1 = 4.0 and C2 /C3 = 4.0 while peak to
peak values of the MHOMA are increased to 5.1 and 5.1, respectively. Therefore, the discrimination between peaks to choose the
best of them as a true event becomes easier and more reliable. The
confusion between genuine events and noise fluctuations can be
occurred when SNR value is low or due to existence of heart valvular
disorders, arrhythmia, and hemodynamic abnormalities. In these
situations, existing events might not be strong enough relative to
noisy oscillations, thus, the detection algorithm might falsely mark
the other peaks as the real events incidence locations.
In the next step, the obtained trend MHOMA is normalized. Normalization of the MHOMA is an algebraic linear transformation
mapping. By application of this mapping to a signal, its first and
second moments are transferred to zero and unity, respectively,
i.e.,
MHOMNA =
[MHOMA − mean(MHOMA )]
std(MHOMA )
std[MHOMNA ] = 1.0
mean[MHOMNA ] = 0;
(17)
where mean(·) and std(·) are mean value and standard deviation
operators, respectively. The aim of normalization is to diminish
subject dependency of the signal MHOMA if so, once applied thresholds and decision criteria are tuned, their performance is preserved
during application of the method to any subject [22,23].
2.5. Adaptive smoothing filtering (ASF)
2.5.1. General structure and programming guidelines
In this section, a new signal smoothing procedure based on
simple mathematics is presented. The aim of the proposed smoothing method is to fade fast fluctuations of a hypothetical signal
x(t) without having any a priori knowledge about the frequency
band or statistical specifications associated with the fluctuations.
After adaptive smoothing of a signal, application of first or second
order derivative-based tests would be possible to find maximum
(minimum) or zero slope locations helping fulfillment of the segmentation process. In signals including fast fluctuations, detection
of locations with maximum (minimum) or zero slopes might turn
in a difficult faulty task with large delineation errors. In order to cast
x0 (t) = x(t)
(18)
where c is the smoothing step, ssn is the smoothing stage number and mean(·) is the mean value operator. For small integer
smoothing step c, the rate of smoothing is moderate while for large
c, for each k increment, substantial changes will appear in the xk (t)
signal. The pseudo code associated with Eq. (18) can be written as
follows:
initialize;
m = length(x);
xs = x; c = 2; ssn = 10;
for k = 1: ssn,
for t = 1: m,
LeftEdge = max(1, t-c);
RightEdge = min(m, t+c);
xs(t) = mean(xs(LeftEdge:RightEdge));
end;
end;
Comprehensive assessments indicate that implementation of
the ASF can improve the edge detection accuracy in the presence of strong impulsive and measurement noises of the recorded
ABP signal approximately 4–20 ms validated by the cardiologist
(CHECK#1). In Fig. 5 some instances of ASF application to the DWT
obtained from an arbitrary PSDB record given ssn = 2, 10, 25 and
50 are depicted. As it is seen, by application of ASF with larger
ssn, the sudden and deep fluctuations of the input signal will be
compensated as much as a smooth and differentiable curve is
obtained.
2.5.2. ASF method evaluation versus higher-order statistical
(HOS) techniques
In this study, to evaluate performance quality of the proposed
denoising algorithm, its performance was compared with the HOSbased filtering methods.
2.5.2.1. Principal components analyzed denoising technique. In this
method, after normalization of each statistical feature, using a linear orthonormal projection method called principal components
analysis (PCA), the dimensionality of the original feature space was
reduced to one. By reduction of the feature space dimensionality
to unity using the PCA algorithm, components having least correlation with the principal and determinant components could be
eliminated from further calculations. Therefore, the accuracy of
the detection–delineation method was improved. In the proposed
PCA-based method, the major point is that by application of the
linear orthonormal projection approach, the asymmetric dispersion shape of the original data in the feature space was mapped to
a symmetric distribution. Consequently, discrimination between
the noise and signal using a threshold-based method can be
fulfilled.
2.5.2.2. Independent components analyzed denoising technique. As
another method, after normalization of each element of feature space, using a blind source separation (BSS) method called
independent components analysis (ICA), the original dependent
feature space was mapped to a new linearly independent space.
By this method, using a sequential method called projection pursuit, the extraction of the most leptokurtotic elements of the
feature space became possible and the detection–delineation decision statistic was generated. Because of the fact that the noise and
high-frequency components decrease the Kurtosis quality of the
employed features, by choosing the most leptokurtotic elements,
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
b
4
Original
Smoothed
3
Normalized Signals
Normalized Signals
a
2
1
0
-1
-2
500
4
Original
Smoothed
3
2
1
0
-1
-2
1000
1500
2000
2500
3000
3500
500
1000
Sample Number
4
Original
Smoothed
3
2
1
0
-1
d
4
3
2000
2500
3000
3500
Original
Smoothed
2
1
0
-1
-2
-2
500
1500
Sample Number
Normalized Signals
Normalized Signals
c
157
1000
1500
2000
2500
3000
500
3500
1000
1500
2000
2500
3000
3500
Sample Number
Sample Number
Fig. 5. Illustration of the smoothing stage number (ssn) effect on weakening the fast fluctuations of a DWT scale obtained by application of à trous algorithm for (a) ssn = 2,
(b) ssn = 10, (c) ssn = 25 and (d) ssn = 50.
the effect of noise and other random phenomena are eliminated
during generation of the detection–delineation decision statistic
and for that reason, better accuracy can be expected.
For fulfillment of this stage, a sliding widow was moved on
the original signal and the excerpted segments in four successive
slides of the window were used as the feature vectors. The successive excerpted segments were put into a matrix in a row-wise
fashion. Then, ASF, PCA and ICA algorithms were applied to this feature space. By application of the PCA technique, the number of the
feature matrix rows was reduced to 1. Also, by application of the
projection pursuit technique, the most leptokurtotic element of the
feature space was calculated. The estimation qualities of the PCA,
ICA and ASF denoising methods for an artificially contaminated
noisy ABP trend were illustrated via Fig. 6.
4
Estimation of Error std
3
2.5
2
1.5
1
0.5
0
-20
-15
-10
-5
0
5
10
15
20
25
2.6. Detection of ABP events via ˛-level Neyman–Pearson based
classifier
2.6.1. Design of false alarm probability (FAP) – bounded classifier
In this study, in order to detect ABP pulses via segmentation of the corrected MHOM metric (Section 2.3), the ␣-level
Neyman–Pearson classifier which is of controlled false alarm
probability (FAP) is designed by implementation of the Gaussian
(normal) stochastic structure.
In order to cast the detection of ABP pulses into a probabilistic framework, assume that the observation set {Z = z|Z ∈ ˝,
˝ = 0 ∪ 1 } consists of two states of nature Hypo. 0 and Hypo.
1, i.e., the samples of observation (in this study observation set
is MHOM metric) Z = z are distributed according to two probability density functions (pdf) p0 (z) and p1 (z). The structure of the
hypotheses Hypo. 0 and Hypo. 1 is shown via Eq. (19), in which
A0 and A1 are the parameters of the hypothesis test problem and
N is a stationary stochastic process with the corresponding model
and parameters.
ASF: ssn=10, c=1
ASF: ssn=10, c=2
ASF: ssn=10, c=3
ASF: ssn=20, c=1
ASF: ssn=20, c=2
ASF: ssn=20, c=3
PCA Denoising
ICA Denoising
3.5
According to this figure, although all examined methods could
successfully follow the original trend, however, the ASF estimation
was accompanied by lower estimation bias.
The denoising results associated with the PCA, ICA and ASF
methods are shown in Fig. 6. As can be seen from this figure, the
estimation error of the ASF algorithm approaches to zero indicating
that the ASF method is an unbiased estimator, while the PCA or ICA
estimation techniques are accompanied by a constant bias by increment of the SNR value. This observation indicates that beside the
higher estimation error versos ASF algorithm, the PCA- or ICA-based
denoising techniques are biased reference signal estimators.
30
SNR Value <dB>
Fig. 6. Comparison between several denoising techniques based on ASF, PCA and
ICA algorithms. In this assessment, ASF with various ssn and window lengths (c) was
applied to artificially contaminated ABP Signal.
Hypo.0 : Z = N + A0
,
Hypo.1 : Z = N + A1
A1 > A0
(19)
In this study, it is supposed that N is a Gaussian zero mean stationary process with standard deviation . Therefore, the pdf of the
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
158
observation under states of nature Hypo. 0 and Hypo. 1 (p0 (z) and
p1 (z), respectively) can be calculated from Eq. (19) as follows:
⎧
⎨ p0 (z) = √ 1
−1
(z − A0 )2
2 2
2
⎩ p1 (z) = √ 1 exp −1 (z − A1 )2
2 2
2
exp
(20)
In order to construct a FAP-controlled Neyman–Pearson classifier,
the likelihood ratio (LR) which is the ratio between Hypo. 1 and
Hypo. 0 pdfs, should be determined as follows:
√1
2
√1
2
p1 (z)
(z) =
=
p0 (z)
exp
exp
= exp
−1
−1
(z − A1 )2
2 2
−1
(z
2 2
2 2
− A0 )2
(z − A1 )2 +
+1
(z − A0 )2
2 2
(21)
By taking log(·) (natural logarithm) from the obtained LR (the
resulted ratio is consequently called log-LR), and by solving the
inequality log(L(z)) > for finding the equivalent inequality for z,
the following result is obtained:
+1
[(z − A0 )2 − (z − A1 )2 ]
2 2
log(L(z)) =
log(L(z))
+1
[(z − A0 )2 − (z − A1 )2 ] > ⇒ z > ′
2 2
> ⇒
(22)
According to Eq. (22), if z > ′ , where ′ is a threshold determined
according to FAP value and the other test parameters, then z will
belong to the alternative class Hypo. 1. The FAP of this test is the
integration of pdf p0 (z) over the 1 (z > ′ ) region as follows:
FAP = P0 (
1) =
+∞
′
dz = 1 − ˚
p0 (z)dz =
′ − A
0
+∞
′
√
1
2
exp
−1
2 2
(z − A0 )2
(23)
where ˚(x) is the unit variance and zero mean Gaussian cumulative probability distribution function with the following integral
definition:
1
˚(x) = √
2
x
2
e(−1)/(2)t dt
(24)
−∞
By equaling the latest term of Eq. (23) to the given FAP value ˛ and
by solving it to find the threshold ′ as a function of A0 , and ˛, the
following result is obtained:
1−˚
′ − A
0
=˛
= A0 + ˚−1 (1 − ˛)
′
(25)
A0 , according to Fig. 7, first, histogram (discrete probability density function) of the MHOM metric is estimated using an existing
simple method (Parzen approximation or kernel method) [34] and
then the maximum value of the obtained histogram is determined
and assigned as A0 . In other words, A0 is the maximum distribution (or equivalently A0 is the value of the MHOM metric baseline).
In order to identify the first corner of the decision statistic histogram, at the first step, easily the absolute value of the histogram
is determined. Due to the fact that the obtained decision statistic is
strictly non-negative quantity (because all constitutive elements of
the MHOM measure are strictly non-negative), thus if the counting
window (Parzen window) lies in the histogram negative region, the
number of samples in this window will be exactly zero. Correspondingly, the value of histogram will be zero (see Fig. 7a). According to
this conclusion, a backward local search starting from the absolute
maximum location of the histogram to find the first zero value to
be assigned as the histogram first corner, is sufficient. As shown
in Fig. 7b, is the distance between A0 value and the first corner of histogram taking place just before A0 . As a matter of fact,
histogram estimation can be fulfilled recursively (on-line) or cumulatively (off-line). The accuracy of the second method is higher than
the on-line methods but in the expense of heavier computational
burden as well as missing real-time implementation of the method.
More details relating to this section are omitted and left to be seen
in [36].
2.7. The MIT-BIH Polysomnographic Database (PSDB) [39]
The MIT-BIH Polysomnographic Database is a collection of
recordings of multiple physiologic signals during sleep. Subjects
were monitored in Boston’s BIH Sleep Laboratory for evaluation of
chronic obstructive sleep apnea syndrome, and to test the effects
of constant positive airway pressure (CPAP), a standard therapeutic
intervention that usually prevents or substantially reduces airway obstruction in these subjects. The database contains over 80 h
worth of four-, six-, and seven-channel polysomnographic recordings, each with an ECG signal annotated beat-by-beat, and EEG
and respiration signals annotated with respect to sleep stages and
apnea.
In Fig. 8, the overall block-diagram of the MHOM-based ABP
events detection–segmentation algorithm is presented. According
to Section 2, the algorithm consists of pre-processing, unifying,
decision statistic preparing, classifier regulating, pulses detecting
and finally segmenting sub-routines.
3. Computer implementation of the proposed ABP events
detection–delineation algorithm and performance
assessment
(26)
where ˛ represents the level (false-alarm probability) of the
binary Neyman–Pearson test and is chosen as 0.005 ≤ ˛ ≤ 0.05 [36].
The structure of Eq. (26) is similar to the conventional threshold = + used frequently as the empirical separation level.
It should be noted that although Eqs. (20)–(26) were derived
based on simplifying assumptions (independent samples, identical
distribution, etc.); however, similar to derivation of Kalman filtering equations, classifier operation depends only on the first and
second-order moments of the signal. Consequently, it can be easily
implemented in actual cases and a higher performance would be
resulted from the algorithm in practical applications [34–36].
2.6.2. Parametric identification of the Neyman–Pearson classifier
parameters via histogram analysis
The parameters A0 and are the test parameters and should be
determined properly to achieve acceptable results. In order to find
3.1. Detection and delineation of the ABP events
3.1.1. Detection of the end-diastolic end-systolic (EDES) regions
According to Eq. (26) of Section 2.6, after determination of the
test parameters A0 and using histogram estimation and given
FAP ˛, the threshold ′ is determined and is used for segmentation of the MHOM metric into EDES (Hypo. 1) and Baseline (Hypo.
0) regions. The decision function ı(n) is then obtained with the
following structure:
ı(n) =
′
EDES
MHOM(n) ≥ ′
Baseline MHOM(n) < ′
(27)
= A0 + ˚−1 (1 − ˛)
To this end, first all samples of the MHOM signal is analyzed
via decision rule presented by Eq. (27) and the edges of each
obtained rectangular pulse is determined using appropriate cal-
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
159
4
4
3.5
3.5
3
3
Normalized MHOM
MHOM Metric Value
a
2.5
2
1.5
2.5
2
1.5
1
1
0.5
A0 Line
0
0.4
Baseline Standard
Deviation (σ )
0.5
0
0.2
0
68.5
69
69.5
Discrete PDF
70
70.5
71
71.5
Time <sec>
b
0.3
Discrete PDF
0.25
0.2
0.15
0.1
σ
0.05
A0
0
0
1
2
3
4
5
6
7
MHOM Metric Value
Fig. 7. Determination of ␣-level Neyman–Pearson classifier parameters via estimation of the MHOM metric histogram (discrete probability density function). The depicted
histogram is obtained from an arbitrary PSDB record. (a) Simultaneous illustration of normalized MHOM metric and associated discrete PDF, (b) parameterization of the
discrete PDF for finding and A0 values.
culations and are assigned as the primary edges of the detected
events, these points are then transferred to the scale 25 and the
index for extremum values are specified between edges. If the minimum is occurred before the maximum, a dominant local minimum
value would exist in the original signal and vice versa. In this way,
the morphology of the detected event can be realized. In Fig. 9,
three EDES regions segmented from the ABP signal by means f the
MHOM-based detection algorithm are shown. As it can be seen, the
obtained metric has a static baseline by which the EDES region can
be segmented.
3.1.2. Detected EDES edges location enhancement
To meet edges detection accuracy enhancement goal, first,
the absolute maximum value of the DS between each detected
onset–offset pair is determined and its location is set as two forward
and backward local searches origin. Also, in each onset–offset interval, the maximum value of DS absolute slope is determined. Next,
by conducting two backward and forward local searches from the
appointed origin, new locations with following satisfied conditions
are identified:
• Point slope equal to or lower than (1/15–1/10) of the determined
maximum slope.
• Do not be far from the initial detected locations more than
1.15 × (offset-peak) interval for offset location and 1.15 × (peakonset) interval for onset location.
• The elevation of new edge position equal to or lower than the first
edge height on DS metric.
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
160
Original ABP
Signal
Digital Band-pass FIR Filter
(0.4 ~ 45) Hz
Discrete Wavelet Transform
Matlab
Scale 2λ
λ=1,…,8
Superposition of Obtained DWT Scales 2λ (λ=5,6)
Detection-Segmentation
Parameters Unification for BP
Acquisition Systems with any
Sampling Frequency
Resampling the Cumulative DWT into Frequency 1 kHz
Application of Sliding Window (WL samples)
C++/MEX
Calculation of Mean, Variance, Skewness and Kurtosis Features
MHOM(k) = Cmean(k)+Cvrc(k)+Cskw(k)+Ckts(k)
Application of Sliding Window (WL/2 samples)
Nonlinear (exponential) Amplification
Histogram
MHOM
Metric
Adaptive Smoothing
Classifier Parameters
Application of Sliding Window (WL/2 Samples)
FAP
α- Level Binary Neyman-Pearson Classifier
EDES Region
Detection
Local Search to Find Two
Absolute Extrema
Matlab
SBP, DBP
Segment Excerption Between Two
End-Systolic End-Diastolic Locations
Median Morphology of the 5
Previous Dicrotic Pulses
ASF
ssn=5, c=1
Local Search to Find All
Absolute Extrema
Dicrotic Notch
Pressure
Fig. 8. General block diagram of the proposed C++/MEX ABP events (SBP, DBP and dicrotic pressures) detection–delineation algorithm via FAP-bounded Neyman–Pearson
segmentation of the MHOM metric.
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
BP
MHOM
100
50
Threshold=3.238
0
2000
2500
3000
b
Normalized Signals
Normalized Signals
a
BP
MHOM
100
50
0
3500
Threshold=3.571
6500
7000
Sample Number
7500
8000
8500
Sample Number
c
BP
MHOM
60
Normalized Signals
161
40
20
0
Threshold=3.2487
1500
2000
2500
3000
3500
4000
Sample Number
Fig. 9. EDES region detection by identification of each MHOM pulse edges and by transferring them to the original ABP to detect the SBP and the DBP incidence locations.
These figures are selected arbitrarily from the PSDB records.
3.5
EDES Region
ABP
MHOM (WL=20)
Normalized Signals
3
MHOM (WL=30)
2.5
MHOM (WL=40)
2
MHOM (WL=50)
1.5
1
0.5
0
400
Threshold=0.3118
500
600
700
800
900
Sample Number
Fig. 10. The effect of sliding window length on the detection of EDES regions. As it
can be seen, the window with larger sample size, causes wider estimation for the
EDES regions, (the unit of WL is ms).
3.1.3. The effect of several window length values on the
delineation performance
According to Fig. 10, the larger window length causes greater
delineation error. On the other hand, if the window length is chosen
small, an oscillatory decision statistic will be obtained. Therefore,
for choosing the window length, a trade-off between delineation
error and high-frequency content of the metric exists.
determined. For specifying proper location of the dicrotic notch
presence, each successive minimum–maximum pair pressure difference, their heights relative to the former SBP and the latter DBP
(surrounding SBP–DBP locations) are compared with the median
of at least five previous beats corresponding parameters Eventually, the best minimum–maximum pair is chosen as the dicrotic
notch. If the algorithm cannot find any local minimum–maximum
pair, it looks for the smoothly varying samples locations and if
was successful in finding this location, the dicrotic edges are determined based on the median of at least five beats occurring before
the present beat. Otherwise, the absence of the dicrotic notch is
announced by the algorithm. In Fig. 11, some examples illustrating
performance of the algorithm in the detection and delineation of
dicrotic notch pressure for some cases such as very weak dicrotic
notches, some random “M” systolic patterns, superfluous notches
resembling actual dicrotic notches, absent dicrotic notches and
severe DBP variations, are shown. In Fig. 12, the long-duration
detected SBP, DBP, calculated MAP and dicrotic pressure versus
beat number for three arbitrary cases of the PSDB are depicted.
In this figure, some examples of blood pressure trends obtained
from PSDB records are shown. As an important application of the
presented algorithm, its probable implementations in ICU monitoring and in the prediction of cardiogenic shock [16,18] can be
mentioned.
3.3. Validation of the proposed ABP pulses detection–delineation
algorithm
3.2. Detection and delineation of the dicrotic pressure
After segmentation of the EDES regions and detection of the
SBP and DBP incidence locations, the segment between each pair
SBP–DBP pair is excerpted from the filtered ABP waveform. Afterwards the ASF with ssn = 5–10 is implemented to the obtained
segment and all local maximum and minimum locations are
In this section, for performance validation of the proposed algorithm, it was applied to ABP signals of all 18 subjects of the MIT-BIH
Polysomnographic Database (PSDB) sampled with 250 Hz rate [39]
(all physionet [31] records can be converted to MAT-files using the
WFDB Software package [38]) and the corresponding end-systolic
and end-diastolic locations of the ABP waveform were extracted.
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
162
a
c
Normalized Trends <Arbitrary Unit>
e
g
i
ABP
DBP
SBP
DP
DN
ABP
b
d
ABP
DBP
SBP
DP
DN
ABP
DBP
DBP
SBP
SBP
DP
DP
DN
DN
ABP
DBP
SBP
DP
DN
ABP
DBP
SBP
DP
DN
ABP
DBP
SBP
DP
DN
f
h
j
ABP
DBP
SBP
DP
DN
ABP
DBP
SBP
DP
DN
ABP
DBP
SBP
DP
DN
Fig. 11. Several abnormal ABP cases parameterized and identified by the proposed MHOM-based algorithm rendering the robustness of the algorithm versus several possible
conditions (ABP: arterial blood pressure, DBP: diastolic blood pressure, SBP: systolic blood pressure, DN: dicrotic notch, DP: dicrotic peak).
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
k
m
o
q
s
ABP
DBP
SBP
DP
DN
ABP
DBP
SBP
DP
DN
ABP
DBP
SBP
DP
DN
ABP
DBP
SBP
DP
DN
ABP
DBP
SBP
DP
DN
l
n
p
r
163
ABP
DBP
SBP
DP
DN
ABP
DBP
SBP
DP
DN
ABP
DBP
SBP
DP
DN
ABP
DBP
SBP
DP
DN
t
Fig. 11. (continued)
ABP
DBP
SBP
DP
DN
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
164
a
DBP
SBP
DP
MAP
150
100
DBP
SBP
DP
MAP
200
ABP <mmHg>
ABP <mmHg>
200
b
150
100
50
50
0
0
0
1000
2000
3000
4000
5000
6000
7000
0
0.5
1
Beat Number
c
2
104
DBP
SBP
DP
MAP
250
200
ABP <mmHg>
1.5
Beat Number
150
100
50
0
0
0.5
1
1.5
2
2.5
3
104
Beat Number
Fig. 12. Detected SBP, DBP, MAP and dicrotic pressure sequences via application of the MHOM-based method. All signals were chosen from the PSDB.
Although the performance of the QRS detection algorithms can
easily be verified using the standard databases such as the MITBIH arrhythmia database [37]; however, in the area of the blood
pressure events detection–delineation, validation of a proposed
algorithm for the segmentation of the ABP pulses, has turned to a
difficult problem due to the lack of an appropriate gold standard as
an accepted universal reference. To overcome the abovementioned
problem, one way is to employ the multi-parameter databases
including annotated ECG trends parallel to the recorded ABP waveform. Because, each triggered QRS complex is accompanied by a
EDES region (see Fig. 13), the QRS annotations can be used to validate indirectly the ABP EDES detector (CHECK#0). But the problem
for validation of the dicrotic notch presence still remains as a hindrance. In this case the only way is to resort to the votes of a
Table 2
Application of the MHOM-based algorithm to the MIT-BIH Polysomnographic Database and obtained results.
MIT-BIH record
slp01a
slp01b
slp02a
slp02b
slp03
slp04
slp14
slp16
slp32
slp37
slp41
slp45
slp48
slp59
slp60
slp61
slp66
slp67x
Total # of beats
7806
11,467
16,145
11,317
24,917
27,029
22,920
27,604
21,718
30,611
25,884
27,686
15,904
16,901
25,017
25,482
15,775
5374
Total # of records
Total # of cardiac cycles
FP
FN
1
6
3
10
19
0
23
3
2
7
15
5
1
0
8
53
30
0
3
5
15
21
6
24
41
21
17
121
80
15
22
7
9
26
56
8
18
359,557
FP + FN
4
11
18
31
25
24
64
24
19
128
95
20
23
7
17
79
86
8
Error (%)
0.051
0.096
0.111
0.274
0.100
0.089
0.279
0.087
0.087
0.418
0.367
0.072
0.145
0.041
0.068
0.310
0.545
0.149
Cumulative Se%
Cumulative P+ %
Se (%)
P+ (%)
99.96
99.96
99.91
99.81
99.98
99.91
99.82
99.92
99.92
99.60
99.69
99.95
99.86
99.96
99.96
99.90
99.65
99.85
99.99
99.95
99.98
99.91
99.92
100
99.90
99.99
99.99
99.98
99.94
99.98
99.99
100
99.97
99.79
99.81
100
99.86
99.95
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
a
6
ECG
ABP
QRS
4
Normalized Signals
Normalized Signals
b
ECG
ABP
QRS
7
165
5
4
3
2
1
2
0
-2
-4
0
-6
-1
5.87
5.872
5.874
5.876
5.878
5.88
5
x 10
Sample Number
2.4518
2.452
c
2.4524
2.4526
6
x 10
ECG
ABP
QRS
10
Normalized Signals
2.4522
Sample Number
5
0
-5
-10
1.37
1.375
1.38
1.385
1.39
1.395
Sample Number
1.4
5
x 10
Fig. 13. Utilization of the QRS annotations as a reference to validate the EDES locations detection accuracy; (a), (b) two subjects with no failure in the blood pressure
acquisition system, (c) a record in which some temporary failures in the blood pressure measurement occur. These intervals are excluded from the algorithm error rate and
accuracy calculations.
cardiologist (CHECK#1) and available medical residents (CHECK#2)
[16,18,19,22,23]. The results of the algorithm applied to blood pressure waveforms are shown in Table 2. In challenging cases, results
were delivered to the cardiologist and accordingly the detection
algorithm was re-validated (CHECK#1). In cases of ABP with very
abnormal morphologies, the results were also checked by some residents (CHECK#2). In order to detect EDES regions, the presented
algorithm was applied to all PSDB signals and the results of SBP and
DBP detection were compared with annotation files by a computer
program (CHECK#0). Due to the fact that a universally accepted
gold standard exists in the area of QRS detection, EDES validation
procedure was significantly simplified to the step CHECK#0. Using
the CHECK#1 and CHECK#2 approaches, it was realized that the
presented algorithm has acceptable performance in the determination of incidence and edges of dicrotic notches. As the final step, in
each case, the trends of SBP, DBP and dicrotic pressure are obtained
and plotted. If the instantaneous pressure is remarkably less than
the mean value, probable false positive (FP) error may exist.
Table 3
Detection–delineation performance evaluation of the MHOM-based ABP characterization algorithm versus artificially contaminated ABP signal.
Test #
SNR (dB)
No. of beats
TP
FP
FN
Error (%)
Se (%)
P+ (%)
Delineation
error SBP (ms)
Delineation
error DP (ms)
Delineation
error DN (ms)
Delineation
error DBP (ms)
1
2
3
4
5
6
7
8
9
10
11
12
13
−35
−30
−25
−20
−15
−10
−6
−2
2
6
10
20
30
1730
1730
1730
1730
1730
1730
1730
1730
1730
1730
1730
1730
1730
876
1473
1728
1730
1730
1730
1730
1730
1730
1730
1730
1730
1730
528
403
1
0
0
0
0
0
0
0
0
0
0
854
257
2
0
0
0
0
0
0
0
0
0
0
79.88
38.15
0.17
0
0
0
0
0
0
0
0
0
0
50.64
85.14
99.88
100
100
100
100
100
100
100
100
100
100
62.39
78.52
99.94
100
100
100
100
100
100
100
100
100
100
Failed
Failed
10.05
7.53
5.47
3.17
2.24
1.74
1.25
0.94
0.82
0.40
0.31
Failed
Failed
7.96
5.09
3.38
2.37
1.41
1.10
0.87
0.68
0.48
0.25
0.09
Failed
Failed
11.24
8.47
6.07
4.09
2.95
2.02
1.49
1.11
0.78
0.49
0.28
Failed
Failed
9.28
6.40
4.63
2.86
2.05
1.45
1.09
0.83
0.67
0.33
0.19
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
166
Table 4
Performance characteristics of some proposed ABP detection–delineation algorithms to be compared with each other (NA: not applicable).
Development
environment
Conventional Hilbert transform
[18]
Modified hilbert transform [18]
Conventional discrete wavelet
transform [23]
DWT-based area curve length
method [23]
This study
Matlab
Matlab
Matlab
Speed
samples (s)
Detection/delineation Noise
robustness
level (dB)
Dependency of
parameters to
sampling frequency
Maximum
delineation
error (ms) (RMS)
Se/P+ (%)
52,710
Yes/no (only SBP)
0
Yes
NA
96.10/91.36
43,830
47,934
Yes/no (only SBP)
Yes/yes (SBP–DBP)
0
−10
Yes
Yes
NA
18.13
99.80/99.86
99.80/99.82
C++/MEX (Matlab)
101,701
Yes/yes (all)
−15
Yes
6.17
99.81/99.81
C++/MEX (Matlab)
189,943
Yes/yes (all)
−30
No
4.36
99.86/99.95
a
DBP
SBP
DP
DN
10
8
6
4
2
0
-20
-10
0
10
20
30
Normalized Signal with SNR -15 dB
Delineation Error <msec>
Detection algorithm
1
0
-1
-2
-3
SNR <dB>
3.4. Robustness assessment of the presented ABP
detection–delineation algorithm versus measurement noise
For robustness assessment of the proposed algorithm, examinations were categorized into two sections as follow.
3.4.1.1. Performance assessment in the presence of abnormal ABP
waveforms
In this case, for rendering performance of the proposed algorithm versus several different cases, the behavior of the method was
scrupulously checked in abnormal situations. For instance, about
20 ABP waveforms and corresponding effects of the algorithm on
them including some explanations and interpretations associated
with each case can be seen in Fig. 11.
3.4.1.2. Performance assessment of the method versus
high-frequency noise
In this stage, the algorithm delineation performance versus
noise level was evaluated. Toward this aim, first the clean signal was
parameterized and identified by the algorithm. Then, the results
were validated by cardiologist to remove probable occurred errors
either in detection or in delineation. Final results were saved in a
m-file script for further robustness assessments. After generation
of the reference annotation file for the selected signal, the clean
trend was embedded in the appropriate artificial noise levels for
obtaining several SNR values. For each SNR value, the algorithm
was applied to the corresponding noisy signal and its performance was evaluated by comparing its outputs with the reference
1.678
1.68
1.682
1.684
1.686
1.688
5
x 10
Sample Number
b
Normalized Signal with SNR -25 dB
On the other hand, if any instantaneous pressure is significantly
greater than the mean value, the false negative error may probably
exist.
1.676
3
ABP
DBP
SBP
DP
DN
2
1
0
-1
-2
-3
1.676
1.678
1.68
1.682
1.684
1.686
1.688
1.69
5
x 10
Sample Number
c
Normalized Signal with SNR -30 dB
Fig. 14. Delineation performance assessment of the proposed algorithm versus SNR
value. The evaluation includes. separate measurement of DBP, SBP, DN and DP localization errors in several SNR values.
ABP
DBP
SBP
DP
DN
2
ABP
DBP
SBP
DP
DN
2
1
0
-1
-2
-3
-4
1.676
1.678
1.68
1.682
1.684
1.686
Sample Number
1.688
1.69
5
x 10
Fig. 15. Graphical illustration of ABP signal with several SNR values, (a) SNR = −15
(dB), (b) −25 (dB) and (c) −30 (dB).
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
annotation file. In Table 3, the quantitative performance evaluation of the algorithm versus several SNR values is depicted. Also, in
Fig. 14 variation of the delineation error (in terms of root mean
square value) versus SNR level is shown. In Fig. 15, a segment
excerpted from original signal and signal with minimum marginal
SNR value is depicted. This figure shows graphically the worst
case that can be parameterized and identified by the proposed
MHOM-based algorithm. After conducting numerous simulations,
it was concluded that for negative SNR values, the performance of
the algorithm gradually decreases and finally for SNR = −30 dB or
lower, the performance of the algorithm is suddenly decreased and
algorithm would rapidly miss its accuracy until complete failure.
Comparison between MHOM-based method and the previous ones
[16,18,23] demonstrates that the method endures in lower SNR
values due to the implementation of appropriate pre-processings
including the proposed static-baseline MHOM metric.
In Table 4, some specifications of the ABP events
detection–delineation are shown. According to this table, the
presented detection–delineation algorithm possesses the best performance characteristics namely as speed of processing, accuracy
and robustness against noise, motion artifacts and circumferential
disturbances.
167
4. Conclusion
In this study, the ABP signal events detection–delineation algorithm was described. To meet this end, first, the original ABP signal
was pre-processed by application of à trous DWT and then, a
fixed sample size sliding window was moved on the appropriately selected scale and in each slid, a newly defined MHOM metric
was generated to be used as the detection DS. After application of
an adaptive-nonlinear transformation for making the DS baseline
static, the histogram parameters of the enhanced DS were used
to regulate the ␣-level Neyman–Pearson tester for FAP-bounded
delineation of the ABP events. The presented algorithm was applied
to more than 359,000 beats of the PSDB and the EDES locations
of the ABP signal as well as the dicrotic notch pressure were
extracted and values of Se = 99.86% and P+ = 99.95% were obtained
for the detection of all ABP events. High robustness against measurement noises, acceptable detection–delineation accuracy of the
ABP events in the presence of severe heart valvular and arrhythmic
dysfunctions within a tolerable computational burden (processing time) and having no parameters dependency to the acquisition
sampling frequency can be mentioned as the merits of the proposed
algorithm.
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
168
Appendix A. The C++/MEX framework for generating the
detection–delineation statistic
//
//
//
//
//
//
//
//
//
High-Speed C++/MEX Framework for Calculation of the Detection-Delineation
Decision Metric. The Approximate Speed is 666,000 Samples/sec.
Author: Mohammad Reza Homaeinezhad.
Co-Authors: Mohammad Aghaee, Hamid Najjaran Toosi.
CardioVascular Research Group- CVRG, Department of Mechanical Engineering,
K.N. Toosi University of Technology.
This Code is an Open-Source Program Released Freely under the Terms of
GNU Public Licence, (http://www.physionet.org/).
Release Date: March, 2010.
#include
#include
#include
#include
#include
#include
"mex.h"
<stdio.h>
<cmath>
<vector>
<windows.h>
<process.h>
using namespace std;
#define DIFF_LEN 5
void PreProcess(void *Data);
struct Params
{
double *py;
int mm;
double *MHOM;
};
double Fs, WL;
void mexFunction( int nlhs, mxArray *plhs[],int nrhs, const mxArray *prhs[])
{
if (nrhs != 10 && nlhs != 1);
// Sampling Frequency <Input # 7>
Fs = mxGetScalar(prhs[6]);
// Sliding Window Length <Input # 8>
WL = mxGetScalar(prhs[7]);
// Channel # 1 DWT <Input # 4>
// Calculation Thread # 1
int mm = mxGetM(prhs[3]);
HANDLE threads[3];
Params param1;
param1.py = mxGetPr(prhs[3]);
param1.mm = mm;
param1.MHOM = new double[mm];
threads[0] =(HANDLE )_beginthread(PreProcess, 0, ¶m1);
// Channel # 2 DWT <Input # 5>
// Calculation Thread # 2
Params param2;
param2.py = mxGetPr(prhs[4]);
param2.mm = mm;
param2.MHOM = new double[mm];
threads[1] =(HANDLE )_beginthread(PreProcess, 0, ¶m2);
// Channel # 3 DWT <Input # 6>
// Calculation Thread # 3
Params param3;
param3.py = mxGetPr(prhs[5]);
param3.mm = mm;
param3.MHOM = new double[mm];
threads[2] =(HANDLE )_beginthread(PreProcess, 0, ¶m3);
// Channel # 1 Original Signal <Input # 1>
double *x001 = mxGetPr(prhs[0]);
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
// Channel # 2 Original Signal <Input # 2>
double *x002 = mxGetPr(prhs[1]);
// Channel # 3 Original Signal <Input # 3>
double *x003 = mxGetPr(prhs[2]);
for(int i = 0; i < mm; i++)
x001[i] += x002[i] + x003[i];
WaitForMultipleObjects(3, threads, true, INFINITE);
for (int i = 0; i < mm; i++)
param1.MHOM[i] = sqrt(param1.MHOM[i] * param1.MHOM[i] + param2.MHOM[i] * param2.MHOM[i]
+ param3.MHOM[i] * param3.MHOM[i]);
mm = mxGetM(prhs[0]);
plhs[0] = mxCreateDoubleMatrix(mm,1,mxREAL);
double *out = mxGetPr(plhs[0]);
for (int i = 0; i < mm; i++)
out[i] = param1.MHOM[i];
delete [] param1.MHOM;
delete [] param2.MHOM;
delete [] param3.MHOM;
}
//=========================================================
inline double MeanValue (const double *V, int L)
{
int j;
double Sum = 0;
for(j = 0 ; j < L ; j++)
Sum += V[j];
return Sum/L;
}
//=========================================================
inline double StdValue (const double *V, int L , double MV)
{
int j;
double Sum = 0;
for(j = 0 ; j < L ; j++)
Sum = Sum + (V[j] - MV) * (V[j] - MV);
return sqrt( Sum / L );
}
//=========================================================
inline double SkewValue (const double *V, int L , double MV)
{
int j;
double Sum = 0;
for(j = 0 ; j < L ; j++)
Sum = Sum + (V[j] - MV) * (V[j] - MV) * (V[j] - MV);
return Sum / L;
}
//=========================================================
inline double KurtoValue (const double *V, int L , double MV)
{
int j;
double Sum = 0;
for(j = 0 ; j < L ; j++)
Sum = Sum + (V[j] - MV) * (V[j] - MV) * (V[j] - MV) * (V[j] - MV);
return Sum / L;
}
//=========================================================
void Normalizer (const double *V, int L , double *rv)
{
double MeanVal = MeanValue(V , L);
double StdVal = StdValue(V , L , MeanVal);
for(int j = 0 ; j < L ; j++)
rv[j] = V[j] / StdVal;
}
169
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
170
//=========================================================
void Differentiator(const double *v, int l, double *rv)
{
for (int i = 0; i < l - 1; i++)
rv[i] = v[i+1] - v[i];
}
//=========================================================
void AbsoluteArray(const double *v, int l, double *rv)
{
for (int i = 0; i < l; i++)
rv[i] = abs(v[i]);
}
//=========================================================
double SumMer(const double *v, int l)
{
double rv = 0;
for (int i = 0; i < l; i++)
rv += v[i];
return rv;
}
//=========================================================
void PreProcess(void *Data)
{
Params params= *(Params *)Data;
int mm = params.mm;
int LeftEdge, RightEdge;
double
double
double
double
double
*Measure1 =
*Measure2 =
*Measure3 =
*Measure4 =
*Sigy = new
new double[mm];
new double[mm];
new double[mm];
new double[mm];
double[WL * 2 + 1];
for (int ii = 0 ; ii < mm ; ii++)
{
LeftEdge = ii - WL;
RightEdge = ii + WL;
if (LeftEdge < 0 )
LeftEdge = 0;
if (RightEdge >= mm )
RightEdge = mm - 1;
int counter = 0;
for (int jj = LeftEdge ; jj <= RightEdge ; jj++)
{
Sigy[jj - LeftEdge] = abs(params.py[jj]);
++counter;
}
counter--;
double MeanVal
Measure1[ii] =
Measure2[ii] =
Measure3[ii] =
Measure4[ii] =
}
Normalizer
Normalizer
Normalizer
Normalizer
= MeanValue (Sigy, counter);
MeanVal;
StdValue (Sigy, counter, MeanVal);
abs(SkewValue (Sigy, counter, MeanVal));
KurtoValue (Sigy, counter, MeanVal);
(Measure1,
(Measure2,
(Measure3,
(Measure4,
mm
mm
mm
mm
,
,
,
,
Measure1);
Measure2);
Measure3);
Measure4);
double *MeasureT = Measure1;
double *MHOMSeg = Sigy;
for (int jj = 0 ; jj < mm ; jj++)
MeasureT[jj] += Measure2[jj] + Measure3[jj] + Measure4[jj];
Normalizer(MeasureT, mm , MeasureT);
for (int ii = 0 ; ii < mm ; ii++)
{
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
171
LeftEdge = ii - WL / 2;
RightEdge = ii + WL /2;
if (LeftEdge < 0 )
LeftEdge = 0;
if (RightEdge >= mm )
RightEdge = mm - 1;
int counter = 0;
for (int jj = LeftEdge ; jj <= RightEdge ; jj++ )
{
MHOMSeg[jj - LeftEdge] = MeasureT[jj];
++counter;
}
counter--;
double MedVar = MeanValue(MHOMSeg, counter);
Differentiator(MHOMSeg, counter, MHOMSeg);
AbsoluteArray(MHOMSeg, counter - 1, MHOMSeg);
double mhom = (MedVar * exp( MeanValue(MHOMSeg, counter - 1) ));
if (mhom > 1)
params.MHOM[ii] = 10;
else
params.MHOM[ii] = mhom;
}
Normalizer(params.MHOM, mm , params.MHOM);
delete
delete
delete
delete
delete
[]Measure1;
[]Measure2;
[]Measure3;
[]Measure4;
[]Sigy;
}
// End of Program
References
[1] P. Libby, R.O. Bonow, D.L. Mann, D.P. Zipes, Braunwald’s Heart Disease: A Textbook of Cardiovascular Medicine , 8th ed., Saunders, 2007.
[2] T.H. Schindler, E.U. Nitzsche, H.R. Schelbert, M. Olschewski, J. Sayre, M. Mix,
I. Brink, X.-L. Zhang, M. Kreissl, N. Magosaki, H. Just, U. Solzbach, Positron
emission tomography-measured abnormal responses of myocardial blood flow
to sympathetic stimulation are associated with the risk of developing cardiovascular events , Journal of American College of Cardiology 45 (2005)
1505–1512.
[3] M.P.S. Chawla, H.K. Verma, V. Kumar, A new statistical PCA–ICA algorithm for
location of R-peaks in ECG , International Journal of Cardiology 129 (2008)
146–148.
[4] M.P.S. Chawla, A comparative analysis of principal component and independent component techniques for electrocardiograms , Neural Computing and
Applications 18 (2009) 539–556.
[5] M. Milanesi, N. Martini, N. Vanello, V. Positano, M.F. Santarelli, L. Landini, Independent component analysis applied to the removal of motion artifacts from
electrocardiographic signals , Medical, Biological, Engineering and Computing
46 (2008) 251–261.
[6] T. He, G. Clifford, L. Tarassenko, Application of independent component analysis in removing artefacts from the electrocardiogram , Neural Computing and
Applications 15 (2006) 105–116.
[7] V.S. Chouhan, S.S. Mehta, N.S. Lingayat, Delineation of QRS-complex, P and Twave in 12-lead ECG , IJCSNS International Journal of Computer Science and
Network Security 8 (4) (2008) 185–190.
[8] S.S. Mehta, N.S. Lingayat, Application of support vector machine for the detection of P- and T-waves in 12-lead electrocardiogram , Computer Methods and
Programs in Biomedicine 93 (2009) 46–60.
[9] S.S. Mehta, N.S. Lingayat, Identification of QRS complexes in 12-lead
electrocardiogram , Expert Systems with Applications 36 (2009) 820–
828.
[10] S.S. Mehta, N.S. Lingayat, SVM based QRS detection in electrocardiogram using
signal entropy , IETE Journal of Research 54 (3) (2008) 231–240.
[11] S.S. Mehta, N.S. Lingayat, Support vector machine for cardiac beat detection in
single lead electrocardiogram , IAENG International Journal of Applied Mathematics 36 (2) (2007).
[12] S.S. Mehta, D.A. Shete, N.S. Lingayat, V.S. Chouhan, K-means algorithm for the
detection and delineation of QRS-complexes in electrocardiogram , IRBM 31
(1) (2010) 48–54.
[13] S.S. Mehta, N.S. Lingayat, Development of entropy based algorithm for cardiac beat detection in 12-lead electrocardiogram , Signal Processing 87 (2007)
3190–3201.
[14] O. Sayadi, M.B. Shamsollahi, A model-based Bayesian framework for ECG beat
segmentation , Physiological Measurement 30 (2009) 335–352.
[15] M. Arzeno Natalia, Z.-D. Deng, C.-S. Poon, Analysis of first-derivative based
QRS detection algorithms , IEEE Transactions on Biomedical Engineering 55
(February (2)) (2008) 478–484.
[16] A. Ghaffari, M.R. Homaeinezhad, M. Akraminia, M. Atarod, A methodology for
prediction of acute hypotensive episodes in ICU via a risk scoring model including analysis of ST-segment variations , Cardiovascular Engineering 10 (1) (2010)
12–29.
[17] D. Benitez, P.A. Gaydecki, A. Zaidi, A.P. Fitzpatrick, The use of the Hilbert transform in ECG signal analysis , Computers in Biology and Medicine 31 (2001)
399–406.
[18] A. Ghaffari, M.R. Homaeinezhad, M. Atarod, M. Akraminia, Parallel processing of ECG and blood pressure waveforms for detection of acute hypotensive
episodes: a simulation study using a risk scoring model , Computer Methods in
Biomechanics and Biomedical Engineering 13 (2) (2009) 197–213.
[19] A. Ghaffari, M.R. Homaeinezhad, M. Khazraee, M. Daevaeiha, Segmentation
of Holter ECG Waves via Analysis of a Discrete Wavelet-Derived Multiple
Skewness–Kurtosis Based Metric, Annals of Biomedical Engineering, vol. 38(4),
Springer Publishing, 2010, pp. 1497–1510.
[20] M. Mitra, S. Mitra, A software based approach for detection of QRS vector of
ECG signal , IFMBE Proceedings 15 (2007) 348–351.
[21] J.P. Martinez, R. Almeida, S. Olmos, A.P. Rocha, P. Laguna, A wavelet-based ECG
delineator: evaluation on standard databases , IEEE Transactions on Biomedical
Engineering 51 (4) (2004) 570–581.
[22] A. Ghaffari, M.R. Homaeinezhad, M. Akraminia, M. Atarod, M. Daevaeiha, A
robust wavelet-based multi-lead electrocardiogram delineation algorithm ,
Medical Engineering & Physics 31 (10) (2009) 1219–1227.
[23] A. Ghaffari, M.R. Homaeinezhad, M. Akraminia, M. Daevaeiha, Finding events
of electrocardiogram and arterial blood pressure signals via discrete wavelet
transform with modified scales , Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine 224 (1) (2010)
27–42.
[24] Anonelli, W. Ohley, R. Khamlach, Dicrotic notch detection using wavelet transform analysis , in: Proc. 16th Annu. Int. Conf. IEEE Engineering in Medicine and
Biology Society, vol. 2, 1994, pp. 1216–1217.
[25] M. Aboy, J. McNames, T. Thong, D. Tsunami, M.S. Ellenby, B. Goldstein, An
automatic beat detection algorithm for pressure signals , IEEE Transactions on
Biomedical Engineering 52 (10) (2005) 1662–1670.
[26] P.F. Kinias, M. Norusis, A real time pressure algorithm , Computers in Biology
and Medicine 11 (1981) 211.
[27] M. Aboy, J. McNames, B. Goldstein, Automatic detection algorithm of intracranial pressure waveform components , in: Proc. 23th Int. Conf. IEEE Engineering
in Medicine and Biology Society, vol. 3, 2001, pp. 2231–2234.
[28] M. Aboy, C. Crespo, J. McNames, B. Goldstein, Automatic detection algorithm
for physiologic pressure signal components , in: Proc. 24th Int. Conf. IEEE Engi-
172
[29]
[30]
[31]
[32]
[33]
[34]
M.R. Homaeinezhad et al. / Biomedical Signal Processing and Control 7 (2012) 151–172
neering in Medicine and Biology Society and Biomedical Engineering Society,
vol. 1, 2002, pp. 196–197.
N. Kannathal, C.M. Lim, U.R. Acharya, P.K. Sadasivan, Cardiac state diagnosis
using adaptive neuro fuzzy technique , Medical Engineering & Physics 28 (2006)
809–815.
G. de Lannoy, B. Frenay, M. Verleysen, J. Delbeke, Supervised ECG delineation
using the wavelet transform and hidden Markov models , The Proceedings of
IFMBE 22 (2008) 22–25.
Physionet.org, http://www.physionet.org/physiobank/database/.
N.J. Fliege, Multirate Digital Signal Processing: Multirate Systems – Filter
Banks–Wavelets , John Wiley & Sons, 1994.
S. Mallet, A Wavelet Tour of Signal Processing , Academic Press, 1999.
D.C. Montgomery, G.C. Runger, Applied Statistics and Probability for Engineers
, third ed., John Wiley & Sons, 2003.
[35] S.M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory ,
Prentice-Hall Inc., 1979.
[36] H. Vincent Poor, An Introduction to Signal Detection and Estimation , in: second
ed., Springer-Verlag, 1994, Chapters 1–3.
[37] G.B. Moody, R.G. Mark, The MIT-BIH arrhythmia database on CD-Rom and
Software for it , in: The Proceeding of Computers in Cardiology, 1990, pp.
185–188.
[38] G.B. Moody, WFDB Applications Guide , tenth ed., Harvard-MIT Division of Health Sciences and Technology, 2006, http://www.physionet.org/
physiotools/wag/.
[39] A.L. Goldberger, L.A.N. Amaral, L. Glass, J.M. Hausdorff, P.Ch. Ivanov, R.G. Mark,
J.E. Mietus, G.B. Moody, C.K. Peng, H.E. Stanley, PhysioBank, PhysioToolkit, and
PhysioNet: components of a new research resource for complex physiologic
signals , Circulation 101 (2000) e215–e220.
Reza Rahmani