ecgsyn.c

ecg心电信号产生程序

/* 09 Oct 2003 "ecgsyn.c"                                                     */
/*                                                                            */
/* Copyright (c)2003 by Patrick McSharry & Gari Clifford, All Rights Reserved */
/* See IEEE Transactions On Biomedical Engineering, 50(3),289-294, March 2003.*/
/* Contact P. McSharry (patrick AT mcsharry DOT net) or                       */
/* G. Clifford (gari AT mit DOT edu)                                          */
/*                                                                            */      
/*   This program is free software; you can redistribute it and/or modify     */ 
/*   it under the terms of the GNU General Public License as published by     */ 
/*   the Free Software Foundation; either version 2 of the License, or        */ 
/*   (at your option) any later version.                                      */ 
/*                                                                            */ 
/*   This program is distributed in the hope that it will be useful,          */ 
/*   but WITHOUT ANY WARRANTY; without even the implied warranty of           */ 
/*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the            */ 
/*   GNU General Public License for more details.                             */
/*                                                                            */      
/*   You should have received a copy of the GNU General Public License        */
/*   along with this program; if not, write to the Free Software Foundation   */
/*   Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA            */ 
/*                                                                            */ 
/*  ecgsyn.m and its dependents are freely availble from Physionet -          */ 
/*  http://www.physionet.org/ - please report any bugs to the authors above.  */

#include <stdio.h>   
#include <math.h>  
#include <stdlib.h> 
#include "opt.h"
#define PI (2.0*asin(1.0))
#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 OFFSET 1
#define ARG1 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)

/*--------------------------------------------------------------------------*/
/*    DEFINE PARAMETERS AS GLOBAL VARIABLES                                 */
/*--------------------------------------------------------------------------*/

char outfile[100]="ecgsyn.dat";/*  Output data file                   */ 
int N = 256;                   /*  Number of heart beats              */
int sfecg = 256;               /*  ECG sampling frequency             */
int sf = 256;                  /*  Internal sampling frequency        */
double Anoise = 0.0;           /*  Amplitude of additive uniform noise*/
double hrmean = 60.0;          /*  Heart rate mean                    */
double hrstd = 1.0;            /*  Heart rate std                     */
double flo = 0.1;              /*  Low frequency                      */
double fhi = 0.25;             /*  High frequency                     */
double flostd = 0.01;          /*  Low frequency std                  */
double fhistd = 0.01;          /*  High frequency std                 */
double lfhfratio = 0.5;        /*  LF/HF ratio                        */

int Necg = 0;                  /*  Number of ECG outputs              */
int mstate = 3;                /*  System state space dimension       */
double xinitial = 1.0;         /*  Initial x co-ordinate value        */
double yinitial = 0.0;         /*  Initial y co-ordinate value        */
double zinitial = 0.04;        /*  Initial z co-ordinate value        */
int seed = 1;                  /*  Seed                               */
long rseed; 
double h; 
double *rr,*rrpc;
double *ti,*ai,*bi;

/* prototypes for externally defined functions */
void dfour1(double data[], unsigned long nn, int isign);
float ran1(long *idum);

/*---------------------------------------------------------------------------*/
/*      ALLOCATE MEMORY FOR VECTOR                                           */
/*---------------------------------------------------------------------------*/

double *mallocVect(long n0, long nx)
{
        double *vect;
 
        vect=(double *)malloc((size_t) ((nx-n0+1+OFFSET)*sizeof(double)));
        if (!vect){
	  printf("Memory allocation failure in mallocVect");
	}
        return vect-n0+OFFSET;
}

/*---------------------------------------------------------------------------*/
/*      FREE MEMORY FOR MALLOCVECT                                           */
/*---------------------------------------------------------------------------*/
 
void freeVect(double *vect, long n0, long nx)
{
        free((ARG1) (vect+n0-OFFSET));
}

/*---------------------------------------------------------------------------*/
/*      MEAN CALCULATOR                                                      */
/*---------------------------------------------------------------------------*/
 
double mean(double *x, int n)
/* n-by-1 vector, calculate mean */
{
        int j;
        double add;
 
        add = 0.0;
        for(j=1;j<=n;j++)  add += x[j];
 
        return (add/n);
}


/*---------------------------------------------------------------------------*/
/*      STANDARD DEVIATION CALCULATOR                                        */
/*---------------------------------------------------------------------------*/

double stdev(double *x, int n)
/* n-by-1 vector, calculate standard deviation */
{
        int j;
        double add,mean,diff,total;

        add = 0.0;
        for(j=1;j<=n;j++)  add += x[j];
        mean = add/n;

        total = 0.0;
        for(j=1;j<=n;j++)  
        {
           diff = x[j] - mean;
           total += diff*diff;
        } 
  
        return (sqrt(total/(n-1)));
}

/*--------------------------------------------------------------------------*/
/*    WRITE VECTOR IN A FILE                                                */
/*--------------------------------------------------------------------------*/

void vecfile(char filename[], double *x, int n)
{
   int i;
   FILE *fp;
  
   fp = fopen(filename,"w");
   for(i=1;i<=n;i++)  fprintf(fp,"%e\n",x[i]);
   fclose(fp);
}

/*--------------------------------------------------------------------------*/
/*    INTERP                                                                */
/*--------------------------------------------------------------------------*/

void interp(double *y, double *x, int n, int r)
{
   int i,j;
   double a;

   for(i=1;i<=n-1;i++)
   {
      for(j=1;j<=r;j++) 
      {
         a = (j-1)*1.0/r;
         y[(i-1)*r+j] = (1.0-a)*x[i] + a*x[i+1];
      }
   }
}


/*--------------------------------------------------------------------------*/
/*    GENERATE RR PROCESS                                                   */
/*--------------------------------------------------------------------------*/

void rrprocess(double *rr, double flo, double fhi, 
double flostd, double fhistd, double lfhfratio,  
double hrmean, double hrstd, double sf, int n)
{
   int i,j;
   double c1,c2,w1,w2,sig1,sig2,rrmean,rrstd,xstd,ratio;
   double df,dw1,dw2,*w,*Hw,*Sw,*ph0,*ph,*SwC;

   w = mallocVect(1,n);
   Hw = mallocVect(1,n);
   Sw = mallocVect(1,n);
   ph0 = mallocVect(1,n/2-1);
   ph = mallocVect(1,n);
   SwC = mallocVect(1,2*n);


   w1 = 2.0*PI*flo;
   w2 = 2.0*PI*fhi;
   c1 = 2.0*PI*flostd;
   c2 = 2.0*PI*fhistd;
   sig2 = 1.0;
   sig1 = lfhfratio;
   rrmean = 60.0/hrmean;
   rrstd = 60.0*hrstd/(hrmean*hrmean);

   df = sf/n;
   for(i=1;i<=n;i++) w[i] = (i-1)*2.0*PI*df;
   for(i=1;i<=n;i++) 
   {
      dw1 = w[i]-w1;
      dw2 = w[i]-w2;
      Hw[i] = sig1*exp(-dw1*dw1/(2.0*c1*c1))/sqrt(2*PI*c1*c1) 
            + sig2*exp(-dw2*dw2/(2.0*c2*c2))/sqrt(2*PI*c2*c2); 
   }
   for(i=1;i<=n/2;i++) Sw[i] = (sf/2.0)*sqrt(Hw[i]);
   for(i=n/2+1;i<=n;i++) Sw[i] = (sf/2.0)*sqrt(Hw[n-i+1]);


   /* randomise the phases */
   for(i=1;i<=n/2-1;i++) ph0[i] = 2.0*PI*ran1(&rseed);
   ph[1] = 0.0;
   for(i=1;i<=n/2-1;i++) ph[i+1] = ph0[i];
   ph[n/2+1] = 0.0;
   for(i=1;i<=n/2-1;i++) ph[n-i+1] = - ph0[i]; 


   /* make complex spectrum */
   for(i=1;i<=n;i++) SwC[2*i-1] = Sw[i]*cos(ph[i]);
   for(i=1;i<=n;i++) SwC[2*i] = Sw[i]*sin(ph[i]);

   /* calculate inverse fft */
   dfour1(SwC,n,-1);

   /* extract real part */
   for(i=1;i<=n;i++) rr[i] = (1.0/n)*SwC[2*i-1];

   xstd = stdev(rr,n);
   ratio = rrstd/xstd; 

   for(i=1;i<=n;i++) rr[i] *= ratio;
   for(i=1;i<=n;i++) rr[i] += rrmean;

   freeVect(w,1,n);
   freeVect(Hw,1,n);
   freeVect(Sw,1,n);
   freeVect(ph0,1,n/2-1);
   freeVect(ph,1,n);
   freeVect(SwC,1,2*n);
}

/*--------------------------------------------------------------------------*/
/*    THE ANGULAR FREQUENCY                                                 */
/*--------------------------------------------------------------------------*/

double angfreq(double t)
{
   int i;
  
   i = 1 + (int)floor(t/h);
  
   return 2.0*PI/rrpc[i];
}

/*--------------------------------------------------------------------------*/
/*    THE EXACT NONLINEAR DERIVATIVES                                       */
/*--------------------------------------------------------------------------*/

void derivspqrst(double t0,double x[],double dxdt[])
{
   int i,k;
   double a0,w0,r0,x0,y0,z0;
   double t,dt,dt2,*xi,*yi,zbase;
 
   k = 5; 
   xi = mallocVect(1,k);
   yi = mallocVect(1,k); 
  
   w0 = angfreq(t0);
   r0 = 1.0; x0 = 0.0;  y0 = 0.0;  z0 = 0.0;
   a0 = 1.0 - sqrt((x[1]-x0)*(x[1]-x0) + (x[2]-y0)*(x[2]-y0))/r0;

   for(i=1;i<=k;i++) xi[i] = cos(ti[i]);
   for(i=1;i<=k;i++) yi[i] = sin(ti[i]);   

   zbase = 0.005*sin(2.0*PI*fhi*t0);

   t = atan2(x[2],x[1]);
   dxdt[1] = a0*(x[1] - x0) - w0*(x[2] - y0);
   dxdt[2] = a0*(x[2] - y0) + w0*(x[1] - x0); 
   dxdt[3] = 0.0;  
   for(i=1;i<=k;i++)  
   {
      dt = fmod(t-ti[i],2.0*PI);
      dt2 = dt*dt;
      dxdt[3] += -ai[i]*dt*exp(-0.5*dt2/(bi[i]*bi[i])); 
   }
   dxdt[3] += -1.0*(x[3] - zbase);
  
   freeVect(xi,1,k);
   freeVect(yi,1,k);
}

/*--------------------------------------------------------------------------*/
/*    RUNGA-KUTTA FOURTH ORDER INTEGRATION                                  */
/*--------------------------------------------------------------------------*/

void drk4(double y[], int n, double x, double h, double yout[], 
          void (*derivs)(double, double [], double []))
{
        int i;
        double xh,hh,h6,*dydx,*dym,*dyt,*yt;
 
        dydx=mallocVect(1,n);
        dym=mallocVect(1,n);
        dyt=mallocVect(1,n);
        yt=mallocVect(1,n);

        hh=h*0.5;
        h6=h/6.0;
        xh=x+hh;
        (*derivs)(x,y,dydx);
        for (i=1;i<=n;i++) yt[i]=y[i]+hh*dydx[i];
        (*derivs)(xh,yt,dyt);
        for (i=1;i<=n;i++) yt[i]=y[i]+hh*dyt[i];
        (*derivs)(xh,yt,dym);
        for (i=1;i<=n;i++) {
                yt[i]=y[i]+h*dym[i];
                dym[i] += dyt[i];
        }
        (*derivs)(x+h,yt,dyt);
        for (i=1;i<=n;i++)
                yout[i]=y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);