/* file: rrgen.c
     Entry number  : 161
     Authors       : Miguel A. García-González, Juan Ramos-Castro
     Organization  : Instrumentation and Bioengineering Division, 
				Electronic Engineering Department,
				Universitat Politècnica de Catalunya
				Barcelona, Spain
     email address : magarcia at eel.upc.es

   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 rr32 and rr34 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 PI 3.14159263
#define DIM 1024

void simcirc(int *, long, float *);
void sese(int, float , float *);
void moving(float *,int,int, float *);
void escalado(float *,float, float );
void fillmat(float [][5],float *,int , int);
float mean(float *,int);
float stdv(float *,int);
void prepro(float *,int);
float generate(void);
float bbv(float,float,float,float,float,float);
float sqsig(float, float);
void square(float *, float , int , float *);
float randn(void);
float randu(void); 

/**********************************/
/* GLOBAL VARIABLES             ***/
float lf1;
float f1;
float sd1;
float xrr[1024];
long tgmax;
/**********************************/


/* This function is called once, before generate() is called.  See main()
	for a description of its arguments.
*/

/***********************************************************************/






void initialize(long seed, long tmax)
{
	float mrr,ddn;
	int mat5[5]={33,10,5,2,1};
/* Initialize global variables */
	srand((unsigned int) seed);
   mrr=0.75+0.3*randu(); /* Mean rr */
   sd1=0.07*mrr+0.04*(randu()); /* generated signal standard deviation if mean rr equals 1 second */
   lf1=0.5+randu()/2.0; /* ratio related to lf/hf index */
   f1=2.0+18.0*randu(); /* ratio related to determistic no deterministic data */
	ddn=0.5*mrr+0.05*randu();
	tgmax=tmax;
	simcirc(mat5 , tmax, xrr);
	escalado(xrr,mrr,ddn);
	return;
}


/* This function is called once per RR interval.  It should return
 the length of the next (simulated) RR interval in seconds.
*/

float generate(void)
{
	static float rr,temp,dfr,fmod,bfp,bfd,bfm,dut,ph,ts,sdf;
	static float t,fr=0.15,kf,ss[512],ss2[512],bar[512];
	int index,k;
	static int lresp,kresp;
	index = (int) floor(t/( (long) tgmax)* ((int) DIM));
	if(index>(DIM-1)) index=DIM-1;
	if(t>=kf){
		kresp=0;
		if (fr > 0.35) fr=fr-0.01*randu();
		if (fr < 0.18) {fr=fr+0.03*randu();}
		else {fr=fr+0.02*(randu()-0.5);};
		lresp= (int) floor(250.0+100.0*randu());
		dfr=0.01*(1.0+fabs(xrr[index]-1))*randu();
		kf=t+xrr[index]*lresp;
      /*Fractal signal used for simulate complex behaviour */
		sese(512, 1.3+0.1*randu(),ss);
		moving(ss,512,100,ss2);
      fmod=0.000*randu();
      bfp=0.12*(1+0.3*(randu()-0.5));
      bfd=0.003*randu();
      bfm=0.07*randu();
      dut=80.0+10.0*randu();
      ph=PI/2.0*(1+0.1*randu());
      for (k=0;k<512;k++) {
			 ts=(t+xrr[index]*k);
	  bar[k]=sqsig(2.0*PI*bfp*(1.0+bfd/2.0*(sin(2.0*PI*bfm*ts)+1.0))*ts+ph,dut);
      }
		moving(bar,512,1,bar);
		prepro(bar,512);
		temp=mean(bar,512);
		for (k=0;k<512;k++) bar[k]=bar[k]-temp;
		sdf=0.05*randu()+0.1;
	}
    /* Simulate non chaotic data */
	temp=bbv(sd1,t,xrr[index],fr,dfr,fmod)*(3*ss2[kresp]+1);
	temp=temp+(pow(xrr[index],-1.5)+0.03*randu())*bar[kresp]*sd1/lf1;
	rr=temp*pow(xrr[index],1.5)+ss[kresp]*f1*sd1*sd1*pow(xrr[index],-1)/3+0.01*randn();
	rr=rr*pow(xrr[index],2)+xrr[index];
	if ((rr > 2) || (rr < 0.5)) rr=xrr[index];
	t += rr;
	kresp++;
   return (rr);
}


/********************* Aqui pongo mi codigo ***************************/
float randn(void)  /* Generador de numeros aleatorios con dist. gaussiana */
{
	float x1,x2,z;
	x1=((float) rand())/RAND_MAX;
	x2=((float) rand())/RAND_MAX;
	x1=x1+1e-10;
	z=sqrt(-2*log(x1))*cos(2*PI*x2);
	return(z);
}

float randu(void)   /* Generador de numeros aleatorios con dist. uniforme (0-1) */
{
	return((float) rand()/RAND_MAX);

}

float mean( float *vector, int dim) /* Calcula la media del vector */
{
	int k;
	float med=0;
	for(k=0;k<dim;k++)
	{
		med=med+vector[k];
	}
	med=med/dim;
	return(med);
}

float stdv(float *vector, int dim)  /* Calcula la desv. estandar ajustada vector*/
{
	int k;
	float var=0.0;

	for(k=0;k<dim;k++)
	{
		var=var+vector[k]*vector[k];
	}
	var=sqrt((var/dim-pow(mean(vector,dim),2.0))*dim/(dim-1));
	return(var);
}

void square(float *ti, float duty, int dim, float *x)
{
	int k;
	float w,temp;
	w=2.0*PI*duty/100.0;
	for(k=0;k<dim;k++){
		temp=-2.0*PI*floor(ti[k]/(2.0*PI));
		x[k]=ti[k]+temp;
		if(x[k]<w) {x[k]=1.0;}
		else {x[k]=0.0;}
	}

}

float sqsig(float ti, float duty)
{
	float w,temp;

	w=2.0*PI*duty/100.0;
	temp=-2.0*PI*floor(ti/(2.0*PI));
	temp=ti+temp;
	if(temp<w) {return(1.0);}
	else {return(0.0);}
}

void fillmat(float mat[][5],float *x,int col, int dim)
{
	int k;
	for(k=0;k<dim;k++) mat[k][col]=x[k];
}

void moving(float *v,int dim,int wl, float *vf)
{
	int k,j;
	for(k=wl;k<(dim-wl);k++) {
		vf[k]=0;
		for(j=-wl;j<=wl;j++) vf[k]=vf[k]+v[k+j];
		vf[k]=vf[k]/(wl*2.0+1.0);
	}
	for(k=0;k<wl;k++) {
		vf[k]=vf[k+wl];
		vf[dim-k-1]=vf[dim-k-wl];
	}
}

void prepro(float *v, int dim)
{
	int k;
	float ma,mi;
	mi=1.0e10;
	ma=-1.0e10;
	for(k=0;k<dim;k++) {
		if(v[k]<mi) mi=v[k];
		if(v[k]>ma) ma=v[k];
	}
	for(k=0;k<dim;k++) v[k]=(v[k]-mi)/(ma-mi)-0.5;
}

void escalado(float *v,float mrr, float dnd)
{
	float m;
	int k;
	float vf[1024];
	prepro(v,1024);
	for(k=0;k<1024;k++) v[k]=(v[k]+1.5)*dnd;
	m=mean(v,1024);
	for(k=0;k<1024;k++){
	 vf[k]=v[k]-m+mrr;
	 }
	moving(vf,1024,2,v);
	for (k=1;k<1023;k++) {
		if (fabs(v[k+1]-v[k])>0.05) v[k]=(v[k-1]+v[k+1])/2.0;
		}
	v[1023]=v[1021];
	v[1022]=v[1020];
	
	return;
}

void simcirc(int dcd[5], long tmax, float *x)
/* Genera el ritmo circadiano */
{
	float ct[1024][5];
	float ti[1024],xf[1024],ph,pes[5];
	int k,zc[5],zcp[5],pos[5],nzc=0;
	int i,n1,n2,j;

	ph=1.5*PI+PI*randu()*0.1;
	for(k=0;k<1024;k++) ti[k]=(PI*2.0*k*tmax/(3600.0*24.0))/1023+ph;
	square(ti,dcd[0]+randu(),1024,x);
	fillmat(ct,x,0,1024);
	nzc=0;
	for(k=0;k<1023;k++){
		if(x[k]!=x[k+1]) {
			zcp[nzc]=1023-k;

			nzc++;
		}
	}
	for(k=0;k<nzc;k++){
		zc[k]=zcp[nzc-k-1];

		}
	for(i=1;i<5;i++) {
		ph=(randu()-0.5)*PI/4;
		for(k=0;k<1024;k++) ti[k]=(PI*2.0*k*tmax/(3600.0*24.0))/1023.0+ph+(k+1)*PI/4;
		square(ti,dcd[i]+randn(),1024,x);
		fillmat(ct,x,i,1024);
		pes[i]=0.5*(randu()-0.5);
	}
	pes[0]=1.0;
	for(k=0;k<1024;k++) {
		x[k]=0;
		for(i=0;i<5;i++) x[k]=x[k]+ct[1023-k][i]*pes[i];
	}

	n1= (int) ceil(((43.0-7.11)*randu())/2*tmax/3600.0/24.0);
	n2= (int) floor((13.11*randu())/2*tmax/3600.0/24.0);
	zc[nzc]=1023;
	i=zc[0];
	for(k=0;k<nzc;k++){
		pos[k]= (int) ceil((i+zc[k+1])/2);
		i=zc[k+1];
	}
	i=0;
	for(j=0;j<nzc;j++){

		for(k=0;k<(pos[j]-i);k++) ti[k]=x[k+i];
		if(j%2 == 0) moving(ti,pos[j]-i,n1,xf);
		else moving(ti,pos[j]-i,n2,xf);
		for(k=0;k<(pos[j]-i);k++) x[k+i]=xf[k];
		i=pos[j];
		}
	sese(1024,1.5,ti);
	moving(ti,1024,40,xf);
	n1= (int) (3*randu()+2);
	moving(x,1024,n1,ti);
	for(k=0;k<1024;k++) x[k]=ti[k]+xf[k]*1.5;
	return;


}

void sese(int n, float d, float *x)
/* Genera serie fractal de longitud n y dimension d, la salida esta en x */
{
	static float h,b,ti[1024],q[1024],p[1024],a[1024],c[1024],temp,temp2;
	int k,n2,j;

	n2=n/2;
	h=2-d;
	b=1+2*h;
	for(k=0;k<n2;k++){
		ti[k]=k+1;
		if(k==0) p[k]=0;
		else p[k]=pow(k,-b/2)*randn();
		q[k]=2*PI*randu();
		a[k]=p[k]*cos(q[k]);
		c[k]=p[k]*sin(q[k]);
	}
	for(k=0;k<n;k++){
		temp=0;
		for(j=0;j<n2;j++){
			temp=temp+a[j]*cos(2*PI*ti[j]*k/n)+c[j]*sin(2*PI*ti[j]*k/n);
		}
		x[k]=temp;
	}
	for(k=0;k<(n-1);k++){
		x[k]=x[k+1]-x[k];
	}
	x[n]=x[n-1];
	temp=mean(x,n);
	temp2=stdv(x,n);
	for(k=0;k<n;k++) x[k]=(x[k]-temp)/temp2;
	return;

}

/* Beat to beat deterministic variatio due to respiration and baroreceptor */
float bbv(float gsd1,float t,float rm,float fr, float dfr, float fmod)
{
/* gsd1: generated signal standard deviation if mean rr equals 1 second */
/* rm:  mean rr */
/* blf1: ratio related to lf/hf index */
/* bf1: ratio related to determistic no deterministic data */
/* t: simulation time */
/* fr: mean respiratory frequency */


   float resp; /*  */


   resp=(1.0+sin(2.0*PI*fr/rm*(1.0+dfr*sin(2.0*PI*fmod*t))*t))/2.0-0.5;
	resp=resp*(pow(rm,1.5)+0.001*randn())*gsd1;

	return(resp);

}



/***********************************************************************/


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