/* file: rrgen.c Entry number : 201 Authors : P. E. McSharry, G. D. Clifford Organization : Dept Maths & Dept Engineering, University of Oxford, UK email address : mcsharry at maths.ox.ac.uk, gari at robots.ox.ac.uk This program should compile without errors or warnings using: gcc -Wall rrgen.c -lm See http://www.physionet.org/challenge/2002/ for further information on the CinC Challenge 2002. This program was used to generate series rr19 and rr31 of the challenge dataset. */ /* NOTE - some algorithms from numerical recipes in C are used. If you have not done so, please purchase a set of disks from the appropriate source before copying/running this code */ #include #include #include #define SPS 128 /* sampling frequency (determines quantization of RR intervals; see main()) */ #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr #define MIN(a,b) (a < b ? a : b) #define MAX(a,b) (a > b ? a : b) #define PI (2.0*asin(1.0)) #define NR_END 1 #define FREE_ARG char* #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) long rseed; #define AWAKE 1.0 #define ASLEEP 0.0 #define YES 1 #define NO 0 int first_pass = 0; float last_rr = 0.0; /* intialise the first rr interval to be zero */ float last_Rt = 0.0; /* intialise the first R-peak time stamp to be zero */ float *trpeaks_new; long N; /* number of seconds */ long count = 1; /* global counter */ /*---------------------------------------------------------------------------*/ /* modify beat for ectopy */ /*---------------------------------------------------------------------------*/ float modify_ectopic(float rr2, float rr1) { /* shorten beat by 80% +/- 10% randomly chosen */ float new_rr; float delta; float rrint; delta = ((float)rand()/RAND_MAX); delta = 0.2*(delta-0.5); /* length of last rr interval */ rrint = rr2 - rr1; /* shorten it by 0.7 +/- a bit */ rrint = (0.7 - delta)*rrint; /* work out new timing */ new_rr = rr1 + rrint; return(new_rr); } /*---------------------------------------------------------------------------*/ /* generate an extra beat (artefact) */ /*---------------------------------------------------------------------------*/ float make_noise(float rr2, float rr1) { /* shorten beat by 80% +/- 10% randomly chosen */ float new_rr; float delta; float rrint; delta = ((float)rand()/RAND_MAX); delta = ((delta*(SPS-(2.0/SPS))) + (1.0/SPS))/SPS; /* length of last rr interval */ rrint = rr2 - rr1; /* shorten it by 0.5 +/- a bit */ rrint = delta*rrint; /* work out new timing */ new_rr = rr1 + rrint; return(new_rr); } /*--------------------------------------------------------------------------*/ /* CHECK IF ECTOPICS SHOULD BE ADDED */ /*--------------------------------------------------------------------------*/ int check_for_ectopic(void) { /* no state or hr dependency */ int rtn=NO; float rand_no; float Pe=0.0003; /* probability of ectopic (~1 per hour)*/ /* pick a random number between 0 and 1*/ rand_no = ((float)rand()/RAND_MAX); if(rand_no=0;j--) { k=(*idum)/IQ; *idum=IA*(*idum-k*IQ)-IR*k; if (*idum < 0) *idum += IM; if (j < NTAB) iv[j] = *idum; } iy=iv[0]; } k=(*idum)/IQ; *idum=IA*(*idum-k*IQ)-IR*k; if (*idum < 0) *idum += IM; j=iy/NDIV; iy=iv[j]; iv[j] = *idum; if ((temp=AM*iy) > RNMX) return RNMX; else return temp; } /*--------------------------------------------------------------------------*/ /* GAUSSIAN DEVIATES */ /*--------------------------------------------------------------------------*/ float gasdev(long *idum) { float ran1(long *idum); static int iset=0; static float gset; float fac,rsq,v1,v2; if (iset == 0) { do { v1=2.0*ran1(idum)-1.0; v2=2.0*ran1(idum)-1.0; rsq=v1*v1+v2*v2; } while (rsq >= 1.0 || rsq == 0.0); fac=sqrt(-2.0*log(rsq)/rsq); gset=v1*fac; iset=1; return v2*fac; } else { iset=0; return gset; } } /*--------------------------------------------------------------------------*/ /* FFT */ /*--------------------------------------------------------------------------*/ void four1(float data[], unsigned long nn, int isign) { unsigned long n,mmax,m,j,istep,i; double wtemp,wr,wpr,wpi,wi,theta; float tempr,tempi; n=nn << 1; j=1; for (i=1;i i) { SWAP(data[j],data[i]); SWAP(data[j+1],data[i+1]); } m=n >> 1; while (m >= 2 && j > m) { j -= m; m >>= 1; } j += m; } mmax=2; while (n > mmax) { istep=mmax << 1; theta=isign*(6.28318530717959/mmax); wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi=sin(theta); wr=1.0; wi=0.0; for (m=1;m dRRamp) { u = (float)rand()/RAND_MAX; dRR = 0.5*dRRamp*(1.0 + (float)exp(gasdev(&rseed))/10.0)*u/fabs(u); } rrmean = RR0[time] + dRR; rrtrend = RRtrendamp*(2*(float)rand()/RAND_MAX - 1); rrstd = RRstdmin + (RRstdmax - RRstdmin)*(float)rand()/RAND_MAX; lfhfratio = LFHFmin + (LFHFmax-LFHFmin)*(float)rand()/RAND_MAX; /* generate RR over this block */ trblock(tr0, flo, fhi, flostd, fhistd, lfhfratio, rrmean, rrtrend, rrstd, Nb); for(j=1;j<=Nb;j++) trpeaks[i+j] = trpeaks[i] + tr0[j]; i += Nb; /* convert beat time stamp back to nearest second */ time = (long)(trpeaks[i]); rra = tr0[Nb] - tr0[Nb-1]; // printf("time = %ld\n",time); /* for(j=2;j<=Nb;j++) printf("%d %e\n",j,tr0[j]- tr0[j-1]); exit(1); */ /* generate record */ while(time < N) { fhi = 0.2 + 0.1*(float)rand()/RAND_MAX; /* choose a new block time in beats */ Nb = blocklength(); /* choose a transition duration in beats */ Nt = translength(); /* choose rr generation states: rrmean, rrstd, lfhfratio */ dRR = dRRamp*(float)rand()/RAND_MAX; rrmean = RR0[time] + dRR; rrtrend = RRtrendamp*(2.0*(float)rand()/RAND_MAX - 1.0); rrstd = RRstdmin + (RRstdmax - RRstdmin)*(float)rand()/RAND_MAX; lfhfratio = LFHFmin + (LFHFmax - LFHFmin)*(float)rand()/RAND_MAX; /* generate RR over this block */ trblock(tr0, flo, fhi, flostd, fhistd, lfhfratio, rrmean, rrtrend, rrstd, Nb); rrb = tr0[2] - tr0[1]; /* generate RR during transition between blcoks with rr-intervals rra, rrb */ rrtrend = 0.0; trblock(trb, flo, fhi, flostd, fhistd, lfhfratio, rrmean, rrtrend, rrstd, Nt); trtrans(tr1, Nt, rra, rrb, trb); /* fill in transition and new block */ for(j=1;j<=Nt;j++) trpeaks[i+j] = trpeaks[i] + tr1[j]; i += Nt; for(j=1;j<=Nb;j++) trpeaks[i+j] = trpeaks[i] + tr0[j]; i += Nb; /* remember last rr-interval */ rra = rrb; time = (long)(trpeaks[i]); } Nbeats = i; /* check for adding ectopics or artefacts */ /* loop over entire data set */ j=1; trpeaks_new[1] = trpeaks[1]; for(i=2;i<=Nbeats;i++){ /* find out what heart rate is at a particular rpeak time */ hr=60.0/RR0[(long)(trpeaks[i])]; /* check so see if we (randomly) want to make this ectopic */ test = check_for_ectopic(); /* no state dependency */ if (test!=NO){ trpeaks_new[j] = modify_ectopic(trpeaks[i],trpeaks[i-1]); } if (test == NO){ /* just keep the current beat */ test = check_for_noise(state[i], hr); /* hr and state dependent */ if (test!=NO){ /* we have noise, add an extra beat */ trpeaks_new[j] = make_noise(trpeaks[i],trpeaks[i-1]); j++; } /* regardless, keep the current beat */ trpeaks_new[j] = trpeaks[i]; } j++; } free_vector(RR0,1,2*N); free_vector(trpeaks,1,2*N); free_vector(tr0,1,2*N); free_vector(tr1,1,2*N); free_vector(trb,1,2*N); } /* This function is called once per RR interval. It should return the length of the next (simulated) RR interval in seconds. The example code generates samples of a noisy sine wave. */ float generate(void) { float rr; count++; rr= trpeaks_new[count+1]-trpeaks_new[count]; /* trap rounding errors */ if(rr<(2.0/SPS)) rr = 2.0/SPS; return(rr); } int main(int argc, char **argv) { float t = 0.0; /* sum of intervals since the beginning of the simulation, in seconds */ long ts = 0; /* t, in sample intervals */ long tsp = 0; /* previous value of ts */ long tmax = 24*60*60; /* 24 hours, in seconds */ long seed; /* a 32-bit random variable that can be used to initialize the generator */ long atol(); if (argc < 2) { fprintf(stderr, "usage: %s seed [tmax]\n", argv[0]); exit(1); } seed = atol(argv[1]); if (argc > 2) tmax = atol(argv[2]); initialize(seed, tmax); while ((t += generate()) < tmax) { /* add RR interval to running time */ /* calculate and output a quantized RR interval */ ts = (long)(SPS*t + 0.5); printf("%5.3f\n", (ts - tsp)/((float)SPS)); tsp = ts; } exit(0); }