diff -Naur --exclude Makefile --exclude info wfdb-10.5.12/app/gqfuse.c wfdb-10.5.13/app/gqfuse.c --- wfdb-10.5.12/app/gqfuse.c 1969-12-31 19:00:00.000000000 -0500 +++ wfdb-10.5.13/app/gqfuse.c 2012-05-07 03:42:34.000000000 -0400 @@ -0,0 +1,318 @@ +/* file: gqfuse.c G. Moody 6 May 2012 + Last revised: 7 May 2012 +------------------------------------------------------------------------------- +gqfuse: combine QRS annotation files +Copyright (C) 2012 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, 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 (george@mit.edu) or postal mail +(MIT Room E25-505A, Cambridge, MA 02139 USA). For updates to this software, +please visit PhysioNet (http://www.physionet.org/). +_______________________________________________________________________________ + +This program produces a QRS annotation file based on two or more input QRS +annotation files A1, A2, .... Each one-minute segment of the output annotation +file is a copy of the corresponding segment of one of the input annotation +files. In each segment, the program copies the input that best matches a +predicted heart rate. If there are N inputs, the prediction is the median of +N+1 values (the previous prediction and the number of beats marked within the +current segment of each of the N input files). Although this process allows +the input to be switched once per minute, the policy for resolving ties (within +2 beats) favors not switching if the previously chosen input is one of those +belonging to the tie. + +One way to use this program is to combine input annotation files for each +available ECG signal in a record, made using a single detector such as gqrs. +Another is to combine input annotation files made using a variety of QRS +detectors. These ideas can be combined as desired. +*/ + +#include +#include +#include +#include + +/* The `getconf' macro is used to check a line of input (already in `buf') for + the string named by getconf's first argument. If the string is found, the + value following the string (and an optional `:' or `=') is converted using + sscanf and the format specifier supplied as getconf's second argument, and + stored in the variable named by the first argument. */ +#define getconf(a, fmt) if (p=strstr(buf,#a)){sscanf(p,#a "%*[=: \t]" fmt,&a);} + +void help(void); +char *prog_name(char *p); +void *gcalloc(size_t nmemb, size_t size); +void cleanup(int status); + +char *pname; /* name of this program, used in messages */ +char *record = NULL; /* name of input record */ + +int icomp(const void *a, const void *b) +{ + if (*(int *)a > *(int *)b) return -1; + if (*(int *)a < *(int *)b) return 1; + return 0; +} + +main(int argc, char **argv) +{ + char *oaname = "gqf"; + double HR; + FILE *config = NULL; + int a0 = 0, a1 = 0, **count, dn, i, ibest, ichosen = 0, j, *mbuf, + n, niann = 0, nseg, tdn; + WFDB_Anninfo *a; + WFDB_Annotation annot; + WFDB_Frequency spm; + WFDB_Time t0, tf; + + pname = prog_name(argv[0]); + + for (i = 1; i < argc; i++) { + if (*argv[i] == '-') switch (*(argv[i]+1)) { + case 'a': /* input annotator names */ + for (a0 = a1 = ++i; a1 < argc && *argv[a1] != '-'; a1++, i++) + ; + if ((niann = a1 - a0) < 2) { + (void)fprintf(stderr, + "%s: input annotator names must follow -a\n", + pname); + cleanup(1); + } + break; + case 'c': /* configuration file */ + if (++i >= argc) { + (void)fprintf(stderr, + "%s: name of configuration file must follow -c\n", pname); + cleanup(1); + } + if ((config = fopen(argv[i], "rt")) == NULL) { + (void)fprintf(stderr, + "%s: can't read configuration file %s\n", pname, argv[i]); + cleanup(1); + } + break; + case 'h': /* help requested */ + help(); + cleanup(0); + break; + case 'o': /* output annotator name */ + if (++i >= argc) { + (void)fprintf(stderr, + "%s: output annotator name must follow -o\n", + pname); + cleanup(1); + } + oaname = argv[i]; + break; + case 'r': /* record name */ + if (++i >= argc) { + (void)fprintf(stderr, "%s: input record name must follow -r\n", + pname); + cleanup(1); + } + record = argv[i]; + break; + default: + (void)fprintf(stderr, "%s: unrecognized option %s\n", pname, + argv[i]); + cleanup(1); + } + else { + (void)fprintf(stderr, "%s: unrecognized argument %s\n", pname, + argv[i]); + cleanup(1); + } + } + + /* Quit unless record name and at least two annotator names specified. */ + if (record == NULL || niann < 2) { + help(); + cleanup(1); + } + + /* Set up annotator info. */ + a = gcalloc(sizeof(WFDB_Anninfo), niann+1); + for (i = a0, j = 0; i < a1; i++, j++) { + a[j].name = argv[i]; a[j].stat = WFDB_READ; + } + a[j].name = oaname; a[j].stat = WFDB_WRITE; + + /* Read a priori physiologic parameters from the configuration file if + available. They can be adjusted in the configuration file for pediatric, + fetal, or animal ECGs. */ + if (config) { + char buf[256], *p; + + /* Read the configuration file a line at a time. */ + while (fgets(buf, sizeof(buf), config)) { + /* Skip comments (empty lines and lines beginning with `#'). */ + if (buf[0] == '#' || buf[0] == '\n') continue; + /* Set parameters. Each `getconf' below is executed once for + each non-comment line in the configuration file. */ + getconf(HR, "%lf"); + } + fclose(config); + } + + /* If any a priori parameters were not specified in the configuration file, + initialize them here (using values chosen for adult human ECGs). */ + if (HR == 0.0) HR = 75; + + spm = 60.0 * sampfreq(record); /* sample intervals per minute */ + nseg = strtim("e")/spm; /* number of complete 1-minute intervals */ + + /* Set up buffer for median calculation. */ + mbuf = gcalloc(sizeof(int), niann + 1); + mbuf[niann] = HR; /* initial median is HR */ + + /* Set up arrays for HR minute-by-minute time series. */ + count = gcalloc(sizeof(int *), niann + 1); + for (i = 0; i < niann; i++) + count[i] = gcalloc(sizeof(int), nseg+1); + + /* Pass 1: Load the HR arrays. */ + if (annopen(record, a, niann) < 0) cleanup(2); + for (i = 0; i < niann; i++) { + for (tf = spm, j = 0; j < nseg+1; tf += spm, j++) { + n = 0; + while (getann(i, &annot) == 0 && annot.time < tf) + if (isqrs(annot.anntyp)) n++; + count[i][j] = n; + } + } + wfdbquit(); + + /* Pass 2: Select the input annotations to be copied, and copy them. */ + if (annopen(record, a, niann+1) < 0) cleanup(2); + for (j = 0, t0 = 0, tf = spm; j < nseg+1; j++, t0 = tf, tf += spm) { + for (i = 0; i < niann; i++) + mbuf[i] = count[i][j]; + /* Note that the last element, mbuf[niann], is the previous median. */ + qsort(mbuf, niann+1, sizeof(int), icomp); + /* If mbuf has an even number of members, the median is the mean of + the two central values. */ + if (niann & 1) mbuf[niann] = (mbuf[(niann-1)/2] + mbuf[(niann+1)/2])/2; + /* Otherwise it's the single central value. */ + else mbuf[niann] = mbuf[niann/2]; + /* Find the input that best matches the prediction (the median). */ + for (i = 0, dn = 9999; i < niann; i++) { + tdn = mbuf[niann] - count[i][j]; + if (tdn < 0) tdn = -tdn; + if (tdn < dn) { dn = tdn; ibest = i; } + } + /* If the best match is not the previously chosen input, avoid switching + the input unless it's more than 2 beats further from the median. */ + if (ibest != ichosen) { + tdn = mbuf[niann] - count[ichosen][j]; + if (tdn < 0) tdn = -tdn; + if (tdn < dn + 3) ibest = ichosen; + } + ichosen = ibest; + /* Copy the chosen annotations for this one-minute segment. */ + n = 0; + while (getann(ichosen, &annot) == 0 && n <= count[ichosen][j]) { + if (annot.time >= t0 && isqrs(annot.anntyp)) { + putann(0, &annot); + n++; + } + } + ungetann(ichosen, &annot); + } + cleanup(0); +} + +/* prog_name() extracts this program's name from argv[0], for use in error and + warning messages. */ + +char *prog_name(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); +} + +/* help() prints a (very) concise summary of how to use this program. + A more detailed summary is in the man page (gqpost.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:", + " -a ANNOTATOR read annotations from the specified ANNOTATOR (default: qrs)", + " -c FILE initialize parameters from the specified configuration FILE", + " -f TIME begin processing at specified time", + " -h print this usage summary", + " -H read multifrequency signals in high resolution mode", + " -m THRESH set interpolated event acceptance threshold to THRESH", + " (default: 1)", + " -o ANNOTATOR write annotations to the specified ANNOTATOR (default: gqp)", + " -t TIME stop processing at specified time", + "If too many true beats are rejected, decrease THRESH; if too many false", + "detections are accepted, increase THRESH.", + "Note that the output is a complete copy of the input (with rejected events", + "flagged as ARFCT). The -f and -t options only limit the interval during", + "which events may be rejected.", +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]); +} + +/* gcalloc() is a wrapper for calloc() that handles errors and maintains + a list of allocated buffers for automatic on-exit deallocation via + cleanup(). */ + +size_t gcn = 0, gcnmax = 0; +void **gclist = NULL; + +void *gcalloc(size_t nmemb, size_t size) +{ + void *p = calloc(nmemb, size); + + if ((p == NULL) || + ((gcn >= gcnmax) && + (gclist = realloc(gclist, (gcnmax += 32)*(sizeof(void *)))) == NULL)) + cleanup(3); /* could not allocate requested memory */ + return (gclist[gcn++] = p); +} + +void cleanup(int status) /* close files and free allocated memory */ +{ + if (status == 3) + fprintf(stderr, "%s: insufficient memory (%d)\n", pname, status); + if (record) wfdbquit(); + while (gcn-- > 0) + if (gclist[gcn]) free(gclist[gcn]); + exit(status); +} diff -Naur --exclude Makefile --exclude info wfdb-10.5.12/app/gqpost.c wfdb-10.5.13/app/gqpost.c --- wfdb-10.5.12/app/gqpost.c 1969-12-31 19:00:00.000000000 -0500 +++ wfdb-10.5.13/app/gqpost.c 2012-05-06 11:56:50.000000000 -0400 @@ -0,0 +1,332 @@ +/* file: gqpost.c G. Moody 16 November 2006 + Last revised: 6 May 2012 +------------------------------------------------------------------------------- +gqpost: A post-processor for gqrs +Copyright (C) 2006-2012 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, 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 (george@mit.edu) or postal mail +(MIT Room E25-505A, Cambridge, MA 02139 USA). For updates to this software, +please visit PhysioNet (http://www.physionet.org/). +_______________________________________________________________________________ + +The companion to this program, gqrs, is a QRS detector that is optimized for +sensitivity. To reduce the number of false positives, gqpost can read the +output annotation file produced by gqrs and can look for interpolated beat +detections in it. Such events are identified by looking for pairs or longer +groups of inter-beat intervals that sum to values close to the predicted +interbeat interval (in other words, an interpolated event is one that occurs +between a pair of beats that occur at their expected times). Although +true beats can be interpolated, such beats are relatively rare in comparison +with false detections that are almost always interpolated, so a high proportion +of interpolated events is false, and a simple strategy for improving positive +predictivity is to reject all interpolated events. To avoid rejecting +(at least some of) the true interpolated beats, gqpost also looks at the size +of each interpolated event as recorded by gqrs in the annotation 'num' field, +and it rejects only those events that are smaller than a threshold value +(which may be adjusted using the -m option). + +gqpost is most effective in the context of regular rhythms. During highly +irregular rhythms, such as atrial fibrillation, the predicted time of the +next beat is highly uncertain, and gqpost rejects very few if any events in +such cases. + +gqpost is very fast since it operates entirely on the input annotation +file, without reference to the raw signals. The output of gqpost is a +complete copy of its input, except that the anntyp fields of any rejected +events occurring between the start and end times (specified on the command +line using the -f and -t options) are set to ARFCT. gqpost can accept its +own output as input, which allows one to vary the rejection threshold if +necessary by use of multiple post-processing passes. +*/ + +#include +#include +#include +#include +#define PBLEN 12 /* size of predictor array */ + +/* The `getconf' macro is used to check a line of input (already in `buf') for + the string named by getconf's first argument. If the string is found, the + value following the string (and an optional `:' or `=') is converted using + sscanf and the format specifier supplied as getconf's second argument, and + stored in the variable named by the first argument. */ +#define getconf(a, fmt) if (p=strstr(buf,#a)){sscanf(p,#a "%*[=: \t]" fmt,&a);} + +void help(void); +char *prog_name(char *p); +void *gcalloc(size_t nmemb, size_t size); +void cleanup(int status); + +char *pname; /* name of this program, used in messages */ +char *record = NULL; /* name of input record */ + +main(int argc, char **argv) +{ + FILE *config = NULL; + int dt, i, meanrr, meanrrd, minrrd, state = 0, thresh = 20; + static double HR, RR, RRdelta; + WFDB_Anninfo a[2]; + WFDB_Annotation annot; + WFDB_Frequency sps; + WFDB_Time from = (WFDB_Time)0, to = (WFDB_Time)0; + + a[0].name = "qrs"; a[0].stat = WFDB_READ; + a[1].name = "gqp"; a[1].stat = WFDB_WRITE; + pname = prog_name(argv[0]); + + for (i = 1; i < argc; i++) { + if (*argv[i] == '-') switch (*(argv[i]+1)) { + case 'a': /* input annotator name */ + if (++i >= argc) { + (void)fprintf(stderr, + "%s: input annotator name must follow -a\n", + pname); + cleanup(1); + } + a[0].name = argv[i]; + break; + case 'c': /* configuration file */ + if (++i >= argc) { + (void)fprintf(stderr, + "%s: name of configuration file must follow -c\n", pname); + cleanup(1); + } + if ((config = fopen(argv[i], "rt")) == NULL) { + (void)fprintf(stderr, + "%s: can't read configuration file %s\n", pname, argv[i]); + cleanup(1); + } + break; + case 'f': /* starting time */ + if (++i >= argc) { + (void)fprintf(stderr, "%s: time must follow -f\n", pname); + cleanup(1); + } + from = i; + break; + case 'h': /* help requested */ + help(); + cleanup(0); + break; + case 'm': /* threshold */ + if (++i >= argc) { + (void)fprintf(stderr, "%s: threshold must follow -m\n", pname); + cleanup(1); + } + thresh = 10*atof(argv[i]) + 0.5; + break; + case 'o': /* output annotator name */ + if (++i >= argc) { + (void)fprintf(stderr, + "%s: output annotator name must follow -o\n", + pname); + cleanup(1); + } + a[1].name = argv[i]; + break; + case 'r': /* record name */ + if (++i >= argc) { + (void)fprintf(stderr, "%s: input record name must follow -r\n", + pname); + cleanup(1); + } + record = argv[i]; + break; + case 't': /* end time */ + if (++i >= argc) { + (void)fprintf(stderr, "%s: time must follow -t\n", pname); + cleanup(1); + } + to = i; + break; + default: + (void)fprintf(stderr, "%s: unrecognized option %s\n", pname, + argv[i]); + cleanup(1); + } + else { + (void)fprintf(stderr, "%s: unrecognized argument %s\n", pname, + argv[i]); + cleanup(1); + } + } + if (record == NULL) { + help(); + cleanup(1); + } + + /* Read a priori physiologic parameters from the configuration file if + available. They can be adjusted in the configuration file for pediatric, + fetal, or animal ECGs. */ + if (config) { + char buf[256], *p; + + /* Read the configuration file a line at a time. */ + while (fgets(buf, sizeof(buf), config)) { + /* Skip comments (empty lines and lines beginning with `#'). */ + if (buf[0] == '#' || buf[0] == '\n') continue; + /* Set parameters. Each `getconf' below is executed once for + each non-comment line in the configuration file. */ + getconf(HR, "%lf"); + getconf(RR, "%lf"); + getconf(RRdelta, "%lf"); + } + fclose(config); + } + + /* If any a priori parameters were not specified in the configuration file, + initialize them here (using values chosen for adult human ECGs). */ + if (HR != 0.0) RR = 60.0/HR; + if (RR == 0.0) RR = 0.8; + if (RRdelta == 0.0) RRdelta = RR/4; + + sps = sampfreq(record); + if (annopen(record, a, 2) < 0) cleanup(2); + meanrr = RR * sps; + meanrrd = RRdelta * sps; + if ((minrrd = meanrrd/2) < 4) minrrd = 4; + if ((dt = meanrr/80) < 1) dt = 1; + + while (getann(0, &annot) == 0) { + switch (state) { + case 0: + if (annot.time < from) continue; + else state++; + /* fall through to case 1 */ + case 1: + if (to > (WFDB_Time)0 && annot.time > to) state++; + else { /* reject sub-threshold interpolated events */ + static WFDB_Time tpred, tpredmin, tprev; + static int rr, rrd; + + if (annot.time < tpredmin && + isqrs(annot.anntyp) && + annot.num < thresh) { /* event is small and early */ + WFDB_Annotation nextann; + if (getann(0, &nextann) == 0) { /* peek ahead */ + if (nextann.time <= tpred) + annot.anntyp = ARFCT; /* reject event */ + ungetann(0, &nextann); + } + } + if (isqrs(annot.anntyp)) { + rr = annot.time - tprev; + if ((rrd = rr - meanrr) < 0) rrd = -rrd; + if (rr > meanrr + dt) meanrr += dt; + else if (rr > meanrr - dt) meanrr = rr; + else meanrr -= dt; + if (rrd > meanrrd + dt) meanrrd += dt; + else if (rrd > meanrrd - dt) meanrrd = rrd; + else meanrrd -= dt; + if (meanrrd < minrrd) meanrrd = minrrd; + tpred = annot.time + meanrr + minrrd; + tpredmin = annot.time + meanrr - meanrrd; + tprev = annot.time; + } + } + break; + case 2: + /* do nothing */ + break; + } + putann(0, &annot); + } + cleanup(0); +} + + +/* prog_name() extracts this program's name from argv[0], for use in error and + warning messages. */ + +char *prog_name(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); +} + +/* help() prints a (very) concise summary of how to use this program. + A more detailed summary is in the man page (gqpost.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:", + " -a ANNOTATOR read annotations from the specified ANNOTATOR (default: qrs)", + " -c FILE initialize parameters from the specified configuration FILE", + " -f TIME begin processing at specified time", + " -h print this usage summary", + " -H read multifrequency signals in high resolution mode", + " -m THRESH set interpolated event acceptance threshold to THRESH", + " (default: 1)", + " -o ANNOTATOR write annotations to the specified ANNOTATOR (default: gqp)", + " -t TIME stop processing at specified time", + "If too many true beats are rejected, decrease THRESH; if too many false", + "detections are accepted, increase THRESH.", + "Note that the output is a complete copy of the input (with rejected events", + "flagged as ARFCT). The -f and -t options only limit the interval during", + "which events may be rejected.", +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]); +} + +/* gcalloc() is a wrapper for calloc() that handles errors and maintains + a list of allocated buffers for automatic on-exit deallocation via + cleanup(). */ + +size_t gcn = 0, gcnmax = 0; +void **gclist = NULL; + +void *gcalloc(size_t nmemb, size_t size) +{ + void *p = calloc(nmemb, size); + + if ((p == NULL) || + ((gcn >= gcnmax) && + (gclist = realloc(gclist, (gcnmax += 32)*(sizeof(void *)))) == NULL)) + cleanup(3); /* could not allocate requested memory */ + return (gclist[gcn++] = p); +} + +void cleanup(int status) /* close files and free allocated memory */ +{ + if (status == 3) + fprintf(stderr, "%s: insufficient memory (%d)\n", pname, status); + if (record) wfdbquit(); + while (gcn-- > 0) + if (gclist[gcn]) free(gclist[gcn]); + exit(status); +} diff -Naur --exclude Makefile --exclude info wfdb-10.5.12/app/gqrs.c wfdb-10.5.13/app/gqrs.c --- wfdb-10.5.12/app/gqrs.c 1969-12-31 19:00:00.000000000 -0500 +++ wfdb-10.5.13/app/gqrs.c 2012-05-06 13:56:09.000000000 -0400 @@ -0,0 +1,755 @@ +/* file: gqrs.c G. Moody 16 November 2006 + Last revised: 23 April 2012 +------------------------------------------------------------------------------- +gqrs: A QRS detector +Copyright (C) 2006-2012 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, 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 (george@mit.edu) or postal mail +(MIT Room E25-505A, Cambridge, MA 02139 USA). For updates to this software, +please visit PhysioNet (http://www.physionet.org/). +_______________________________________________________________________________ +*/ + +#include +#include +#include +#include + +#define BUFLN 32768 /* must be a power of 2, see qf() */ +#define NPEAKS 64 /* number of peaks buffered (per signal) */ + +/* The `getconf' macro is used by gqrs_init() (below) to check a line of input + (already in `buf', defined in gqrs_init) for the string named by getconf's + first argument. If the string is found, the value following the string + (and an optional `:' or `=') is converted using sscanf and the format + specifier supplied as getconf's second argument, and stored in the variable + named by the first argument. */ +#define getconf(a, fmt) if (p=strstr(buf,#a)){sscanf(p,#a "%*[=: \t]" fmt,&a);} + +/* The 'peak' structure contains information about a local maximum in the + filtered signal (qfv). Peaks are stored in circular buffers of peak + structures, one buffer per input signal. The time of a flat-topped peak + is the time of the first sample that has the maximum value. A peak is + secondary if there is a larger peak within its neighborhood (time +- rrmin), + of if it has been identified as a T-wave associated with a previous primary + peak. A peak is primary if it is largest in its neighborhood, or if the + only larger peaks are secondary. */ +struct peak { + struct peak *prev, *next; /* pointers to neighbors (in time) */ + WFDB_Time time; /* time of local maximum of qfv */ + int amp; /* value of qfv at time of peak */ + short type; /* 1: primary, 2: secondary, 0: not determined */ +} *peaks, *cpeak; + +/* Prototypes of functions defined below. The definitions of these functions + follow that of main(), in the order shown below. */ +WFDB_Sample sm(WFDB_Time t); +void qf(void); +void addpeak(WFDB_Time t, int peak_amplitude); +int peaktype(struct peak *p); +void gqrs_init(WFDB_Time from, WFDB_Time to); +void rewind_gqrs(void); +struct peak *find_missing(struct peak *previous_peak, struct peak *next_peak); +void gqrs(WFDB_Time from, WFDB_Time to); +void help(void); +char *prog_name(char *p); +void *gcalloc(size_t nmemb, size_t size); +void cleanup(int status); + +char auxbuf[1+255+1]; /* 'aux' string buffer for annotations */ +char *pname; /* name of this program, used in messages */ +char *record = NULL; /* name of input record */ +FILE *config = NULL; /* configuration file, if any */ +int countdown = -1; /* if > 0, ticks remaining (see gqrs()) */ +int debug; /* if non-zero, generate debugging output */ +int minutes = 0; /* minutes elapsed in the current hour */ +int nsig; /* number of signals */ +int pthr; /* peak-detection threshold */ +int qthr; /* QRS-detection threshold */ +int pthmin, qthmin; /* minimum values for pthr and qthr */ +int rrdev; /* mean absolute deviation of RR from rrmean */ +int rrinc; /* maximum incremental change in rrmean */ +int rrmean; /* mean RR interval, in sample intervals */ +int rrmax; /* maximum likely RR interval */ +int rrmin; /* minimum RR interval, in sample intervals */ +int rtmax; /* maximum RT interval, in sample intervals */ +int rtmin; /* minimum RT interval, in sample intervals */ +int rtmean; /* mean RT interval, in sample intervals */ +int tamean; /* mean T-wave amplitude in qfv */ +double thresh = 1.0; /* normalized detection threshold */ +long v1; /* integral of dv in qf() */ +long v1norm; /* normalization for v1 */ +WFDB_Annotation annot; /* most recently written annotation */ +WFDB_Sample *v; /* latest sample of each input signal */ +WFDB_Siginfo *si; /* characteristics of each input signal */ +WFDB_Signal sig; /* signal number of signal to be analyzed */ +WFDB_Time next_minute; /* sample number of start of next minute */ +WFDB_Time spm; /* sample intervals per minute */ +WFDB_Time sps; /* sample intervals per second */ +WFDB_Time t; /* time of the current sample being analyzed */ +WFDB_Time t0; /* time of the start of the analysis period */ +WFDB_Time tf; /* time of the end of the analysis period */ +WFDB_Time tf_learn; /* time of the end of the learning period */ + +/* Current operating mode of the detector */ +enum { LEARNING, RUNNING, CLEANUP } state = LEARNING; + +int dt, dt2, dt3, dt4, smdt; /* time intervals set by gqrs_init() depending + on the sampling frequency and used by + filter functions sm() and qf(); units + are sample intervals */ +WFDB_Time smt = 0, smt0; /* current and starting time for sm() */ +long *qfv, *smv; /* filter buffers allocated by gqrs_init() */ + +/* The q() and s() macros can be used for fast lookup of filter outputs. */ + +#define q(T) (qfv[(T)&(BUFLN-1)]) +#define s(T) (smv[(T)&(BUFLN-1)]) + +main(int argc, char **argv) +{ + char *p; + int gvmode = 0, i, isiglist = 0, j, nisig, s; + WFDB_Anninfo a; + + pname = prog_name(argv[0]); + a.name = "qrs"; a.stat = WFDB_WRITE; + + for (i = 1; i < argc; i++) { + if (*argv[i] == '-') switch (*(argv[i]+1)) { + case 'c': /* configuration file */ + if (++i >= argc) { + (void)fprintf(stderr, + "%s: name of configuration file must follow -c\n", pname); + cleanup(1); + } + if ((config = fopen(argv[i], "rt")) == NULL) { + (void)fprintf(stderr, + "%s: can't read configuration file %s\n", pname, argv[i]); + cleanup(1); + } + break; + case 'd': /* record debugging info in the annotation file */ + debug = 1; + break; + case 'f': /* starting time */ + if (++i >= argc) { + (void)fprintf(stderr, "%s: time must follow -f\n", pname); + cleanup(1); + } + t0 = i; + break; + case 'h': /* help requested */ + help(); + cleanup(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); + cleanup(1); + } + thresh = atof(argv[i]); + break; + case 'o': /* write output annotations as specified annotator */ + if (++i >= argc) { + (void)fprintf(stderr,"%s: annotator name must follow -o\n", + pname); + cleanup(1); + } + a.name = argv[i]; + break; + case 'r': /* record name */ + if (++i >= argc) { + (void)fprintf(stderr, "%s: input record name must follow -r\n", + pname); + cleanup(1); + } + record = argv[i]; + break; + case 's': /* signal name or number follows */ + if (++i >= argc) { + (void)fprintf(stderr, "%s: signal name or # must follow -s\n", + pname); + cleanup(1); + } + sig = i; + break; + case 't': /* end time */ + if (++i >= argc) { + (void)fprintf(stderr, "%s: time must follow -t\n",pname); + cleanup(1); + } + tf = i; + break; + default: + (void)fprintf(stderr, "%s: unrecognized option %s\n", pname, + argv[i]); + cleanup(1); + } + else { + (void)fprintf(stderr, "%s: unrecognized argument %s\n", pname, + argv[i]); + cleanup(1); + } + } + + if (record == NULL) { + help(); + cleanup(1); + } + + if (gvmode == 0 && (p = getenv("WFDBGVMODE"))) + gvmode = atoi(p); + setgvmode(gvmode|WFDB_GVPAD); + + if ((nsig = isigopen(record, NULL, 0)) < 1) cleanup(2); + si = (WFDB_Siginfo *)gcalloc((size_t)nsig, sizeof(WFDB_Siginfo)); + v = (WFDB_Sample *)gcalloc((size_t)nsig, sizeof(WFDB_Sample)); + if ((nsig = wfdbinit(record, &a, 1, si, nsig)) < 1) cleanup(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"); + cleanup(3); + } + if (t0 > 0L && (t0 = strtim(argv[t0])) < 0L) + t0 = -t0; + if (tf > 0L) { + if ((tf = strtim(argv[tf])) < 0L) + tf = -tf; + tf += sps; /* make sure that the last beat before 'to' is marked */ + } + else + tf = strtim("e"); + t = t0; + sps = strtim("1"); + spm = strtim("1:0"); + next_minute = t + spm; + + if (sig > 0) { + int i = findsig(argv[sig]); + if (i < 0) { + (void)fprintf(stderr, + "%s: (warning) no signal %d in record %s\n", + pname, s, record); + cleanup(4); + } + sig = i; + } + + /* Record the argument list in a NOTE annotation. */ + for (i = j = 0, p = auxbuf+1; i < argc; i++) { + int len; + j += len = strlen(argv[i]); + if (j < 255) { + strncpy(p, argv[i], len); + p += len; + *p++ = ' '; + j++; + } + else { + j -= len; + break; + } + } + *(--p) = '\0'; + auxbuf[0] = j; + annot.subtyp = annot.chan = annot.num = 0; + annot.aux = auxbuf; + annot.anntyp = NOTE; + annot.time = (WFDB_Time)0; + putann(0, &annot); + + gqrs_init(t0, tf); /* initialize variables for gqrs() */ + if (spm >= BUFLN) { + if (tf - t0 > BUFLN) tf_learn = t0 + BUFLN - dt4; + else tf_learn = tf - dt4; + } + else { + if (tf - t0 > spm) tf_learn = t0 + spm - dt4; + else tf_learn = tf - dt4; + } + + state = LEARNING; + gqrs(t0, tf_learn); + + rewind_gqrs(); + state = RUNNING; + t = t0 - dt4; + gqrs(t0, tf); /* run the detector and collect output */ + + printf(" %s\n", timstr(strtim("i"))); + cleanup(0); +} + +/* sm() implements a trapezoidal low pass (smoothing) filter (with a gain of + 4*smdt) applied to input signal sig before the QRS matched filter qf(). + Before attempting to 'rewind' by more than BUFLN-smdt samples, reset smt + and smt0. + */ + +WFDB_Sample sm(WFDB_Time t) +{ + while (t > smt) { + int i; + + if (++smt > smt0) { /* fast update by summing first differences */ + s(smt) = s(smt-1) + + sample(sig, smt+smdt) + sample(sig, smt+smdt-1) - + sample(sig, smt-smdt) - sample(sig, smt-smdt-1); + } + else { /* get initial value by full convolution */ + int j, v; + + for (j = 1, v = sample(sig, smt); j < smdt; j++) + v += sample(sig, smt+j) + sample(sig, smt-j); + s(smt) = (v << 1) + sample(sig, smt+j) + sample(sig, smt-j) + - si[sig].adczero * (smdt << 2); /* FIXME: needed? */ + } + } + return (s(t)); +} + +void qf() /* evaluate the QRS detector filter for the next sample */ +{ + long dv, dv1, dv2, v0; + + dv2 = sm(t+dt4);/* do this first, to ensure that all of the other + smoothed values needed below are in the buffer */ + dv2 -= s(t-dt4); + dv1 = s(t+dt) - s(t-dt); + dv = (dv1 << 1); + dv -= s(t+dt2) - s(t-dt2); + dv <<= 1; + dv += dv1; + dv -= s(t+dt3) - s(t-dt3); + dv <<= 1; + dv += dv2; + v1 += dv; + v0 = v1 / v1norm; /* scaling is needed to avoid overflow */ + q(t) = v0 * v0; +} + +void addpeak(WFDB_Time t, int peak_amplitude) +{ + struct peak *p = cpeak->next; + + p->time = t; + p->amp = peak_amplitude; + p->type = 0; + cpeak = p; + (p->next)->amp = 0; +} + +/* peaktype() returns 1 if p is the most prominent peak in its neighborhood, 2 + otherwise. The neighborhood consists of all other peaks within rrmin. + Normally, "most prominent" is equivalent to "largest in amplitude", but this + is not always true. For example, consider three consecutive peaks a, b, c + such that a and b share a neighborhood, b and c share a neighborhood, but a + and c do not; and suppose that amp(a) > amp(b) > amp(c). In this case, if + there are no other peaks, a is the most prominent peak in the (a, b) + neighborhood. Since b is thus identified as a non-prominent peak, c becomes + the most prominent peak in the (b, c) neighborhood. This is necessary to + permit detection of low-amplitude beats that closely precede or follow beats + with large secondary peaks (as, for example, in R-on-T PVCs). +*/ + +int peaktype(struct peak *p) +{ + if (p->type) + return (p->type); + else { + int a = p->amp; + struct peak *pp; + WFDB_Time t0 = p->time - rrmin, t1 = p->time + rrmin; + + if (t0 < 0) t0 = 0; + for (pp = p->prev; t0 < pp->time && pp->time < (pp->next)->time; + pp = pp->prev) { + if (pp->amp == 0) break; + if (a < pp->amp && peaktype(pp) == 1) + return (p->type = 2); + } + for (pp = p->next; pp->time < t1 && pp->time > (pp->prev)->time; + pp = pp->next) { + if (pp->amp == 0) break; + if (a < pp->amp && peaktype(pp) == 1) + return (p->type = 2); + } + return (p->type = 1); + } +} + +/* rewind_gqrs resets the sample pointers and annotation fields to their + initial values. */ +void rewind_gqrs(void) +{ + int i; + struct peak *p; + + countdown = -1; + (void)sample(0,t); + annot.time = (WFDB_Time)0; + annot.anntyp = NORMAL; + annot.subtyp = annot.chan = annot.num = 0; + annot.aux = NULL; + for (i = 0, p = peaks; i < NPEAKS; i++, p++) + p->type = p->amp = p->time = 0; +} + +/* gqrs_init() is intended to be called once only, to initialize variables for + the QRS detection function gqrs(). */ + +void gqrs_init(WFDB_Time from, WFDB_Time to) +{ + int i, dv; + static double HR, RR, RRdelta, RRmin, RRmax, QS, QT, RTmin, RTmax, + QRSa, QRSamin; + + /* Allocate workspace. The allocator, gcalloc, is defined below; it + includes an error handler and on-exit deallocation so that this + code can be kept simpler. */ + qfv = (long *)gcalloc((size_t)BUFLN, sizeof(long)); + smv = (long *)gcalloc((size_t)BUFLN, sizeof(long)); + peaks = (struct peak *)gcalloc((size_t)NPEAKS, sizeof(struct peak)); + + /* Gather peak structures into circular buffers */ + for (i = 0; i < NPEAKS; i++) { + peaks[i].next = &peaks[i+1]; + peaks[i].prev = &peaks[i-1]; + } + peaks[0].prev = &peaks[NPEAKS-1]; + cpeak = peaks[NPEAKS-1].next = &peaks[0]; + + /* Read a priori physiologic parameters from the configuration file if + available. They can be adjusted in the configuration file for pediatric, + fetal, or animal ECGs. */ + if (config) { + char buf[256], *p; + + /* Read the configuration file a line at a time. */ + while (fgets(buf, sizeof(buf), config)) { + /* Skip comments (empty lines and lines beginning with `#'). */ + if (buf[0] == '#' || buf[0] == '\n') continue; + /* Set parameters. Each `getconf' below is executed once for + each non-comment line in the configuration file. */ + getconf(HR, "%lf"); + getconf(RR, "%lf"); + getconf(RRdelta, "%lf"); + getconf(RRmin, "%lf"); + getconf(RRmax, "%lf"); + getconf(QS, "%lf"); + getconf(QT, "%lf"); + getconf(RTmin, "%lf"); + getconf(RTmax, "%lf"); + getconf(QRSa, "%lf"); + getconf(QRSamin, "%lf"); + } + fclose(config); + } + + /* If any a priori parameters were not specified in the configuration file, + initialize them here (using values chosen for adult human ECGs). */ + if (HR != 0.0) RR = 60.0/HR; + if (RR == 0.0) RR = 0.8; + if (RRdelta == 0.0) RRdelta = RR/4; + if (RRmin == 0.0) RRmin = RR/4; + if (RRmax == 0.0) RRmax = 3*RR; + if (QS == 0.0) QS = 0.07; + if (QT == 0.0) QT = 5*QS; + if (RTmin == 0.0) RTmin = 3*QS; + if (RTmax == 0.0) RTmax = 5*QS; + if (QRSa == 0.0) QRSa = 750; + if (QRSamin == 0.0) QRSamin = QRSa/5; + + /* Initialize gqrs's adaptive parameters based on the a priori parameters. + During its learning period, gqrs will adjust them based on the observed + input; after the learning period, gqrs continues to adjust these + parameters, but with slower rates of change than during the learning + period. */ + rrmean = RR * sps; + rrdev = RRdelta * sps; + rrmin = RRmin * sps; + rrmax = RRmax * sps; + if ((rrinc = rrmean/40) < 1) rrinc = 1; + if ((dt = QS * sps / 4) < 1) { + dt = 1; + fprintf(stderr, "%s (warning): sampling rate may be too low\n", pname); + } + rtmin = RTmin * sps; /* minimum RTpeak interval */ + rtmean = 0.75 * QT * sps; /* expected RTpeak interval, about 75% of QT */ + rtmax = RTmax * sps; /* maximum RTpeak interval */ + + dv = muvadu(sig, (int)(QRSamin)); + pthr = (thresh * dv * dv) / 6; + qthr = pthr << 1; + pthmin = pthr >> 2; + qthmin = (pthmin << 2)/3; + tamean = qthr; /* initial value for mean T-wave amplitude */ + + /* Filter constants and thresholds. */ + dt2 = 2*dt; + dt3 = 3*dt; + dt4 = 4*dt; + smdt = dt; + v1norm = smdt * dt * 64; + smt = t0; + smt0 = t0 + smdt; + t = t0 - dt4; + for (i = 0; i < BUFLN; i++) + qfv[i] = smv[i] = 0; +} + +/* find_missing() is invoked by gqrs() whenever it is suspected that a + low-amplitude beat may have been missed between two consecutive detected + beats at r and p. The primary peak closest to the expected time of + the missing beat, if any, is returned as the suggested missing beat. */ + +struct peak *find_missing(struct peak *r, struct peak *p) +{ + int rrerr, rrtmp, minrrerr; + struct peak *q, *s = NULL; + + if (r == NULL || p == NULL) return (NULL); + minrrerr = p->time - r->time; + for (q = r->next; q->time < p->time; q = q->next) { + if (peaktype(q) == 1) { + rrtmp = q->time - r->time; + rrerr = rrtmp - rrmean; + if (rrerr < 0) rrerr = -rrerr; + if (rrerr < minrrerr) { + minrrerr = rrerr; + s = q; + } + } + } + return (s); +} + +/* gqrs() is the main QRS detection function. It attempts to find all + beats between the time limits specified by its arguments, and to label + them using annotations of type NORMAL (gqrs() does not attempt to + differentiate normal and ectopic beats). */ + +void gqrs(WFDB_Time from, WFDB_Time to) +{ + int c, i, qamp, q0, q1 = 0, q2 = 0, rr, rrd, rt, rtd, rtdmin; + struct peak *p, *q = NULL, *r = NULL, *tw; + WFDB_Time last_peak = from, last_qrs = from; + + while (t <= to + sps) { + if (countdown < 0 && sample_valid()) + qf(); + else if (countdown < 0) { + countdown = strtim("1"); + state = CLEANUP; + } + else if (countdown-- <= 0) + break; + + q0 = q(t); q1 = q(t-1); q2 = q(t-2); + if (q1 > pthr && q2 < q1 && q1 >= q0 && t > dt4) { + addpeak(t-1, q1); + last_peak = t-1; + for (p = cpeak->next; p->time < t - rtmax; p = p->next) { + if (p->time >= annot.time + rrmin && peaktype(p) == 1) { + if (p->amp > qthr) { + rr = p->time - annot.time; + if (rr > rrmean + 2 * rrdev && + rr > 2 * (rrmean - rrdev) && + (q = find_missing(r, p))) { + p = q; + rr = p->time - annot.time; + annot.subtyp = 1; + } + if ((rrd = rr - rrmean) < 0) rrd = -rrd; + rrdev += (rrd - rrdev) >> 3; + if (rrd > rrinc) rrd = rrinc; + if (rr > rrmean) rrmean += rrd; + else rrmean -= rrd; + if (p->amp > qthr * 4) qthr++; + else if (p->amp < qthr) qthr--; + if (qthr > pthr * 20) qthr = pthr * 20; + last_qrs = p->time; + if (state == RUNNING) { + int qsize; + + annot.time = p->time - dt2; + annot.anntyp = NORMAL; + annot.chan = sig; + qsize = p->amp * 10.0 / qthr; + if (qsize > 127) qsize = 127; + annot.num = qsize; + putann(0, &annot); + annot.time += dt2; + } + /* look for this beat's T-wave */ + tw = NULL; rtdmin = rtmean; + for (q = p->next; q->time > annot.time; q = q->next) { + rt = q->time - annot.time - dt2; + if (rt < rtmin) continue; + if (rt > rtmax) break; + if ((rtd = rt - rtmean) < 0) rtd = -rtd; + if (rtd < rtdmin) { + rtdmin = rtd; + tw = q; + } + } + if (tw) { + static WFDB_Annotation tann; + + tann.time = tw->time - dt2; + if (debug && state == RUNNING) { + tann.anntyp = TWAVE; + tann.chan = sig+1; + tann.num = rtdmin; + tann.subtyp = (tann.time > annot.time + rtmean); + tann.aux = NULL; + putann(0, &tann); + } + rt = tann.time - annot.time; + if ((rtmean += (rt - rtmean) >> 4) > rtmax) + rtmean = rtmax; + else if (rtmean < rtmin) + rtmean = rrmin; + tw->type = 2; /* mark T-wave as secondary */ + } + r = p; q = NULL; qamp = 0; annot.subtyp = 0; + } + else if (t - last_qrs > rrmax && qthr > qthmin) + qthr -= (qthr >> 4); + } + } + } + else if (t - last_peak > rrmax && pthr > pthmin) + pthr -= (pthr >> 4); + + if (++t >= next_minute) { + next_minute += spm; + (void)fprintf(stderr, "."); + (void)fflush(stderr); + if (++minutes >= 60) { + (void)fprintf(stderr, " %s\n", timstr(-t)); + minutes = 0; + } + } + } + + if (state == LEARNING) + return; + + /* Mark the last beat or two. */ + for (p = cpeak->next; p->time < (p->next)->time; p = p->next) { + // if (to > 0 && p->time > to - sps) + // break; + if (p->time >= annot.time + rrmin && p->time < tf && peaktype(p) == 1) { + annot.anntyp = NORMAL; + annot.chan = sig; + annot.time = p->time; + putann(0, &annot); + } + } +} + +/* prog_name() extracts this program's name from argv[0], for use in error and + warning messages. */ + +char *prog_name(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); +} + +/* help() prints a (very) concise summary of how to use this program. + A more detailed summary is in the man page (gqrs.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:", + " -c FILE initialize parameters from the specified configuration FILE", + " -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: 1.00)", + " -o RECORD save filtered signals in a new RECORD", + " -s SIGNAL analyze specified SIGNAL (default: 0)", + " (Note: SIGNAL may be specified by number or name.)", + " -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]); +} + +/* gcalloc() is a wrapper for calloc() that handles errors and maintains + a list of allocated buffers for automatic on-exit deallocation via + gfreeall(). */ + +size_t gcn = 0, gcnmax = 0; +void **gclist = NULL; + +void *gcalloc(size_t nmemb, size_t size) +{ + void *p = calloc(nmemb, size), **q = NULL; + + if ((p == NULL) || + ((gcn >= gcnmax) && + (q = realloc(gclist, (gcnmax += 32)*(sizeof(void *)))) == NULL)) { + fprintf(stderr, "%s: insufficient memory\n", pname); + cleanup(3); + } + if (q) gclist = q; + return (gclist[gcn++] = p); +} + +void gfreeall() /* free memory allocated using gcalloc() */ +{ + while (gcn-- > 0) + if (gclist[gcn]) free(gclist[gcn]); + free(gclist); +} + +void cleanup(int status) /* close files and free allocated memory */ +{ + if (record) wfdbquit(); + gfreeall(); + exit(status); +} diff -Naur --exclude Makefile --exclude info wfdb-10.5.12/app/gqrs.conf wfdb-10.5.13/app/gqrs.conf --- wfdb-10.5.12/app/gqrs.conf 1969-12-31 19:00:00.000000000 -0500 +++ wfdb-10.5.13/app/gqrs.conf 2012-05-06 13:27:36.000000000 -0400 @@ -0,0 +1,79 @@ +# file: gqrs.conf G. Moody 18 October 2007 +# +# Configuration file for `gqrs' and `gqpost' applications + +# To use this file, use the -c option of gqrs and gqpost: +# gqrs -c gqrs.conf ... +# gqpost -c gqrs.conf ... + +# This file contains parameter definitions, which can be provided in any of +# the following forms: +# parameter = value +# parameter: value +# parameter value +# Neither parameter names nor values may contain embedded whitespace. + +# The order of parameter definitions in this file is not significant, except +# that a later definition overrides any earlier definition of the same +# parameter. Two or more parameter definitions may appear on a single line if +# separated by whitespace. Empty lines, and lines beginning with `#', are +# ignored. + +# The parameter values given below have been chosen for use with adult human +# ECGs. In general, pediatric and small animal ECGs have more rapid heart +# rates, and QRS complexes that are shorter and lower in amplitude. + +HR 75 # Typical heart rate, in beats per minute +# RR 0.8 # Typical inter-beat interval, in seconds +# Note that HR overrides RR if both are set -- in this case, RR = 60/HR + +RRdelta 0.2 # Typical difference between successive RR intervals in seconds +RRmin 0.28 # Minimum RR interval ("refractory period"), in seconds +RRmax 2.4 # Maximum RR interval, in seconds; thresholds will be adjusted + # if no peaks are detected within this interval +QS 0.07 # Typical QRS duration, in seconds +QT 0.35 # Typical QT interval, in seconds +RTmin 0.25 # Minimum interval between R and T peaks, in seconds +RTmax 0.33 # Maximum interval between R and T peaks, in seconds + +QRSa 750 # Typical QRS peak-to-peak amplitude, in microvolts +QRSamin 130 # Minimum QRS peak-to-peak amplitude, in microvolts (* see note) + +# For infants, these settings may yield better results than the defaults for +# adults given above: +# HR 120 +# QS 0.05 +# QRSa 500 + +# Settings for other animals + +# By choosing appropriate values for HR, QS, and QRSa, gqrs and gqpost can be +# used to analyze ECGs from a wide variety of animals. In general, gqrs and +# gqpost can adapt as needed if these parameters are within a factor of 2 of +# their ideal values. + +# For mammals, the typical heart rate is about 240 bpm, divided by the fourth +# root of the body mass in kilograms: +# HR ~= 240 bpm / (M^(1/4)) +# (Schmidt-Nielsen K. Animal Physiology: Adaptation and Environment. Cambridge +# University Press, 1997). + +# The typical QRS duration is proportional to the cube root of the mass of +# the heart, although there is evidence that for the largest mammals (elephants +# and whales), the normal QRS duration does not exceed 200 ms, suggesting that +# the His-Purkinjie system may be denser in these species (Meijler et al., +# JACC 1992(Aug):475-479). + +# QRS amplitudes vary considerably depending on electrode type and location; +# set QRSa to match the typical peak-to-peak amplitude for your recordings. + +# The other settings (RR, RRdelta, RRmin, RRmax, QT, RTmax, and QRSamin) may be +# set to zero, in which case gqrs and gqpost determine reasonable values for +# them as needed from HR, QS, and QRSa. + +# * SAFETY NOTE: Applicable standards do not permit ECG monitors intended for +# human use to detect QRS complexes with peak-to-peak amplitudes below 150 +# microvolts, so QRSamin must not be reduced for human ECG monitoring +# applications. This regulatory requirement is intended to prevent adaptive +# detectors such as gqrs from triggering on low-amplitude noise during +# asystole. diff -Naur --exclude Makefile --exclude info wfdb-10.5.12/app/Makefile.tpl wfdb-10.5.13/app/Makefile.tpl --- wfdb-10.5.12/app/Makefile.tpl 2012-04-06 15:04:13.000000000 -0400 +++ wfdb-10.5.13/app/Makefile.tpl 2012-05-07 03:59:57.000000000 -0400 @@ -1,20 +1,23 @@ # file: Makefile.tpl G. Moody 23 May 2000 -# Last revised: 6 April 2012 +# Last revised: 7 May 2012 # This section of the Makefile should not need to be changed. -CFILES = ann2rr.c bxb.c calsig.c ecgeval.c epicmp.c fir.c hrstats.c ihr.c \ - mfilt.c mrgann.c mxm.c nguess.c nst.c plotstm.c pscgen.c pschart.c psfd.c \ - rdann.c rdsamp.c rr2ann.c rxr.c sampfreq.c sigamp.c sigavg.c signame.c \ - signum.c skewedit.c snip.c sortann.c sqrs.c sqrs125.c sumann.c sumstats.c \ - tach.c time2sec.c wabp.c wfdb-config.c wfdbcat.c wfdbcollate.c wfdbdesc.c \ - wfdbmap.c wfdbsignals.c wfdbtime.c wfdbwhich.c wqrs.c wrann.c wrsamp.c xform.c +CFILES = ann2rr.c bxb.c calsig.c ecgeval.c epicmp.c fir.c gqfuse.c gqpost.c \ + gqrs.c hrstats.c ihr.c mfilt.c mrgann.c mxm.c nguess.c nst.c plotstm.c \ + pscgen.c pschart.c psfd.c rdann.c rdsamp.c rr2ann.c rxr.c sampfreq.c sigamp.c \ + sigavg.c signame.c signum.c skewedit.c snip.c sortann.c sqrs.c sqrs125.c \ + sumann.c sumstats.c tach.c time2sec.c wabp.c wfdb-config.c wfdbcat.c \ + wfdbcollate.c wfdbdesc.c wfdbmap.c wfdbsignals.c wfdbtime.c wfdbwhich.c \ + wqrs.c wrann.c wrsamp.c xform.c +CFFILES = gqrs.conf HFILES = signal-colors.h -XFILES = ann2rr bxb calsig ecgeval epicmp fir hrstats ihr \ - mfilt mrgann mxm nguess nst plotstm pscgen pschart psfd \ - rdann rdsamp rr2ann rxr sampfreq sigamp sigavg signame \ - signum skewedit snip sortann sqrs sqrs125 sumann sumstats \ - tach time2sec wabp wfdb-config wfdbcat wfdbcollate wfdbdesc \ - wfdbmap wfdbsignals wfdbtime wfdbwhich wqrs wrann wrsamp xform +XFILES = ann2rr bxb calsig ecgeval epicmp fir gqfuse gqpost \ + gqrs hrstats ihr mfilt mrgann mxm nguess nst plotstm \ + pscgen pschart psfd rdann rdsamp rr2ann rxr sampfreq sigamp \ + sigavg signame signum skewedit snip sortann sqrs sqrs125 \ + sumann sumstats tach time2sec wabp wfdb-config wfdbcat \ + wfdbcollate wfdbdesc wfdbmap wfdbsignals wfdbtime wfdbwhich \ + wqrs wrann wrsamp xform SCRIPTS = cshsetwfdb setwfdb PSFILES = pschart.pro psfd.pro 12lead.pro MFILES = Makefile @@ -67,7 +70,7 @@ # `make listing': print a listing of WFDB applications sources listing: - $(PRINT) README $(MFILES) $(CFILES) $(HFILES) $(PSFILES) + $(PRINT) README $(MFILES) $(CFILES) $(HFILES) $(CFFILES) $(PSFILES) # Rules for compiling applications that require non-standard options diff -Naur --exclude Makefile --exclude info wfdb-10.5.12/doc/wag-src/gqrs.1 wfdb-10.5.13/doc/wag-src/gqrs.1 --- wfdb-10.5.12/doc/wag-src/gqrs.1 1969-12-31 19:00:00.000000000 -0500 +++ wfdb-10.5.13/doc/wag-src/gqrs.1 2012-05-05 19:35:56.000000000 -0400 @@ -0,0 +1,104 @@ +.TH GQRS 1 "5 May 2012" "WFDB 10.5.13" "WFDB Applications Guide" +.SH NAME +gqrs, gqpost \- QRS detector and post-processor +.SH SYNOPSIS +\fBgqrs -r\fR \fIrecord\fR [ \fIoptions\fR ... ] +.br +\fBgqpost -r\fR \fIrecord\fR [ \fIoptions\fR ... ] +.SH DESCRIPTION +.PP +\fBgqrs\fR attempts to locate QRS complexes in an ECG signal in the +specified \fIrecord\fR. The detector algorithm is new and as yet +unpublished. The output of \fBgqrs\fR is an annotation file (with +annotator name \fBqrs\fR) in which all detected beats are labelled +normal ("N"). +.PP +As a QRS detector for research, \fBgqrs\fR has been optimized for sensitivity. +\fBgqpost\fR can post-process \fBgqrs\fR's output annotation file to improve +positive predictivity, at a cost of reduced sensitivity. It does this by +copying its input annotation file, changing N annotations into artifact +("!") annotations if they are likely to be erroneous. +.PP +A shared configuration file can be used to describe some of the expected +characteristics of the ECG signal. This is unnecessary when processing +adult human ECGs, but an appropriately constructed configuration file allows +\fBgqrs\fR to analyze fetal, pediatric, and animal ECGs. +.PP +\fIOptions\fR include: +.TP +\fB-a\fR \fIannotator\fR +[gqpost only] Read annotations from the specified \fIannotator\fR (default: +qrs). +.TP +\fB-c\fR \fIfile\fR +Initialize parameters based on the specified (text) configuration +\fIfile\fR. See the example configuration file, \fIgqrs.conf\fR, for details. +.TP +\fB-f\fR \fItime\fR +Begin at the specified \fItime\fR in \fIrecord\fR (default: the beginning of +\fIrecord\fR). +.TP +\fB-h\fR +Print a usage summary. +.TP +\fB-H\fR +Read the signal files in high-resolution mode (default: standard mode). +.TP +\fB-m\fR \fIthreshold\fR +Specify the \fIthreshold\fR (default: 1.0) for detection [qqrs] or +acceptance [gqpost]. Use higher values to reduce false detections, or lower +values to reduce the number of missed beats. +.TP +\fB-n\fR \fIname\fR +[gqrs only] Save the filtered signals in a new record with the specified record +\fIname\fR. +.TP +\fB-o\fR \fIname\fR +[gqpost only] write annotations to an annotation file with the specified +annotator \fIname\fR. +.TP +\fB-s\fR \fIsignal\fR +[gqrs only] Specify the \fIsignal\fR to be used for QRS detection (default: 0). Note that signals may be specified by number or name. +.TP +\fB-t\fR \fItime\fR +Process until the specified \fItime\fR in \fIrecord\fR (default: the end of the +\fIrecord\fR). +.PP +Note that \fBgqpost\fR always copies its entire input annotation file. The +\fB-f\fR and \fB-t\fR options, if present, only define the interval during which +\fBgqpost\fR may change annotations. Since \fBgqpost\fR can reprocess its own +output, this feature allows multiple passes using different threshold values +and processing intervals, if necessary. +.SH ENVIRONMENT +.PP +It may be necessary to set and export the shell variable \fBWFDB\fR (see +\fBsetwfdb\fR(1)). +.SH EXAMPLES +.PP +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: +.br + \fBgqrs -r 100 -f 5:0 -t 10:35 -s 1\fR +.br +The output annotations may be read using (for example): +.br + \fBrdann -a qrs -r 100\fR +.PP +To evaluate the performance of this program, run it on the entire record, by: +.br + \fBgqrs -r 100\fR +.br +and then compare its output with the reference annotations by: +.br + \fBbxb -r 100 -a atr qrs\fR +.SH SEE ALSO +\fBbxb\fR(1), \fBecgpuwave\fR(1), \fBrdann\fR(1), \fBsetwfdb\fR(1), +\fBsqrs\fR(1), \fBwqrs\fR(1), \fBxform\fR(1) +.SH AUTHORS +George B. Moody (george@mit.edu). +.SH SOURCES +http://www.physionet.org/physiotools/wfdb/app/gqrs.c +.br +http://www.physionet.org/physiotools/wfdb/app/gqpost.c +.br +http://www.physionet.org/physiotools/wfdb/app/gqrs.conf diff -Naur --exclude Makefile --exclude info wfdb-10.5.12/doc/wag-src/header.5 wfdb-10.5.13/doc/wag-src/header.5 --- wfdb-10.5.12/doc/wag-src/header.5 2010-06-29 10:33:47.000000000 -0400 +++ wfdb-10.5.13/doc/wag-src/header.5 2012-05-04 12:30:19.000000000 -0400 @@ -287,7 +287,6 @@ the initial `#' comment character) are referred to as `info strings'. There must be no whitespace preceding the initial `#' in any line that is to be recognized by \fIgetinfo\fP. -form ``record \fIxxx\fP, signal \fIn\fP''. .SS "Multi-segment records" .LP Each non-empty, non-comment line following the record line in the diff -Naur --exclude Makefile --exclude info wfdb-10.5.12/doc/wpg-src/wpg0.tex wfdb-10.5.13/doc/wpg-src/wpg0.tex --- wfdb-10.5.12/doc/wpg-src/wpg0.tex 2012-04-25 08:45:16.000000000 -0400 +++ wfdb-10.5.13/doc/wpg-src/wpg0.tex 2012-05-13 15:10:53.000000000 -0400 @@ -9049,6 +9049,14 @@ @unnumberedsec WFDB 10.5 +@unnumberedsubsec Changes in version 10.5.13 (13 May 2012) + +Versions of the WFDB library up to 10.5.10 ignored embedded empty lines +within the @code{info} sections of @file{.hea} files, but versions 10.5.10 and +10.5.11 treat them as markers of the end of the @code{info} section. This +version restores the previous treatment of embedded empty lines. Thanks +to Justin Leo Cheang Loong for reporting this issue. + @unnumberedsubsec Changes in version 10.5.12 (25 April 2012) When called with a NULL argument, getinfo() sometimes behaves differently in diff -Naur --exclude Makefile --exclude info wfdb-10.5.12/lib/signal.c wfdb-10.5.13/lib/signal.c --- wfdb-10.5.12/lib/signal.c 2012-04-24 16:20:19.000000000 -0400 +++ wfdb-10.5.13/lib/signal.c 2012-05-04 13:24:59.000000000 -0400 @@ -1,5 +1,5 @@ /* file: signal.c G. Moody 13 April 1989 - Last revised: 24 April 2012 wfdblib 10.5.12 + Last revised: 4 May 2012 wfdblib 10.5.13 WFDB library functions for signals _______________________________________________________________________________ @@ -3070,16 +3070,20 @@ if (*buf != '#') break; /* skip initial comments, if any */ while (wfdb_fgets(buf, 256, ifile)) if (*buf == '#') break; /* skip header content */ - while (*buf == '#') { /* read and save info */ - p = buf + strlen(buf) - 1; - if (*p == '\n') *p-- = '\0'; - if (*p == '\r') *p = '\0'; - if (ninfo >= nimax) { - nimax += 16; - SREALLOC(pinfo, nimax, sizeof(char *)); + while (*buf) { /* read and save info */ + if (*buf == '#') { /* skip anything that isn't info */ + p = buf + strlen(buf) - 1; + if (*p == '\n') *p-- = '\0'; + if (*p == '\r') *p = '\0'; + if (ninfo >= nimax) { + int j = nimax; + nimax += 16; + SREALLOC(pinfo, nimax, sizeof(char *)); + memset(pinfo + j, 0, (size_t)(16*sizeof(char *))); + } + SSTRCPY(pinfo[ninfo], buf+1); + ninfo++; } - SSTRCPY(pinfo[ninfo], buf+1); - ninfo++; if (wfdb_fgets(buf, 256, ifile) == NULL) break; } wfdb_fclose(ifile); @@ -3093,8 +3097,10 @@ if (*p == '\n') *p-- = '\0'; if (*p == '\r') *p = '\0'; if (ninfo >= nimax) { + int j = nimax; nimax += 16; SREALLOC(pinfo, nimax, sizeof(char *)); + memset(pinfo + j, 0, (size_t)(16*sizeof(char *))); } SSTRCPY(pinfo[ninfo], buf+1); ninfo++; diff -Naur --exclude Makefile --exclude info wfdb-10.5.12/lib/wfdb.h wfdb-10.5.13/lib/wfdb.h --- wfdb-10.5.12/lib/wfdb.h 2012-04-25 08:19:35.000000000 -0400 +++ wfdb-10.5.13/lib/wfdb.h 2012-05-13 15:11:24.000000000 -0400 @@ -32,7 +32,7 @@ /* WFDB library version. */ #define WFDB_MAJOR 10 #define WFDB_MINOR 5 -#define WFDB_RELEASE 12 +#define WFDB_RELEASE 13 #define WFDB_NETFILES 1 /* if 1, library includes code for HTTP, FTP clients */ #define WFDB_NETFILES_LIBCURL 1 diff -Naur --exclude Makefile --exclude info wfdb-10.5.12/MANIFEST wfdb-10.5.13/MANIFEST --- wfdb-10.5.12/MANIFEST 2012-04-06 15:44:35.000000000 -0400 +++ wfdb-10.5.13/MANIFEST 2012-05-07 04:02:50.000000000 -0400 @@ -7,6 +7,10 @@ app/ecgeval.c app/epicmp.c app/fir.c +app/gqfuse.c +app/gqpost.c +app/gqrs.c +app/gqrs.conf app/hrstats.c app/ihr.c app/Makefile @@ -283,6 +287,7 @@ doc/wag-src/fixag.sed doc/wag-src/fixag.sh doc/wag-src/getpagenos.c +doc/wag-src/gqrs.1 doc/wag-src/header.5 doc/wag-src/hrfft.1 doc/wag-src/hrstats.1 diff -Naur --exclude Makefile --exclude info wfdb-10.5.12/NEWS wfdb-10.5.13/NEWS --- wfdb-10.5.12/NEWS 2012-04-25 08:54:40.548019711 -0400 +++ wfdb-10.5.13/NEWS 2012-05-13 16:03:06.861389753 -0400 @@ -1,3 +1,16 @@ +10.5.13 (13 May 2012): + A flexible high-sensitivity QRS detector for research (app/gqrs.c), and + a companion post-processor (app/gqpost.c) for improving its positive + predictivity with minimal loss of sensitivity, are included in this + release. The algorithms employed in these applications are previously + unpublished. + + Versions of the WFDB library up to 10.5.10 ignored embedded empty lines + within the 'info' sections of '.hea' files, but versions 10.5.10 and + 10.5.11 treat them as markers of the end of the 'info' section. This + version restores the previous treatment of embedded empty lines. Thanks + to Justin Leo Cheang Loong for reporting this issue. + 10.5.12 (24 April 2012): When called with a NULL argument, getinfo() sometimes behaves