ecgcalc.java

egc 心电信号检测的源程序

/*
 * EcgCalc.java
 *
 * See EcgLicense.txt for License terms.
 */

/**
 *
 * @author  Mauricio Villarroel ([e-mail])
 */

public class EcgCalc {

    /** Creates a new instance of EcgCalc */
    public EcgCalc(EcgParam parameters, EcgLogWindow logOb) {
        paramOb = parameters;
        ecgLog = logOb;

        /* variables for static function ranq() */
        iy=0;
        iv = new long[NTAB];     
    }
    
    public boolean calculateEcg(){
        boolean RetValue;

        ecgLog.println("Starting to calculate ECG....");
        
        RetValue = dorun();

        ecgLog.println("Finished calculating ECG table data.\n");

        return(RetValue);        
    }

    public int getEcgResultNumRows(){
        return ecgResultNumRows;
    }

    public double getEcgResultTime(int index){
        return ecgResultTime[index];
    }

    public double getEcgResultVoltage(int index){
        return ecgResultVoltage[index];
    }

    public int getEcgResultPeak(int index){
        return ecgResultPeak[index];
    }

    /* C defines */
    private final double PI = 2.0*Math.asin(1.0);
    private final int NR_END = 1;
    private final int IA = 16807;
    private final long IM = 2147483647;
    private final double AM = (1.0/IM);
    private final long IQ = 127773;
    private final int IR = 2836;
    private final int NTAB = 32;
    private final double NDIV = (1+(IM-1)/NTAB);
    private final double EPS = 1.2e-7;
    private final double RNMX = (1.0-EPS);

    /*****************************************************************************
     *    DEFINE PARAMETERS AS GLOBAL VARIABLES                                 *
     *****************************************************************************/
    //private String outfile ="ecgsyn.dat";
    // Order of extrema: [P Q R S T]
    private double[] ti = new double[6];  /* ti converted in radians             */
    private double[] ai = new double[6];  /* new calculated a                    */
    private double[] bi = new double[6];  /* new calculated b                    */        

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

    /*
     * Variables for static function ran1()
     */
    private long iy;
    private long[] iv;
    
    /*
     * ECG Result Variables
     */
    /* Result Vectors*/
    private double[] ecgResultTime;
    private double[] ecgResultVoltage;
    private int[] ecgResultPeak;
    
    private int ecgResultNumRows;

    /* Object Variables */
    private EcgParam paramOb;
    private EcgLogWindow ecgLog;    

    /*--------------------------------------------------------------------------*/
    /*    UNIFORM DEVIATES                                                      */
    /*--------------------------------------------------------------------------*/

    private double ran1(){

        int j;
        long k;
        double temp;
        boolean flg;

        if(iy == 0)
            flg = false;
        else
            flg = true;

        if((rseed <= 0) || !flg){
            if (-(rseed) < 1)
                rseed = 1;
            else
                rseed = -rseed;

            for (j=NTAB+7; j>=0; j--) {
                    k=(rseed)/IQ;
                    rseed=IA*(rseed-k*IQ)-IR*k;
                    if (rseed < 0)
                        rseed += IM;
                    if (j < NTAB)
                        iv[j] = rseed;
            }
            iy=iv[0];
        }

        k=(rseed)/IQ;
        rseed=IA*(rseed-k*IQ)-IR*k;
        if (rseed< 0)
            rseed += IM;

        j = (int)(iy/NDIV);
        iy=iv[j];
        iv[j] = rseed;

        if ((temp=AM*iy) > RNMX)
            return RNMX;
        else
            return temp;
    }

    /*
     * FFT
     */
    private void ifft(double[] data, long nn, int isign){

        long n, mmax, m, istep, i, j;
        double wtemp,wr,wpr,wpi,wi,theta;
        double tempr,tempi;
        double swap;

        n=nn << 1;
        j=1;
        for (i=1; i< n; i+=2) {
                if (j > i) {
                    //SWAP(data[j],data[i]);
                    swap = data[(int) j];
                    data[(int)j] = data[(int)i];
                    data[(int)i] = swap;
                    //SWAP(data[j+1],data[i+1]);
                    swap = data[(int)j+1];
                    data[(int)j+1] = data[(int)i+1];
                    data[(int)i+1] = swap;
                }
                m=n >> 1;
                while (m >= 2 && j > m) {
                        j -= m;
                        m >>= 1;
                }
                j += m;
        }
        mmax=2;
        while (n > mmax) {
                istep=mmax << 1;
                theta=isign*(6.28318530717959/mmax);
                wtemp=Math.sin(0.5*theta);
                wpr = -2.0*wtemp*wtemp;
                wpi=Math.sin(theta);
                wr=1.0;
                wi=0.0;
                for (m=1; m<mmax; m+=2) {
                        for (i=m; i<=n; i+=istep) {
                                j= i + mmax;
                                tempr=wr * data[(int)j] - wi * data[(int)j+1];
                                tempi=wr * data[(int)j+1] + wi * data[(int)j];
                                data[(int)j]= data[(int)i] - tempr;
                                data[(int)j+1]= data[(int)i+1] - tempi;
                                data[(int)i] += tempr;
                                data[(int)i+1] += tempi;
                        }
                        wr=(wtemp=wr)*wpr-wi*wpi+wr;
                        wi=wi*wpr+wtemp*wpi+wi;
                }
                mmax=istep;
        }
    }

    /*
     * STANDARD DEVIATION CALCULATOR
     */
    /* n-by-1 vector, calculate standard deviation */
    private double stdev(double[] x, int n){
            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 (Math.sqrt(total/((double)n-1)));
    }        

    /*
     * THE ANGULAR FREQUENCY
     */
    private double angfreq(double t){
       int i = 1 + (int)Math.floor(t/h);
       return(2.0*PI/rrpc[i]);
    }

    /*--------------------------------------------------------------------------*/
    /*    THE EXACT NONLINEAR DERIVATIVES                                       */
    /*--------------------------------------------------------------------------*/
    private void derivspqrst(double t0,double[] x, double[] dxdt){

        int i,k;
        double a0,w0,r0,x0,y0,z0;
        double t,dt,dt2,zbase;
        double[] xi, yi;

        k = 5; 
        xi = new double[k + 1];
        yi = new double[k + 1];
        w0 = angfreq(t0);
        r0 = 1.0; x0 = 0.0;  y0 = 0.0;  z0 = 0.0;
        a0 = 1.0 - Math.sqrt((x[1]-x0)*(x[1]-x0) + (x[2]-y0)*(x[2]-y0))/r0;

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


        zbase = 0.005* Math.sin(2.0*PI*paramOb.getFHi()*t0);

        t = Math.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 = Math.IEEEremainder(t-ti[i], 2.0*PI);
          dt2 = dt*dt;
          dxdt[3] += -ai[i] * dt * Math.exp(-0.5*dt2/(bi[i]*bi[i])); 
        }
        dxdt[3] += -1.0*(x[3] - zbase);
    }

    /*
     * RUNGA-KUTTA FOURTH ORDER INTEGRATION
     */
    private void Rk4(double[] y, int n, double x, double h, double[] yout){
        int i;
        double xh,hh,h6;
        double[] dydx, dym, dyt, yt;

        dydx=  new double[n + 1];
        dym =  new double[n + 1];
        dyt =  new double[n + 1];
        yt  =  new double[n + 1];

        hh= h * 0.5;
        h6= h/6.0;
        xh= x + hh;

        derivspqrst(x,y,dydx);
        for (i=1; i<=n; i++)
            yt[i]=y[i]+hh*dydx[i];

        derivspqrst(xh,yt,dyt);
        for (i=1; i<=n; i++)
            yt[i]=y[i] + hh * dyt[i];

        derivspqrst(xh,yt,dym);
        for (i=1; i<=n; i++){
                yt[i]=y[i] + h * dym[i];
                dym[i] += dyt[i];
        }

        derivspqrst(x+h,yt,dyt);
        for (i=1; i<=n; i++)
                yout[i]=y[i] + h6 * (dydx[i]+dyt[i]+2.0*dym[i]);
    }

    /*
     * GENERATE RR PROCESS
     */
    private 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;
        double[] w, Hw, Sw, ph0, ph, SwC;

        w =  new double[n+1];
        Hw = new double[n+1];
        Sw = new double[n+1];
        ph0= new double[(int)(n/2-1 +1)];
        ph = new double[n+1];
        SwC= new double[(2*n)+1];

        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/(double)n;
        for(i=1; i<=n; i++)
           w[i] = (i-1)*2.0*PI*df;