17
Project: Detection of the QRS Complex
Adapted from
• “Laboratory Exercises and Projects on Biomedical Signal Analysis” by Professor Rangaraj Rangayyan, University of Calgary.
• Rangayyan, Rangaraj M. Biomedical signal analysis. Vol. 33. John Wiley & Sons, 2015.
• Pan, Jiapu, and Willis J. Tompkins. “A real-time QRS detection algo- rithm.” IEEE transactions on biomedical engineering 3 (1985): 230-236.
• Kim, Jinkwon, and Hangsik Shin. “Simple and robust realtime qrs de- tection algorithm based on spatiotemporal characteristic of the qrs com- plex.” PloS one 11.3 (2016): e0150144.
17.1 Submission
• Due 29 Nov 2020.
• Group-based, maximum of 4 students per group (1-4).
• Each group must submit: one report via Canvas: Matlab code via Matlab Grader Platform (by only one representative member).
17.2 Objectives
• To detect QRS complexes in ECG signals using the Pan–Tompkins algo-
rithm.
• To measure parameters of the ECG for analysis of the heart rate and cardiac rhythm.
17.3 Background
17.3.1 The electrocardiogram (ECG)
The ECG is the electrical manifestation of the contractile activity of the heart and can be recorded fairly easily with surface electrodes on the limbs or chest. The rhythm of the heart in terms of beats per minute (bpm) may be easily estimated by counting the readily identifiable waves. More important
Signals and Systems in Biomedical Engineering Semester A 2020/21
is the fact that the ECG waveshape is altered by cardiovascular diseases and abnormalities such as myocardial ischemia and infarction, ventricular hypertrophy, and conduction problems. You are referred to ecg_note.pdf f2o8r backINgTRroOuDUnCdTIOreNaTdOiBnIOgM.EDICAL SIGNALS
P Chirarattananon
Figure 17.1: Summary of the parts and waves of a cardiac cycle as seen in an ECG signal.
R
QS
80 ms
AV node fires SA node fires
P ECG
T
160 ms
T: ventricular repolarization
80 ms
P: atrial contraction
ventricular relaxation
or diastole
up to next beat
80 ms
100 ms
QRS: ventricular depolarization, contraction,
or systole
Figure 1.26
full compensatory pause.)
PQ segment: isoelectric (AV delay)
ST segment:
isoelectric,
plateau of action potential of ventricular myocyte
Summary of the parts and waves of a cardiac cycle as seen in an ECG signal.
A cardiac cycle is represented in a period of the repetitive ECG signal
as the series of waves labeled P, QRS, and T. If we view the cardiac cycle
as a series of events, we have the following epochs in an ECG waveform as sbhiogwemnininy,Fwigh.er1e7e.1v:ery second pulse from the SA node is replaced by a PVC with a
P wave: The P wave is the epoch related to the event of atrial contraction.
The QRS waveshape is affected by conduction disorders; for example, bundle-
The atria do not possess any specialized conduction system
branch block causes a widened and possibly jagged QRS. Figure 1.28 shows the ECG
signal of a paastitehnet wveitnhtricglhetsbduon;dales-bsruacnhc,hcbolnoctrka.ctOiobnseorvfeththeeawtriidaelr-mthuasnc-lneosrmal
QRS complex, which also displays a waveshape that is significantly different from
takes place in a slow squeezing manner, with the excitation
the normal QRS waves. Ventricular hypertrophy (enlargement) could also cause a
stimulus being propagated by the muscle cells themselves. For
wider-than-normal QRS.
this reason, the P wave is a slow waveform, with a duration of
The ST segment, which is normally isoelectric (flat and in line with the PQ seg-
about 80 ms. The P wave amplitude is much smaller (about 0.1 − ment) may be elevated or depressed due to myocardial ischemia (reduced blood sup-
0.2 mV ) than that of the QRS because the atria are smaller than ply to a part of the heart muscles due to a block in the coronary arteries) or due to
myocardial tinhfearvcetinotnri(cdleas.d myocardial tissue incapable of contraction due to total lack of blood supply). Many other diseases cause specific changes in the ECG wave-
PQ segment: The PQ segment, of about 80 ms duration, is a “nonevent”; shape: the ECG is an important signal that is useful in heart-rate (rhythm) monitoring
and the diagnosis of cardiovascular diseases.
however, it is important in recognizing the baseline as the inter-
val is usually isoelectric.
QRS wave: The specialized system of Purkinje fibers stimulate contrac-
tion of ventricular muscles in a rapid sequence from the apex upwards. The almost simultaneous contraction of the entire ven- tricular musculature results in a sharp and tall QRS complex of
City University of Hong Kong
142
Signals and Systems in Biomedical Engineering Semester A 2020/21
about 1 mV amplitude and 80 − 100 ms duration. The event of ventricular contraction is represented by the QRS epoch.
ST segment: The normally flat (isoelectric) ST segment is related to the plateau in the action potential of the left ventricular muscle cells. The duration of the plateau in the action potential is about 200 ms; the ST segment duration is usually about 100 − 120 ms.
As in the case of the PQ segment, the ST segment may also be called a nonevent.
P Chirarattananon
T wave:
The T wave appears in a normal ECG signal as a discrete wave separated from the QRS by an isoelectric ST segment. However, it relates to the last phase of the action potential of ventricular muscle cells, when the potential returns from the plateau of the depolarized state to the resting potential through the process of repolarization. The T wave is elusive, being low in amplitude (0.1 − 0.3 mV ) and being a slow wave extending over 120 − 160 ms. It is almost absent in many ECG recordings.
17.3.2 Motivation
An electrocardiograph (ECG) is a graphical representation of the electrical
8 Laboratory Exercise: Detection of the QRS and Parameteri-
activity of a heart over time, and it is the most basic examination method
• To measure parameters of the ECG for analysis of the heart rate and cardiac rhythm. detecting the QRS complex. With the results of the QRS complex detection,
other fiducial points in the ECG, such as P, Q, S, or T, can also be detected to provide further information.
zation of the ECG
that can be used for cases of heart disease. However, it may be necessary to record and analyze long-term ECGs since the symptoms of some types of
8à Objectives
heart disease are generally intermittent. Thus, research related to automatic ECG processing algorithms has been actively conducted over the last several
• To detect QRS complexes in ECG signals using the Pan–Tompkins algorithm. decades. The first step of an automatic ECG processing algorithm involves
8à2 Background
The algorithm developed by Pan and Tompkins identifies QRS complexes based on analysis of
17.3.3 QRS Complex Detection
the slope, amplitude, and width of the QRS. The various stages of the algorithm are shown in
Figure 2.
Figure 2: Block diagram of the Pan–Tompkins algorithm for the detection of QRS complexes in
The algorithm developed by Pan and Tompkins identifies QRS complexes Figure 17.2: Block di- ECG signals.
based on analysis of the slope, amplitude, and width of the QRS. The vari- agram of the Pan–
Bandpass filter
Differentiator
The bandpass filter, formed using lowpass and highpass filters, reduces noise in the ECG
ous stages of the algorithm are shown in Fig. 17.2. Tompkins algorithm
signal. Noise such as muscle noise, 60 Hz interference, and baseline drift are removed by bandpass The bandpass filter, formed using lowpass and highpass filters, reduces for the detection of filtering. The signal is then passed through a differentiator to provide a large response at the high
noise in the ECG signal. Noise such as muscle noise, 60 Hz interference, QRS complexes in ECG slopes that distinguish QRS complexes from low-frequency ECG components such as the P and
and baseline drift are removed by bandpass filtering. The signal is then signals.
T waves.
passed through a differentiator to provide a large response at the high
The next operation is the squaring operation, which emphasizes the higher values expected
slopes that distinguish QRS complexes from low-frequency ECG compo-
due to QRS complexes and suppresses smaller values related to the P and T waves, as well as noise
nents such as the P and T waves.
in the output of the preceding stage. The squared signal is then passed through a moving-window
The next operation is the squaring operation, which emphasizes the
integrator of window length N = 30 samples (for the sampling frequency of fs = 200 Hz). The highervaelxupeescteexdpercetseudltdisueatsoinQglReSsmcoomotphlepxeaskarnedlastuedpptoretshsesQsmRSalcleormvpalle-xforeachECGcycle. The output of the moving-window integrator may be used to detect QRS complexes, measure RR
intervals, and determine the duration of the QRS complex (see Figure 3).
City University of Hong Kong 143 See Section 4.3.2 of the textbook for details.
8à3 Exercises to be worked by hand
Squaring operation
Moving-window integrator
Signals and Systems in Biomedical Engineering Semester A 2020/21
ues related to the P and T waves, as well as noise in the output of the pre- ceding stage. The squared signal is then passed through a moving-window integrator of window length N = 30 samples (for the sampling frequency of
fs = 200 Hz). The expected result is a single smooth peak related to the QRS complex for each ECG cycle. The output of the moving-window integrator may be used to detect QRS complexes, measure RR intervals (time between consecutive R peaks), and determine the duration of the QRS complex (see Fig. 17.3).
P Chirarattananon
R
QS
W – QS
T
Figure 17.3: The re- lationship of a QRS complex to the moving- window integrator output. (Top) Schematic ECG signal. (Bottom) Output of the moving- window integrator. QS: QRS complex width. W: width of the integrator window, given as N/fs s.
P
QS
QS
ure 3: The relationship of a QRS complex to the moving-window integrator output. Schematic ECG signal. (b) Output of the moving-window integrator. QS: QRS complex
2. Another digithaaltficoltmerpoisestphecPifiane–dToinmptkerinmssalogofrithsmim. Tphuelsfoellroewsipnognfislesaasnd func- tions are provided through the Matlab Grader platform.
– ECG data: ECG3.dat, EChG4(.nda)t=, ECδG(5n.)da−t δ(n − 1). 2
Instruction: Detection of the QRS wave
th. W: width of the integrator window, given as N/fs s.
1. Develop a Matlab program to perform the various filtering procedures
– Matlab function: pan_tompkins_thresholding_fn.m
Derive its transfer function as well as the magnitude and phase parts of its frequency re-
2. The provided template contains the script that loads the raw ECG data sponse. Explain the nature and effects of the filter.
recorded at 200 Hz. The data is then normalized and named necg, ready
to be used as input to the Pan–Tompkins algorithm in Fig. 17.1. You are
3. The two filters given above are used in series to filter a signal. Derive and plot the impulse
to add your own code to this file after the provided lines.
response of the combined filter.
3. The first step in the Pan–Tompkins algorithm is the band-pass filtering,
Derive the transfer function as well as the magnitude and phase parts of the frequency
which is formed using lowpass and highpass filters. Hence, we will start
response of the combined filter. Explain the nature and effects of the combined filter.
by applying a low-pass filter to the signal necg. The transfer function of
the low-pass filter is
Detection of the QRS wave
1 1−z−62
Hlp (z) = 32 (1 − z−1)2 , (17.1)
1. Develop a Matlab program to perform the various filtering procedures that compose the
where z = ejω. With the sampling rate being 200 Hz, the filter has a
Pan–Tompkins algorithm. Use the filter command for each step; see Section 4.3.2 of the
low cutoff frequency of 11 Hz. The filter provides an attenuation greater
textbook. Process each of the three ECG signals that you acquired in Lab 1 using your
program. PCirteypUanrieveprsloittysoof fHtohneg Ksiogngal being processed over the full duration (or a relevant 144 portion of about 10 − 15 seconds) and over two cardiac cycles, before and after each stage
of the Pan–Tompkins method.
You may also use ECG3àdat, ECG4àdat, ECG5àdat, ECG6àdat, and ECGSàm.
g )
d
4
Signals and Systems in Biomedical Engineering Semester A 2020/21
P Chirarattananon
4.
5.
6.
7.
than 35 dB at 60 Hz and effectively suppresses power-line interference, if present.
After the low-pass filter, the high-pass filter is applie. Together, they form a band-pass filter. A suitable transfer function of the high-pass filter for the Pan-Tompkin algoirthm is
Hhp(z)=z−16− 1 1−z−32. 32 (1−z−1)
The highpass filter has a cutoff frequency of 5 Hz.
The derivative operator used by Pan and Tompkins is specified as
y[n] = 1 (2x[n]+x[n−1]−x[n−3]−2x[n−4]), 8
(17.2)
(17.3)
8.
which approximately is an operator d/dt for up to 30 Hz. The derivative procedure suppresses the low-frequency components of the P and T waves, and provides a large gain to the high-frequency components arising from the high slopes of the QRS complex.
The squaring operation is a nonlinear operation that mathematically described as y [n] = x2 [n]. This makes the result positive and emphasizes large differences resulting from QRS complexes; the small differences arising from P and T waves are suppressed.
Next, the Pan–Tompkins algorithm performs smoothing of the output of the preceding operations through a moving-window integration filter as
1 N−1 y[n]=N ∑x[n−i]
i=0
The choice of the window width N is to be made with the following considerations: too large a value will result in the outputs due to the QRS and T waves being merged, whereas too small a value could yield multiple peaks for a single QRS. A window width of N = 30 samples was found to be suitable for the sample rate of 200 Hz.
Fig. 17.3 illustates the the output of the integrator and its relationship to the QRS width. The shape of the signal in the region of interest is approximately trapezoidal. The width of the trapezoid is W + QS.
To finally obtain the QRS width, you may use the provided thresholding function:
[yout, n_i, n_f] = pan_tompkins_thresholding_fn(x_in,thredhold,min_len) Using the output from the moving average filter as the input of this x_in function, you must also provide the threshold value for recognizing the trapezoids (see Fig. 17.4 for an example). The value of min_len could be
safely set as N or 30. The output of the function are
-yout The envelop of the recognized trapezoids, as shown in Fig. 17.4.
The width of each envelop is expected to be QS+W-QS+QS=W+QS as
shown in Fig. 17.3.
-n_i Vector of time indices indicating the beginning of the trapezoid
(blue circles in Fig. 17.4)
n_f Vector of time indices indicating the end of the trapezoid (red circles
in Fig. 17.4).
City University of Hong Kong
145
Signals and Systems in Biomedical Engineering Semester A 2020/21
P Chirarattananon
8 10-4
7 6 5 4 3 2 1 0
-1
1850 1900
1950 2000 2050 index n
2100 2150
x_in, signal after the moving average filter.
yout, signal after applying The thresholding function
threshold value n_i n_f
Example Matlab Code
Suppose you have an array of x [n] data in Matlab variable xn. To directly
compute yn, y [n], that corresponds to the application of the transfer func-
% discrete−time signal sampling
Figure 17.4: The re- lationship of a QRS complex to the moving- window integrator output. (Top) Schematic ECG signal. (Bottom) Output of the moving- window integrator. QS: QRS complex width. W: widthoftheintegrator window, given as N/fs s.
tion H ejω = 1 1−z−1 to x [n] in the frequency domain as Y ejω =
jω 2 1−3z−2
H e X e . The following Matlab code could be used.
jω
fs = 200; frequency
z = tf(’z’,1/fs); transfer functions .
% creates a ’z’ object for
H = (1/2)*(1−z^−1)/(1−3*z^−2); % desired transfer function .
yn = lsim(H,xn) ; % apply the transfer function ( in the frequency domain) to signal xn in the time domain .
Project Report
You project report should be limited to the maximum of 10 pages (exclud- ing the appendix if needed). The project report should include standard components (course title, lab members, etc). In addition, you project report should include
• Name of the group representative who submits the final, complete code on the Matlab Grader platform (please ignore the grading function of the Grader platform when you submit your code, multiple submissions before the submission deadline are allowed ). Clearly comment your code.
• The results from each step, including comments and observations for the output from each step of the algorithm.
• The plots of output signals (with proper axis labels) after each step shown in Fig. 17.2.
City University of Hong Kong
146
ECG after MA and thresholding
Signals and Systems in Biomedical Engineering Semester A 2020/21
• The difference equations corresponding to the transfer functions of Eqs. (17.1) and (17.1).
• Transfer function of the difference equation (17.3).
• Describe and Include steps in your program to compute automatically the following parameters for each ECG signal (ECG3.dat, ECG4.dat, ECG5.dat) from the final results:
– Total number of beats detected in each signal and the heart rate in beats per minute.
– Average RR interval (time separation between consecutive R peaks) and standard deviation of RR intervals of each signal (in ms).
– Average QRS width computed over all the beats in each signal (in ms). These numbers should be tabulated.
P Chirarattananon
City University of Hong Kong
147