Electrocardiograms QRS Peak and Heart Rate detection Using DWT in Matlab
December 17, 2021
- Machine Learning
The QRS combines three deflections (Q, R, and S) seen on a typical ECG. It corresponds to the depolarization of the right and left ventricles of the human heart and contraction of the large ventricular muscles.
In numerical and functional analysis, a discrete wavelet transform (DWT) is any wavelet transform in which the wavelets are discretely sampled.
The discrete wavelet transform has many engineering, mathematics, and computer science applications. Most notably, it is used for signal coding to represent a discrete signal in a more redundant form, often as preconditioning for data compression.
sym4 wavelet resembles the QRS, suitable for QRS detection. Therefore, this process can help to diagnose various heart diseases. This tutorial will look at how to obtain the peak and rate of detection of these ECGs using the ECG database. From this method, we can get the heart rate.
To follow through this tutorial, you’ll need:
Table of content
- Table of content
- The QRS complex
- ECG database on PhysioNet
- Use of symlet4 wavelet for ecg signal analysis
- Proposed DWT based QRS detection
- Matlab code to get QRS peak and heart rate from ecg signals
The QRS complex
As we said earlier, it is a combination of three deflections (Q, R, and S) seen on a typical ecg signal:
P is the first deflection.
Q is the first negative deflection to the baseline.
R is the highest positive deflection to the baseline.
S the second negative deflection to the baseline.
The amplitude of a normal QRS is 5 to 30mm, and the duration is 0.06 to 0.12 seconds. The width, amplitude, and shape of the QRS complex help diagnose ventricular arrhythmias, conduction abnormalities, ventricular hypertrophy, myocardial infarction, electrolyte rearrangements, and other diseases state.
Note that the QRS complex does not always have all three QRS. It can have various shapes, as shown below:
ECG database on PhysioNet
For this tutorial, we use signals from MIT-BIH arrhythmia, and the ECG-ID database downloaded from PhysioNet. Each ecg signal on PhysioNet has the following three files:
- *.atr: Reference Annotation.
- *.dat: Datafile(signal).
- *.hea: Header file.
However Matlab cannot read such files, we therefore have to convert our ecg to a
.mat file. To do that, we use the PhysioNet ATM. The interface of the ATM bank is as shown below:
You can select your database in the input by clicking on the dropdown arrow to choose your database. Note that all the PhysioNet ecg databases are available here:
You can select the record, signals, annotation, output length, time format, and data format since they all have options. When you reach the toolbox section, you also select your options, when you choose
plot waveforms, you will have the plots of the waveform as shown below:
Since we need to read it in Matlab, we export it. To do that, we select the
export signal as .mat and then download it on the toolbox:
Since we only need the signal, we download the
Use of symlet4 wavelet for ecg signal analysis
sym4 wavelet is similar to the QRS complex. That is why it’s preferred for QRS detection. To make this clear, look at the image of extracted QRS complex and dilated
sym4 wavelet and make a comparison:
As you can see, the QRS complex of the ecg is quite similar to the
sym4 wavelet in shape. That’s why
sym4 wavelets are always preferred for the ecg signal analysis.
Proposed DWT based QRS detection
Below are the essential ecg signals, and if we look at them carefully, we can locate the labeled areas with a particular frequency contribution.
f1: Represents the high-frequency noise and has some frequency f1.
f2: It is the QRS that has the frequency contribution of f2.
f3: Slow varying content of the ecg and have a frequency contribution f3.
The relationship between these three frequencies will be
Our objective to preserve all the R-peaks and eliminate all the other frequencies. To make it clear, we say that we want to eliminate
f3 but preserve
This is known as bandpass filtering. You achieve it with the help of the wavelet transform. Wavelet transform groups signals of the same frequency bands. Therefore, You can implement bandpass filtering by eliminating some frequency bands.
This bandpass filtering can be achieved by eliminating wavelet coefficients of some lower scale (high frequencies) and higher scales (lower frequency) of ecg signals. For this purpose, an undecimated wavelet transform is used to get wavelet coefficients.
What is an undecimated wavelet transform?
Well, in a normal
mra wavelet, transform signals are downsampled to two after every decomposition level, by which its size reduces at every decomposition level.
Therefore, in an undecimated wavelet, the signal length remains the same. A 4-level decomposition of an ecg signal using
sym4 is shown in the figure below:
The first plot is the ecg signal. The
d's are the detailed coefficients at every level of the ecg signal.
a4 is the approximate coefficients at level 4.
We will obtain the bandpass filtering by removing the co-efficient
a4 since it will not be considered—similarly, we eliminate
The reason why we don’t consider it is because it is an approximated coefficient. It carries all the low-frequency details.
d2 are not considered because they contain details of the signal’s high frequency.
d4 are considered to reconstruct or achieve the signal the bandpass is filtering.
We get the following signals by considering only
d4 and taking the inverse wavelet transform.
With the help of a standard peak detection algorithm, we can locate these R-peaks. Also, you find the number of total R-peaks for a given time interval to find the heart rate.
For example, suppose we have a 10-second ecg signal and the total number of R-peaks have some values, then we can find the number of R-peaks in a minute, representing the beat per minute which is the heart rate.
Matlab code to get QRS peak and heart rate from ecg signals
The first step is to input our signal. The user should input the signal, so Matlab should ask for it. For Matlab to allow the user to select the signals from the folder, we use the
uigetfile function. This function takes into consideration the path and the file name:
%program to get QRS peaks and heart rate from ecg signal [filename, pathname]=uigetfile('*.*', 'select your ecg signal'); filewithpath = strcat(pathname, filename);
Next, we need the sampling frequency of the signal. These sampling frequencies are defined in the database. We use the
input function since the user defines the sampling frequency. This function reads the user’s input. After this, the data is loaded using the
Fs = input('Enter sampling rate:'); ecg = load(filename); %reading ecg signal
Afterwards, we normalize the amplitude. It is done by dividing the ecg value by the gain. This gain value is given in the database too. We also get the length of the signal using function
length, this function takes in the signal as the input.
This length helps in determining the time taken by the signal:
ecgsig = (ecg.val)./200; %Normalize gain t = 1:length(ecgsig); %No. of samples tx = t./Fs; %Getting time vector
Next, we need to compute the undecimated wavelet transform of the 4-level using
sym4. This ensures that the length of the signal remains the same. To compute this, we use the
modwt function. This function takes
ecg signal and the
sym4 level 4.
wt = modwt(ecgsig,4, 'sym4'); %4-level undecimated DWT using sym4 wtrec = zeros(size(wt));
As explained earlier, our wavelet transform has 5 rows, that is, $a_n, d_4, d_3, d_2$, and $d_1$. We don’t need the approximated and high-frequency coefficients $d_1$ and $d_2$. So we extract the $d_3$ and $d_4$, which are the 3rd and 4th rows.
We use the rows to extract as the argument for the undecimated DWT:
wtrec(3:4, :)= wt(3:4,:); %extracting only d3 and d4 coefficients
By performing the inverse discrete wavelet transform (IDWT), we will have a signal that has only the r-peaks well preserved. Inverse DWT returns the signal to the original form after performing the DWT.
In Matlab, we use
imodwt function to do the IDWT with the arguments as the signals with the extracted parts
y=imodwt(wtrec, 'sym4'); %IDWT with only d3 and d4. y=abs(y).^2; %magnitude square
We then find the average of the signal. It will be used as the threshold when finding the signal’s peak. Finding the average is done by using the
avg = mean(y); %getting average of y^2 as threshold %Finding peaks [Rpeaks, locs] = findpeaks(y,t, 'MinPeakHeight', 8*avg, 'MinPeakDistance', 50);
find peaks is a variable available in the signal processing toolbox to find the peaks. We have the minimum peak distance as
50 to avoid false detection if the peaks are close to each other. It could happen due to improper filtering.
locs give the location of the R-peaks.
Let’s find the location of the R-peaks in consideration of the length of the signal. It represents the number of beats. This number of beats is then converted to beats per minute:
nohb = length(locs); %No. of beats timelimit = length(ecgsig)/Fs; %getting the time function of the signal hbpermin = (nohb*60)/timelimit; %Getting Beat per minute. disp(strcat('HeartRate= ', num2str(hbpermin))) %displaying the heartrate
Plot the normal ecg signal against time so we could be able to see the difference:
%displaying ecg signal and detected R-peaks subplot(211) plot(tx, ecgsig); xlim([0, timelimit]) grid on xlabel('seconds') ylabel('ECG signal')
Additionally plot the filtered signal along with the detected peaks:
subplot(212) plot(t,y) grid on xlim([0, length(ecgsig)]); hold on plot(locs, Rpeaks, 'ro') xlabel('samples') title(strcat('R peaks found and heartrate: ', num2str(hbpermin))) % displayes the heartrate.
When we execute our program, we will have the following as our output:
ECG, QRS, and heart rate detection are easier using the discrete wavelet transform. Matlab is the best software for wavelet analysis. As we have seen, these transforms are already done and exist in in-built form. Therefore, it makes it easy to perform operations.
Also, Matlab has other built-in functions that help analyze the signal apart from having the transform in in-built form. These functions include
length for getting the length. Furthermore, the database for the ecg signal is compatible with Matlab since it gives options to download Matlab files.
I hope you found this tutorial helpful.
Peer Review Contributions by: Monica Masae