algorithmlp.java,v

完整的模式识别库

				    data_points_d.color_dset3);
       }
       
       if(set4_d.size() > 0 )
       {	
	   display_set4.removeAllElements();
	   iset4.removeAllElements();
	   mset4.removeAllElements();

	   interpol(set4_d, display_set4);
	   average4 = mean(display_set4, mset4);
	   interpol(mset4, iset4);
   
	   output_panel_d.addOutput(set4_d, (Classify.PTYPE_INPUT * 4), 
				    data_points_d.color_dset4);
       }

	pro_box_d.setProgressCurr(1);
	output_panel_d.repaint();
	return true;
    }
 
    /**
     *
     * Calculates the interpolated points for the data inputs   
     *
     * @@param  v  input data points
     * @@param  iset  interpolated data points
     *
     */
    public void interpol(Vector<MyPoint> v , Vector<MyPoint> iset)
    {	
	double sample;
	double d ;
	MyPoint r = new MyPoint();	
	
	double[] y2 = new double[v.size()];  
	double[] y  = new double[v.size()];
	double[] x  = new double[v.size()];
       
       for (int k = 0; k < v.size(); k++)
       {
	   x[k] = ((MyPoint)v.elementAt(k)).x;
	   y[k] = ((MyPoint)v.elementAt(k)).y;
       }

       // Reference:
       // "Numerical Recipe in C", Cambridge University Press, pp 113-116
       
       // total time interval of the data points calculated
       //
       d = ((MyPoint)v.lastElement()).x - ((MyPoint)v.firstElement()).x ; 
       
       // time step calculation
       //
       sample = d / scale;
       spline(x, y, y2, v.size());
       
       // iset holds new data points - interpolated values
       // First value added 
       //
       iset.addElement(new MyPoint((MyPoint)v.firstElement()));
       
       // Increment to next Time value
       //
      	r.x = sample + ((MyPoint)v.firstElement()).x;

	while (r.x <=  ((MyPoint)v.lastElement()).x)
	{
	    for(int j = 1; j < v.size(); j++)
	    {
		// the New Interpolated amplitude will be calculated
		// for the two points within which the time intervals 
		// lies, Hence, the condition check
		//
		if((r.x <= ((MyPoint)v.elementAt(j)).x ) && 
		   (r.x >=  ((MyPoint)v.elementAt(j - 1)).x)) 
		{
		    r.y = 0 ;
		    splint ((MyPoint)v.elementAt(j - 1), 
			    (MyPoint)v.elementAt(j), r, y2, (j - 1));

		    // New value appended to the array
		    //
		    iset.addElement(new MyPoint(r));
		    break; 
		}
	    }	    
   	    r.x = r.x + sample;
	}
	// sometimes the last element of the user entered 
	// values is missed out hence the condition to 
	// cross-check and put it in
	//
	if (((MyPoint)iset.lastElement()).x != ((MyPoint)v.lastElement()).x )
	    iset.addElement( new MyPoint((MyPoint)v.lastElement()));
    }

    /**
     *
     * Actually interpolates the points
     *
     * @@param x    array containing the x coordinates of datapoints
     * @@param y    array containing the y coordinates of datapoints
     * @@param y2   array containing the interpolated y coordinates
     * @@param size the size of the array to be interpolated
     * 
     */
    public void spline (double[] x, double[] y, double[] y2, int size )
    { 
	double[] u = new double[size];
	double p; 
	double sig; 
	
	y2[0] = 0;
	u[0] = 0;
	 
	for( int i = 1; i < size - 1; i++ )
	{ 
	    sig = ( x[i] - x[i - 1] ) / ( x[i + 1] - x[i - 1] ) ;
	    p = sig * y2[i - 1] + 2 ;
	    y2[i] = ( sig - 1) / p ;
	    u[i] = ( y[i + 1] - y[i] ) / ( x[i + 1] - x[i] ) ;
	    u[i] = ( 6 * u[i] / ( x[i + 1] - x[i - 1] ) 
		     - sig * u[i - 1] ) / p ;
	}
	y2[size - 1] = 0 ;
	
	// This is the back substitution loop of the tridiagonal algorithm
	//
	for(int i = size - 2; i >= 0; i-- )     
	    y2[i] = y2[i] * y2[i + 1] + u[i] ;
    }

    /**
     * Interpolates for a point between the two known points
     * using Cubic Interpolation
     *
     * @@param   u1 start point for the interpolation
     * @@param   u2 end point for the interpolation
     * @@param   r  returning point, basically the interpolated point
     * @@param   y2 array used for reassigning of r
     * @@param   i  the sample number
     *
     */
    public void splint ( MyPoint u1, MyPoint u2, MyPoint r, 
			 double[] y2, int i )
    {
	double h;
	double a;
	double b;
 
	h = u2.x - u1.x ;
	a = ( u2.x - r.x ) / h ;
	b = ( r.x - u1.x ) / h ;
	r.y = a * u1.y + b * u2.y + ((a * a * a - a) * y2[i] 
	                             + ( b * b * b - b ) * y2[i + 1] ) 
	                                                 * ( h * h ) / 6 ; 
    }


    /**
     *
     * Calculates the mean and the zero-mean data points  
     *
     * @@param    v orginal datapoints
     * @@param    mv zero mean datapoints 
     *
     * @@return   The average in a double
     * 
     */
    public double mean(Vector<MyPoint> v, Vector<MyPoint> mv)
    {
        double average = 0;
	MyPoint p      = new MyPoint();

       for(int i = 0; i < v.size(); i++ )
	   average = average + ((MyPoint)v.elementAt(i)).y;
       average = average / v.size();
       
       for(int i = 0; i < v.size(); i++ )
       {
	   p.y = ((MyPoint)v.elementAt(i)).y - average;
	   p.x =  ((MyPoint)v.elementAt(i)).x; 
	   mv.addElement(new MyPoint(p));
       }        
       return average ;
    }  
   
    /**
    * Computes the Auto Correlation and  Linear Coefficients and 
    * displays the coefficients.
    *
    * @@return    True
    *
    */
    boolean step2()
    {
	autoCorrelation();
	lpcCoefficient();
	final_estimate();

	// need a subroutine here to calculate the actual error:
	// 1. Get first two points on the estimated waveform.
	// 2. Check the point given by the user for the x-coordinate,
	//    if the x-coordinate of the user point is same as 
	//    either of two points [from one], 
	// 2a. YES. note the difference in the y-coordinate, this is the error.
	// 2b. Error Energy is square of this term.
	// 2c. NO. check if the user point is outside these two points
	// 2c.1. Yes, pick the next point in the estimated waveform, 
	//       and repeat from step 2 again.
	// 2c.2. No, note the x-coordinate of the user point, 
	//       find the y-coordinate lying on the straight line between
	//       the two estimated points.
	// 2c.2a. Find the difference in the y-coordinates calculate in
	//        step 2c.2 and the user given. This is the error.
	// 2c.3 The Error Energy is the square of this term.
	// 3. Continue with step 2 
	//    till you have covered all the user defined points.
	if(set1_d.size() > 0)
	{
	    actual_err_1 = (double) 0;
	    actual_err_1 = actual_error(y_estimate1, set1_d); 
	}
	if(set2_d.size() > 0)
	{
	    actual_err_2 = (double) 0;
	    actual_err_2 = actual_error(y_estimate2, set2_d); 
	}    
	if(set3_d.size() > 0)
	{
	    actual_err_3 = (double) 0;
	     actual_err_3 = actual_error(y_estimate3, set3_d); 
	}
	if(set4_d.size() > 0)
	{
	    actual_err_4 = (double) 0;
	    actual_err_4 = actual_error(y_estimate4, set4_d); 
	}
	
	step2_display();
	
	return true;
    }    

    /**
    *
    * Computes the autocorrelation coeffient from the data sets   
    *
    */
    public void autoCorrelation()
    {
	if (iset1.size() > 0 )
	{
	    auto_co_1 = new double[lporder + 1];
	    autocorrelate(iset1, auto_co_1);
	}	

	if (iset2.size() > 0)
        {
	    auto_co_2 = new double[lporder + 1];
	    autocorrelate(iset2, auto_co_2);
	}	

	if (iset3.size() > 0)
        {
	    auto_co_3 = new double[lporder + 1];
	    autocorrelate(iset3, auto_co_3);
	}	

	if (iset4.size() > 0)
        {
	    auto_co_4 = new double[lporder + 1];
	    autocorrelate(iset4, auto_co_4);
	}	
    }

    /**
     *
     * Actaully computes the autocorrelation coefficients
     *
     * @@param  v Vector of datapoints
     * @@param  autoCoeff_co array of autocorrelation coefficients
     *
     */
    public void autocorrelate ( Vector v, double[] autoCoeff_co) 
    {
	int size;
	size                 = v.size();
      	double[] amplitude_y = new double[size];

      	// to store the amplitude for calculation autocoefficient
      	//
    	for (int i = 0; i < size; i++)
	    amplitude_y[i] = ((MyPoint)v.elementAt(i)).y;

      	for (int j = 0; j <= lporder; j++)
      	{
      	    autoCoeff_co[j] = 0;
      	    if ( j < size)
      		for (int k = 0; k < size - j; k++)
		    autoCoeff_co[j] = autoCoeff_co[j] 
			            + amplitude_y[k] * amplitude_y[k + j];
	}

   	// normalize the autocorrelation cofficient
      	//
     	for (int i = lporder; i >= 0; i--) 
      	    autoCoeff_co[i] = autoCoeff_co[i] / autoCoeff_co[0];	  
    }

    /**
     *
     * Computes the Linear Prediction coefficient from the data sets   
     *
     */
    public void lpcCoefficient()
    {
	
	// Rabiner, Schafer, "Digital Processing of Speech signals", 
	// Bell Laboratories, Inc. 1978, 
	// section 8.3.2 Durbin's recursive solution for 
	// the autocorrelation equations pp 411-412
	
	if (set1_d.size() > 0)
	{
	    final_lpc_1 = new double [ lporder + 1];
	    ref_coeff_1 = new double [ lporder + 1];
	    estimate_err_1 = calculate_lpc( auto_co_1, 
					    final_lpc_1, ref_coeff_1);
	}

	if (set2_d.size() > 0)
	{
	    final_lpc_2 = new double [ lporder + 1];
	    ref_coeff_2 = new double [ lporder + 1];
	    estimate_err_2 = calculate_lpc( auto_co_2, 
					    final_lpc_2, ref_coeff_2);
	}

	if (set3_d.size() > 0)
	{
	    final_lpc_3 = new double [ lporder + 1];
	    ref_coeff_3 = new double [ lporder + 1];
	    estimate_err_3 = calculate_lpc( auto_co_3, 
					    final_lpc_3, ref_coeff_3);
	}

	if (set4_d.size() > 0)
	{
	    final_lpc_4 = new double [ lporder + 1];
	    ref_coeff_4 = new double [ lporder + 1];
	    estimate_err_4 = calculate_lpc( auto_co_4, 
					    final_lpc_4, ref_coeff_4);
	}
    }

    /**
     *
     * Actually calculate the LP coefficient and the Residual Error 
     * Energy, and Reflection Coefficients
     *
     * @@param  auto_coeff array of auto correlation coefficients
     * @@param  lpc        array of linear prediction coefficients
     * @@param  rc_reg     array of reflection coefficients
     * 
     * @@return residual error energy in a double
     *
     */
    public double calculate_lpc (double[] auto_coeff, 
				 double[] lpc, double[] rc_reg)
    {
	int j_bckwd       = 0;
	long middle_index = 0;

	// Reference:
	// J.D. Markel and A.H. Gray, "Linear Prediction of Speech,"
	// Springer-Verlag Berlin Heidelberg, New York, USA, pp. 219,
	// 1980.
 
	// initialization
	//
	lpc[0] = 1;
	double err_energy = auto_coeff[0];

	// if the energy of the signal is zero we set all the predictor
	// coefficients to zero and exit
	//
	if (err_energy == (double)0)
	{
	    for (int i = 0; i < lporder + 1 ; i++)
		lpc[i] = 0;
	}
	else
        {
	    // do the first step manually
	    //
	    double sum = 0;
	    double tmp = 0;
	    
	    sum         = -auto_coeff[1] / err_energy;
	    rc_reg[1]   = sum;
	    lpc[1]      = rc_reg[1];
	    tmp         = 1 - rc_reg[1] * rc_reg[1];
	    err_energy *= tmp;
	    // recursion
	    //
	    for (int i = 2; i <= lporder; i++)
	    {
		sum = 0;
		for (int j = 1; j < i; j++)
		    sum += lpc[j] * auto_coeff[i - j];
 
		rc_reg[i]    = -(auto_coeff[i] + sum) / err_energy;	    
		lpc[i]       = rc_reg[i];
		j_bckwd      = i - 1;
		middle_index = i / 2;
		for (int j = 1; j <= middle_index; j++)
		{
		    sum          = lpc[j_bckwd] + rc_reg[i] * lpc[j];
		    lpc[j]       = lpc[j] + rc_reg[i] * lpc[i - j];
		    lpc[j_bckwd] = sum;
		    j_bckwd--;
		}
		// compute new error
		//
		tmp = 1.0 - rc_reg[i] * rc_reg[i];
		err_energy *= tmp;
	    }
	}
	return err_energy;
    }

    /**
     *
     * Calculates the estimated points for the data inputs