/* file: rrgen.c
Entry number : 171
Authors : Dragan Gamberger, Ivan Maric, Tomislav Smuc,
Gordan Bosanac, Nikola Bogunovic, Goran Krstacic
Organization : Rudjer Boskovic Institute,
Institute for Cardiovascular Prevention and Rehabilitation
Zagreb, Croatia
email address : dragan.gamberger at irb.hr
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 rr27 and rr43 of the challenge
dataset.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define SPS 128 /* sampling frequency (determines quantization of
RR intervals; see main()) */
#define STACK 20 /* stack daepth for mean value computation */
#define def_load 0.35 /* default load during the day */
float val[STACK]; /* previous interval values stack*/
float ip[STACK]; /* integral power values stack */
float val_mean; /* mean value of interval values in stack */
float ip_mean; /* mean value of integral power in stack */
float rr_high; /* expected interval value for load 0.0*/
float rr_low; /* expected interval value for load 1.0*/
float abs_max; /* absolut maximal interval value */
float abs_min; /* absolut minimal interval value */
float heart_load; /* present heart load 0.0 - 1.0 */
float int_pow; /* integral power (needed by the body) */
float int_pow_change; /* change of power */
float mem_prev; /* previous disorder */
int aux_counter; /* auxilliary counter */
/* load variables */
float load0; /* long term load */
float load00; /* intended long term load */
float load1; /* permanent load, steps */
int load1c; /* permanent load counter */
float load2; /* permanent load, continuous */
float load2p; /* positive step */
float load2m; /* negative step */
int load2cp; /* positive step counter */
int load2cm; /* negative step counter */
float load3; /* long term load, continuous */
float load3p; /* positive step */
float load3m; /* negative step */
int load3cp; /* positive step counter */
int load3cm; /* negative step counter */
float load4; /* medium term load, continuous */
float load4p; /* positive step */
float load4m; /* negative step */
int load4cp; /* positive step counter */
int load4cm; /* negative step counter */
float load5; /* medium term load, steps */
int load5c; /* counter */
float load6; /* short term load, steps */
int load6c; /* counter */
/* global time constants */
float start_ss,per_ss1,per_ss2; /* basic sleep period */
float incr1,incr2; /* increments for sleep */
float start_wk1,per_wk1; /* basic work period */
float start_wk2,per_wk2; /* intensive work period */
/* main random generating routine */
float randn(void)
{
return(((float)rand())/((float)RAND_MAX));
}
float distri(xa,xbl,xbh,xcl,xch)
float xa,xbl,xbh,xcl,xch;
{
float x,y;
double base,exp,res;
x=randn();
if(x>.5) { /* increase interval */
x=2.0*(x-0.5);
exp=x; /* x = 0,1 */
base=2.0; /* select distribution slope 2=flat, 25 steep */
if(xch/xbh > 2.65) base=4.0;
if(xch/xbh > 3.30) base=8.0;
if(xch/xbh > 4.10) base=12.0;
if(xch/xbh > 4.70) base=16.0;
if(xch/xbh > 5.50) base=25.0;
res=pow(base,exp);
y=xa+xch*((float) res -1.0)/((float) base -1.0);
}
else { /* decrease interval */
exp=2.0*x; /* x = 0,1 */
base=2.0; /* select distribution slope 2=flat, 25 steep */
if(xcl/xbl > 2.65) base=4.0; /* both xbl xcl MUST be negative */
if(xcl/xbl > 3.30) base=8.0;
if(xcl/xbl > 4.10) base=12.0;
if(xcl/xbl > 4.70) base=16.0;
if(xcl/xbl > 5.50) base=25.0;
res=pow(base,exp);
y=xa+xcl*((float) res -1.0)/((float) base -1.0); /* xcl is negative */
}
return(y);
}
/* This function is called once, before generate() is called. See main()
for a description of its arguments.
*/
void initialize(long seed, long tmax)
{
int j,jj; /* auxilliary counters */
srand((unsigned int)seed); /* initialize random generator */
for (j=0;;) {
j =(int) (randn()*3000) -2000; /* do not use few first random numbers */
if(j>0) break;
}
for(jj=0;jj<j+10;jj++) randn();
/* select and compute constants */
rr_low = 0.48 + 0.05*randn(); /* select rr at maximal load */
rr_high = rr_low * 1.8 + .05*randn(); /* select rr at minimal load */
for(;;) { /* prevent this rr to be too long */
if(rr_high<1.150) break;
rr_high-=0.050;
}
abs_min = rr_low * .7; /* select rr at minimal load */
for(;;) { /* prevent this rr to be too short */
if(abs_min>0.35 && abs_min<rr_low-.05) break;
if(abs_min<0.35) abs_min+=.015;
else abs_min-=.015;
}
abs_max = rr_high * 1.52; /* select theoretically maximal rr */
for(;;) { /* prevent this rr to be too long */
if(abs_max<1.550) break;
abs_max-=.05;
}
/* initialize previous rr intervals */
for(j=0;j<STACK;j++) val[j]=(rr_low+rr_high)/2.0;
for(j=0;j<STACK;j++) ip[j]=0.0;
val_mean=(rr_low+rr_high)/2.0;
ip_mean=0.0;
int_pow=0; /* integral power is in stable position */
load0=def_load; /* long term load default value*/
load1=0; /* permanent load, steps */
load1c=0; /* counter */
load2=0; /* permanent load, continuous */
load2cp=0; /* positive counter */
load2cm=0; /* negative counter */
load3=0; /* long term load, continuous */
load3cp=0; /* positive counter */
load3cm=0; /* negative counter */
load4=0; /* medium term load, continuous */
load4cp=0; /* positive counter */
load4cm=0; /* negative counter */
load5=0; /* medium term load, steps */
load5c=0; /* counter */
load6=0; /* short term load, steps */
load6c=0; /* counter */
mem_prev=0.0; /* no previous memory */
/* set global time constants */
start_ss=45000+30000*randn(); /* go to sleep whenever you like */
per_ss2=20000+10000*randn(); /* duration about 7 hours */
/* ajust so that you wake at least 1 hour before end */
if(start_ss + per_ss2 >82000) start_ss = 82000- per_ss2;
per_ss1=per_ss2*(0.7+0.2*randn()); /* duration of the first sleep part */
per_ss2 -= per_ss1; /* duratio of the second part */
incr1 = (def_load*100)/ per_ss1;
incr2 = (def_load*100)/ per_ss2;
start_wk1=7000+1000*randn(); /* start increased load during the day */
per_wk1=start_ss-start_wk1-10000-1000*randn();
start_wk2=start_wk1+10000+1000*randn(); /* start second increased load during the day */
per_wk2=start_wk1+per_wk1-start_wk2-10000-1000*randn();
aux_counter=0;
return;
}
/* 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; /* new interval value */
static float t;
float xa,xbl,xbh,xcl,xch; /* distribution constants */
int j; /* auxilliary counter */
float tt; /* auxilliary time */
float corr1,corr2; /* correction values for rr*/
float change; /* statistical instabilities of rr */
/* compute load */
/* long term */
tt=t;
for(;;) { /* ajust for more than 24 hours */
if(tt<=86400) break;
else tt-=86400;
}
/* medium load is default value during the day*/
if(tt<start_ss || tt>start_ss+per_ss1+per_ss2) load0=def_load + 0.05*sin(tt*7.26/100000);
else {
if(aux_counter>0) aux_counter--;
else {
aux_counter=100;
if(tt>start_ss && tt<start_ss+per_ss1) load0-=incr1*2.0*randn(); /* first sleep part */
if(tt>start_ss+per_ss1 && tt<start_ss+per_ss1+per_ss2) load0+=incr2*2.0*randn(); /* second sleep part */
if(load0<0) load0=0;
if(load0>def_load) load0=def_load + 0.05*sin(tt*7.26/100000);
}
}
if(tt>start_wk1 && tt<start_wk1+per_wk1) load0+=0.1;
if(tt>start_wk2 && tt<start_wk2+per_wk2) load0+=0.1;
if(tt>start_ss-10000 && tt<start_ss-3000) load0-=0.1;
/* permanent load changes, steps */
if(load1c>0) load1c--;
if(load1c==0) {
load1c= (int) (1000*randn()); /* duration of the new offset */
if(load0<def_load-0.05) load1=0.1*randn(); /* its magnitude */
else load1=0.1*randn();
}
/* permanent load changes, continuous */
if(load2cp>0) { /* first increase load */
load2-=load2p;
load2cp--;
}
else if(load2cm>0) { /* after that decrese it */
load2+=load2m;
load2cm--;
}
if(load2cp==0 && load2cm==0) { /* determine new period */
load2=0;
load2cp= (int) (10+100*randn()); /* up and down times */
load2cm= (int) (load2cp*randn());
load2cp -= load2cm;
if(load0<def_load-0.05) load2p= 0.1*randn(); /* peak value */
else load2p= 0.1*randn();
load2m= load2p/((float) (load2cm)); /* decrease value */
load2p= load2p/((float) (load2cp)); /* increase value */
}
/* long term load during the work, continuous */
if(load3cp>0) { /* first increase load */
load3-=load3p;
load3cp--;
}
else if(load3cm>0) { /* after that decrease it */
load3+=load3m;
load3cm--;
}
if(load3cp==0 && load3cm==0) { /* determine new period */
load3=0;
if(load0>def_load-0.05 && randn()<0.3) if(randn()<0.001) {
load3=0.1;
load3cp= (int) (4500+4000*randn()); /* total up and down times */
load3cm= (int) (load3cp*randn()); /* decrease time */
load3cp -= load3cm; /* increase time */
load3p= 0.1*randn()+0.1; /* peak value */
load3m= load3p/((float) (load3cm)); /* decrease value */
load3p= load3p/((float) (load3cp)); /* increase value */
}
}
/* medium term load during the work, continuous */
if(load4cp>0) { /* first increase load */
load4-=load4p;
load4cp--;
}
else if(load4cm>0) { /* after that decrese it */
load4+=load4m;
load4cm--;
}
if(load4cp==0 && load4cm==0) { /* determine new period */
load4=0;
if(load0>def_load-0.05 && randn()<0.1) if(randn()<0.001) {
load4=0.1;
load4cp= (int) (4500+4000*randn()); /* total up and down times */
load4cm= (int) (load4cp*randn()); /* decrease time */
load4cp -= load4cm; /* increase time */
load4p= 0.1*randn()+ 0.1; /* peak value */
if(randn()<0.1) load4p=load4p*(-2.0);
if(load0+load1+load2+load3>0.7) load4p= 0.1*randn();
load4m= load4p/((float) (load4cm)); /* decrease value */
load4p= load4p/((float) (load4cp)); /* increase value */
}
}
/* medium term load during the whole day, steps*/
if(load5c>0) {
load5c--;
}
else{
load5=0.0;
if(randn()<0.0005*(load0+load1+load2+load3+load4+0.1)) {
load5c= (int) (2500+5000*randn());
load5=0.15*randn()+0.1;
if(load0+load1+load2+load3+load4>0.7) load5=load5/4.0;
if(load0<def_load-0.05) load5=load5/2.0;
}
}
/* short term load during the whole day, steps*/
if(load6c>0) {
load6c--;
}
else{
load6=0.0;
if(randn()<0.005*(load0+load1+load2+load3+load4+load5+0.1)) {
load6c= (int) (10+50*randn());
load6=0.1*randn()+0.05;
if(load0+load1+load2+load3+load4+load5>0.9) load6=load6/4;
if(val[0]<0.5) load6=0.0;
}
}
heart_load=load0+load1+load2+load3+load4+load5+load6; /* total heart load */
if(heart_load>1.0) heart_load=1.0; /* normalization */
if(heart_load<0.0) heart_load=0.0;
/* compute changes in integral power of the body */
int_pow_change=heart_load+1.0-2*(rr_high-rr_low)/(val[0]+rr_high-2*rr_low);
/* integral power correction, long term stability */
/* ** increse second constant if global oscilations are too high */
/* maximal oscilation depression with C=0.015 */
if(int_pow>0) corr1=-.002;
else corr1=.002;
if(int_pow>0 && int_pow>ip_mean && int_pow_change>0) corr1-=0.015*randn();
if(int_pow<0 && int_pow<ip_mean && int_pow_change<0) corr1+=0.015*randn();
if(int_pow>ip_mean && int_pow_change>0) corr1-=0.015*randn();
if(int_pow<ip_mean && int_pow_change<0) corr1+=0.015*randn();
if(int_pow_change>0) corr1-=0.003;
if(int_pow_change<0) corr1+=0.003;
/* recent history correction, use it only in stable situation */
corr2=0;
if(int_pow>-0.25 && int_pow<0.25 ) { /* stable condition */
corr2+=(val_mean-val[0])*(0.0+1.0*randn())*1.2*val[0];
corr2+=(val_mean-val[1])*(0.0+1.0*randn())*1.4*val[0];
}
/* introduce instabilities */
/* basic algorithm, use distribution to determine rr changes */
xa=0.000;
xbl=-0.010*val[0];
xbh=0.010*val[0];
xcl=-0.020*(1.0+val[0]);
xch=0.030*(1.0+val[0]);
change=distri(xa,xbl,xbh,xcl,xch); /* use distribution to compute difference */
/* asimetric instabilities */
if(randn()<0.005 && val[0]>0.7) change+=0.4*randn()*val[0];
/*predicted interval value is previous + corrections + instabilities */
rr=val[0]+corr1+corr2+change;
/* suprises */
if(mem_prev>0.001 || mem_prev<-0.001) {
if(randn()<0.5) {
rr-=mem_prev*(0.5+0.5*randn());
mem_prev=0.0;
}
}
else if(randn() < 0.002) {
if(val[0]> 0.8 && randn()<0.1) mem_prev=(-1.0)*(0.1+0.1*randn());
else{
if(val[0]<1.0 && val[0]>0.6 && randn()<0.5) rr+=.2;
else mem_prev=0.05+.1*randn();
}
rr+=mem_prev;
}
/* prevent this rr to be too short or too long*/
for(;;) { /* prevent this rr to be too short */
if(rr>abs_min) break;
if(int_pow>1.0) rr+=0.05*randn();
else rr+=0.25*randn();
}
for(;;) { /* prevent this rr to be too long */
if(rr<abs_max) break;
rr-=0.5*randn();
}
int_pow+=(int_pow_change); /* update integral power */
val_mean+=(rr-val[9])/10.0; /* compute means */
ip_mean+=(int_pow-ip[9])/10.0;
/* save in stack */
for(j=STACK-1;j>0;j--) val[j]=val[j-1];
val[0]=rr;
for(j=STACK-1;j>0;j--) ip[j]=ip[j-1];
ip[0]=int_pow;
t += rr;
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);
}