/* Levinson-Durbin algorithm */
/* return 1 if ok, otherwise 0 */


int levdur(real acorr[], real  coeff[], int order)
{
  /* Local variables */
  static int minc, minc2;
  static float s;
  static int ib, mh;
  static float at;
  static int ip;
  static float rc[20];
  static float alpha;
  static float tmp;

  /* Parameter adjustments */
  --acorr;
  --coeff;

/* .......CHECK FOR ZERO SIGNAL OR ILLEGAL ZERO-LAG AUTOCORRELATION */
  if ((acorr[1] <= (float)0.)||(acorr[order+1] == 0))
    return 0;


/* .......START DURBIN's RECURSION */
  rc[1] = -acorr[2] / acorr[1];
  coeff[1] = (float)1.;
  coeff[2] = rc[1];
  alpha = acorr[1] + acorr[2] * rc[1];
  if (alpha <= (float)0.)
    return 0;

/* Outer loop is recursive */
  for (minc = 2; minc <= order; ++minc) {
    minc2 = minc + 2;
    s = acorr[minc + 1];
    for (ip = 2; ip <= minc; ++ip) {
      s += acorr[minc2 - ip] * coeff[ip];
    }
    rc[minc] = -s / alpha;
    mh = minc / 2 + 1;

/* Inner loop is parallelizable */
    for (ip = 2; ip <= mh; ++ip) {
      ib = minc2 - ip;
      at = rc[minc] * coeff[ib];
      at += coeff[ip];
      tmp = rc[minc] * coeff[ip];
      coeff[ib] += tmp;
      coeff[ip] = at;
    }
    coeff[minc + 1] = rc[minc];
    alpha += rc[minc] * s;

/* ...........IF RESIDUAL ENERGY LESS THAN ZERO (DUE TO ILL-CONDITIONING), */
/* ...........RETURN WITHOUT UPDATING FILTER COEFFICIENTS (USE OLD ONES). */
    if (alpha <= (float)0.)
      return 0;
  }
  return 1;
}