function peakloc = PeakDetection2(data,fs,varargin)
%
% peaks = PeakDetection2(x,fs,wlen,fp1,fp2,th,flag),
% R-peak detector based on Pan-Tompkins method.
%
% inputs:
% x: vector of input data
% fs: sampling rate
% wlen: moving average window length (default = 150ms)
% fp1: lower cut-off frequency (default = 10Hz)
% fp2: upper cut-off frequency (default = 33.3Hz)
% th: detection threshold (default = 0.2)
% flag: search for positive (flag=1) or negative (flag=0) peaks. By default
% the maximum absolute value of the signal, determines the peak sign.
%
% output:
% peaks: vector of R-peak impulse train
%
%
% Open Source ECG Toolbox, version 2.0, March 2008
% Released under the GNU General Public License
% Copyright (C) 2008 Reza Sameni
% Sharif University of Technology, Tehran, Iran -- GIPSA-Lab, INPG, Grenoble, France
% reza.sameni@gmail.com
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 2 of the License, or (at your
% option) any later version.
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
% Public License for more details.
%
% See also ECGtask_QRS_detection
%
% Adapted by Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% to ecg-kit toolbox for Matlab.
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate : 21/4/2015
% Copyright 2008-2015
%
if(nargin>2 && ~isempty(varargin{1})),
winlen = varargin{1};
else
winlen = .150; % 150ms
end
if(nargin>3 && ~isempty(varargin{2})),
fp1 = varargin{2};
else
fp1 = 10;
end
if(nargin>4 && ~isempty(varargin{3})),
fp2 = varargin{3};
else
fp2 = 33.3;
end
if(nargin>5 && ~isempty(varargin{4})),
thr = varargin{4};
else
thr = 0.2;
end
if(nargin>6 && ~isempty(varargin{5})),
flag = varargin{5};
else
flag = abs(max(data))>abs(min(data));
end
N = length(data);
data = data(:);
L1 = round(fs/fp2); % first zero of the LP filter is placed at f = 33.3Hz;
L2 = round(fs/fp1); % first zero of the HP filter is placed at f = 3Hz;
fprintf(1, '.');
% x0 = data - Median(data',N,round(fs*winlen/3));
q = ceil(fs/150);
N_padded = (ceil(N / q) * q);
data = [data; zeros(N_padded-N,1)];
data = data - resample(MedianFilt(resample(data,1,q), round(fs*winlen/3/q)), q, 1);
data = data(1:N);
% x0 = data;
fprintf(1, '.');
% LP filter
x = filter([1 zeros(1,L1-1) -1],[L1 -L1],data);
x = filter([1 zeros(1,L1-1) -1],[L1 -L1],x);
x = [x(L1:end);zeros(L1-1,1) + x(end)]; % lag compensation
fprintf(1, '.');
% HP filter
y = filter([L2-1 -L2 zeros(1,L2-2) 1],[L2 -L2],x);
% differentiation
z = diff([y(1) ; y]);
% squaring
w = z.^2;
fprintf(1, '.');
% moving average
L3 = round(fs*winlen);
v = filter([1 zeros(1,L3-1) -1],[L3 -L3],w);
v = [v(round(L3/2):end);zeros(round(L3/2)-1,1) + v(end)]; % lag compensation
vmax = max(v);
p = v > (thr*vmax);
fprintf(1, '.');
% edge detection
rising = find(diff([0 ; p])==1); % rising edges
falling = find(diff([p ; 0])==-1); % falling edges
if( length(rising) == length(falling)-1 ),
rising = [1 ; rising];
elseif( length(rising) == length(falling)+1 ),
falling = [falling ; N];
end
fprintf(1, '.');
peakloc = zeros(length(rising),1);
width = zeros(length(rising),1);
Nrising = max(1, length(rising)/20);
j = 0;
if(flag)
for i=1:length(rising)
if( j > Nrising)
fprintf(1, '.');
j = 0;
else
j = j + 1;
end
[val mx] = max( data(rising(i):falling(i)) );
peakloc(i) = mx - 1 + rising(i);
width(i) = falling(i) - rising(i);
end
else
for i=1:length(rising)
if( j > Nrising)
fprintf(1, '.');
j = 0;
else
j = j + 1;
end
[val mn] = min( data(rising(i):falling(i)) );
peakloc(i) = mn - 1 + rising(i);
width(i) = falling(i) - rising(i);
end
end
fprintf(1, '\n');
% peaks = zeros(1,N);
% peaks(peakloc) = 1;