/* file: sqrs125.c G. Moody 27 October 1990
Last revised: 22 March 2018
-------------------------------------------------------------------------------
sqrs125: Single-channel QRS detector for data sampled at 100 - 150 Hz
Copyright (C) 1990-2010 George B. Moody
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, see <http://www.gnu.org/licenses/>.
You may contact the author by e-mail (wfdb@physionet.org) or postal mail
(MIT Room E25-505A, Cambridge, MA 02139 USA). For updates to this software,
please visit PhysioNet (http://www.physionet.org/).
_______________________________________________________________________________
The detector algorithm is based on example 10 in the WFDB Programmer's
Guide, which in turn is based on a Pascal program written by W.A.H. Engelse
and C. Zeelenberg, "A single scan algorithm for QRS-detection and feature
extraction", Computers in Cardiology 6:37-42 (1979). `sqrs' does not include
the feature extraction capability of the Pascal program. The output of sqrs125
is an annotation file (with annotator name `qrs') in which all detected beats
are labelled normal; the annotation file may also contain `artifact'
annotations at locations which `sqrs125' believes are noise-corrupted.
The filter calculations and default thresholds used by this version have been
optimized for signals sampled at 125 Hz, and may be used for signals sampled
at rates between 100 and 150 Hz. (The original program, sqrs.c, was intended
for use with signals sampled at 250 Hz).
sqrs125 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). For best results on adult human ECGs, use `xform'
to resample the input signal at 125 Hz if a significantly different sampling
frequency was used originally, or use the original program if the sampling
frequency is closer to 250 Hz than to 125 Hz. For other ECGs, it may be
necessary to experiment with the input sampling frequency and the time
constants indicated below.
This program is provided as an example only, and is not intended for any
clinical application. At the time the algorithm was originally published,
its performance was typical of state-of-the-art QRS detectors. Recent designs,
particularly those which can analyze two or more input signals, may exhibit
significantly better performance.
Usage:
sqrs125 -r RECORD [ OPTIONS ]
where RECORD is the record name, and OPTIONS may include:
-f TIME to specify the starting TIME (default: the beginning of
the record)
-m THRESHOLD to specify the detection THRESHOLD (default: 250);
use higher values to reduce false detections, or lower
values to reduce the number of missed beats
-s SIGNAL to specify the SIGNAL to be used for QRS detection
(default: 0)
-t TIME to specify the ending TIME (default: the end of the
record)
For example, to mark QRS complexes in record 100 beginning 5 minutes from the
start, ending 10 minutes and 35 seconds from the start, and using signal 1, use
the command:
sqrs125 -r 100 -f 5:0 -t 10:35 -s 1
The output may be read using (for example):
rdann -a qrs -r 100
To evaluate the performance of this program, run it on the entire record, by:
sqrs125 -r 100
and then compare its output with the reference annotations by:
bxb -r 100 -a atr qrs
*/
#include <stdio.h>
#include <wfdb/wfdb.h>
#include <wfdb/ecgcodes.h>
#define abs(A) ((A) >= 0 ? (A) : -(A))
char *pname;
main(argc, argv)
int argc;
char *argv[];
{
char *p, *record = NULL, *prog_name();
int filter, i, minutes = 0, nsig, time = 0,
slopecrit, sign, maxslope = 0, nslope = 0,
qtime, maxtime, t0, t1, t2, t3, t4, t5,
ms160, ms200, s2, scmax, scmin = 250, signal = -1, *v;
long from = 0L, next_minute, now, spm, to = 0L;
WFDB_Anninfo a;
WFDB_Annotation annot;
static int gvmode = 0;
static WFDB_Siginfo *s;
void help();
pname = prog_name(argv[0]);
for (i = 1; i < argc; i++) {
if (*argv[i] == '-') switch (*(argv[i]+1)) {
case 'f': /* starting time */
if (++i >= argc) {
(void)fprintf(stderr, "%s: time must follow -f\n", pname);
exit(1);
}
from = i;
break;
case 'h': /* help requested */
help();
exit(0);
break;
case 'H': /* operate in WFDB_HIGHRES mode */
gvmode = WFDB_HIGHRES;
break;
case 'm': /* threshold */
if (++i >= argc) {
(void)fprintf(stderr, "%s: threshold must follow -m\n", pname);
exit(1);
}
scmin = atoi(argv[i]);
break;
case 'r': /* record name */
if (++i >= argc) {
(void)fprintf(stderr, "%s: input record name must follow -r\n",
pname);
exit(1);
}
record = argv[i];
break;
case 's': /* signal */
if (++i >= argc) {
(void)fprintf(stderr,
"%s: signal number or name must follow -s\n",
pname);
exit(1);
}
signal = i;
break;
case 't': /* end time */
if (++i >= argc) {
(void)fprintf(stderr, "%s: time must follow -t\n",pname);
exit(1);
}
to = i;
break;
default:
(void)fprintf(stderr, "%s: unrecognized option %s\n", pname,
argv[i]);
exit(1);
}
else {
(void)fprintf(stderr, "%s: unrecognized argument %s\n", pname,
argv[i]);
exit(1);
}
}
if (record == NULL) {
help();
exit(1);
}
if (gvmode == 0 && (p = getenv("WFDBGVMODE")))
gvmode = atoi(p);
setgvmode(gvmode|WFDB_GVPAD);
if ((nsig = isigopen(record, NULL, 0)) < 1) exit(2);
if ((s = malloc(nsig * sizeof(WFDB_Siginfo))) == NULL ||
(v = malloc(nsig * sizeof(WFDB_Sample))) == NULL) {
(void)fprintf(stderr, "%s: insufficient memory\n", pname);
exit(2);
}
if ((nsig = isigopen(record, s, nsig)) < 1) exit(2);
if (sampfreq((char *)NULL) < 50.) {
(void)fprintf(stderr, "%s: sampling frequency (%g Hz) is too low%s",
pname, sampfreq((char *)NULL),
gvmode & WFDB_HIGHRES ? "\n" : ", try -H option\n");
exit(3);
}
if (sampfreq((char *)NULL) < 120. || sampfreq((char *)NULL) > 130.)
setifreq(125.);
else if (gvmode & WFDB_HIGHRES)
setafreq(sampfreq(NULL));
a.name = "qrs"; a.stat = WFDB_WRITE;
if (annopen(record, &a, 1) < 0) exit(2);
if (from > 0L) {
if ((from = strtim(argv[from])) < 0L)
from = -from;
if (isigsettime(from) < 0)
exit(2);
}
if (to > 0L) {
if ((to = strtim(argv[to])) < 0L)
to = -to;
}
spm = strtim("1:0");
next_minute = from + spm;
if (signal >= 0) signal = findsig(argv[signal]);
if (signal < 0 || signal >= nsig) signal = 0;
scmin = muvadu((unsigned)signal, scmin);
if (scmin < 1) scmin = muvadu((unsigned)signal, 250);
slopecrit = scmax = 10 * scmin;
now = from;
/* These time constants may need adjustment for pediatric or small
mammal ECGs. */
ms160 = strtim("0.16"); ms200 = strtim("0.2"); s2 = strtim("2");
annot.subtyp = annot.chan = annot.num = 0; annot.aux = NULL;
(void)getvec(v);
t5 = t4 = t3 = t2 = t1 = v[signal];
do {
filter = (t0 = v[signal]) + 2*t1 + t2 - t3 - 2*t4 - t5;
if (time % s2 == 0) {
if (nslope == 0) {
slopecrit -= slopecrit >> 4;
if (slopecrit < scmin) slopecrit = scmin;
}
else if (nslope >= 5) {
slopecrit += slopecrit >> 4;
if (slopecrit > scmax) slopecrit = scmax;
}
}
if (nslope == 0 && abs(filter) > slopecrit) {
nslope = 1; maxtime = ms160;
sign = (filter > 0) ? 1 : -1;
qtime = time;
}
if (nslope != 0) {
if (filter * sign < -slopecrit) {
sign = -sign;
maxtime = (++nslope > 4) ? ms200 : ms160;
}
else if (filter * sign > slopecrit &&
abs(filter) > maxslope)
maxslope = abs(filter);
if (maxtime-- < 0) {
if (2 <= nslope && nslope <= 4) {
slopecrit += ((maxslope>>2) - slopecrit) >> 3;
if (slopecrit < scmin) slopecrit = scmin;
else if (slopecrit > scmax) slopecrit = scmax;
annot.time = now - (time - qtime) - 4;
annot.anntyp = NORMAL; (void)putann(0, &annot);
time = 0;
}
else if (nslope >= 5) {
annot.time = now - (time - qtime) - 4;
annot.anntyp = ARFCT; (void)putann(0, &annot);
}
nslope = 0;
}
}
t5 = t4; t4 = t3; t3 = t2; t2 = t1; t1 = t0; time++; now++;
if (now >= next_minute) {
next_minute += spm;
(void)fprintf(stderr, ".");
(void)fflush(stderr);
if (++minutes >= 60) {
(void)fprintf(stderr, " %s\n", timstr(now));
minutes = 0;
}
}
} while (getvec(v) > 0 && (to == 0L || now <= to));
if (minutes) (void)fprintf(stderr, " %s\n", timstr(now));
wfdbquit();
exit(0); /*NOTREACHED*/
}
char *prog_name(s)
char *s;
{
char *p = s + strlen(s);
#ifdef MSDOS
while (p >= s && *p != '\\' && *p != ':') {
if (*p == '.')
*p = '\0'; /* strip off extension */
if ('A' <= *p && *p <= 'Z')
*p += 'a' - 'A'; /* convert to lower case */
p--;
}
#else
while (p >= s && *p != '/')
p--;
#endif
return (p+1);
}
static char *help_strings[] = {
"usage: %s -r RECORD [OPTIONS ...]\n",
"where RECORD is the name of the record to be analyzed, and OPTIONS may",
"include any of:",
" -f TIME begin at specified time",
" -h print this usage summary",
" -H read multifrequency signals in high resolution mode",
" -m THRESH set detector threshold to THRESH (default: 250)",
" -s SIGNAL analyze specified signal (default: 0)",
" -t TIME stop at specified time",
"If too many beats are missed, decrease THRESH; if there are too many extra",
"detections, increase THRESH.",
NULL
};
void help()
{
int i;
(void)fprintf(stderr, help_strings[0], pname);
for (i = 1; help_strings[i] != NULL; i++)
(void)fprintf(stderr, "%s\n", help_strings[i]);
}