/* 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 #include #include #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]); 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 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 (patrick@mcsharry.net) or G. Clifford (gari@mit.edu)\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 */ }