refl_05.cc

这是一个从音频信号里提取特征参量的程序

  // apply dynamic range threshold
  //
  VectorFloat temp_vector;
  float value;
  
  // add the dynamic range threshold
  //
  if (dyn_range_d < (float)0.0) {
    float scale_factor =
      pow(10.0, ((float)dyn_range_d / Integral::DB_POW_SCALE_FACTOR));
    long num_rows = cor_a.getNumRows();
    temp_vector.setLength(num_rows);
    for (long i = 0; i < num_rows; i++) {
      value = cor_a.getValue(i, i);
      temp_vector(i) = value;
      value *= 1 + scale_factor;
      const_cast<MatrixFloat&>(cor_a).setValue(i, i, value);
    }
  }
  else {
    return Error::handle(name(), L"computeCovarCholesky", ERR_DYNRANGE,
			 __FILE__, __LINE__);
  }
  
  // compute the order
  //
  long refl_order = cor_a.getNumRows() - 1;

  // declare some scratch variables
  //
  MatrixFloat b(refl_order + (long)1, refl_order + (long)1);
  VectorFloat beta(refl_order + 1);
    
  // create space
  //
  if (!pred_coef.setLength(refl_order + 1)) {
    return Error::handle(name(), L"computeCovarCholesky", ERR,
			 __FILE__, __LINE__, Error::WARNING);
  }
  if (!output_a.setLength(refl_order)) {
    return Error::handle(name(), L"computeCovarCholesky", ERR,
			 __FILE__, __LINE__, Error::WARNING);
  }

  // declare local variables
  //
  long m  = refl_order;
  long mf = refl_order;
  long minc;
  long ip;
  float gam;
  long j;
  long jp;
  float s;

  // Note: the equation numbers indicated here are from the COVAR
  // algorithm in following refence:
  //
  //  J.D. Markel and A.H. Gray, "Linear Prediction of Speech,"
  //  Springer-Verlag Berlin Heidelberg, New York, USA, pp. 52, 1980.
  
  // initialization
  //
  b.setValue(1, 1, 1.0);                            // Eq. 3.40
  err_energy_a = cor_a.getValue(1 - 1, 1 - 1);      // Eq. 3.41
  beta(1) = cor_a.getValue(2 - 1, 2 - 1);           // Eq. 3.41
  output_a(0) = -(float)cor_a.getValue(2 - 1, 1 - 1) / 
                    cor_a.getValue(2 - 1, 2 - 1);   // Eq. 3.42
  pred_coef(0) = 1;                               // Eq. 3.43
  pred_coef(1) = output_a((long)0);            // Eq. 3.43
  err_energy_a = err_energy_a - output_a(0) *
    output_a(0) * beta(1);                       // Eq. 3.44

  // iteration
  //
  for (minc = 2; minc <= mf; minc++) {              // do 400 ...
    m = minc - 1;
    b.setValue(minc, minc, 1);                      // Eq. 3.51a

    for (ip = 1; ip <= m; ip++) {                   // do 200 ...
      if ((float)beta(ip) < 0) {                   // if 600, 200, 130
	return Error::handle(name(), L"computeCovarCholesky", ERR_BETA,
			     __FILE__, __LINE__, Error::WARNING);
      }

      // Eq. 3.50 and 3.51b
      //
      else if ((float)beta(ip) > 0) {
        gam = 0;
        for (j = 1; j <= ip; j++) {
          gam += cor_a.getValue(m + 2 - 1, j + 1 - 1) * b.getValue(ip, j);
	}
        gam /= (float)beta(ip);
        for (jp = 1; jp <= ip; jp++) {
	  value = b.getValue(minc, jp) - gam * b.getValue(ip, jp);
	  b.setValue(minc, jp, value);
	}
      }
    }

    // Eq. 3.52
    //
    beta(minc) = 0;
    for (j = 1; j <= minc; j++) {
      beta(minc) += cor_a.getValue(m + 2 - 1, j + 1 - 1) * b.getValue(minc, j);
    }
    m++;
    
    if ((float)beta(m) < 0) {                      // If 600,360,260
      return Error::handle(name(), L"computeCovarCholesky", ERR_BETA,
			   __FILE__, __LINE__, Error::WARNING);
    }

    // Eq. 3.55
    //
    else if ((float)beta(m) > 0) {
      s = 0;
      for (ip = 1; ip <= m; ip++) {
        s += cor_a.getValue(m + 1 - 1, ip - 1) * pred_coef(ip - 1);
      }

      // Eq. 3.56
      //
      output_a(m - 1) = -s / beta(m);
      for (ip = 2; ip <= m; ip++) {
        pred_coef(ip - 1) += output_a(m - 1) * b.getValue(m, ip - 1);
      }
      pred_coef(m) = output_a(m - 1);
    }

    // if beta(m) goes to zero, output zero reflection
    // coefficients. this is not part of Markel and Gray's algorithm.
    //
    else if ((double)beta(m) == 0) {
      err_energy_a.assign(0.0);
      output_a.assign(0.0);      
      break;
    }
    
    // Eq. 3.57
    //
    err_energy_a -= output_a(m - 1) * output_a(m - 1) *
      beta(m);
    
    if (err_energy_a <= (float)0) {
      err_energy_a = 0;
      return Error::handle(name(), L"computeCovarCholesky", ERR_ENERGY,
			   __FILE__, __LINE__, Error::WARNING);
    }
  }

  // resume the values before applying the dynamic threshold
  //
  if (dyn_range_d < (float)0.0) {
    long num_rows =  temp_vector.length();
    for (long i = 0; i < num_rows; i++) {
      const_cast<MatrixFloat&>(cor_a).setValue(i, i, temp_vector(i));
    }
  }
  
  // exit gracefully
  //
  return true;
}

// method: computeLatticeBurg
//
// arguments:
//  VectorFloat& output: (output) reflection coefficients
//  Float& err_energy: (output) error energy  
//  const VectorFloat& input: (input) input data
//
// return: a boolean value indicating status
//
// this method computes the reflection coefficients using the burg
// analysis method. it also outputs the error energy.
//
boolean Reflection::computeLatticeBurg(VectorFloat& output_a,
				       Float& err_energy_a,
				       const VectorFloat& input_a) {

  // declare local variable
  //
  VectorFloat pred_coef;
  
  // resize output vectors
  //
  long pc_order = (long)order_d + 1;

  if (!pred_coef.setLength(pc_order)) {
    return Error::handle(name(), L"computeLatticeBurg", ERR,
			 __FILE__, __LINE__, Error::WARNING);
  }

  // declare local variables:
  //  Float Registers/Accumulators
  //
  float sum, pe, akk, ai, aj; 
  float sn, sd;
  akk = 0;

  // Indices In Autocorrelation
  //
  long i, k;                                   
  long n = input_a.length();

  // find the signal power
  //
  sum = input_a.sumSquare();

  // if the signal power is zero, output zero reflection coefficients
  //
  if (sum == 0) {
    output_a.assign(order_d, 0.0);    
    return true;
  }
  
  // check for divide by zero error
  //
  if (n == 0) {
    return Error::handle(name(), L"computeLatticeBurg", Error::ARG, __FILE__, __LINE__);
  }
  
  // compute error energy
  //
  err_energy_a = sum / n;

  // declare some internal scratch space:
  //  note that this data is declared static so that it can be shared across
  //  all objects.
  //
  float *tmp_r;
  float *tmp_rc;
  float *tmp_a;
  float *tmp_pef;
  float *tmp_per;

  // create new space
  //
  tmp_r = new float[pc_order];
  tmp_rc = new float[pc_order];
  tmp_a = new float[pc_order];
  tmp_pef = new float[input_a.length() + 1];
  tmp_per = new float[input_a.length() + 1];

  // initialize
  //
  float scale_factor = 1.0 +
    pow(10.0, ((float)dyn_range_d / Integral::DB_POW_SCALE_FACTOR));
  pe = sum * scale_factor;
  tmp_r[0] = sum * scale_factor;
  tmp_a[0] = 1.0;

  // main loop
  //
  for (k = 1; k <= order_d; k++) {

    // compute initial prediction errors
    //
    if (k == 1) {
      for (i = 2; i <= n; i++) {
        tmp_pef[i] = (float)input_a(i - 1);
        tmp_per[i] = (float)input_a(i - 2);
      }
    }

    // update prediction errors
    //
    else {
      for (i = n; i >= k + 1; i--) {
        tmp_pef[i] = tmp_pef[i] + akk * tmp_per[i];
        tmp_per[i] = tmp_per[i - 1] + akk * tmp_pef[i - 1];
      }
    }

    // compute new reflection coefficient
    //
    sn = 0;
    sd = 0;
    for (i = k + 1; i <= n; i++) {
      sn = sn - tmp_pef[i] * tmp_per[i];
      sd = sd + tmp_pef[i]*tmp_pef[i] + tmp_per[i]*tmp_per[i];
    }

    // check for divide by zero errors
    //
    if (sd == 0) {
      return Error::handle(name(), L"computeLatticeBurg", Error::ARG, __FILE__, __LINE__);
    }
    
    akk = 2.0 * sn / sd;
    tmp_rc[k] = akk;

    // convert rc[] to a[] using the levinson recursion
    //
    tmp_a[k] = akk;
    for (i = 1; i <= k / 2; i++) {
      ai = tmp_a[i];
      aj = tmp_a[k - i];
      tmp_a[i] = ai + akk * aj;
      tmp_a[k - i] = aj + akk * ai;
    }

    // compute new autocorrelation estimate
    //
    sum = 0;
    for (i = 1; i <= k; i++) {
      sum = sum + tmp_a[i] * tmp_r[k - i];
    }
    tmp_r[k] = sum;

    // update prediction-error power
    //
    pe *= 1 - akk*akk;
    if (pe < 0) {
      break;
    }
  }

  // copy predictor coefficients
  //
  if (!pred_coef.assign(pc_order, tmp_a)) {
    return Error::handle(name(), L"computeLatticeBurg", ERR,
			 __FILE__, __LINE__, Error::WARNING);
  }

  // copy reflection coefficients
  //  
  if (!output_a.assign(order_d, &tmp_rc[1])) {
    return Error::handle(name(), L"computeLatticeBurg", ERR,
			 __FILE__, __LINE__, Error::WARNING);
  }

  // compute the filter gain
  //
  err_energy_a = 1.0;
  for (i = 1; i <= order_d; i++) {
    err_energy_a *= 1 - tmp_rc[i] * tmp_rc[i];
  }
  err_energy_a = sqrt(err_energy_a);

  // clean up old space
  //
  delete [] tmp_r;
  delete [] tmp_rc;
  delete [] tmp_a;
  delete [] tmp_pef;
  delete [] tmp_per;

  // exit gracefully
  //
  if (pe >= (float)0) {
    return true;
  }
  else {
    return Error::handle(name(), L"computeLatticeBurg", ERR_PREDERR,
			 __FILE__, __LINE__, Error::WARNING);
  }
}

// method: computePredictionStepUp
//
// arguments:
//  VectorFloat& refl_coef: (output) reflection coefficients
//  const VectorFloat& pred_coef: (input) prediction coefficients
//
// return: a boolean value indicating status
//
// this method computes the reflection coefficients from the
// prediction coefficents
//
// Reference:
//  J.D. Markel and A.H. Gray, "Linear Prediction of Speech,"
//  Springer-Verlag, New York, pp. 95 - 97 1980.
//
boolean Reflection::computePredictionStepUp(VectorFloat& refl_coef_a,
					    const VectorFloat& pred_coef_a) {

  // declare local variables
  //
  int mp1, mrp1, mr, mm;
  float alpha, d;
  VectorFloat b;

  // set the filter order
  //
  long m = pred_coef_a.length() - 1;
  
  // make a copy of the predictor coefficients
  //
  VectorFloat temp_pred(pred_coef_a);  

  // set the lengths of the vectors
  //
  b.setLength(m);
  refl_coef_a.setLength(m);
  
  mp1 = m + 1;
  alpha = 1.0;

  // compute the inverse of the step-up procedure, that is, given the synthesis
  // filter 1 / A(z), it is possible to determine the corresponding acoustic
  // tube reflection coeficients
  //
  for (int j = 0; j < m; j++) {
    mr = m + 1 - (j + 1);
    mrp1 = mr + 1;
    d = 1.0 - temp_pred(mrp1 - 1) * temp_pred(mrp1 - 1);

    // check for divide by zero error
    //
    if (d == 0) {
      return Error::handle(name(), L"computePredictionStepUp", Error::ARG, __FILE__, __LINE__);
    }
    
    alpha = alpha / d;

    for (int k = 0; k < mr; k++) {
      mm = mr + 2 - (k + 1);
      b(k) = temp_pred(mm - 1);
    }

    for (int k = 0; k < mr; k++) {

    // check for divide by zero error
    //
    if (d == 0) {
      return Error::handle(name(), L"computePredictionStepUp", Error::ARG, __FILE__, __LINE__);
    }
    
      temp_pred(k) = (temp_pred(k) - temp_pred(mrp1 - 1) * b(k)) / d;      
    }
    refl_coef_a(mr - 1) = temp_pred(mrp1 - 1);

    // make sure that the reflection coefficients are not greater that one
    //
    if (abs(refl_coef_a(mr - 1)) >= 1) {
      Error::handle(name(), L"The synthesis filter 1 / A(z) is unstable", \
		    Error::IO, __FILE__, __LINE__, Error::WARNING);      
    }
  }
  
  // exit gracefully
  //
  return true;
}

// method: computeLogAreaKellyLochbaum
//
// arguments:
//  VectorFloat& output: (output) prediction coefficients
//  const VectorFloat& input: (input) log erea ratio coefficients
//
// return: a boolean value indicating status
//
// this method computes the reflection coefficients from the
// log-area-ratio, using LOG_AREA_RATIO algorithm and KELLY_LOCHBAUM
// implementation
//
boolean Reflection::computeLogAreaKellyLochbaum(VectorFloat& output_a,
						const VectorFloat&
						input_a) {

  // declare local variables
  //
  double g;
  long order = input_a.length();
  output_a.setLength(order);
  
  //                             1 - log_area_ratio(i)
  // reflection coefficient(i) = ---------------------
  //                             1 + log_area_ratio(i)
  //
  for (long i = 0; i < order; i++) {
    g = exp(input_a(i));

    // check for divide by zero error
    //
    if ((1 + g) == 0) {
      return Error::handle(name(), L"computeLogAreaKellyLochbaum", Error::ARG, __FILE__, __LINE__);
    }
    output_a(i) = (1 - g) / (1 + g);
  }

  // exit gracefully
  //
  return true;
}