function QRSOnsets = QRSDetection(ECG)
%QRSDETECTION
% This function performs Pan-Tompkins algorithm for QRS complex detection.
% Sampling frequency is assumed to be 200 Hz.
% ECG = electrocardiogram
% QRSOnsets = the locations (indices) where the QRS-comlexes begin
T = 1/200;
% **Bandpass filter**
% Lowpass filter
% The parameters of the lowpass filter are detemined.
b_lowpass = (1/32)*[1 0 0 0 0 0 -2 0 0 0 0 0 1];
a_lowpass = [1 -2 1];
% The signal is filtered with the lowpass filter
ECG_filtered1 = filter(b_lowpass,a_lowpass,ECG);
% **Highpass filter**
% The parameters of the highpass filter are determined.
b_highpass = zeros(1,33);
b_highpass(1) = -1;
b_highpass(17) = 32;
b_highpass(18) = -32;
b_highpass(33) = 1;
b_highpass = (1/32)*b_highpass;
a_highpass = [1 -1];
% The lowpassfiltered signal is filtered with the highpass filter.
ECG_filtered2 = filter(b_highpass,a_highpass,ECG_filtered1);
% **Differentiator**
% The parameters of the derivative filter are determined.
b_diff = (1/8)*[2 1 0 -1 -2];
a_diff = 1;
% The highpass filtered signal is filtered with the derivative filter.
ECG_filtered3 = filter(b_diff,a_diff,ECG_filtered2);
% **Squaring operation**
% The derivative filtered signal is squared.
ECG_filtered4 = ECG_filtered3.^2;
% **Moving-window integrator**
% The squared signal is filtered with a moving-window integrator.
% The size of the window is set.
N = 30;
% The parameters of the moving-average integrator are determined.
b_integ = (1/N)*ones(1,N);
a_integ = 1;
% The squared signal is filtered with the moving-average integrator.
ECG_filtered5 = filter(b_integ,a_integ,ECG_filtered4);
% **Thresholding**
% The amplitude threshold for QRS detection is set.
amplitudeThreshold = 0.01;
% The maximum heart rate is set.
maximumHeartRate = 200; %Bpm
% The blanking interval is calculated based on the maximum heart rate.
blankingInterval = (1/(maximumHeartRate/60))*(1/T);
%detectPeaks-function is applied to the moving-average integrated signal
% for detection of peaks.
[QRSOnsets QRSOffsets] = detectPeaks(ECG_filtered5,blankingInterval,amplitudeThreshold);
function [peakOnsets peakOffsets] = detectPeaks(data,blankingInterval,amplitudeThreshold)
% This function determines the number of peaks in the given data based on
% the amplitudeThreshold and blankingInterval. The results are stored in
% the vectors peakOnsets and peakOffsets (equal length) indicating the
% positions of onsets and offsets of peaks, respectively.
% The vectors are initialized.
peakOnsets = [];
peakOffsets = [];
for i = 1:length(data)-1
% Every position where the threshold is crossed is stored as peak onset
% or offset depending on the direction.
if data(i) <= amplitudeThreshold && data(i+1) > amplitudeThreshold
peakOnsets = [peakOnsets i];
end
if data(i) > amplitudeThreshold && data(i+1) <= amplitudeThreshold
peakOffsets = [peakOffsets i];
end
end
% The last peak onset is left out if the data ends during a peak.
if data(length(data)) > amplitudeThreshold
peakOnsets = peakOnsets(1:length(peakOnsets)-1);
end
% The first peak offset is left out if the data starts during a peak.
if data(1) > amplitudeThreshold
peakOffsets = peakOffsets(2:length(peakOnsets));
end
trueQRS = 1;
for i = 2:length(peakOnsets)
% The peaks occuring too early (during blanking interval) after
% previous peak and peaks with too short duration are left out.
if peakOnsets(i)-peakOnsets(i-1) > blankingInterval || peakOffsets(i)-peakOnsets(i) < 10
trueQRS = [trueQRS i];
end
end
peakOnsets = peakOnsets(trueQRS);
peakOffsets = peakOffsets(trueQRS);
% Time delay correction.
peakOnsets = peakOnsets-30;
peakOffsets = peakOffsets-30;