/* file: hp_scatter.c
/****************************************************************
# This program replaces premature ventricular (V)-beats
# by normal (N) beats, i.e. it interpolates normal-normal (N-N) time intervals.
# These N-N time intervals are averaged over a number of beats,
# typically 4 to 10 beats.
# given by the variable window.
# Then, for each V-beat we print the
# - coupling interval:
# interval between V-beat and previous normal beat, if the previous beat
# was not a normal beat, a zero is printed
# - V-V time interval:
# interval between consecutive V-beats
# - Number of Intervening normal Beats (NIB):
# number of normal beats appearing between consecutive V-beats.
# - corresponding interpolated and averaged N-N interval
# The output is the file file_root.scatter, wich contains
# first the histogram of the interpolated N-N intervals in 2 columns,
# then, in five columns, coupling interval, VV interval, NIB,
# and corresponding N-N interval
#
# The input needs the format:
# column 1: Annotation V=V-beat, N=normal, and others
# column 2: RR-interval
#
# command-line: hp_scatter file_name NN_min NN_max binsize max
# > file_name.scatter
#
# where
# binsize = 1/sampling frequency
# NN_min = minimum of NN-interval range under consideration
# NN_max = maximum of NN-interval range under consideration
# max = maximum number of data points
# These values have to be taken from the data, e.g. by plotting
# the RR-intervals.
#
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# The interpolation of the V-beats by normal beats can be tricky if there
# are many V-beats, especially many consecutive V-beats (runs of v-tach)
# or interpolated beats (i.e. no compensatory pause).
# Note, that we replace each V-beat by a normal beat since we are only
# interested in the value of the time interal between the underlying normal
# beats. This might generate more normal beats than were present originally
# if long runs of ventricular tachycardia are present in the data.
#
# We consider the following cases:
# - isolated V-beats followed by compensatory pause (N-V-N):
# V-beat is replaced by normal beat which is placed half way [(NV+VN)/2]
# between the surrounding N-beats.
# - isolated V-beat without compensatory pause:
# V-beat should be just taken out, while the N-N interval results
# from adding NV+VN. But in order to keep the number of beats unchanged,
# we replace both, the NV-interval and the VN-interval, by NV+VN.
# - two consecutive V-beats blocking two normal beats
# (= with compensatory pause)
# - two consecutive V-beats blocking one normal beat
# (= without compensatory pause)
# - three consecutive V-beats blocking two or three normal beats
# (= with compensatory pause)
# - three consecutive V-beats blocking one normal beat
# (= without compensatory pause)
# The identification of these cases involves thresholds that might have
# to to be adjusted to the data at hand. Also, other cases might be
# possible that is not dealt with presently. The result of the interpolation
# should therefore be checked by printing out RRneu!!!
**********************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define NUM_OF_ARGC 6
#define MAX_CHAR_IN_LINE 4096
#define STRL 80
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
main ( int argc, char *argv[] )
{
FILE *fp, *ofp;
double binsize, NN_min, NN_max;
double *RR, *RRneu, *RRave, *ttime;
double Voldbeat, x;
long *hist_RRave;
char *type;
long i,l,k,jj,j, max, min, steps, window, Voldindex, NIB;
double coupint, comp_pause, Ts, Vtime;
long round(double );
char dummy;
if ( argc!=NUM_OF_ARGC )
{
printf( "Usage: %s RRLIST NNmin NNmax BINSIZE LENGTH\n\n",
argv[0]);
printf( "argc = %d\n", argc);
for (i = 1; i < argc; i++)
printf("argv[%d] = %s\n", i, argv[i]);
printf( "(steps= (NN_max-NN_min)/binsize must be integer \n\n" );
exit(0);
}
NN_min = atof(argv[2]); /* minimum of NN-interval range */
NN_max = atof(argv[3]); /* maximum of NN-interval range */
binsize = atof(argv[4]); /* 1/sampling frequency */
max = atol(argv[5]); /* maximum number of data points */
Voldbeat=0;
min= 1;
steps=round((NN_max-NN_min)/binsize);
window=10;
RR = (double *) malloc( (max+1)*sizeof(double) );
for (i=0; i<=max; i++)
RR[i]=0.0;
RRneu = (double *) malloc( (max+1)*sizeof(double) );
for (i=0; i<=max; i++)
RRneu[i]=0.0;
RRave = (double *) malloc( (max+1)*sizeof(double) );
for (i=0; i<=max; i++)
RRave[i]=0.;
type = (char *) malloc( (max+1)*sizeof(char) );
for (i=0; i<=max; i++)
type[i]=0;
ttime = (double *) malloc( (max+1)*sizeof(double) );
for (i=0; i<=max; i++)
ttime[i]=0.0;
hist_RRave = (long *) malloc( (steps+1)*sizeof(long) );
for (i=0; i<=steps; i++)
hist_RRave[i]=1;
fp = fopen( argv[1], "r" );
i = 1;
/* while ((i<max)&&(fscanf(fp, "%c %lf\n", &type[i], &RR[i]) == 2) )*/
/* this is how is used to be */
while ((i<max)&&(fscanf(fp, "%lf %c\n", &RR[i], &type[i]) == 2) ) {
if (type[i] == 'r' || type[i] == '!')
type[i] = 'V'; /* count R-on-T PVCs and ventricular flutter waves as
V-beats */
i++;
}
max=i;
fclose( fp );
for( i=min; i<=max;i++)
{
ttime[i] = ttime[i-1] + RR[i];
}
/****************************
* replace V-beats:
* this part may have to be adapted to the data at hand!!!!
* note that the number of beats is not changed
****************************/
for (i=min+1;i<max;i++)
{
/*** if not a V-beat, just copy beat ***/
if (type[i]!='V')
{
RRneu[i]=RR[i];
}
/*** if isolated V-beat check if compensatory pause ***/
if ((type[i]=='V')&&(type[i+1]!='V'))
{
x=(RR[i]+RR[i+1])/2;
if (x > 2*RRneu[i-1]/3)
{
RRneu[i]=RRneu[i+1]=x;
i++;
}
else if (x < 2*RRneu[i-1]/3)
{
RRneu[i]=RRneu[i+1]=RR[i]+RR[i+1];
i++;
}
}
/*** if 2 consecutive V-beats check if compensatory pause ***/
if ((type[i]=='V')&&(type[i+1]=='V')&&(type[i+2]!='V'))
{
x=(RR[i]+RR[i+1]+RR[i+2])/3;
if (x > 3*RRneu[i-1]/4)
{
RRneu[i]=RRneu[i+1]=RRneu[i+2]=x;
i+=2;
}
else
{
RRneu[i]=RRneu[i+1]=RRneu[i+2]=(RR[i]+RR[i+1]+RR[i+2])/2;
i+=2;
}
}
/*** if 3 consecutive V-beats check if compensatory pause ***/
if ((type[i]=='V')&&(type[i+1]=='V')&&(type[i+2]=='V')&&(type[i+3]!='V'))
{
x=(RR[i]+RR[i+1]+RR[i+2]+RR[i+3])/4;
if (x > 3*RRneu[i-1]/4)
{
RRneu[i]=RRneu[i+1]=RRneu[i+2]=RRneu[i+3]=x;
i+=3;
}
else
{
x=(RR[i]+RR[i+1]+RR[i+2]+RR[i+3])/3;
if (x > 3*RRneu[i-1]/4)
{
RRneu[i]=RRneu[i+1]=RRneu[i+2]=RRneu[i+3]=x;
i+=3;
}
else
{
x=(RR[i]+RR[i+1]+RR[i+2]+RR[i+3])/2;
RRneu[i]=RRneu[i+1]=RRneu[i+2]=RRneu[i+3]=x;
i+=3;
}
}
}
/*** if more then 3 consecutive v-beats or 2 v-beats
*** and then not a normal beat fill with interpolated
*** NN-interval of previous beat until normal beat appears ***/
else if ((type[i]=='v')&&(type[i+1]=='v'))
{
while (type[i]!='n')
{
RRneu[i]=RRneu[i-1];
i++;
}
RRneu[i]=RRneu[i-1];
}
}
/*** print reconstructed normal sinus beats on file: **/
ofp = fopen("NN_intervals", "w" );
for (i=1;i<max;i++)
fprintf(ofp,"%d %f %f\n", i, RRneu[i], RR[i]);
/*** calculate average and RR_ave histogram: ***/
for (i=min;i<min+window/2;i++)
RRave[i]=RRneu[i];
for (i=min+1;i<=min+window;i++)
RRave[window/2+min]+=RRneu[i];
for (i=window+min+1;i<=max;i++)
RRave[i-window/2]=RRave[i-window/2-1]-RRneu[i-window]+RRneu[i];
for (i=window+min;i<=max;i++)
{
RRave[i]= RRave[i]/window;
if ((RRave[i]<=NN_max)&&(RRave[i]>=NN_min))
{
hist_RRave[round((RRave[i]-NN_min)/binsize)]++;
}
}
for (i=0;i<=steps;i++)
printf("%f %d\n", i*binsize+NN_min, hist_RRave[i]);
/* calculate and print calculate and print vv-intervals, NIB and coupling intervals:*/
Voldindex=1;
for (i=min;i<max;i++)
{
if (type[i]=='V')
{
if (NIB!=0) /* V-V interval is not a coupling interval! */
coupint=RR[Voldindex];
else
coupint=0;
comp_pause=RR[Voldindex+1];
NIB=i-Voldindex-1;
Vtime=ttime[i]-ttime[Voldindex];
/* for interpolated beats you may want to subtract the NIB by one
* to match the NIB triplets in parasystole: */
if (((coupint+comp_pause)/2 < 2*RR[Voldindex-1]/3)&&(NIB>0))
NIB=NIB-1;
/* you can use the averaged RR-intervals or the non-averaged RR-intervals:*/
/* Ts=RRneu[Voldindex];*/
Ts=RRave[Voldindex];
if (Ts != 0)
{
printf("%f %f %d %f %d\n", coupint, Vtime, NIB, Ts, Voldindex);
}
Voldindex=i;
}
}
free(RR);
free(type);
free(RRneu);
free(RRave);
free(ttime);
free(hist_RRave);
}
int round(double x)
{
if (ceil(x)-floor(x) < 0.5)
return floor(x);
else
return ceil(x);
}