/* file: rrgen.c
Entry number : 191
Author(s) : D.C. Lin
Organization : Ryerson University
email address : derlin at acs.ryerson.ca
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 rr10 and rr37 of the challenge
dataset.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define SPS 128
/* >>>>>>>>>>>>>>>>>>>> MY define's <<<<<<<<<<<<<<<<<<<< */
/* N is the maximum number of RR intervals; it must be enough to generate a
sequence of up to 48 hours */
#define N 131072
#define J 17
#define max(a,b) ((a>b) ? a:b)
#define min(a,b) ((a>b) ? b:a)
#define dabs(a) ((a>0) ? a:-a)
/* >>>>>>>>>>>>>>>>>>> MY functions <<<<<<<<<<<<<<<<<<<< */
float randn()
{
float R=2,v1,v2;
while (R>1) {
v1=2*drand48()-1;
v2=2*drand48()-1;
R=v1*v1+v2*v2;
}
return(v2*sqrt(-2*log(R)/R));
}
void bsort(int *a, int *b, int n)
{
int i, j;
int temp;
for (i=0;i<n;i++) b[i]=a[i];
for (i=0;i<n-1;i++)
for (j=i+1;j<n;j++)
if (b[i]>b[j]) {
temp=b[j];
b[j]=b[i];
b[i]=temp;
}
}
/* >>>>>>>>>>>>>>>>>>>>>>>>>> <<<<<<<<<<<<<<<<<<<<<<<<<< */
float rr[N]; /* generated RR intervals, in seconds */
void initialize(long seed, long tmax)
{
float alpha=-0.1264, sig0=0.3;
float meanRR, stdRR, per_pvc=0.0025;
int i,j,k,nsg,nsg_1,npt_sg,LLid,LRid,RLid,RRid,il,ir;
float meanr0=0, stdr0=0, mu0, wleft, wright;
float midr;
static float sigmaj[J], r0[N];
float *r = rr;
int ipvc=0, npvc, itmp;
static int id[N];
float jump;
float maxRR, minRR, RRspan;
srand48(seed);
meanRR=1.;
stdRR=.1;
for (i=0;i<J;i++) sigmaj[i]=sig0*pow(2,alpha*i);
mu0=pow((double) meanRR, (double) 1.0/J);
for (i=0;i<N;i++) r0[i]=mu0;
for (j=2;j<J+2;j++) {
nsg=pow(2,j-1);
npt_sg=N/nsg;
nsg_1=nsg/2+1;
for (k=1;k<nsg_1;k++) {
LLid=2*(k-1)*npt_sg; LRid=LLid+npt_sg;
RLid=(2*k-1)*npt_sg; RRid=RLid+npt_sg;
wleft=randn()*sigmaj[j-2]+mu0;
wright=randn()*sigmaj[j-2]+mu0;
for (il=LLid;il<LRid;il++) r0[il] *= wleft;
for (ir=RLid;ir<RRid;ir++) r0[ir] *= wright;
}
}
/* ADJUST r0 TO meanRR AND stdRR */
for (i=0;i<N;i++) meanr0 += r0[i]; meanr0 /= N;
for (i=0;i<N;i++) r[i]=r0[i]*meanRR/meanr0;
for (i=0;i<N;i++) r[i]=r[i] - meanRR;
for (i=0;i<N;i++) stdr0 += r[i]*r[i]; stdr0=sqrt(stdr0/N);
for (i=0;i<N;i++) r[i] = r[i]*stdRR/stdr0 + meanRR;
npvc=(int)(N*per_pvc);
while (ipvc<npvc) {
itmp=(int)(drand48()*N);
if (ipvc==0)
id[ipvc++]=itmp;
else
for (i=0;i<ipvc;i++)
if (itmp != id[i]) {
id[ipvc]=itmp;
ipvc += 1;
break;
}
}
bsort(id,id,npvc);
for (i=0;i<npvc;i++) {
j=id[i];
if (j > 15) {
midr=(r[j-1]+r[j])*.5;
maxRR=-100;
minRR=1000;
for (k=j-15;k<j+15;k++) {
maxRR = max(maxRR,r[k]);
minRR = min(minRR,r[k]);
}
RRspan=maxRR-minRR;
if (max(dabs(midr-r[j-1]),dabs(midr-r[j]))/midr < 0.2) {
jump=dabs(drand48()*RRspan);
r[j-1]=midr+jump+drand48()*midr*0.1;
r[j] =midr-jump+drand48()*midr*0.1;
}
}
}
/* for (i=0;i<N;i++) r[i] *= 128; */
}
float generate()
{
static int i;
if (i >= N) {
fprintf(stderr, "rr array is too short -- increase N and try again\n");
exit(1);
}
return (rr[i++]);
}
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);
}