/* file: pnnlist.c Joe E. Mietus 19 November 2002
Last revised: 23 February 2003
-------------------------------------------------------------------------------
pnnlist: derive pNNx statistics from an annotation interval list
Copyright (C) 2003 Joe E. Mietus
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 (joe at physionet dot org) or postal mail
(BIDMC Room KB-28, 300 Brookline Ave., Boston, MA 02215 USA). For updates to
this software, please visit PhysioNet (http://www.physionet.org/).
_______________________________________________________________________________
A detailed description of the application of this algorithm can be
found in: JE Mietus, C-K Peng, I Henry, R L Goldsmith and AL Goldberger.
"The pNNx files: re-examining a widely used heart rate variability measure."
Heart 2002; vol. 88, pp. 378-380.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define ABS(A) ((A) < 0 ? -(A) : (A))
#define ROFFERR 1e-9
char *pname, *prog_name();
main(argc, argv)
int argc;
char *argv[];
{
char ann, lastann[2];
int pflag, sflag, dblcmp();
double rr, lastrr, inc, nninc, mult, *pnn = NULL;
long i, maxdat, sum, tot, totn, totp;
void bsd_qsort(), help();
pname = prog_name(argv[0]);
inc = pflag = sflag = 0;
mult = 1000;
for (i = 0; ++i < argc && *argv[i] == '-'; ) {
switch(argv[i][1]) {
case 'i':
if (i < argc-1)
inc = atof(argv[++i]);
else {
fprintf(stderr, "%s: increment must follow -i\n", pname);
exit(1);
}
break;
case 'p':
pflag = 1;
mult = 100;
break;
case 's':
sflag = 1;
break;
default:
help();
exit(1);
}
}
if (inc < 0) {
fprintf(stderr, "%s -i inc : increment must be greater than 0\n",
pname);
exit(1);
}
else if (inc > 0) {
if (pflag)
inc /= 100;
else
inc /= 1000;
}
if (scanf("%lf %c", &rr, &ann) != 2) {
fprintf(stderr, "%s : improperly formatted data\n\n", pname);
help();
exit(2);
}
i = maxdat = 0;
lastann[0] = 'X';
lastann[1] = ann;
while (scanf("%lf %c", &rr, &ann) == 2) {
if (lastann[0] == 'N' && lastann[1] == 'N' && ann == 'N') {
if (i >= maxdat) {
double *s;
maxdat += 50000; /* allow the pnn array to grow (the
increment is arbitrary) */
if ((s = realloc(pnn, maxdat * sizeof(double))) == NULL) {
fprintf(stderr,
"%s: truncating input after %ld increments\n",
pname, i);
break;
}
pnn = s;
}
if (pflag)
pnn[i] = (rr-lastrr)/lastrr;
else
pnn[i] = rr-lastrr;
if (!sflag)
pnn[i] = ABS(pnn[i]);
i++;
}
lastrr = rr;
lastann[0] = lastann[1];
lastann[1] = ann;
}
tot = i;
if (tot <= 0) {
fprintf(stderr, "%s: no consecutive NN intervals in input\n", pname);
exit(3);
}
bsd_qsort(pnn, tot, sizeof(double), dblcmp);
if (sflag) { /* signed distributions */
for (totn = totp = i = 0; i < tot; i++) {
if (pnn[i] <= 0)
totn++;
if (pnn[i] >= 0)
totp++;
}
if (inc == 0) {
printf("%g\t0\n", mult*pnn[0]);
for (sum = i = 1; i < tot; i++) {
if (pnn[i] < 0) {
if (pnn[i] - pnn[i-1] > ROFFERR)
printf("%g\t%g\n", mult*pnn[i], 100*(double)sum/totn);
sum++;
}
else if (pnn[i] == 0 && pnn[i-1] != 0) {
printf("%g\t%g\n", mult*pnn[i], 100*(double)sum/totn);
sum = 1;
}
else if (pnn[i] == 0 && pnn[i-1] == 0)
sum++;
else {
if (pnn[i] - pnn[i-1] > ROFFERR)
printf("%g\t%g\n", mult*pnn[i-1],
100*(1-(double)sum/totp));
sum++;
}
}
printf("%g\t%g\n", mult*pnn[i-1], 100*(1-(double)sum/totp));
}
else { /* inc > 0 */
for (nninc = 0; nninc > pnn[0]; nninc -= inc)
;
printf("%g\t0\n", mult*nninc);
nninc += inc;
for (sum = i = 0; i < tot; i++) {
if (pnn[i] < 0) {
if (pnn[i] - nninc > ROFFERR) {
printf("%g\t%g\n", mult*nninc, 100*(double)sum/totn);
nninc += inc;
while (pnn[i] - nninc > ROFFERR) {
printf("%g\t%g\n", mult*nninc,
100*(double)sum/totn);
nninc += inc;
}
}
sum++;
}
else if (pnn[i] == 0 && pnn[i-1] != 0) {
printf("%g\t%g\n", mult*nninc, 100*(double)sum/totn);
sum = 1;
}
else if (pnn[i] == 0 && pnn[i-1] == 0)
sum++;
else {
if (pnn[i] - nninc > ROFFERR) {
printf("%g\t%g\n", mult*nninc,
100*(1-(double)sum/totp));
nninc += inc;
while (pnn[i] - nninc > ROFFERR) {
printf("%g\t%g\n", mult*nninc,
100*(1-(double)sum/totp));
nninc += inc;
}
}
sum++;
}
}
printf("%g\t%g\n", mult*nninc, 100*(1-(double)sum/totp));
}
}
else { /* unsigned distributions */
if (inc == 0) {
for (sum = i = 1; i < tot; i++) {
if (pnn[i] - pnn[i-1] > ROFFERR)
printf("%g\t%g\n", mult*pnn[i-1], 100*(1-(double)sum/tot));
sum++;
}
printf("%g\t%g\n", mult*pnn[i-1], 100*(1-(double)sum/tot));
}
else { /* inc > 0 */
for (nninc = 0, sum = i = 0; i < tot; i++) {
if (pnn[i] - nninc > ROFFERR) {
printf("%g\t%g\n", mult*nninc, 100*(1-(double)sum/tot));
nninc += inc;
while (pnn[i] - nninc > ROFFERR) {
printf("%g\t%g\n", mult*nninc,
100*(1-(double)sum/tot));
nninc += inc;
}
}
sum++;
}
printf("%g\t%g\n", mult*nninc, 100*(1-(double)sum/tot));
}
}
if (pnn) free(pnn);
exit(0);
}
int dblcmp(y1, y2)
double *y1, *y2;
{
if (*y1 < *y2)
return(-1);
else if (*y1 > *y2)
return(1);
else
return(0);
}
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 [OPTIONS ...]\n",
" Takes as standard input an annotation interval list, containing intervals",
" in seconds and the (beat and non-beat) annotations that terminate each",
" interval; and outputs on standard output each unique NN increment (x) in",
" milliseconds, and the percentage of NN interval increments (pNNx) greater",
" than x.",
" ",
"OPTIONS may include:",
" -h print this usage summary",
" -i INC compute and output pNNx for x = 0, INC, 2*INC, ... milliseconds",
" -p compute and output increments as percentage of initial intervals",
" -s compute and output separate distributions of positive and negative",
" intervals",
" ",
"For further information, see http://www.physionet.org/physiotools/pNNx/ .",
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]);
}