ecgsyn.c

ecg心电信号产生程序

        freeVect(dydx,1,n);
        freeVect(dym,1,n);
        freeVect(dyt,1,n);
        freeVect(yt,1,n);
}

/*--------------------------------------------------------------------------*/
/*    DETECT PEAKS                                                          */
/*--------------------------------------------------------------------------*/

double detectpeaks(double *ipeak, double *x, double *y, double *z, int n)
{
   int i,j,j1,j2,jmin,jmax,d;
   double thetap1,thetap2,thetap3,thetap4,thetap5;
   double theta1,theta2,d1,d2,zmin,zmax;
   
   /* use globally defined angles for PQRST */
   thetap1 = ti[1];
   thetap2 = ti[2];
   thetap3 = ti[3];
   thetap4 = ti[4];
   thetap5 = ti[5];

   for(i=1;i<=n;i++) ipeak[i] = 0.0;
   theta1 = atan2(y[1],x[1]);
   for(i=1;i<n;i++)
   {
      theta2 = atan2(y[i+1],x[i+1]);
      if( (theta1 <= thetap1) && (thetap1 <= theta2) )  
      {
	d1 = thetap1 - theta1;
        d2 = theta2 - thetap1;
        if(d1 < d2)  ipeak[i] = 1.0;
        else         ipeak[i+1] = 1.0;
      }
      else if( (theta1 <= thetap2) && (thetap2 <= theta2) )  
      {
	d1 = thetap2 - theta1;
        d2 = theta2 - thetap2;
        if(d1 < d2)  ipeak[i] = 2.0;
        else         ipeak[i+1] = 2.0;
      }
      else if( (theta1 <= thetap3) && (thetap3 <= theta2) )  
      {
	d1 = thetap3 - theta1;
        d2 = theta2 - thetap3;
        if(d1 < d2)  ipeak[i] = 3.0;
        else         ipeak[i+1] = 3.0;
      }
      else if( (theta1 <= thetap4) && (thetap4 <= theta2) )  
      {
	d1 = thetap4 - theta1;
        d2 = theta2 - thetap4;
        if(d1 < d2)  ipeak[i] = 4.0;
        else         ipeak[i+1] = 4.0;
      }
      else if( (theta1 <= thetap5) && (thetap5 <= theta2) )  
      {
	d1 = thetap5 - theta1;
        d2 = theta2 - thetap5;
        if(d1 < d2)  ipeak[i] = 5.0;
        else         ipeak[i+1] = 5.0;
      }
      theta1 = theta2; 
   }

   /* correct the peaks */
   d = (int)ceil(sfecg/64);
   for(i=1;i<=n;i++)
   { 
     if( ipeak[i]==1 || ipeak[i]==3 || ipeak[i]==5 )
     {
        j1 = MAX(1,i-d);
        j2 = MIN(n,i+d);
        jmax = j1;
        zmax = z[j1];
        for(j=j1+1;j<=j2;j++)
	{ 
	   if(z[j] > zmax) 
           {
	      jmax = j;
              zmax = z[j];
	   }
	}
        if(jmax != i)
	{
           ipeak[jmax] = ipeak[i];
           ipeak[i] = 0;
	}
     }
     else if( ipeak[i]==2 || ipeak[i]==4 )
     {
        j1 = MAX(1,i-d);
        j2 = MIN(n,i+d);
        jmin = j1;
        zmin = z[j1];
        for(j=j1+1;j<=j2;j++)
	{ 
	   if(z[j] < zmin) 
           {
	      jmin = j;
              zmin = z[j];
	   }
	}
        if(jmin != i)
	{
           ipeak[jmin] = ipeak[i];
           ipeak[i] = 0;
	}
     }
   }

}

/*--------------------------------------------------------------------------*/
/*      MAIN PROGRAM                                                         */
/*---------------------------------------------------------------------------*/

main(argc,argv)
int     argc;
char    **argv;
{

    /* First step is to register the options */

    optregister(outfile,CSTRING,'O',"Name of output data file");  
    optregister(N,INT,'n',"Approximate number of heart beats");    
    optregister(sfecg,INT,'s',"ECG sampling frequency [Hz]");   
    optregister(sf,INT,'S',"Internal Sampling frequency [Hz]"); 
    optregister(Anoise,DOUBLE,'a',"Amplitude of additive uniform noise [mV]");
    optregister(hrmean,DOUBLE,'h',"Heart rate mean [bpm]");
    optregister(hrstd,DOUBLE,'H',"Heart rate standard deviation [bpm]");
    optregister(flo,DOUBLE,'f',"Low frequency [Hz]");
    optregister(fhi,DOUBLE,'F',"High frequency [Hz]");
    optregister(flostd,DOUBLE,'v',"Low frequency standard deviation [Hz]");
    optregister(fhistd,DOUBLE,'V',"High frequency standard deviation [Hz]");
    optregister(lfhfratio,DOUBLE,'q',"LF/HF ratio");
    optregister(seed,INT,'R',"Seed");    
    opt_title_set("ECGSYN: A program for generating a realistic synthetic ECG\n" 
     "Copyright (c) 2003 by Patrick McSharry & Gari Clifford. All rights reserved.\n");

    getopts(argc,argv);

    dorun();
}



/*--------------------------------------------------------------------------*/
/*    DORUN PART OF PROGRAM                                                 */
/*--------------------------------------------------------------------------*/

int dorun()
{
   int i,j,k,q,Nrr,Nt,Nts;
   double *x,tstep,tecg,rrmean,qd,hrfact,hrfact2;
   double *xt,*yt,*zt,*xts,*yts,*zts;
   double timev,*ipeak,zmin,zmax,zrange;
   FILE *fp;
   void (*derivs)(double, double [], double []);

   /* perform some checks on input values */
   q = (int)rint(sf/sfecg);
   qd = (double)sf/(double)sfecg;
   if(q != qd) {
     printf("Internal sampling frequency must be an integer multiple of the \n"); 
     printf("ECG sampling frequency!\n"); 
     printf("Your current choices are:\n");
     printf("ECG sampling frequency: %d Hertz\n",sfecg);
     printf("Internal sampling frequency: %d Hertz\n",sf);
     exit(1);}


   /* declare and initialise the state vector */
   x=mallocVect(1,mstate);
   x[1] = xinitial; 
   x[2] = yinitial;
   x[3] = zinitial;

   /* declare and define the ECG morphology vectors (PQRST extrema parameters) */
   ti=mallocVect(1,5);
   ai=mallocVect(1,5);
   bi=mallocVect(1,5);
   /* P            Q            R           S           T        */
   ti[1]=-60.0; ti[2]=-15.0; ti[3]=0.0;  ti[4]=15.0; ti[5]=90.0;
   ai[1]=1.2;   ai[2]=-5.0;  ai[3]=30.0; ai[4]=-7.5; ai[5]=0.75;
   bi[1]=0.25;  bi[2]=0.1;   bi[3]=0.1;  bi[4]=0.1;  bi[5]=0.4;

   /* convert angles from degrees to radians */
   for(i=1;i<=5;i++) ti[i] *= PI/180.0;

   /* adjust extrema parameters for mean heart rate */
   hrfact = sqrt(hrmean/60.0);
   hrfact2 = sqrt(hrfact);
   for(i=1;i<=5;i++) bi[i] *= hrfact;
   ti[1]*=hrfact2;  ti[2]*=hrfact; ti[3]*=1.0; ti[4]*=hrfact; ti[5]*=1.0;


   /* calculate time scales */
   h = 1.0/sf;
   tstep = 1.0/sfecg;

   printf("ECGSYN: A program for generating a realistic synthetic ECG\n" 
    "Copyright (c) 2003 by Patrick McSharry & Gari Clifford. All rights reserved.\n"
    "See IEEE Transactions On Biomedical Engineering, 50(3), 289-294, March 2003.\n"
    "Contact P. McSharry ([e-mail]) or G. Clifford ([e-mail])\n"); 

   printf("Approximate number of heart beats: %d\n",N);
   printf("ECG sampling frequency: %d Hertz\n",sfecg);
   printf("Internal sampling frequency: %d Hertz\n",sf);
   printf("Amplitude of additive uniformly distributed noise: %g mV\n",Anoise);
   printf("Heart rate mean: %g beats per minute\n",hrmean);
   printf("Heart rate std: %g beats per minute\n",hrstd);
   printf("Low frequency: %g Hertz\n",flo);
   printf("High frequency std: %g Hertz\n",fhistd);
   printf("Low frequency std: %g Hertz\n",flostd);
   printf("High frequency: %g Hertz\n",fhi);
   printf("LF/HF ratio: %g\n",lfhfratio);

   /* initialise seed */
   rseed = -seed;  


   /* select the derivs to use */
   derivs = derivspqrst;

   /* calculate length of RR time series */
   rrmean = (60/hrmean);
   Nrr = (int)pow(2.0, ceil(log10(N*rrmean*sf)/log10(2.0)));	 
   printf("Using %d = 2^%d samples for calculating RR intervals\n",
           Nrr,(int)(log10(1.0*Nrr)/log10(2.0))); 


   /* create rrprocess with required spectrum */
   rr = mallocVect(1,Nrr);
   rrprocess(rr, flo, fhi, flostd, fhistd, lfhfratio, hrmean, hrstd, sf, Nrr); 
   vecfile("rr.dat",rr,Nrr);

   /* create piecewise constant rr */
   rrpc = mallocVect(1,2*Nrr);
   tecg = 0.0;
   i = 1;
   j = 1;
   while(i <= Nrr)
   {  
      tecg += rr[j];
      j = (int)rint(tecg/h);
      for(k=i;k<=j;k++) rrpc[k] = rr[i];
      i = j+1;
   }
   Nt = j;
   vecfile("rrpc.dat",rrpc,Nt);

   printf("Printing ECG signal to file: %s\n",outfile);

   /* integrate dynamical system using fourth order Runge-Kutta*/
   xt = mallocVect(1,Nt);
   yt = mallocVect(1,Nt);
   zt = mallocVect(1,Nt);

   timev = 0.0;
   for(i=1;i<=Nt;i++)
   {
      xt[i] = x[1];
      yt[i] = x[2];
      zt[i] = x[3];
      drk4(x, mstate, timev, h, x, derivs);
      timev += h;
   }


   /* downsample to ECG sampling frequency */
   xts = mallocVect(1,Nt);
   yts = mallocVect(1,Nt);
   zts = mallocVect(1,Nt);

   j=0;
   for(i=1;i<=Nt;i+=q)
   { 
      j++;
      xts[j] = xt[i];
      yts[j] = yt[i];
      zts[j] = zt[i];
   }
   Nts = j;


   /* do peak detection using angle */
   ipeak = mallocVect(1,Nts);
   detectpeaks(ipeak, xts, yts, zts, Nts);
 
   /* scale signal to lie between -0.4 and 1.2 mV */
   zmin = zts[1];
   zmax = zts[1];
   for(i=2;i<=Nts;i++)
   {
     if(zts[i] < zmin)       zmin = zts[i];
     else if(zts[i] > zmax)  zmax = zts[i];
   }
   zrange = zmax-zmin;
   for(i=1;i<=Nts;i++) zts[i] = (zts[i]-zmin)*(1.6)/zrange - 0.4;

   /* include additive uniformly distributed measurement noise */
   for(i=1;i<=Nts;i++) zts[i] += Anoise*(2.0*ran1(&rseed) - 1.0);    

   /* output ECG file */
   fp = fopen(outfile,"w");
   for(i=1;i<=Nts;i++) fprintf(fp,"%f %f %d\n",(i-1)*tstep,zts[i],(int)ipeak[i]);
   fclose(fp);


   printf("Finished ECG output\n");

freeVect(x,1,mstate);
freeVect(rr,1,Nrr);
freeVect(rrpc,1,2*Nrr);
freeVect(ti,1,5);
freeVect(ai,1,5);
freeVect(bi,1,5);
freeVect(xt,1,Nt);
freeVect(yt,1,Nt);
freeVect(zt,1,Nt);
freeVect(xts,1,Nt);
freeVect(yts,1,Nt);
freeVect(zts,1,Nt);
freeVect(ipeak,1,Nts);

/* END OF DORUN */
}