/* file: parsescp.c	G. Moody and E. Moody	 10 January 2000
			Last revised:		 11 December 2017
-------------------------------------------------------------------------------
parsescp: parse an SCP-ECG file (read from the standard input)
Copyright (C) 2000-2014 George B. Moody and Edna S. 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/).
_______________________________________________________________________________

'parsescp' converts SCP-ECG files produced by SpaceLabs/Burdick ECG carts into
more easily usable formats.  It was written in 2000 and has been used to convert
about a million ECGs collected at Boston's Beth Israel Deaconess Medical Center
since then.

'parsescp' can also be used to create deidentified SCP-ECG files, although
it does not perform this function by default;  see the comments below.

Compiling and installing 'parsescp'
===================================

Prerequisites:
  * gcc (the free GNU C Compiler, available for all popular platforms;
      other ANSI/ISO C compilers should also work but you are on your own)
  * parsescp.c
  * (optional) the WFDB library (included in the free WFDB software package,
      from http://physionet.org/physiotools/wfdb.shtml)

To compile 'parsescp' using the GNU C compiler (gcc), use one of these two
methods:
 (1) If you have installed the WFDB Software Package, use this command:
        gcc -o parsescp -O parsescp.c -lwfdb
 (2) Otherwise, use this command:
        gcc -o parsescp -DNOWFDB -O parsescp.c

To install 'parsescp', compile it first, then copy the executable 'parsescp'
generated by gcc into a directory in your PATH.

Using 'parsescp' to convert an SCP-ECG file
===========================================

Run 'parsescp' within a directory containing the SCP-ECG file to be converted,
which can have any name.  For example, if the file is named 'RECORD.scp', use
a command such as:
       parsescp -o RECORD <RECORD.scp
The argument following '-o' doesn't need to match the name of the input file
as in this example, but this choice may reduce opportunities for confusion.

A variety of options that affect the output may also be given in the command
line:
   -a    anonymize and copy stdin to stdout (suppresses other output)
   -f    force parsing of a file with CRC errors
   -h    print brief help
   -l    low-pass filter (smooth) the output waveforms
   -v    option to print a (very) detailed analysis of RECORD.scp
   -w    produce a WFDB-format record
   -z    zero-mean the input

Additional options for debugging only:
   -b    show baseline (residuals);  suppress templates
   -s N  shift templates by N samples before adding
   -S N  skip parsing of section N
   -t    show templates; suppress baseline
   -x    show hex dumps (implies -v)

Unless the -a option is used, this program produces at least these three files:
  RECORD.des    (text) description (age, sex, recording bandwidth,
                       measurements, diagnoses)
  RECORD.ecg    (binary) reconstructed ECGs
  RECORD.key    (text) patient's name and ID (medical record number)

If invoked with the -v option, parsescp also produces:
  RECORD.txt    (text) reconstructed ECGs
With -v, parsescp also writes a (very) detailed analysis of the contents of the
SCP-ECG input on the standard output.

If parsescp was compiled with the WFDB library, and if it is invoked with the
-w option, it also produces a pair of PhysioBank-compatible output files:
  RECORD.dat    (binary) signal file containing 12 continuous leads
  RECORD.hea    (text) header file describing RECORD.dat

Supported SCP versions
======================

This program was written using AAMI SCP-1999 (Standard communications
protocol for computer-assisted electrocardiography, 25 October 1999 draft)
as a reference for SCP format.  It has been tested only with SCP records
produced by SpaceLabs/Burdick ECG carts (these produce second-difference
encoded data with reference beat subtraction using a single reference beat,
Huffman encoded using the SCP standard Huffman table).  Amplitude (unencoded)
data and first-difference encoded data should be readable using this program,
but these formats have not been tested.  Use of custom Huffman tables is
recognized but not otherwise supported.  Use of multiple reference beats is
recognized but not otherwise supported.

ECG signals in Spacelabs/Burdick SCP-ECG files 
==============================================

Spacelabs/Burdick ECG carts of the type for which this program was designed
record 2 of the 3 Einthoven leads and all 6 precordial leads simultaneously for
10 seconds, at 500 samples per second per lead, with 16-bit precision over a
range of +/-32.767 mV.  Thus the sampling interval is 2 ms, and the amplitude
resolution is 5 microvolts (5000 nanovolts) per ADC unit.

Note that although the SCP-ECG standard specifies how to record the sampling
frequency and amplitude resolution in SCP-ECG files, the Spacelabs/Burdick
carts don't do this, so parsescp assumes the sampling frequency and resolution
above.  parsescp will need modification in order to convert SCP-ECGs with other
sampling frequencies or resolutions correctly.

ECG signals in parsescp's output files
======================================

This program derives the third Einthoven lead and the three augmented leads
using the standard relationships among the leads:
    III = II - I
    aVR = -(I + II)/2
    aVL = II/2 - III
    aVF = I/2 + III

In all of its output formats, parsescp represents the samples of each signal as
a sequence of unscaled integers, exactly as they appear in the original SCP-ECG
input file.  Thus, in the .ecg, .txt, and .dat output files, the unit of
amplitude is equivalent to 5 microvolts (5000 nanovolts), as in the SCP input.
If the recording is shorter than 10 seconds, or if a signal is missing and
cannot be reconstructed from the relationships above, each missing sample is
assigned a special value (WFDB_INVALID_SAMPLE, or -32768).

The -l and -z options modify the input values as noted above;  if neither
option is used, the output sample values are numerically identical to the
input sample values.

The .ecg file contains selected and rearranged segments of the signals in the
commonly-used layout of twelve 2.5 second segments arranged in groups of 4 above
a continuous 10-second lead II.  Each sample is represented as a big-endian
16-bit two's complement signed integer.  The file begins with a 512-byte prolog
containing the record name and recording date and time, which are HIPAA-defined
protected health information (PHI) unless the input SCP-ECG has been
deidentified.  The prolog is followed by four "traces", each representing the
same 10-second interval.  The first three of these traces are made by
concatenating 2.5 second segments (1250 samples) of each of the 12 leads, in
this order:
	(  I aVR V1 V4)
	( II aVL V2 V5)
	(III aVF V3 V6)
The fourth trace is a continuous 10-second segment (5000 samples) of lead II.

The optional .txt and .dat files contain the ECG signals only (no metadata, and
no PHI).  The signals appear in the standard order:  
    I, II, III, aVR, aVF, aVL, V1, V2, V3, V4, V5, V6
The .txt and .dat files contain the full 10 seconds (5000 samples at 500 samples
per second) for each lead (i.e., 12*5000 = 60000 samples in all).  

In the .txt file, each line begins with a sample number (0 to 4999)
and is followed by a sample from each of the 12 leads, in order.  Each
sample is represented as a base 10 numeral, with spaces inserted
between samples so that the columns line up.  Thus the sample numbers
are in column 0, samples of lead I are in column 1, those of lead II
are in column 2, etc.

In the .dat file, the first 24 bytes contain the first sample of each signal,
in the standard order as for the .txt file.  As in the .ecg files, each sample
is represented as a big-endian 16-bit two's complement signed integer.  The
next 24 bytes contain the second sample of each signal, etc.

Other output files
==================

The .des file contains a variety of information extracted from the SCP-ECG
input file, in human-readable form.  It does not contain the ECG signals
themselves, or the patient's name or medical record number.  Note that .des
files made from SCP-ECG files that have not been anonymized will generally
contain HIPAA-defined PHI such as the recording date and the patient's age
(even if over 90).

The .key file contains PHI: the recording date and time, the patient's name,
and the medical record number, if recorded in the input file.

The .hea file, if generated, contains metadata (information about the
corresponding .dat file) only;  it does not contain any PHI, even if the
input was not anonymized.  Age and sex are recorded in the .hea file if
present in the input file, except that ages of 90 and more are recorded
as 90.  The recording date and time are not recorded in the .hea file.

Using parsescp to create deidentified SCP-ECG files
===================================================

The SCP-ECG standard defines how to record a variety of information
that includes elements defined by HIPAA as PHI (protected health information).
These include the patient's name, medical record number, birth day and month,
recording day and month, and (if the age is over 90) birth year and age.

If invoked with the -a option, parsescp reads the input SCP-ECG file and
writes an anonymized (deidentified) version of it to the standard output.
For example:
        parsescp -a <12345678.scp >anonymous.scp
In this case, none of the other output files are produced.

parsescp removes all of the PHI as well as names of physicians and technicians,
names of hospitals or clinics, and room numbers, replacing them with 'xxx'.
It changes all dates to January 1, and if the age is over 90, it resets the
age to 90 and the birth year to 90 years before the recording year.  Finally,
it recalculates the SCP-ECG CRCs so that the output is still a valid
SCP-ECG file.  Note that the original input file is not modified.

Note that parsescp does not deidentify other types of data (including its
own .des and .key files);  it can only deidentify SCP-ECG files.
_______________________________________________________________________________

21 Sep 2000 --  added referring dr to .des file  esm
24 Jan 2001 --  added comments to .des files  esm
25 Jan 2001 --  added -s, -b, -t, -x debugging options  gbm
 9 Dec 2002 --  added GPL notice  gbm
28 Oct 2004 --  added -w option to create WFDB record  gbm
 4 Nov 2004 --  added -a option to create an anonymized copy of the input  gbm
19 Jun 2007 --  added -S debugging option for parsing defective inputs  gbm
21 Nov 2008 --  added -z option to zero-mean the signals  gbm
11 Dec 2008 --  added Institution # and Department #  esm
19 May 2009 --  added ncompare to zflag  esm
 4 May 2010 --  added check for missing comments in section 1  gbm
24 Jun 2010 --  added checks for missing patid and refdr in section 1  gbm
21 Feb 2011 --  added age (if under 90) and sex to .hea file  gbm
 8 Mar 2014 --  corrected comments above  gbm
*/

#include <stdio.h>

#ifndef NOWFDB
#include <wfdb/wfdb.h>
#endif
#ifdef _WINDOWS
#include <io.h>
#include <fcntl.h>
#endif

#define VERSION_STRING	"#ECG/1.0"
#define min(A, B)	(((A) <= (B)) ? (A) : (B))

/* Define SWAP if this program is compiled for a little-endian CPU (e.g., HP)
   but don't do so on a big-endian CPU (e.g., Intel x86) */
#ifndef SWAP
#define fwriteswapped	fwrite
#else
size_t fwriteswapped(void *p, size_t size, size_t nmemb, FILE *stream)
{
    char *q;
    size_t n = nmemb*size;

    if ((q = (char *)malloc(n)) == NULL) {
	fprintf(stderr,
		"insufficient memory (%d bytes needed) in fwriteswapped\n",
		n);
	return (0);
    }
    swab(p, q, n);   /* Copy n bytes from *p to *q, swapping adjacent bytes
			swab() is standard in System V and BSD Unix but is
			not in the ANSI C library, so you may need to provide
			your own implementation if this program is ported to
			a non-Unix environment. */
    n = fwrite(q, size, nmemb, stream);
    free(q);
    return (n);
}
#endif

unsigned short get16(char *p)
{
    unsigned short result;

    result  = *p++ & 0xff;
    result |= *p << 8;
    return (result);
}

void put16(unsigned short val, char *p)
{
    *p++ = val & 0xff;
    *p = val >> 8;
}

unsigned long get32(char *p)
{
    unsigned long result;

    result  =  *p++ & 0xff;
    result |= (*p++ & 0xff) << 8;
    result |= (*p++ & 0xff) << 16;
    result |=  *p << 24;
    return (result);
}

int bflag = 0;		/* 0: normal operation, 1: show baselines only */
int fflag = 0;		/* 0: quit on CRC error, 1: continue if possible */
int tflag = 0;		/* 0: normal operation, 1: show templates only */
int vflag = 0;		/* 0: run quietly, 1: print (lots of) messages */
int xflag = 0;		/* 0: run quietly, 1: print hexdump */

int shift = 0;		/* 0: normal operation; otherwise, shift templates
			   by this many samples before adding to residuals */

unsigned short getcrc(char *p, long length)
{
    unsigned char a, b, crchigh = 0xff, crclow = 0xff;
    unsigned short crc;

    while (length-- > 0L) {
	crchigh = a = *p++ ^ crchigh;
	a >>= 4;		/* top 4 bits of a are now zero */
	a ^= crchigh;		

	/* swap crchigh, crclow */ 
	crchigh = crclow;
	crclow = a;

	/* rotate a left by four bits (swap nibbles) */
	b = (a & 0xf) << 4;	
 	a >>= 4;
	b = a |= b;
	
	/* rotate a left by one more bit */
	if (a & 0x80) a = (a << 1) | 1;
	else a <<= 1;

	crchigh ^= a &= 0x1f;
	a = b & 0xf0;
	crchigh ^= a;		/* crchigh complete */

	/* rotate b left by one bit */
	if (b & 0x80) b = (b << 1) | 1;
	else b <<= 1;

	crclow ^= b &= 0xe0;	/* crclow complete */
    }
    crc = (crchigh << 8) | crclow;
    return (crc);
}

int crcok(char *p, long length, unsigned short crcref, char *description)
{
    unsigned short crc;

    crc = getcrc(p, length);

    if (crc == crcref) {
	if (vflag) printf("%s: correct (%d)\n", description, crc);
	return (1);
    }
    else {
	if (vflag) printf("%s: error (calculated: %d, expected: %d)\n",
	       description, crc, crcref);
	if (fflag)
	    return (0);
	else
	    exit(1);	/* quit on CRC error unless in -f (force) mode */
    }
}

char *pname;	   /* the name of this program (for use in error messages) */
int nleads = -1;	/* number of ECG leads */
int nqrs = -1;		/* number of beats detected */
unsigned long *nsamp;	/* nsamp[i]: number of samples in lead i */
long rblen;		/* number of samples in refbeat (template) */
long rblenms;		/* length of refbeat in milliseconds */
int fcM;		/* index of fiducial within refbeat */
double gref, gres;   /* gains for refbeat and for residuals (nanovolt/unit) */
double sfreq = -1.;	/* sampling frequency (samples/second/lead) */
char recdate[20];	/* date of recording */
char rectime[20];	/* time of recording */
char recname[12];	/* record name */
char patient_name[100];	/* patient's name (LAST, FIRST) */
char referring_dr[100];	/* referring dr */
char comments[200];	/* comments */
char patient_id[30];	/* patient ID (medical record number) */
int age = -1;		/* patient's age */
int sex = -1;		/* 1: male, 2: female, 0: unknown, 9: unspecified */
int RRint = -1;         /* RR interval */
int Pon = -1;           /* P onset */
int Pend = -1;          /* P end */
int QRSon = -1;         /* QRS onset */
int QRSend = -1;        /* QRS end */
int Tend = -1;          /* T end */
short Paxis = 999;         /* P axis */
short QRSaxis = 999;       /* QRS axis */
short Taxis = 999;         /* T axis */
int blfilt = -1;	/* baseline filter -3 dB point */
int lpfilt = -1;	/* low-pass filter -3 dB point */
int filter_bit = -1;	/* see SCP spec */
int tdevid = -1;	/* ID number of acquiring device */
int tinst = -1;         /* institution number */
int tdept = -1;         /* department number */
char *dept = "<missing>";/* location description or code */
char *inst = "<missing>";/* institution description */
int stmnts = -1;        /* number of statements in report */
char *report[1200];      /* report */
int lflag = 0;		/* if non-zero, this program should low-pass filter
			   its output */
int wflag = 0;		/* if non-zero, this program should produce a WFDB
			   format record */
int zflag = 0;		/* if non-zero, remove final transients and zero-mean
			   the signals */
struct subzone {	/* reference beat subtraction zone */
    short type;		/* beat type (0: dominant) */
    long t0;		/* first sample in zone */
    long t1;		/* fiducial */
    long t2;		/* last sample in zone */
} *subz;

short **ecg;		/* ecg[i]: reconstructed ECG for ith lead */
short **refbeat;	/* reference beat (template) */

/* Indices of the leads within ecg[];  e.g. ecg[aVR][] is the array of samples
   of lead aVR, etc. */
short leadI = -1, leadII = -1, leadIII = -1, aVR = -1, aVF = -1, aVL = -1,
    V1 = -1, V2 = -1, V3 = -1, V4 = -1, V5 = -1, V6 = -1;

struct ECGmeas {	/* sets of measurements in section 10 */
    short leadid;	/* Lead ID (index into leadname[]) */
    short vlen;		/* length of the rest of the data for this lead */
    short Pdur;		/* P duration */
    short PRint;	/* PR interval */
    short QRSdur;	/* QRS duration */ 
    short QTint;	/* QT interval */
    short Qdur; 	/* Q duration */
    short Rdur;		/* R duration */
    short Sdur;		/* S duration */
    short Rpdur;	/* R' duration */
    short Spdur;	/* S' duration */
    short Qamp; 	/* Q amplitude */
    short Ramp;		/* R amplitude */
    short Samp;		/* S amplitude */
    short Rpamp;	/* R' amplitude */
    short Spamp;	/* S' amplitude */
    short Jamp;		/* J-point amplitude */
    short Ppamp;	/* P+ amplitude */
    short Pnamp;	/* P- amplitude */
    short Tpamp;	/* T+ amplitude */
    short Tnamp;	/* T- amplitude */
    short STslope;	/* ST slope */
    short Pmorph;	/* P morphology */
    short Tmorph;	/* T morphology */
    short isoI;		/* isoelectric segment at QRS onset */
    short isoK;		/* isoelectric segment at QRS end */
    short tind;		/* time of intrinsicoid deflection */
    short qual;		/* recording quality */
    short st20;		/* ST amplitude at J + 20 ms */
    short st60;		/* ST amplitude at J + 60 ms */
    short st80;		/* ST amplitude at J + 80 ms */
    short strr16;	/* ST amplitude at J + RRmean/16 */
    short strr8;	/* ST amplitude at J + RRmean/8 */
} *ecgmeas;

/* This low-pass filter is implemented as a 5-tap weighted moving average
   FIR filter.  The -3 dB point is at roughly 1/10 the sampling frequency;
   hence it behaves as a 50 Hz low-pass filter given 500 samples/second
   input.

   This doesn't really belong in this program, but it's easily and efficiently
   implemented here.  It exists so that users who prefer to look at low-pass
   filtered output can do so, by using the -l option.  Note that the filter
   is applied only to the .ecg file, not to the text-format samples. */

void lowpass()
{
    int i, j;
    short *in, *out;

    for (i = 0; i < nleads; i++) {
	if (ecg[i] == NULL) continue;
	if ((out = (short *)calloc(nsamp[i], sizeof(short))) == NULL) {
	    fprintf(stderr, "Insufficient memory for filtering\n");
	    break;
	}
	in = ecg[i];
	out[0] = in[0] * 0.67 + in[1] * 0.22 + in[2] * 0.11;
	out[1] = in[0] * 0.33 + in[1] * 0.34 + in[2] * 0.22 + in[3] * 0.11;
	for (j = 2; j < nsamp[i] - 2; j++)
	    out[j] = in[j]*0.34+(in[j-1]+in[j+1])*0.22+(in[j-2]+in[j+2])*0.11;
	out[j] = in[j-2]*0.11 + in[j-1]*0.22 + in[j]*0.34 + in[j+1]*0.33;
	j++;
	out[j] = in[j-2]*0.11 + in[j-1]*0.22 + in[j]*0.67;
	ecg[i] = out;
	free(in);
    }
}

/* ncompare is used by qsort() (in the calculation of the signal medians).
   It returns the difference between the values pointed to by its arguments
   (two samples of a signal). */
int ncompare(const void *a, const void *b)
{
    return (*(const short *)a - *(const short *)b);
}

int write_output()
{
    char dname[15];		/* name of .des file */
    char ename[15];		/* name of .ecg file */
    char kname[15];		/* name of .key file */
    char tname[15];		/* name of .txt file */
    FILE *dfile; /* (text) description file:  measurements, diagnoses, etc. */
    FILE *efile; /* (binary) ECG signal file */
    FILE *kfile; /* (text) key file: patient name, medical record number */
    FILE *tfile; /* (optional, text) ECG signal file (for debugging) */
    int i;
    static char ecgprolog[512];	/* prolog for .ecg file */
    void precgmeas();

    if (recname[0] == '\0') strcpy(recname, "ecg");  /* default record name */
    sprintf(dname, "%s.des", recname);
    sprintf(ename, "%s.ecg", recname);
    sprintf(kname, "%s.key", recname);
    if (vflag) sprintf(tname, "%s.txt", recname);
    if ((dfile = fopen(dname, "wt")) == NULL ||
	(efile = fopen(ename, "wb")) == NULL ||
	(kfile = fopen(kname, "wt")) == NULL ||
	(vflag && ((tfile = fopen(tname, "wt")) == NULL))) {
	if (dfile) fclose(dfile);
	else {
	    fprintf(stderr, "%s: can't create '%s'\n", pname, dname);
	    return (0);
	}
	if (efile) fclose(efile);
	else {
	    fprintf(stderr, "%s: can't create '%s'\n", pname, ename);
	    return (0);
	}
	if (kfile) fclose(kfile);
	else {
	    fprintf(stderr, "%s: can't create '%s'\n", pname, kname);
	    return (0);
	}
	fprintf(stderr, "%s: can't create '%s'\n", pname, tname);
	return (0);
    }

    /* Write the description file first. */
    fprintf(dfile, VERSION_STRING " description\n");
    fprintf(dfile, "Record: %s (%s %s)\n", recname, recdate, rectime);
    fprintf(dfile, "Age: %d\n", age);
    fprintf(dfile, "Sex: ");
    switch (sex) {
    case 1: fprintf(dfile, "M\n"); break;
    case 2: fprintf(dfile, "F\n"); break;
    case 0: fprintf(dfile, "not known\n"); break;
    case 9: fprintf(dfile, "unspecified\n"); break;
    default: fprintf(dfile, "<error>\n"); break; /* shouldn't happen */
    }
    fprintf(dfile, "Referring Dr: %s\n", referring_dr);
    fprintf(dfile, "Comments: %s\n", comments);
    if (tdevid == 0)
	fprintf(dfile, "Cart ID: <missing>\n");
    else
	fprintf(dfile, "Cart ID: %d\n", tdevid);
    fprintf(dfile, "Department: %s\n", dept);
    if (tdept == 0)
      fprintf(dfile, "Department #: <missing>\n");
    else
      fprintf(dfile, "Department #: %s\n", dept);
    fprintf(dfile, "Institution: %s\n", inst);
    if (tinst == 0)
      fprintf(dfile, "Institution #: <missing>\n");
    else
      fprintf(dfile, "Institution #: %d\n", tinst);

    for (i=0; i < stmnts; i++)
      fprintf(dfile, "Report %d: %s\n", i, report[i]);

    /* FIXME -- according to the SCP spec, units for blfilt should be
       .01 Hz, not .0001 Hz;  but files from Spacelabs carts have blfilt
       equal to 500 and clearly are diagnostic (0.05-100 Hz) bandwidth. */
    fprintf(dfile, "Bandwidth: %g-%d Hz\n",	blfilt*0.0001, lpfilt);
    fprintf(dfile, "Filtering:");
    if (filter_bit & 1) fprintf(dfile, " 60 Hz notch");
    if (filter_bit & 2) fprintf(dfile, " 50 Hz notch");
    if (filter_bit & 4) fprintf(dfile, " Artifact filter");
    if (filter_bit & 8) fprintf(dfile, " Baseline filter");
    if (filter_bit & 0xf0) fprintf(dfile, " <undefined>");
    if (filter_bit == 0) fprintf(dfile, " <not specified>");
    fprintf(dfile, "\n\n");

    fprintf(dfile, "RR interval: %d\n", RRint);
    fprintf(dfile, "P onset: %d\n", Pon);
    fprintf(dfile, "P end: %d\n", Pend);
    fprintf(dfile, "QRS onset: %d\n", QRSon);
    fprintf(dfile, "QRS end: %d\n", QRSend);
    fprintf(dfile, "T end: %d\n", Tend);
    fprintf(dfile, "P axis: %d\n", Paxis);
    fprintf(dfile, "QRS axis: %d\n", QRSaxis);
    fprintf(dfile, "T axis: %d\n", Taxis);

    precgmeas(dfile);
    fclose(dfile);

    /* Generate the key file.  The patient's name and ID number should
       appear only in this file (and in the debugging output, if
       enabled using the -v option). */
    fprintf(kfile, VERSION_STRING " key\n");
    fprintf(kfile, "Record: %s (%s %s)\n", recname, recdate, rectime);
    fprintf(kfile, "Name: %s\n", patient_name);
    fprintf(kfile, "MRN: %s\n", patient_id);
    fclose(kfile);

    /* If the -z option was selected, then for each signal:
         1. Calculate the median amplitude of the first 5000 samples in the
	    signal
	 2. If its final sample is zero, replace that sample with the previous
	    one.  This step eliminates transients that occur occasionally.
	 3. Subtract the calculated median amplitude from all of the samples.
            This step is meant to make it as likely as possible that the signal
            will be visible and neatly centered on a chart.

       If the result after subtracting the median falls outside of the
       +/- 10 millivolt range (+/- 2000 A/D units), it is assumed to be
       invalid, and the sample is stored with the (unique) value
       WFDB_INVALID_SAMPLE (-32768).  Software designed to interpret or
       display the data should avoid displaying samples with this value,
       or using them in calculations;  pschart (used by printchart) obeys
       this rule.
    */

    if (zflag) {
	int i, imax, s;
	static short buf[5000], m, v;

	for (s = 0; s < nleads; s++) {
	    imax = (nsamp[s] > 5000) ? 5000 : nsamp[s];
	    for (i = 0; i < imax; i++)
		buf[i] = ecg[s][i];
	    qsort(buf, imax, sizeof(short), ncompare);
	    m = buf[i/2];
	    if (ecg[s][nsamp[s]-1] == 0)
		ecg[s][nsamp[s]-1] = ecg[s][nsamp[s]-2];
	    for (i = 0; i < nsamp[0]; i++) {
		v = ecg[s][i] -= m;
		if (v < -2000 || v > 2000)
		    ecg[s][i] = WFDB_INVALID_SAMPLE;
	    }
	}
    }

    /* Write the WFDB-format record if requested (-w option). */
    if (wflag)
	write_wfdb_record();

    /* Write the ECG in text form if requested (-v option). */
    if (vflag) {
	int i;

	if (vflag > 1)
	    fprintf(tfile,
		    "samp#      I     II    III    aVR    aVF    aVL"
		         "     V1     V2     V3     V4     V5     V6");
        for (i = 0; i < nsamp[0]; i++) {
	    fprintf(tfile, "%5d", i);
	    fprintf(tfile, "%7d", ecg[  leadI][i]);
	    fprintf(tfile, "%7d", ecg[ leadII][i]);
	    fprintf(tfile, "%7d", ecg[leadIII][i]);
	    fprintf(tfile, "%7d", ecg[    aVR][i]);
	    fprintf(tfile, "%7d", ecg[    aVF][i]);
	    fprintf(tfile, "%7d", ecg[    aVL][i]);
	    fprintf(tfile, "%7d", ecg[     V1][i]);
	    fprintf(tfile, "%7d", ecg[     V2][i]);
	    fprintf(tfile, "%7d", ecg[     V3][i]);
	    fprintf(tfile, "%7d", ecg[     V4][i]);
	    fprintf(tfile, "%7d", ecg[     V5][i]);
	    fprintf(tfile, "%7d", ecg[     V6][i]);
	    fprintf(tfile, "\n");
	}
	fclose(tfile);
    }

    /* Finally, output the ECG in vector (binary) form.
    ** FIXME ** This assumes that the sampling frequency is 500 Hz, that
    nsamp[leadII] >= 5000, and that nsamp[all other leads] >= 1250  */
    if (lflag == 0) {
	sprintf(ecgprolog, VERSION_STRING " signals\n"
		"Record: %s (%s %s)\n", recname, recdate, rectime);
    }
    else {
	sprintf(ecgprolog, VERSION_STRING " signals\n"
		"Record: %s (%s %s)\nLow-pass filtered\n",
		recname, recdate, rectime);
	lowpass();
    }
    fwrite(ecgprolog, 1, sizeof(ecgprolog), efile);
    fwriteswapped(ecg[  leadI]     , 2, 1250, efile);
    fwriteswapped(ecg[    aVR]+1250, 2, 1250, efile);
    fwriteswapped(ecg[     V1]+2500, 2, 1250, efile);
    fwriteswapped(ecg[     V4]+3750, 2, 1250, efile);
    fwriteswapped(ecg[ leadII]     , 2, 1250, efile);
    fwriteswapped(ecg[    aVL]+1250, 2, 1250, efile);
    fwriteswapped(ecg[     V2]+2500, 2, 1250, efile);
    fwriteswapped(ecg[     V5]+3750, 2, 1250, efile);
    fwriteswapped(ecg[leadIII]     , 2, 1250, efile);
    fwriteswapped(ecg[    aVF]+1250, 2, 1250, efile);
    fwriteswapped(ecg[     V3]+2500, 2, 1250, efile);
    fwriteswapped(ecg[     V6]+3750, 2, 1250, efile);
    fwriteswapped(ecg[ leadII]     , 2, 5000, efile);
    fclose(efile);
 
    return (1);
}


int write_wfdb_record()
{
#ifdef NOWFDB
    fprintf(stderr, "Warning: this version of %s was not compiled with"
	    " WFDB support\n", pname);
    return (0);
#else
    char *filename = malloc(strlen(recname) + 5), ptime[42];
    int i;
    static char *leadname[12] = { "I", "II", "III", "aVR", "aVF", "aVL",
				  "V1", "V2", "V3", "V4", "V5", "V6" };
    static char info[20];
    static WFDB_Sample v[12];
    static WFDB_Siginfo s[12];

    if (filename == NULL) {
	fprintf(stderr, "insufficient memory in write_wfdb_record\n");
	return (0);
    }
    sprintf(filename, "%s.dat", recname);

    if (newheader(recname) < 0)
	return(0);

    setsampfreq(500.);	/* FIXME: we should calculate this! */
    for (i = 0; i < 12; i++) {
	s[i].fname = filename;
	s[i].desc = leadname[i];
	s[i].units = "mV";
	s[i].fmt = 16;
	s[i].gain = 200;	/* FIXME: we should calculate this! */
	s[i].adcres = 16;	/* FIXME: we should calculate this! */
    }
    if (osigfopen(s, 12) < 12)
	return (0);
    sprintf(ptime, "%s %s\n", recdate, rectime);
    setbasetime(ptime);

    /* write the ECG samples */
    for (i = 0; i < nsamp[0]; i++) {
	v[ 0] = ecg[  leadI][i];
	v[ 1] = ecg[ leadII][i];
	v[ 2] = ecg[leadIII][i];
	v[ 3] = ecg[    aVR][i];
	v[ 4] = ecg[    aVF][i];
	v[ 5] = ecg[    aVL][i];
	v[ 6] = ecg[     V1][i];
	v[ 7] = ecg[     V2][i];
	v[ 8] = ecg[     V3][i];
	v[ 9] = ecg[     V4][i];
	v[10] = ecg[     V5][i];
	v[11] = ecg[     V6][i];
	if (putvec(v) < 0) break;
    }
    (void)newheader(recname);
    sprintf(info, "<age>: %d <sex>: %c", (age < 90) ? age : 90,
	    ((sex == 1) ? 'M' : (sex == 2) ? 'F' : '?'));
    putinfo(info);
    wfdbquit();
    return (1);
#endif
}

void help()
{
    fprintf(stderr, "usage: %s -o RECORD [OPTIONS ...] <SCPFILE\n",
	    pname);
    fprintf(stderr, "  where RECORD is the name (<= 11 alphanumeric"
		    " characters) of this case,\n"
		    "  SCPFILE is the name of a file containing ECG"
		    " data in SCP-ECG format, and\n"
                    "  OPTIONS may include:\n"
	    	    "   -a    anonymize and copy stdin"
	            " to stdout (suppresses other output).\n"
		    "   -f    force parsing of a file with CRC"
		    " errors.\n"
	            "   -h    print brief help\n"
		    "   -l    low-pass filter (smooth) the"
		    " output waveforms.\n"
		    "   -v    option to print a (very) detailed analysis"
		    " of SCPFILE.\n"
#ifndef NOWFDB
		    "   -w    produce a WFDB-format record\n"
#endif
	            "   -z    zero-mean the input\n\n"
	            "  Additional options for debugging only:\n"
	            "   -b    show baseline (residuals);  suppress templates\n"
	            "   -s N  shift templates by N samples before adding\n"
	            "   -S N  skip parsing of section N\n"
                    "   -t    show templates; suppress baseline\n"
	    	    "   -x    show hex dumps (implies -v)\n"
	    );
    return;
}

int aflag = 0;
int skip[12];

main(int argc, char **argv)
{
    unsigned short crc, sec_id;
    unsigned long bytesread, length, sec_len;
    unsigned char desc[20], header[6], *data, *p;
    int i, j, k, m;
    FILE *ofile, *vfile;
    char *prog_name();

    pname = prog_name(argv[0]);

    for (i = 1; i < argc; i++) {
	if (argv[i][0] == '-') switch (argv[i][1]) {
	case 'a':	/* write an anonymized version of stdin to stdout */
	    aflag = 1;
	    xflag = vflag = 0;	/* incompatible with -v and -x */
	    break;
	case 'b':	/* show baseline (residuals);  suppress templates */
	    bflag = 1; tflag = 0;
	    break;
	case 'f':	/* force parsing even if CRC error(s) found */
	    fflag = 1;
	    break;
        case 'h':	/* print a usage summary */
	    help();
	    exit(1);
	case 'l':	/* low-pass filter the output waveforms */
	    lflag = 1;
	    break;
	case 'o':	/* record name follows */
	    if (++i < argc && strlen(argv[i]) <= 11)
		strcpy(recname, argv[i]);
	    else {
		fprintf(stderr, "%s: a record name of 11 or fewer alphanumeric"
			"characters must\n follow -o\n", pname);
		exit(1);
	    }
	    break;
	case 's':	/* shift templates by N samples */
	    if (++i < argc)
		shift = atoi(argv[i]);
	    else {
		fprintf(stderr,
		   "%s: number of samples to shift templates must follow -s\n",
			pname);
		exit(1);
	    }
	    break;
	case 'S':       /* skip a section (may be repeated) */
	    if (++i < argc) {
		if (argv[i][0] == '-')
		    skip[0] = skip[1] = skip[2] = skip[3] = skip[4] = skip[5] =
		      skip[6] = skip[7] = skip[8] = skip[9] = skip[10] =
			skip[11] = 1;
		else {
		    int index = atoi(argv[i]);
		    if (0 <= index && index <= 11)
			skip[index] = 1 - skip[index];
		}
	    }
	    else {
		fprintf(stderr,
		   "%s: number of section to skip, or '-', must follow -S\n",
			pname);
		exit(1);
	    }
	    break;
	case 't':	/* show templates;  suppress baseline */
	    tflag = 1; bflag = 0;
	    break;
	case 'v':	/* print a detailed analysis of the record contents */
	    vflag++;
	    aflag = 0;	/* incompatible with -a */
	    break;
	case 'w':	/* create a WFDB-format record */
	    wflag++;
	    break;
	case 'x':	/* print a hex dump (implies -v) */
	    xflag = vflag = 1;
	    aflag = 0;
	    break;
        case 'z':	/* suppress final transients and zero-mean the ECGs */
	    zflag++;
	    break;
	default:
	    fprintf(stderr, "%s: unrecognized option '%s'\n", pname, argv[i]);
	    help();
	    exit(1);
	}
	else {
	    fprintf(stderr,"%s: unrecognized argument '%s'\n", pname, argv[i]);
	    help();
	    exit(1);
	}
    }

#ifdef _WINDOWS
    _setmode(_fileno(stdin), _O_BINARY);
    if (aflag) _setmode(_fileno(stdout), _O_BINARY);
#endif

    bytesread = fread(header, 1, 6, stdin);
    if (bytesread != 6) {
	fprintf(stderr, "%s: input too short (only %lu bytes)\n",
		pname, bytesread);
	exit(1);
    }	
    crc = get16(header);
    length = get32(header+2);
    if (vflag) printf("CRC = %d, expected length = %ld bytes\n", crc, length);
    if ((data = (unsigned char *)malloc(length)) == NULL) {
	fprintf(stderr, "%s: insufficient memory (%ld bytes needed)\n",
		pname, length);
	exit(2);
    }
    memcpy(data, header, 6);
    i = 1;
    while (!ferror(stdin) && bytesread < length && i > 0)
	bytesread += i = fread(data+bytesread, 1, length-bytesread, stdin);
    if (bytesread < length) {
	fprintf(stderr, "%s: input too short (%lu byte%s missing)\n",
		pname, length-bytesread, (length-bytesread == 1) ? "" : "s");
	exit(3);
    }
    if (vflag) printf("Record length: %ld bytes\n", length);
    crcok(data+2, length-2, crc, "Top-level CRC");

    /* Now check the sections. */
    p = data + 6;	/* skip the record header -- p points to section 0 */

    while (p < data + length) {	/* true if there is at least 1 more section */
	crc = get16(p);
	sec_id = get16(p+2);
	sec_len = get32(p+4);
	if (vflag) printf("\nSection %d length: %ld bytes\n", sec_id, sec_len);
	if (sec_len < 8) {
	    if (vflag) printf(
    " Warning: section length too short (must be at least 8 bytes)\n"
    "  Remaining data (%lu bytes) following short section will not be read\n",
			      (unsigned long) (data+length-(p+8)));
	    sec_len = 8;
	    //	    break;	/* don't attempt to read any further */
	}
	if (sec_len > data + length - p) {
            if (vflag) printf(
	       " Warning: section length exceeds amount of remaining data\n"
	       "  This section will not be parsed\n");
	    break;
	}
	sprintf(desc, " Section %d CRC", sec_id);
	crcok(p+2, sec_len-2, crc, desc);
	if (vflag) printf(" Section version number: %d\n", p[8]);
	if (vflag) printf(" Protocol version number: %d\n", p[9]);
	if (vflag) printf(" Contents: ");
	switch (sec_id) {
	case 0: if (vflag) printf(" Pointers to data areas in the record\n");
	    if (skip[0]) printf("  [skipped]\n");
	    else section0(p, sec_len);
	    break;
	case 1: if (vflag) printf(" Header information - Patient data/ECG acquisition"
		       " data\n");
	    if (skip[1]) printf("  [skipped]\n");
	    else section1(p, sec_len);
	    break;
	case 2: if (vflag) printf(" Huffman tables used in encoding of ECG data\n");
	    if (skip[2]) printf("  [skipped]\n");
	    else section2(p, sec_len);
	    break;
	case 3: if (vflag) printf(" ECG lead definition\n");
	    if (skip[3]) printf("  [skipped]\n");
	    else section3(p, sec_len);
	    break;
	case 4: if (vflag) printf(" QRS locations\n");
	    if (skip[4]) printf("  [skipped]\n");
	    else section4(p, sec_len);
	    break;
	case 5: if (vflag) printf(" Encoded reference beat data\n");
	    if (skip[5]) printf("  [skipped]\n");
	    else section5(p, sec_len);
	    break;
	case 6: if (vflag) printf(" Residual signal after reference beat subtraction (if"
		       " reference beats\n"
		       "  are stored);  otherwise, encoded rhythm data\n");
	    if (skip[6]) printf("  [skipped]\n");
	    else section6(p, sec_len);
	    break;
	case 7: if (vflag) printf(" Global measurements\n");
	    if (skip[7]) printf("  [skipped]\n");
	    else section7(p, sec_len);
	    break;
	case 8: if (vflag) printf(" Textual diagnosis from the interpretive device\n");
	    if (skip[8]) printf("  [skipped]\n");
	    else section8(p, sec_len);
	    break;
	case 9: if (vflag) printf(" Manufacturer specific diagnostic and overreading"
		       " data\n");
	    if (skip[9]) printf("  [skipped]\n");
	    else section9(p, sec_len);
	    break;
	case 10:if (vflag) printf(" Lead measurement results\n");
	    if (skip[10]) printf("  [skipped]\n");
	    else section10(p, sec_len);
	    break;
	case 11:if (vflag) printf(" Universal statement codes resulting from the"
		       " interpretation\n");
	    if (skip[11]) printf("  [skipped]\n");
	    else section11(p, sec_len);
	    break;
	default:if (vflag) printf(" Manufacturer-defined (non-standard) section %d\n", sec_id);
	     break;
	}
	p += sec_len;
    }

    /* Low-pass filter the residuals. */
    for (i = 0; i < nleads; i++) {
	short *out, *in;
	int v;

	if (ecg[i] == NULL) continue;
	if ((out = (short *)calloc(nsamp[i], sizeof(short))) == NULL) {
	    fprintf(stderr, "Insufficient memory for filtering\n");
	    break;
	}
	in = ecg[i];
	if (tflag == 0) {  /* true unless -t debugging option was given */
	    out[0] = in[0];
	    for (j = 1; j < nsamp[i] - 1; j++) {
		v = (in[j-1] + in[j] + in[j+1]);
		if (v > 0) out[j] = (v + 2)/3;
		else if (v < 0) out[j] = (v - 2)/3;
		else out[j] = 0;
	    }
	    out[j] = in[j];
	    for (k = 0; k < nqrs; k++)
		if (subz[k].type == 0)
		  for (j = subz[k].t0 - 1; j <= subz[k].t2 + 1; j++)
			out[j] = in[j];	/* restore unfiltered samples within
					   protected zones and adjacent samples
					   on either end of protected zones */
	}
	ecg[i] = out;
	free(in);
    }

    /* Add the reference beats back into the residuals to obtain the raw ECG */
    if (bflag == 0) {	/* true unless -b debugging option was given */
	for (k = 0; k < nqrs; k++)
	    for (i = 0; subz[k].type == 0 && i < nleads; i++) {
		if (subz[k].t2 - subz[k].t0 + 1 > rblen ||
		    (m = fcM - subz[k].t1 + subz[k].t0 - 1) < 0) {
		    fprintf(stderr,
			"error: subtraction zone is larger than template\n");
		    exit(1);
		}
		for (j = subz[k].t0-1+shift; j <= subz[k].t2-1+shift; j++, m++)
		    ecg[i][j] += refbeat[i][m];
	    }
    }

    /* Now calculate the missing leads.  Normally the input contains only two
       of the three limb leads (I, II, and III) and the missing one (usually
       III) must be calculated from the other two, using the relationship:
		lead III = lead II - lead I				    */
    if (leadIII < 0) {	/* we need to calculate lead III */
	if (leadI >= 0 && leadII >= 0) {
	    if (vflag) printf(" Calculating lead III\n");
	    leadIII = nleads++;
	    nsamp[leadIII] = min(nsamp[leadI], nsamp[leadII]);
	    for (j = 0; j < nsamp[leadIII]; j++)
		ecg[leadIII][j] = ecg[leadII][j] - ecg[leadI][j];
	}
	else if (leadI < 0) {
	    if (leadII < 0)
		fprintf(stderr, "Leads I, II, and III are missing\n");
	    else
		fprintf(stderr, "Leads I and III are missing\n");
	}
	else
	    fprintf(stderr, "Leads II and III are missing\n");
    }
    else if (leadII < 0) { /* we need to calculate lead II */
	if (leadI >= 0) {
	    if (vflag) printf(" Calculating lead II\n");
	    leadII = nleads++;
	    nsamp[leadII] = min(nsamp[leadI], nsamp[leadIII]);
	    for (j = 0; j < nsamp[0]; j++)
		ecg[leadII][j] = ecg[leadI][j] + ecg[leadIII][j];
	}
	else
	    fprintf(stderr, "Leads I and II are missing\n");
    }
    else if (leadI < 0) { /* we need to calculate lead I */
	if (vflag) printf(" Calculating lead I\n");
	leadI = nleads++;
	nsamp[leadI] = min(nsamp[leadII], nsamp[leadIII]);
	for (j = 0; j < nsamp[0]; j++)
	    ecg[leadI][j] = ecg[leadII][j] - ecg[leadIII][j];
    }

    /* We need both lead I and lead II to derive aVR. */
    if (aVR < 0) {
	if (leadI >= 0 && leadII >= 0) {
	    if (vflag) printf(" Calculating lead aVR\n");
	    aVR = nleads++;
	    nsamp[aVR] = min(nsamp[leadI], nsamp[leadII]);
	    for (j = 0; j < nsamp[0]; j++)
		ecg[aVR][j] = - (ecg[leadI][j] + ecg[leadII][j])/2;
	}
	else
	    fprintf(stderr, "Lead aVR cannot be derived\n");
    }
    
    /* We need both lead II and lead III to derive aVL. */
    if (aVL < 0) {
	if (leadII >= 0 && leadIII >= 0) {
	    if (vflag) printf(" Calculating lead aVL\n");
	    aVL = nleads++;
	    nsamp[aVL] = min(nsamp[leadII], nsamp[leadIII]);
	    for (j = 0; j < nsamp[0]; j++)
		ecg[aVL][j] = ecg[leadII][j]/2 - ecg[leadIII][j];
	}
	else
	    fprintf(stderr, "Lead aVL cannot be derived\n");
    }

    /* We need both lead I and lead III to derive aVF. */
    if (aVF < 0) {
	if (leadI >= 0 && leadII >= 0) {
	    if (vflag) printf(" Calculating lead aVF\n");
	    aVF = nleads++;
	    nsamp[aVF] = min(nsamp[leadI], nsamp[leadIII]);
	    for (j = 0; j < nsamp[0]; j++)
		ecg[aVF][j] = ecg[leadI][j]/2 + ecg[leadIII][j];
	}
	else
	    fprintf(stderr, "Lead aVF cannot be derived\n");
    }

    if (V1 < 0) fprintf(stderr, "Lead V1 is missing\n");
    if (V2 < 0) fprintf(stderr, "Lead V2 is missing\n");
    if (V3 < 0) fprintf(stderr, "Lead V3 is missing\n");
    if (V4 < 0) fprintf(stderr, "Lead V4 is missing\n");
    if (V5 < 0) fprintf(stderr, "Lead V5 is missing\n");
    if (V6 < 0) fprintf(stderr, "Lead V6 is missing\n");

    if (nleads > 0 && (leadI < 0 || leadII < 0 || leadIII < 0 ||
		       aVR < 0 || aVF < 0 || aVL < 0 ||
		       V1 < 0 || V2 < 0 || V3 < 0 ||
		       V4 < 0 || V5 < 0 || V6 < 0)) {
	/* Allocate an array of zero samples to use for the missing lead(s). */
	if (ecg[nleads] == NULL &&
	    (ecg[nleads] = (short *)calloc(nsamp[i], sizeof(short))) == NULL) {
	    fprintf(stderr, "Insufficient memory\n");
	    exit(1);
	}
	if (leadI < 0) leadI = nleads;
	if (leadII < 0) leadII = nleads;
	if (leadIII < 0) leadIII = nleads;
	if (aVR < 0) aVR = nleads;
	if (aVF < 0) aVF = nleads;
	if (aVL < 0) aVL = nleads;
	if (V1 < 0) V1 = nleads;
	if (V2 < 0) V2 = nleads;
	if (V3 < 0) V3 = nleads;
	if (V4 < 0) V4 = nleads;
	if (V5 < 0) V5 = nleads;
	if (V6 < 0) V6 = nleads;
    }

    if (aflag) {
	/* Recalculate the top-level CRC. */
	crc = getcrc(data+2, bytesread-2);
	put16(crc, data);
	fwrite(data, 1, bytesread, stdout);
    }
    else if (nleads > 0) {
	write_output();
    }

    /* Release all allocated memory. */
    if (refbeat) {
	for (i = 0; i < nleads && refbeat[i] != NULL; i++)
	    free(refbeat[i]);
	free(refbeat);
    }
    if (ecg) {
	for (i = 0; i <= nleads && ecg[i] != NULL; i++)
	    free(ecg[i]);
	free(ecg);
    }
    if (nsamp) free(nsamp);
    if (subz) free(subz);
    if (ecgmeas) free(ecgmeas);
    if (data) free(data);

    exit(0);
}

void hexdump(unsigned char *p, long len)
{
    int i;

    printf("Hex dump of data area:\n\n");
    for (i = 0; i < len; i++) {
      if (i % 16 == 0) printf("\n%4d:  ", i);
      printf(" %2x", p[i]);
    }
    printf("\n\n");
}

int section0(unsigned char *p, long len)
{
    if (strncmp("SCPECG", p+10, 6) != 0)
	if (vflag) printf(" Warning: SCPECG identifier is missing in header\n");
    p += 16; len -= 16;	/* move to data area */
    if (len < 10) {
	if (vflag) printf(" Error: section 0 contains no data\n");
	return (0);
    }

    while (len > 0) {
	unsigned short sec_id;
	unsigned long sec_len, sec_index;

	sec_id = get16(p);
	sec_len = get32(p+2);
	sec_index = get32(p+6);
	if (vflag) printf("  Section %2d:%6ld bytes beginning at byte %6ld\n",
	       sec_id, sec_len, sec_index);
	len -= 10;
	p += 10;
    }
    if (len < 0) {
	if (vflag) printf(" Error: section 0 overlaps the next section\n");
	return (0);
    }
    return (1);
}

char *month[] = { "January", "February", "March", "April", "May", "June",
    "July", "August", "September", "October", "November", "December" };

void censor(char *p, unsigned short len, char c)
{
    while (--len)
	*p++ = c;
    *p = '\0';
}

int section1(unsigned char *p, long len)
{
    FILE *pfile;
    unsigned short vlen;
    char format[16], *tpatid = NULL, *trefdr = NULL, *tcomments = NULL, *tdate = NULL, *ttime = NULL, *pyoa = NULL, *pyob = NULL;
    char *firstname = NULL, *lastname = NULL;
    enum { text, date, time, b1, b2, b2b1, bn, mixed } vtype;
    unsigned char *p_save = p;
    long len_save = len;

    p += 16; len -= 16; /* move to data area */
    if (len < 3) {
	if (vflag) printf(" Error: section 1 contains no data\n");
	return (0);
    }

    while (len > 2) {
	if (*p == 255) break;	/* tag 255 is the terminator */
	vlen = get16(p+1);
	if (vlen > len) {
	    if (vflag) printf(" Error: tagged field overlaps next section\n");
	    return (0);
	}
	switch (*p) {
	case  0: if (vflag) printf("\n  Last name\n");
	    vtype = text; lastname = p+3;
	    if (aflag) censor(lastname, vlen, 'x');
	    break;
	case  1: if (vflag) printf("\n  First name\n");
	    vtype = text; firstname = p+3;
	    if (aflag) censor(firstname, vlen, 'x');
	    break;
	case  2: if (vflag) printf("\n  Patient ID number\n");
	    vtype = text; tpatid = p+3;
	    if (aflag) censor(tpatid, vlen, '0');
	    break;
	case  3: if (vflag) printf("\n  Second last name\n");
	    vtype = text;
	    if (aflag) censor(p+3, vlen, 'x');
	    break;
	case  4: if (vflag) printf("\n  Age\n");
	    vtype = b2b1; age = get16(p+3); /* FIXME */
	    if (aflag && age > 90)
		put16(90, p+3);	/* censor age if over 90 */
	    break;
	case  5: if (vflag) printf("\n  Date of birth\n");
	    vtype = date;
	    if (aflag) {
		pyob = p+3;	/* save pointer to year of birth */
		put16(257, p+5); /* censor month and day (257 => 1 Jan) */
	    }
	    break;
	case  6: if (vflag) printf("\n  Height\n");
	    vtype = b2b1; break;
	case  7: if (vflag) printf("\n  Weight\n");
	    vtype = b2b1; break;
	case  8: if (vflag) printf("\n  Sex\n");
	    vtype = b1; sex = *(p+3); break;
        case  9: if (vflag) printf("\n  Race\n");
	    vtype = b1; break;
	case 10: if (vflag) printf("\n  Drug\n");
	    vtype = mixed; break;
	case 11: if (vflag) printf("\n  Systolic blood pressure\n");
	    vtype = b2; break;
	case 12: if (vflag) printf("\n  Diastolic blood pressure\n");
	    vtype = b2; break;
	case 13: if (vflag) printf("\n  Diagnosis or referral indication\n");
	    vtype = text; break;
	case 14: if (vflag) printf("\n  ID of the acquiring device\n");
	    tdevid = get16(p+7);
	    tinst = get16(p+3);  /* institution # */
	    tdept = get16(p+5);  /* department # */
	    if (vflag) {
		char *q;

		printf("   Institution number: %d\n", get16(p+3));
		printf("   Department number: %d\n", get16(p+5));
		printf("   Device ID: %d\n", get16(p+7));
		printf("   Device type: %d (%s)\n", *(p+9),
			(*(p+9) == 0 ? "cart" :
			 (*(p+9) == 1 ? "system or host" : "<error>")));
		printf("   Manufacturer code: %d\n", *(p+10));
		printf("   Text model description: [%d]", *(p+11));
		printf("   SCP-ECG protocol revision number: %g\n",
		       (double)(*(p+17))/10.0);
		printf("   SCP-ECG protocol compatibility level: %d\n",
		       *(p+18));
		printf("   Language support code: %d\n", *(p+19));
		printf("   Device capabilities: %d\n", *(p+20));
		printf("   AC mains frequency environment: %d ", *(p+21));
		switch (*(p+21)) {
		case 0: printf("(unspecified)\n"); break;
		case 1: printf("(50 Hz)\n"); break;
		case 2: printf("(60 Hz)\n"); break;
		default: printf("<error>"); break;
		}
		q = p+39;
		printf("   Analyzing program revision number: [%s]\n", q);
		q += strlen(q) + 1;
		printf("   Serial number of acquisition device: [%s]\n", q);
		q += strlen(q) + 1;
		printf("   Acquisition device system software ID: [%s]\n", q);
		q += strlen(q) + 1;
		printf("   Acquisition device SCP implementation ID: [%s]\n",
		       q);
		q += strlen(q) + 1;
		printf("   Manufacturer of acquisition device: [%s]\n", q);
	    }
	    vtype = mixed; break;
	case 15: if (vflag) printf("\n  ID of the analyzing device\n");
	    vtype = mixed; break;
	case 16: if (vflag) printf("\n  Acquiring institution\n");
	    vtype = text; inst=p+3;
	    if (aflag) censor(p+3, vlen, 'x');
	    break;
	case 17: if (vflag) printf("\n  Analyzing institution\n");
	    vtype = text;
	    if (aflag) censor(p+3, vlen, 'x');
	    break;
	case 18: if (vflag) printf("\n  Acquiring department\n");
	    vtype = text;
	    if (aflag) censor(p+3, vlen, 'x');
	    dept = p+3; break;
	case 19: if (vflag) printf("\n  Analyzing department\n");
	    vtype = text;
	    if (aflag) censor(p+3, vlen, 'x');
	    break;
	case 20: if (vflag) printf("\n  Referring physician\n");
	    vtype = text;
    	    if (aflag) censor(p+3, vlen, 'x');
	    trefdr = p+3; break;
	case 21: if (vflag) printf("\n  Latest confirming physician\n");
	    vtype = text;
	    if (aflag) censor(p+3, vlen, 'x');
	    break;
	case 22: if (vflag) printf("\n  Technician\n");
	    vtype = text;
	    if (aflag) censor(p+3, vlen, 'x');
	    break;
	case 23: if (vflag) printf("\n  Room\n");
	    vtype = text;
	    if (aflag) censor(p+3, vlen, 'x');
	    break;
	case 24: if (vflag) printf("\n  Stat code (urgency)\n");
	    vtype = b1; break;
	case 25: if (vflag) printf("\n  Date of acquisition\n"); tdate = p+3;
	    if (aflag) {
		pyoa = p+3;	/* save pointer to year of acquisition */
		put16(257, p+5); /* censor month and day (257 => 1 Jan) */
	    }
	    vtype = date; break;
	case 26: if (vflag) printf("\n  Time of acquisition\n"); ttime = p+3;
	    vtype = time; break;
	case 27: if (vflag) printf("\n  Baseline filter\n");
	    vtype = b2; blfilt = get16(p+3); break;
	case 28: if (vflag) printf("\n  Low-pass filter\n");
	    vtype = b2; lpfilt = get16(p+3); break;
	case 29: if (vflag) printf("\n  Filter bit map\n");
	    vtype = b1; filter_bit = *(p+3); break;
	case 30: if (vflag) printf("\n  Free text field (comments)\n"); tcomments = p+3;
	    vtype = text; break;
	case 31: if (vflag) printf("\n  Sequence number\n");
	    vtype = text; break;
	case 32: if (vflag) printf("\n  Medical history codes\n");
	    vtype = bn; break;
	case 33: if (vflag) printf("\n  Electrode configuration code\n");
	    vtype = b2; break;
	case 34: if (vflag) printf("\n  Time zone\n");
	    vtype = mixed; break;
	case 35: if (vflag) printf("\n  Free text medical history\n");
	    vtype = text; break;
	default: if (vflag) printf("\n  <Undefined tag>\n");
	    vtype = mixed; break;
	}
	if (vflag) printf("   Tag%3d (%2d bytes): ", *p, vlen);
	p += 3; len -= vlen+3;
	if (vlen == 0) { if (vflag) printf("<not defined>\n"); }
	else switch (vtype) {
	case text:
	    if (vlen > 0) {
		sprintf(format, "[%%%ds]", vlen-1);
		if (vflag) printf(format, p);
		p += vlen;
	    }
	    else
		if (vflag) printf("<undefined>");
	    break;
	case date:
	    if (vlen != 4) {
		if (vflag) printf("  Error: incorrect date format (%d bytes)\n", vlen);
	    }
	    else {
		if (vflag) printf("%4d/%02d/%02d", get16(p), *(p+2), *(p+3));
	    }
	    p += vlen;
	    break;
	case time:
	    if (vlen != 3) {
		if (vflag) printf("  Error: incorrect time format (%d bytes)\n", vlen);
	    }
	    else 
		if (vflag) printf("%02d:%02d:%02d", *p, *(p+1), *(p+2));
	    p += vlen;
	    break;
	case b1:
	    if (vlen != 1) {
		if (vflag) printf("  Error: incorrect format (%d bytes instead of 1)\n",
		       vlen);
	    }
	    else
		if (vflag) printf("%d", *p);
	    p += vlen;
	    break;
	case b2:
	    if (vlen != 2) {
		if (vflag) printf("  Error: incorrect format (%d bytes instead of 2)\n",
		       vlen);
	    }
	    else
		if (vflag) printf("%d", get16(p));
	    p += vlen;
	    break;
	case b2b1:
	    if (vlen != 3) {
		if (vflag) printf("  Error: incorrect format (%d bytes instead of 3)\n",
		       vlen);
	    }
	    else
		if (vflag) printf("%d (code %d)", get16(p), *(p+2));
	    p += vlen;
	    break;
	case bn:
	    while (vlen-- > 0) {
		if (vflag) printf("%d, ", *p);
		p++;
	    }
	    break;
	default:
	    if (vflag) printf("oops!  undefined variable type\n");
	case mixed:
	    if (vflag) {
		while (vlen-- > 0) {
		    if (isprint(*p)) putchar(*p);
		    else if (vflag) printf("\\%03o", *p);
		    p++;
		}
		printf("\n");
	    }
	    else
		p += vlen;
	}
    }

    sprintf(recdate, "%d %s %d",*(tdate+3), month[*(tdate+2)-1], get16(tdate));
    sprintf(rectime, "%02d:%02d:%02d", *ttime, *(ttime+1), *(ttime+2));
    if (lastname == NULL) lastname = "<last name missing>";
    if (firstname == NULL) firstname = "<first name missing>";
    if (strlen(lastname) + strlen(firstname) < sizeof(patient_name) - 3)
	sprintf(patient_name, "%s, %s", lastname, firstname);
    else {
	fprintf(stderr, "Patient name (%s, %s) is too long\n", lastname,
		firstname);
	strcpy(patient_name, "Name not available -- too long");
    }
    if (tpatid == NULL) tpatid = "<patient ID missing>";
    strncpy(patient_id, tpatid, sizeof(patient_id));
    if (trefdr == NULL) trefdr = "<referring physician missing>";
    strncpy(referring_dr, trefdr, sizeof(referring_dr));
    if (tcomments == NULL) tcomments = "<comments missing>";
    strncpy(comments, tcomments, sizeof(comments));

    if (*p != 255)
	if (vflag) printf(" Error: header terminator (tag 255) is missing\n");
    if (tpatid == NULL)
	if (vflag) printf(" Error: patient ID number (tag 2) is missing\n");
    if (tdevid == -1 || tinst == -1 || tdept == -1)
	if (vflag) printf(" Error: device ID number (tag 14) is missing\n");
    if (tdate == NULL)
	if (vflag) printf(" Error: date of acquisition (tag 25) is missing\n");
    if (ttime == NULL)
	if (vflag) printf(" Error: time of acquisition (tag 26) is missing\n");

    if (aflag && pyob != NULL) {
	int yoa = 2004, yob;
	unsigned short crc;

	if (pyoa) yoa = get16(pyoa);
	yob = get16(pyob);
	if (yoa - yob > 90) yob = yoa - 90;	/* censor year of birth if >90
						   years before acquisition */

	/* Recalculate the CRC for this section. */
	crc = getcrc(p_save+2, len_save-2);
	put16(crc, p_save);
    }

    if (*p != 255 || tpatid == NULL || tdevid == -1 ||
	tdate == NULL || ttime == NULL)
	return (0);
    return (1);
}

int section2(unsigned char *p, long len)
{
    unsigned short ntables;

    p += 16; len -= 16; /* move to data area */
    if (len < 2) {
	if (vflag) printf("  Error: section 2 contains no data\n");
	return (0);
    }

    ntables = get16(p);
    if (ntables == 19999) {
	if (vflag) printf("  Using standard Huffman table only\n");
    }
    else {
	if (vflag) printf("  Number of Huffman tables defined: %d\n", get16(p));
	if (vflag) printf("  Number of code structures in table 1: %d\n", get16(p+2));
	/* ** more to do here ** */
    }
    return (1);
}

char *leadname[] = {		/* standard lead names */
    "<Unspecified lead>", "I", "II", "V1", "V2",		 /*  0 -  4 */
    "V3", "V4", "V5", "V6", "V7",				 /*  5 -  9 */
    "V2R", "V3R", "V4R", "V5R", "V6R",			 	 /* 10 - 14 */
    "V7R", "X", "Y", "Z", "CC5",				 /* 15 - 19 */
    "CM5", "Left Arm", "Right Arm", "Left Leg", "I (Frank)",	 /* 20 - 24 */
    "E", "C", "A", "M", "F",				 	 /* 25 - 29 */
    "H", "I-cal", "II-cal", "V1-cal", "V2-cal",			 /* 30 - 34 */
    "V3-cal", "V4-cal", "V5-cal", "V6-cal", "V7-cal",		 /* 35 - 39 */
    "V2R-cal", "V3R-cal", "V4R-cal", "V5R-cal", "V6R-cal",	 /* 40 - 44 */
    "V7R-cal", "X-cal", "Y-cal", "Z-cal", "CC5-cal",		 /* 45 - 49 */
    "CM5-cal", "Left Arm-cal", "Right Arm-cal", "Left Leg-cal","I-cal (Frank)",
    								 /* 50 - 54 */
    "E-cal", "C-cal", "A-cal", "M-cal", "F-cal",		 /* 55 - 59 */
    "H-cal", "III", "aVR", "aVL", "aVF",			 /* 60 - 64 */
    "-aVR", "V8", "V9", "V8R", "V9R",			 	 /* 65 - 69 */
    "D (Nehb - Dorsal)", "A (Nehb - Anterior)", "J (Nehb - Inferior",
    "Defibrillator lead: anterior-lateral",
    "External pacing lead: anterior-posterior",		 	 /* 70 - 74 */
    "A1 (Auxiliary unipolar lead 1)", "A2 (Auxiliary unipolar lead 2)",
    "A3 (Auxiliary unipolar lead 3)", "A4 (Auxiliary unipolar lead 4)",
    "V8-cal",							 /* 75 - 79 */
    "V9-cal", "V8R-cal", "V9R-cal", "D-cal (Nehb - Dorsal)",
    "A-cal (Nehb - Anterior)",					 /* 80 - 84 */
    "J-cal (Nehb - Inferior)"					 /* 85 */
};

int section3(unsigned char *p, long len)
{
    int i = 0, nsmax = 0;

    p += 16; len -= 16; /* move to data area */

    if (len < 2) {
	if (vflag) printf("  Error: section3 contains no data\n");
	return (0);
    }

    nleads = *p++;
    if (vflag) printf("  Number of leads: %d\n", nleads);
    len--;

    if ((nsamp=(unsigned long *)malloc((nleads+4)*sizeof(unsigned long))) ==
	NULL){
	if (vflag) printf("  Error: too many (%d) leads\n", nleads);
	return (0);
    }

    if (vflag) printf("  Flag byte: %xH\n", *p); len--;
    if (*p & 1) {
	if (vflag)
	    printf("   Reference beat subtraction used for compression\n");
    }
    else
	if (vflag)
	    printf("   Reference beat subtraction not used for compression\n");
    if (*p & 4) {
	if (vflag) printf("   All leads recorded simultaneously\n");
    }
    else
	if (vflag) printf("   %d leads recorded simultaneously\n", (*p >>3) & 0xf);
    p++;

    while (len > 8) {
	int j = *(p+8);
	if (j < sizeof(leadname)/sizeof(char *)) {
	    if (vflag) printf("  %s", leadname[j]);
	    /* yes, there are more efficient ways to do this -- but this one
	       won't break */
	    if      (strcmp(  "I", leadname[j]) == 0) leadI   = i;
	    else if (strcmp( "II", leadname[j]) == 0) leadII  = i;
	    else if (strcmp("III", leadname[j]) == 0) leadIII = i;
	    else if (strcmp("aVR", leadname[j]) == 0) aVR     = i;
	    else if (strcmp("aVF", leadname[j]) == 0) aVF     = i;
	    else if (strcmp("aVL", leadname[j]) == 0) aVL     = i;
	    else if (strcmp( "V1", leadname[j]) == 0) V1      = i;
	    else if (strcmp( "V2", leadname[j]) == 0) V2      = i;
	    else if (strcmp( "V3", leadname[j]) == 0) V3      = i;
	    else if (strcmp( "V4", leadname[j]) == 0) V4      = i;
	    else if (strcmp( "V5", leadname[j]) == 0) V5      = i;
	    else if (strcmp( "V6", leadname[j]) == 0) V6      = i;
	}
	else if (*(p+8) < 100) {
	    if (vflag) printf("  <Reserved lead>");
	}
	else
	    if (vflag) printf("  <Manufacturer-defined lead>");
	if (vflag) printf(": ID %d, samples %lu to %lu\n", *(p+8),
	       get32(p), get32(p+4));
	nsamp[i++] = get32(p+4) - get32(p) + 1;
	p += 9; len -= 9;
    }

    /* allocate space for 4 extra signals that will be calculated later,
       and for a NULL pointer to signal the end of the array */
    if ((ecg = (short **)malloc((nleads+5)*sizeof(short *))) == NULL) {
	if (vflag) printf("  Error: too many (%d) leads\n", nleads);
	return (0);
    }

    for (i = 0; i < nleads+4; i++) {
	if (i >= nleads) nsamp[i] = nsmax;
	else if (nsamp[i] > nsmax) nsmax = nsamp[i];
	if (nsamp[i] < 1) nsamp[i] = 1;
	if ((ecg[i] = (short *)calloc(nsamp[i], sizeof(short))) == NULL) {
	    if (vflag) printf("  Error: too many (%d) leads\n", nleads);
	    return (0);
	}
    }
    ecg[nleads+4] = NULL;

    if (len < 0) return (0);
    return (1);
}

int section4(unsigned char *p, long len)
{
    unsigned short i, stat = 1;

    p += 16; len -= 16; /* move to data area */
    if (len < 6) {
	if (vflag) printf("  Error: section 4 contains no data\n");
	return (0);
    }
    rblenms = get16(p);
    fcM = get16(p+2);
    nqrs = get16(p+4);
    if (vflag) {
	printf("  Length of reference beat type 0 data: %ld msec\n", rblenms);
	printf(
          "  Sample # of fiducial relative to start of ref beat type 0: %d\n",
	  fcM);
	printf("  Number of QRS complexes in the entire record: %d\n", nqrs);
    }

    if (nqrs < 1)
	return (0);	/* no beats -- nothing else to do here */

    p += 6; len -= 6;
    if ((subz=(struct subzone *)malloc(sizeof(struct subzone)*nqrs)) == NULL) {
	fprintf(stderr, "  Error: too many (%d) beats\n", nqrs);
	return (0);
    }

    if (len >= nqrs*14) {
	if (vflag) printf("  Reference beat subtraction zones:\n");
	for (i = 0; i < nqrs && len >= 14; i++) {
	    subz[i].type = get16(p);
	    subz[i].t0 = get32(p+2);
	    subz[i].t1 = get32(p+6);
	    subz[i].t2 = get32(p+10);
	    if ((subz[i].type == 0) && (subz[i].t2 == 0))
	      subz[i].type = 9;

	    if (vflag) {
		printf(
		   "   %2d  Type: %2d  start: %7ld"
		   "  fiducial: %7ld  end: %7ld\n",
		   i+1, subz[i].type, subz[i].t0, subz[i].t1, subz[i].t2);
		if (subz[i].type != 0 && (subz[i].t0 != 0 ||
					  subz[i].t2 != 0))
		    printf("  Error: start and end should be zero for beats "
			   "of non-zero types\n");
	    }
	    p += 14; len -= 14;
	}
    }
    else {
	if (vflag)
	    printf("  Error: reference beat subtraction zones are missing\n");
	stat = 0;
    }
    if (len >= nqrs*8) {
	if (vflag) printf("  QRS locations:\n");
	for (i = 0; i < nqrs; i++) {
	    if (vflag) printf("   %2d  start: %7lu  end: %7lu\n",
		   i+1, get32(p), get32(p+4));
	    p += 8; len -= 8;
	}
    }
    else {
	if (vflag) printf("  Error: QRS locations are missing\n");
	stat = 0;
    }
    return (stat);
}

int section5(unsigned char *p, long len)
{
    int i;
    p += 16; len -= 16; /* move to data area */


    if (len <= 0) {
	if (vflag) printf(" Error: section 5 contains no data\n");
	return (0);
    }

    /* Calculate length of template in samples */
    rblen = rblenms * 1000.0 / get16(p+2);

    if ((refbeat = (short **)calloc((nleads+1), sizeof(short *))) == NULL) {
	fprintf(stderr, "  Error: too many (%d) leads\n", nleads);
	return (0);
    }

    for (i = 0; i < nleads; i++) {
	if ((refbeat[i] = (short *)calloc(rblen, sizeof(short))) == NULL) {
	    fprintf(stderr, "  Error: too many (%d) leads\n", nleads);
	    return (0);
	}
    }

    i = Huffdecode(p, len, 5);
    return (i);
}


struct htentry {
    unsigned short pattern;
    unsigned short suffixbits;
    short basevalue;
} htdefaultentries[] = {
    {     0,  0,  0 },
    {     4,  0,  1 },
    {     5,  0, -1 },
    {   0xc,  0,  2 },
    {   0xd,  0, -2 },
    {  0x1c,  0,  3 },
    {  0x1d,  0, -3 },
    {  0x3c,  0,  4 },
    {  0x3d,  0, -4 },
    {  0x7c,  0,  5 },
    {  0x7d,  0, -5 },
    {  0xfc,  0,  6 },
    {  0xfd,  0, -6 },
    { 0x1fc,  0,  7 },
    { 0x1fd,  0, -7 },
    { 0x3fc,  0,  8 },
    { 0x3fd,  0, -8 },
    { 0x3fe,  8,  0 },
    { 0x3ff, 16,  0 }
};

struct htable {
    int npatterns;
    struct htentry *entry;
} htdefault = { 19, htdefaultentries };

struct htable *htab = &htdefault;

/* This function reads bits from the input buffer, beginning 'ibits' bits
   from the MSB of *p, until it recognizes a Huffman-encoded value.  The
   decoded value is then stored into x, and the function returns the number
   of bits it has read. */
unsigned long getbits(short *x, unsigned char *p, unsigned long ibits)
{
    unsigned short bitbuffer = 0L;
    static unsigned char w;
    unsigned int nbits = 0;  /* number of bits read on this call to getbits */
    short i, j, v;

    if (ibits == 0) {	/* first call for this lead -- initialize */
	;
    }

    do {
	bitbuffer <<= 1;	/* make room for another bit */

	if ((ibits & 7) == 0)	/* once every 8 bits ... */
	    w = p[ibits / 8];	/* ... get the next 8 bits */
	/* copy next bit from w into bitbuffer ... */
	if (w & 0x80) bitbuffer |= 1;
	/* ... and discard it from w */
	w <<= 1;

	ibits++;
	nbits++;

	/* check bitbuffer against all of the Huffman table entries */
	for (i = 0; i < htab->npatterns; i++) {
	    if (bitbuffer == htab->entry[i].pattern) {
		/* we have a match */
		if ((j = htab->entry[i].suffixbits) == 0)
		    *x = htab->entry[i].basevalue;
		else if (j < 0) {
		    /* table mode switch */
		    /* ** not implemented ** */
		    ;
		}
		else {	/* j is the number of additional bits we need */
		    /* as above, get 8 more bits if needed */
		    if ((ibits & 7) == 0)
			w = p[ibits / 8];
		    /* Get first (most significant) bit of value.  This
		       bit is 0 if the value is positive or zero, and 1
		       if the value is negative.  If the value is negative,
		       and the number of additional bits is less than the
		       number of bits in a 'short', we need to sign-extend
		       the number by shifting additional ones into it before
		       reading the rest of the value from the bit stream. */
		    if (w & 0x80) {	/* if it's a 1, sign-extend it */
			int jj = j;
			for (v = 1; jj < 8 * sizeof(short); jj++) {
			    v <<= 1;	/* make room for another bit */
			    v |= 1;	/* shift in another 1 */
			}
		    }
		    else v = 0;
		    w <<= 1;
		    ibits++;
		    nbits++;
		    for (j--; j > 0; j--, ibits++, nbits++) {
			v <<= 1;
			if ((ibits & 7) == 0)
			    w = p[ibits / 8];
			if (w & 0x80) v |= 1;
			w <<= 1;
		    }
		    *x = v;
		}
		return (nbits);
	    }
	}
    } while (nbits < sizeof(unsigned long) * 8 - 1);
    if (vflag) printf("Huffman encoding error\n");	
    return (1000000L);
}

int section6(unsigned char *p, long len)
{
    int i;

    p += 16; len -= 16; /* move to data area */

    if (len <= 6) {
	if (vflag) printf(" Error: section 6 contains no data\n");
	return (0);
    }
    i=Huffdecode(p, len, 6);
    return (i);
}

int Huffdecode(unsigned char *p, long len, int section)
{
    int i, comptype = 1, enctype = 2, sampint;
    unsigned long ibits, tbits, *leadlen, tleadlen = 0L;
    double g = get16(p);

    if (section == 5) gref = g / 5000.;
    else if (section == 6) gres = g / 5000.;
    else {
	fprintf(stderr, "oops!  shouldn't call Huffdecode with section = %d\n",
	       section);
	exit(2);
    }
    if (vflag) printf("  Amplitude value multiplier: %g nanovolts per unit\n",
		      g);
    sampint = get16(p+2);
    if (vflag) printf("  Sample time interval: %d microseconds\n", sampint);
    if (sampint > 0) {
	if (sfreq > 0 && sfreq != 1.0e6/sampint)
	    printf("Warning: section 5 and 6 sample intervals don't match");
	    /* FIXME: If this ever happens, we need to redo the reconstruction
	       code in main() */
	sfreq = 1.0e6/sampint;
    }
    else {
	printf("Error: sampling frequency not specified;  assuming 500 Hz\n");
	sfreq = 500;
    }
    if (vflag) printf("  Encoding type: %d ", *(p+3));
    switch (*(p+3)) {
    case 0: if (vflag) printf("(amplitudes)\n"); enctype = 0; break;
    case 1: if (vflag) printf("(first differences)\n"); enctype = 1; break;
    case 2: if (vflag) printf("(second differences)\n"); enctype = 2; break;
    default: if (vflag) printf("<undefined, assuming second differences>\n"); enctype = 2;
	break;
    }
    if (section == 6) {
	if (vflag) printf("  Compression type: %d ", *(p+4));
	switch (*(p+4)) {
	case 0: if (vflag) printf("(non-bimodal)\n"); comptype = 0; break;
	case 1: if (vflag) printf("(bimodal)\n"); comptype = 1; break;
	default: if (vflag) printf("<undefined, assuming bimodal>\n"); comptype = 1;
	    break;
	}
    }
    p += 6; len -= 6;

    if (nleads < 0) {
	if (vflag) printf("  Error: number of leads undefined -- section 3 missing"
	       " or out of order\n");
	return (0);
    }
    if (nleads == 0) {
	if (vflag) printf("  Error: no leads recorded\n");
	return (0);
    }
    if ((leadlen=(unsigned long *)malloc(nleads*sizeof(unsigned long)))==NULL){
	if (vflag) printf("  Error: too many (%d) leads\n", nleads);
	return (0);
    }

    for (i = 0; i < nleads; i++, p += 2, len -= 2) {
	tleadlen += leadlen[i] = get16(p);
	if (vflag) printf("  Lead %d:%5lu bytes\n", i, leadlen[i]);
    }

    if (len < tleadlen) {
	if (vflag) printf("  Error: %lu data byte%s missing from section %d\n",
			  tleadlen - len,
			  (tleadlen - len == 1) ? "" : "s",
			  section);
	free(leadlen);
	return (0);
    }
    else if (len > ((tleadlen + 1) & ~1))
	/* the expression above rounds tleadlen up to an even number */
	if (vflag) printf("  Warning: %lu extra byte%s in section %d\n",
			  len - tleadlen,
			  (len - tleadlen == 1) ? "" : "s",
			  section);

    for (i = 0; i < nleads; i++) {
	char fname[6];
	int explen, j;
	static short dx, x, xx, xxx, xout;

	if (vflag) printf("  Lead %d data\n", i);
	for (j=0; j < 10; j++)
	    if (vflag) printf(" %02x", p[j]);
	if (vflag) printf("\n");
	xxx = xx = 0;
	if (section == 6) explen = nsamp[i];
	else explen = rblen;
	for (j=0, ibits=0, tbits=leadlen[i]*8; j<explen && ibits<tbits; j++ ) {
	    ibits += getbits(&x, p, ibits);
	    if (j == 0) xxx = xx = x;
	    else if (j == 1) { dx = x - xxx;  xxx = x; xx += x; }
	    else { xxx += dx += x; xx += x; }
	    switch (enctype) {
	    case 0: xout = x; break;
	    case 1: xout = xx; break;
	    case 2: xout = xxx; break;
	    }
	    if (section == 5) refbeat[i][j] = xout * gref;
	    else ecg[i][j] = xout * gres;
	    if (vflag > 1) printf("   %5d %5d %5d %5d\n", j, x, xx, xxx);
	}
	if (j < explen)
	    if (vflag) printf("   Error: %ld sample%s missing in lead %d\n",
			      leadlen[i] - j,
			      (leadlen[i] - j == 1) ? "" : "s",
			      i);
	if (ibits < tbits)
	    if (vflag) printf("   Warning: %ld extra bit%s in lead %d\n",
			      tbits - ibits,
			      (tbits - ibits == 1) ? "" : "s",
			      i);
	p += leadlen[i];
	if (vflag) printf("\n");
    }
    free(leadlen);
    return (1);
}

int section7(unsigned char *p, long len)
{
    unsigned char nmb, nps, sec7qrs, *q, *r;
    short i, n, sec7qrs_offset;

    p += 16; len -= 16; /* move to data area */

    if (len <= 0) {
	if (vflag) printf(" Error: section 7 contains no data\n");
	return (0);
    }

    if (xflag && vflag) hexdump(p, len);

    if (len < 20) {                  /* there isn't enough room left for measurements */
	if (vflag) printf(" Error: section 7 is too short (len = %ld)\n", len);
	return (0);
    }

    RRint = get16(p+2);

/*  get P onset, etc from first measurement block */
    q = p + 6;
    Pon = get16(q);
    Pend = get16(q+2);
    QRSon = get16(q+4);
    QRSend = get16(q+6);
    Tend = get16(q+8);
    Paxis = get16(q+10);
    QRSaxis = get16(q+12);
    Taxis = get16(q+14);

    nmb = *p;
    nps = *(p+1);
    sec7qrs_offset = 6 + 16*nmb + 10*nps;
    if (sec7qrs_offset > len) sec7qrs_offset = len;

    if (sec7qrs_offset > len) {
	if (vflag) printf(" Error: section 7 is too short"
			  " (nmb = %d, nps = %d)\n", nmb, nps);
	return (0);
    }

    sec7qrs = get16(p + sec7qrs_offset);

    if (vflag) {
	printf(" QRS complexes measured: %d\n", sec7qrs);
	printf(" Measurement blocks: %d\n", nmb);
	if (nmb == sec7qrs + 1)
	    printf("  (First measurement block contains measurements for"
		   " reference beat type 0;\n"
		 "   Others contain measurements for each individual beat)\n");
	else
	    printf("  (Measurements for each reference beat type)\n");
	printf("  Number of measurement blocks: %d\n", nmb);
	printf("  Number of pacemaker spikes: %d\n", nps);
	printf("  Mean RR interval: %d ms\n", RRint);
	printf("  Mean PP interval: %d ms\n", get16(p+4));
	for (i = 0, q = p+6; i < nmb; i++, q += 16) {
	    printf("   Block %d\n", i);
	    printf("    P onset:   %d ms\n", get16(q));
	    printf("    P end:     %d ms\n", get16(q+2));
	    printf("    QRS onset: %d ms\n", get16(q+4));
	    printf("    QRS  end:  %d ms\n", get16(q+6));
	    printf("    T end:     %d ms\n", get16(q+8));
	    if ((n = get16(q+10)) == 999) printf("    P axis:    undefined\n");
	    else printf("    P axis:    %d degrees\n", n);
	    if ((n = get16(q+12)) == 999) printf("    QRS axis:  undefined\n");
	    else printf("    QRS axis:  %d degrees\n", n);
	    if ((n = get16(q+14)) == 999) printf("    T axis:    undefined\n");
	    else printf("    T axis:    %d degrees\n", n);
	}

	if (nps > 0)
	    printf(" Pacemaker spike measurement data\n");
	for (i = 0, r = q + nps*4; i < nps; i++, q += 4, r += 6) {
	    printf("Spike %d: %u ms, %d microvolts ",
		   i, (unsigned)get16(q), get16(q+2));
	    switch (*r) {
	    case 0: printf("(unknown spike type, "); break;
	    case 1: printf("(triggers neither P nor QRS, "); break;
	    case 2: printf("(triggers a QRS, "); break;
	    case 3: printf("(triggers a P, "); break;
	    default: printf("(spike type %d, ", *r); break;
	    }
	    switch (*(r+1)) {
	    case 0: printf("source unknown, "); break;
	    case 1: printf("internal, "); break;
	    case 2: printf("external, "); break;
	    default: printf("source type %d, ", *(r+1)); break;
	    }
	    if ((n = get16(r+2)) == 0)
		printf(" not linked to a QRS, ");
	    else
		printf(" linked to QRS #%d, ", n);
	    if ((n = get16(r+4)) == 0)
		printf("width unknown)\n");
	    else
		printf("%d microseconds\n)", n);
	}


	/* The next part doesn't seem to be right ... */
	printf(" <!> sec7qrs = <%d>\n", sec7qrs);
        if (sec7qrs > 0)
	    printf(" QRS type information (%d beats)\n", get16(r));
	for (i = 0, q = r+2; i < sec7qrs; i++, q++)
	    printf("  QRS %d: reference beat type %d\n", i, *q);
#if 0
	printf(" Additional global measurements\n");
	printf("  Ventricular rate: %d bpm\n", get16(q));
	printf("  Atrial rate: %d bpm\n", get16(q+2));
	printf("  QTc: %d ms ", get16(q+4));
	switch (*(q+6)) {
	case 0: printf("(unspecified correction formula)\n"); break;
	case 1: printf("(Bazett correction)\n"); break;
	case 2: printf("(Hodges correction)\n"); break;
	default: printf("(correction type %d)\n", *(q+6)); break;
	}
#endif


	printf("%ld bytes remaining in section 7\n", len - (q - p));
	while (q < p + len) {
	    printf(" %5d", get16(q));
	    printf(" [%3d]", *q++);
	}	printf("\n");

    }

    /* ** more to do here ** */
    return (1);
}

int section8(unsigned char *p, long len)
{
    int i, n;

    p += 16; len -= 16; /* move to data area */

    if (len <= 0) {
	if (vflag) printf(" Error: section 8 contains no data\n");
	return (0);
    }

    if (vflag) printf(" Status: %d ", *p);
    switch (*p) {
    case 0: if (vflag) printf("(original report, not confirmed)\n"); break;
    case 1: if (vflag) printf("(confirmed report)\n"); break;
    case 2: if (vflag) printf("(overread report, but not confirmed)\n"); break;
    default: if (vflag) printf("<undefined>\n"); break;
    }

    if (vflag) printf(" Date: %4d/%02d/%02d\n", get16(p+1), *(p+3), *(p+4));
    if (vflag) printf(" Time: %02d:%02d:%02d\n", *(p+5), *(p+6), *(p+7));
    n = *(p+8);
    stmnts = n;
    if (vflag) printf(" Number of statements: %d\n", n);
    p += 9; len -= 9;

    for (i = 0; i < n && len > 2; i++) {
	char format[16];
	int slen;

	report[i] = p+3;
	slen = get16(p+1);
	if (vflag) printf("  Statement %2d: (%d bytes)\n", *p, slen);
	sprintf(format, "   [%%%ds]\n", slen);
	if (vflag) printf(format, p+3);
	p += slen+3; len -= slen+3;
    }

    if (i < n) {
	if (vflag) printf("  Error: data missing in section 8\n");
	return (0);
    }

    return (1);
}

int section9(unsigned char *p, long len)
{
    p += 16; len -= 16; /* move to data area */

    if (len <= 0) {
	if (vflag) printf(" Error: section 9 contains no data\n");
	return (0);
    }

    /* FIXME: more to do here */
    if (vflag) {
	printf("  Data:\n   ");
	while (len-- > 0) {
	    if (isprint(*p)) putchar(*p);
	    else printf("\\%03o", *p);
	    p++;
	}
	printf("\n");
    }

    return (1);
}

int nsecleads;		/* number of leads with measurements in section 10 */

void prval(FILE *ofile, short val)
{
    if (val == 29999)
	fprintf(ofile, "  n/comp");	/* not computed (ever) by analyzer */
    else if (val == 29998)
	fprintf(ofile, "     rej");	/* lead rejected, not computed */
    else if (val == 19999)
	fprintf(ofile, "  absent");	/* correponding wave not present */
    else
	fprintf(ofile, "%8d", val);
}

void precgmeas(FILE *ofile)
{
    int i;

    fprintf(ofile, "Parameter  Units");
    for (i = 0; i < nsecleads; i++) {
	if (0 < ecgmeas[i].leadid && ecgmeas[i].leadid <= 85)
	    fprintf(ofile, "%8s", leadname[ecgmeas[i].leadid]);
	else
	    fprintf(ofile, "  <%4d>", ecgmeas[i].leadid);
    }
    fprintf(ofile, "\nP duration    ms");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].Pdur);
    fprintf(ofile, "\nPR interval   ms");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].PRint);
    fprintf(ofile, "\nQRS duration  ms");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].QRSdur);
    fprintf(ofile, "\nQT interval   ms");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].QTint);
    fprintf(ofile, "\nQ duration    ms");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].Qdur);
    fprintf(ofile, "\nR duration    ms");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].Rdur);
    fprintf(ofile, "\nS duration    ms");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].Sdur);
    fprintf(ofile, "\nR' duration   ms");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].Rpdur);
    fprintf(ofile, "\nS' duration   ms");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].Spdur);
    fprintf(ofile, "\nQ amplitude   uV");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].Qamp);
    fprintf(ofile, "\nR amplitude   uV");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].Ramp);
    fprintf(ofile, "\nS amplitude   uV");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].Samp);
    fprintf(ofile, "\nJ amplitude   uV");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].Jamp);
    fprintf(ofile, "\nP+ amplitude  uV");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].Ppamp);
    fprintf(ofile, "\nP- amplitude  uV");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].Pnamp);
    fprintf(ofile, "\nT+ amplitude  uV");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].Tpamp);
    fprintf(ofile, "\nT- amplitude  uV");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].Tnamp);
    fprintf(ofile, "\nST slope    uV/s");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].STslope);
    fprintf(ofile, "\nP morphology  --");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].Pmorph);
    fprintf(ofile, "\nT morphology  --");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].Tmorph);
    fprintf(ofile, "\nIsoelec (I)   ms");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].isoI);
    fprintf(ofile, "\nIsoelec (K)   ms");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].isoK);
    fprintf(ofile, "\nIntrins defl  ms");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].tind);
    fprintf(ofile, "\nQuality       --");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].qual);
    /* FIXME: SCP spec doesn't specify units of st20, 60, 80, rr16, rr8 --
       assume they are microvolts */
    fprintf(ofile, "\nST (J+20ms)   uV");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].st20);
    fprintf(ofile, "\nST (J+60ms)   uV");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].st60);
    fprintf(ofile, "\nST (J+80ms)   uV");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].st80);
    fprintf(ofile, "\nST (J+RRm/16) uV");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].strr16);
    fprintf(ofile, "\nST (J+RRm/8)  uV");
    for (i = 0; i < nsecleads; i++)
	prval(ofile, ecgmeas[i].strr8);
    fprintf(ofile, "\n");
}

int section10(unsigned char *p, long len)
{
    int i;
    FILE *ofile;

    p += 16; len -= 16; /* move to data area */

    if (len <= 0) {
	if (vflag) printf(" Error: section 10 contains no data\n");
	return (0);
    }

    /* Section 10 data area begins with a 4-byte header.  The first 2 bytes
       specify the number of leads for which measurements are stored. */
    nsecleads = get16(p);
    /* FIXME:  the draft standard doesn't say what the next 2 bytes are --
       so they are ignored here. */

    /* Allocate structures to hold all of the parsed measurements. */
    ecgmeas = (struct ECGmeas *)calloc(nsecleads,  sizeof(struct ECGmeas));
    if (ecgmeas == NULL) {
	fprintf(stderr, "Error: too many (%d) leads\n", nsecleads);
	return (0);
    }
    
    /* Fill the array of structures. */
    for (i = 0, p += 4; i < nsecleads; i++) {
	ecgmeas[i].leadid = get16(p);
	ecgmeas[i].vlen = get16(p+2);
	ecgmeas[i].Pdur = get16(p+4);
	ecgmeas[i].PRint = get16(p+6);
	ecgmeas[i].QRSdur = get16(p+8);
	ecgmeas[i].QTint = get16(p+10);
	ecgmeas[i].Qdur = get16(p+12);
	ecgmeas[i].Rdur = get16(p+14);
        ecgmeas[i].Sdur = get16(p+16);
	ecgmeas[i].Rpdur = get16(p+18);
	ecgmeas[i].Spdur = get16(p+20);
	ecgmeas[i].Qamp = get16(p+22);
	ecgmeas[i].Ramp = get16(p+24);
	ecgmeas[i].Samp = get16(p+26);
	ecgmeas[i].Rpamp = get16(p+28);
	ecgmeas[i].Spamp = get16(p+30);
	ecgmeas[i].Jamp = get16(p+32);
	ecgmeas[i].Ppamp = get16(p+34);
	ecgmeas[i].Pnamp = get16(p+36);
	ecgmeas[i].Tpamp = get16(p+38);
	ecgmeas[i].Tnamp = get16(p+40);
	ecgmeas[i].STslope = get16(p+42);
	ecgmeas[i].Pmorph = get16(p+44);
	ecgmeas[i].Tmorph = get16(p+46);
	ecgmeas[i].isoI = get16(p+48);
	ecgmeas[i].isoK = get16(p+50);
	ecgmeas[i].tind = get16(p+52);
	ecgmeas[i].qual = get16(p+54);
	ecgmeas[i].st20 = get16(p+56);
	ecgmeas[i].st60 = get16(p+58);
	ecgmeas[i].st80 = get16(p+60);
	ecgmeas[i].strr16 = get16(p+62);
	ecgmeas[i].strr8 = get16(p+64);
	p += ecgmeas[i].vlen + 4;
    }

    if (vflag) {
	printf(" Results for %d leads\n", nsecleads);
	precgmeas(stdout);                                   
    }
    return (1);
}

int section11(unsigned char *p, long len)
{
    p += 16; len -= 16; /* move to data area */

    if (len <= 0) {
	if (vflag) printf(" Error: section 11 contains no data\n");
	return (0);
    }

    /* ** more to do here ** */
    return (1);
}

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);
}