function [qrs,jpoints]=wqrsm(data,fs,PWfreq,TmDEF,jflag)
% The WQRS detector rewrited from wqrs.c
% rewrite by Qiao Li, March 6, 2017
% qiaolibme@gmail.com
%
% input:
% data: ECG data, the wqrs.c used the data in raw adus, if the input
% data is in physical units, when the function found the peak-peak
% value < 10, it will multiply the value by WFDB_DEFGAIN (200)
% fs: sampling frequency (default: 125)
% PWfreq: power line (mains) frequency, in Hz (default: 60)
% TmDEF: minimum threshold value (default: 100)
% jflag: annotate J-points (ends of QRS complexes) (default: 0)
% output:
% qrs: QRS fiducial mark in samples
% jpoints:J-points annotation, if jflag==1
%
% Matlab code based on:
% file: wqrs.c Wei Zong 23 October 1998
% Last revised: 9 April 2010 (by G. Moody)
% -----------------------------------------------------------------------------
% wqrs: Single-lead QRS detector based on length transform
% Copyright (C) 1998-2010 Wei Zong
%
% 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.
%
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 59 Temple
% Place - Suite 330, Boston, MA 02111-1307, USA.
%
% You may contact the author by e-mail (wzong@mit.edu) or postal mail
% (MIT Room E25-505, Cambridge, MA 02139, USA). For updates to this software,
% please visit PhysioNet (http://www.physionet.org/).
% ------------------------------------------------------------------------------
%
% This program analyzes an ECG signal, detecting QRS onsets and J-points, using
% a nonlinearly-scaled ECG curve length feature. This version has been optimized
% for ECGs sampled at 125 Hz, but it can analyze ECGs sampled at any frequency
% using on-the-fly resampling provided by the WFDB library.
%
% `wqrs' can process records containing any number of signals, but it uses only
% one signal for QRS detection (signal 0 by default; this can be changed using
% the `-s' option, see below). 'wqrs' has been optimized for adult human ECGs.
% For other ECGs, it may be necessary to experiment with the input sampling
% frequency and the time constants.
% ----------------------------------------------------------------------
% This matlab version has updates to deal with a variable sampling
% frequency but has not been exhanstively tested.
% ----- Add notes here:
%
%
if nargin<5
jflag = 0;
end
if nargin<4
TmDEF = 100;
end
if nargin<3
PWfreq = 60;
end
if nargin<2
fs = 125;
end
clear global wqrsm_Yn wqrsm_Yn1 wqrsm_Yn2 wqrsm_lt_tt wqrsm_lbuf wqrsm_ebuf wqrsm_aet wqrsm_et wqrsm_LPn wqrsm_LP2n wqrsm_lfsc wqrsm_data wqrsm_BUFLN wqrsm_LTwindow
global wqrsm_Yn wqrsm_Yn1 wqrsm_Yn2;
global wqrsm_lt_tt wqrsm_lbuf wqrsm_ebuf;
global wqrsm_aet wqrsm_et;
global wqrsm_LPn wqrsm_LP2n wqrsm_lfsc;
global wqrsm_data;
global wqrsm_BUFLN;
global wqrsm_LTwindow;
wqrsm_BUFLN = 16384; % /* must be a power of 2, see ltsamp() */
EYE_CLS = 0.25; % /* eye-closing period is set to 0.25 sec (250 ms) */
MaxQRSw = 0.13; % /* maximum QRS width (130ms) */
NDP = 2.5; % /* adjust threshold if no QRS found in NDP seconds */
PWFreqDEF = PWfreq; % /* power line (mains) frequency, in Hz (default) */
% TmDEF = 100; % /* minimum threshold value (default) */
WFDB_DEFGAIN = 200.0; % /* default value for gain (adu/physical unit) */
timer_d=0;
gain=WFDB_DEFGAIN;
wqrsm_lfsc = fix(1.25*gain*gain/fs); % /* length function scale constant */
PWFreq = PWFreqDEF; % /* power line (mains) frequency, in Hz */
% test data is physical units (mV) or raw adus units
datatest=data(1:fix(length(data)/fs)*fs);
if length(datatest)>fs
datatest=reshape(datatest,fs,[]);
end
test_ap=median(max(datatest)-min(datatest));
if test_ap<10 % peak-peak < 10 mV, may physical units
data=data*gain;
end
wqrsm_data = data;
wqrsm_lbuf = zeros(wqrsm_BUFLN,1);
wqrsm_ebuf = zeros(wqrsm_BUFLN,1);
wqrsm_ebuf(1:end)=fix(sqrt(wqrsm_lfsc));
wqrsm_lt_tt = 0;
wqrsm_aet = 0;
wqrsm_Yn = 0;
wqrsm_Yn1 = 0;
wqrsm_Yn2 = 0;
qrs = [];
jpoints = [];
Tm = fix(TmDEF / 5.0);
% spm = 60 * fs;
% next_minute = from + spm;
wqrsm_LPn = fix(fs/PWFreq); % /* The LP filter will have a notch at the power line (mains) frequency */
if (wqrsm_LPn > 8)
wqrsm_LPn = 8; % /* avoid filtering too agressively */
end
wqrsm_LP2n = 2 * wqrsm_LPn;
EyeClosing = fix(fs * EYE_CLS); % /* set eye-closing period */
ExpectPeriod = fix(fs * NDP); % /* maximum expected RR interval */
wqrsm_LTwindow = fix(fs * MaxQRSw); % /* length transform window size */
% for i=1:2000
% ltdata(i)=ltsamp(i);
% end
% /* Average the first 8 seconds of the length-transformed samples
% to determine the initial thresholds Ta and T0. The number of samples
% in the average is limited to half of the ltsamp buffer if the sampling
% frequency exceeds about 2 KHz. */
t1 = fs*8;
if t1> fix(wqrsm_BUFLN*0.9)
t1=wqrsm_BUFLN/2;
end
T0=0;
for t=1:t1
T0 = T0 + ltsamp(t);
end
T0 = T0 / t1;
Ta = 3 * T0;
% /* Main loop */
t=1;
learning=1;
while t<length(data)
if (learning==1)
if (t > t1)
learning = 0;
T1 = T0;
t = 1; % /* start over */
else
T1 = 2*T0;
end
end
% /* Compare a length-transformed sample against T1. */
if (ltsamp(t) > T1) % /* found a possible QRS near t */
timer_d = 0; % /* used for counting the time after previous QRS */
maxd = ltsamp(t);
mind = maxd;
for tt = t+1:t + fix(EyeClosing/2)
if (ltsamp(tt) > maxd)
maxd = ltsamp(tt);
end
end
for tt = t-1:-1:t - fix(EyeClosing/2)
if (ltsamp(tt) < mind)
mind = ltsamp(tt);
end
end
if (maxd > mind+10) % /* There is a QRS near tt */
% /* Find the QRS onset (PQ junction) */
onset = fix(maxd/100) + 2;
tpq = t - 5;
for tt = t:-1:t - fix(EyeClosing/2)
if (ltsamp(tt) - ltsamp(tt-1) < onset && ltsamp(tt-1) - ltsamp(tt-2) < onset && ltsamp(tt-2) - ltsamp(tt-3) < onset && ltsamp(tt-3) - ltsamp(tt-4) < onset)
tpq = tt - wqrsm_LP2n; % /* account for phase shift */
break;
end
end
if (learning~=1)
% /* Check that we haven't reached the end of the record. */
if tpq>length(data)
break;
end
% /* Record an annotation at the QRS onset */
qrs = [qrs tpq];
% J-points processing
if (jflag)
tj = t+5;
for tt=t:t + fix(EyeClosing/2)
if ltsamp(tt) > maxd - fix(maxd/10)
tj = tt;
break;
end
end
if tj>length(data)
break;
end
% Record an annotation at the J-point
jpoints = [jpoints tj];
end
end
% /* Adjust thresholds */
Ta = Ta + (maxd - Ta)/10;
T1 = Ta / 3;
% /* Lock out further detections during the eye-closing period */
t = t + EyeClosing;
end
elseif (learning~=1)
% /* Once past the learning period, decrease threshold if no QRS
% was detected recently. */
timer_d = timer_d + 1;
if timer_d > ExpectPeriod && Ta > Tm
Ta = Ta -1;
T1 = Ta / 3;
end
end
t = t+1;
end
end
function lt_data = ltsamp(t)
% /* ltsamp() returns a sample of the length transform of the input at time t.
% Since this program analyzes only one signal, ltsamp() does not have an
% input argument for specifying a signal number; rather, it always filters
% and returns samples from the signal designated by the global variable
% 'sig'. The caller must never "rewind" by more than wqrsm_BUFLN samples (the
% length of ltsamp()'s buffers). */
global wqrsm_Yn wqrsm_Yn1 wqrsm_Yn2;
global wqrsm_lt_tt wqrsm_lbuf wqrsm_ebuf;
global wqrsm_aet wqrsm_et;
global wqrsm_LPn wqrsm_LP2n wqrsm_lfsc;
global wqrsm_data;
global wqrsm_BUFLN;
global wqrsm_LTwindow;
while (t > wqrsm_lt_tt)
wqrsm_Yn2 = wqrsm_Yn1;
wqrsm_Yn1 = wqrsm_Yn;
v0=wqrsm_data(1);
v1=wqrsm_data(1);
v2=wqrsm_data(1);
if wqrsm_lt_tt>0 && wqrsm_lt_tt<=length(wqrsm_data)
v0 = wqrsm_data(wqrsm_lt_tt);
end
if wqrsm_lt_tt-wqrsm_LPn>0 && (wqrsm_lt_tt-wqrsm_LPn)<=length(wqrsm_data)
v1 = wqrsm_data(wqrsm_lt_tt-wqrsm_LPn);
end
if wqrsm_lt_tt-wqrsm_LP2n>0 && (wqrsm_lt_tt-wqrsm_LP2n)<=length(wqrsm_data)
v2 = wqrsm_data(wqrsm_lt_tt-wqrsm_LP2n);
end
if v0~=-32768 && v1~=-32768 && v2~=-32768
wqrsm_Yn = 2*wqrsm_Yn1 - wqrsm_Yn2 + v0 - 2*v1 + v2;
end
dy = fix((wqrsm_Yn - wqrsm_Yn1) / wqrsm_LP2n); % /* lowpass derivative of input */
wqrsm_lt_tt=wqrsm_lt_tt+1;
wqrsm_et = fix(sqrt(wqrsm_lfsc +dy*dy)); % /* length transform */
id = mod(wqrsm_lt_tt,wqrsm_BUFLN);
if id == 0
id = wqrsm_BUFLN;
end
wqrsm_ebuf(id) = wqrsm_et;
id2 = mod(wqrsm_lt_tt-wqrsm_LTwindow,wqrsm_BUFLN);
if id2 == 0
id2 = wqrsm_BUFLN;
end
wqrsm_aet = wqrsm_aet + (wqrsm_et - wqrsm_ebuf(id2));
wqrsm_lbuf(id) = wqrsm_aet;
% /* wqrsm_lbuf contains the average of the length-transformed samples over
% the interval from tt-wqrsm_LTwindow+1 to tt */
end
id3 = mod(t,wqrsm_BUFLN);
if id3 == 0
id3 = wqrsm_BUFLN;
end
lt_data = wqrsm_lbuf(id3);
end