ldcsub.c

Reference Implementation of G.711 standard and other voice codecs

  /* Filter input signal through weighting filter: */
  for (m = 0; m <= 4; ++m)
  {
    xp[(i_1 = m + 10)] = xbuf[(i_2 = m + 60)];
    for (k = 1; k <= 10; ++k)
    {
      xp[(i_1 = m + 10)] = xp[(i_2 = m + 10)] +
	q1[k] * xbuf[(i_4 = m - k + 60)];
    }

    for (k = 1; k <= 10; ++k)
    {
      xp[(i_1 = m + 10)] = xp[(i_2 = m + 10)] -
	q2[(i_3 = k)] * xp[(i_4 = m - k + 10)];
    }
  }

  /* Filter output signal through weighting filter */
  for (m = 0; m <= 4; ++m)
  {
    yp[(i_1 = m + 10)] = ybufp[(i_2 = m + 10)];
    for (k = 1; k <= 10; ++k)
    {
      yp[(i_1 = m + 10)] = yp[(i_2 = m + 10)] +
	q1[(i_3 = k)] * ybufp[(i_4 = m - k + 10)];
    }

    for (k = 1; k <= 10; ++k)
    {
      yp[(i_1 = m + 10)] = yp[(i_2 = m + 10)] -
	q2[(i_3 = k)] * yp[(i_4 = m - k + 10)];
    }
  }

  /* Compute rms of weighted error signal: */
  *d = (double) 0.;
  for (j = 0; j <= 4; ++j)
  {
    /* Computing 2nd power */
    d_1 = yp[(i_1 = j + 10)] - xp[(i_2 = j + 10)];
    *d += d_1 * d_1;
  }

  /* Update state: */
  *d = sqrt(*d / 5);
  for (m = -60; m <= -1; ++m)
  {
    xbuf[(i_1 = m + 60)] = xbuf[(i_2 = m + 65)];
  }

  for (m = -10; m <= -1; ++m)
  {
    ybufp[(i_1 = m + 10)] = ybufp[(i_2 = m + 15)];
  }

  for (m = -10; m <= -1; ++m)
  {
    xp[(i_1 = m + 10)] = xp[(i_2 = m + 15)];
  }

  for (m = -10; m <= -1; ++m)
  {
    yp[(i_1 = m + 10)] = yp[(i_2 = m + 15)];
  }
}
#undef dltbuf
#undef xp
#undef yp
#undef as
#undef ag
#undef ryexp
#undef rxexp
#undef q2
#undef q1
#undef ybufp
#undef rdexp
#undef delta
#undef ybuf
#undef xbuf
/* ............................ End of pwfilt() ........................... */


/* LDCELP decoder */
/* Version 29.11.91 / Ftj */
int             decode_(n, c, y, mode)
  long           *n, *c;
  double         *y;
  long           *mode;
{
  /* System generated locals */
  long            i_1, i_2, i_3, i_4;
  double          d_1;

  /* Builtin functions */
  int             s_stop();
  long            s_rnge();
  double          pow_dd(), d_lg10();

  /* Local variables */
#define ag ((double *)&decod_1 + 5372)
#define as ((double *)&decod_1 + 5445)
#define delta ((double *)&decod_1 + 5326)
#define dltbuf ((double *)&decod_1 + 5337)
#define q1 ((double *)&dist_1 + 81)
#define q2 ((double *)&dist_1 + 92)
#define rdexp ((double *)&decod_1 + 5383)
#define rxexp ((double *)&dist_1 + 133)
#define ryexp ((double *)&decod_1 + 5394)
#define xbuf ((double *)&dist_1 + 144)
#define xp ((double *)&dist_1 + 209)
#define ybuf ((double *)&decod_1 + 5496)
#define ybufp ((double *)&dist_1 + 103)
#define yp ((double *)&dist_1 + 118)
  extern int      durbin_();
  static double   rc[50];
  static double   tmp;
  static long     i, j, k, m;

  /* Test if valid code file */
  if (*c > 1023 || *c < 0)
  {
    fprintf(stderr, "Illegal codeword encountered\n");
    exit(-1);
  }
  if (*n % 4 == 1)
  {
    for (i = 1; i <= 34; ++i)
    {
      decmon_1.deltaw[(i_1 = i - 1)] = decod_1.wg[(i_2 = i - 1)] *
	dltbuf[(i_3 = -i + 34)];
    }

    for (k = 0; k <= 10; ++k)
    {
      tmp = (double) 0.;
      for (i = 21; i <= 24; ++i)
      {
	tmp += decmon_1.deltaw[i - 1] * decmon_1.deltaw[i + k - 1];
      }
      rdexp[k] = decod_1.rhog * rdexp[k] + tmp;
    }

    for (k = 0; k <= 10; ++k)
    {
      tmp = (double) 0.;
      for (i = 1; i <= 20; ++i)
      {
	tmp += decmon_1.deltaw[i - 1] * decmon_1.deltaw[i + k - 1];
      }
      decmon_1.rdelta[k] = rdexp[k] + tmp;
    }

    decmon_1.rdelta[0] = global_1.wnc * decmon_1.rdelta[0];
    if (decmon_1.rdelta[10] == 0.)
    {
      decmon_1.ifailg = 3;
    }
    else
    {
      durbin_(decmon_1.rdelta, &c__10, ag, rc, &decmon_1.ifailg);
    }
    if (decmon_1.ifailg == 0)
    {
      for (k = 1; k <= 10; ++k)
      {
	ag[k] = ag[k] * decod_1.lambdg[k - 1];
      }
    }
  }

  /* Perform logarithmic gain prediction and excitation vector computation: */
  tmp = decod_1.doffst;
  for (k = 1; k <= 10; ++k)
  {
    tmp -= ag[k] * delta[-k + 10];
  }
  if (tmp < decod_1.dltmin)
  {
    tmp = decod_1.dltmin;
  }
  if (tmp > decod_1.dltmax)
  {
    tmp = decod_1.dltmax;
  }
  tmp /= (double) 20.;
  decmon_1.sigma = pow((double) 10., tmp);

  for (j = 0; j <= 4; ++j)
  {
    decmon_1.e[j] = decmon_1.sigma * decod_1.codebk[j][*c];
  }
  tmp = (double) 0.;
  for (j = 0; j <= 4; ++j)
  {
    /* Computing 2nd power */
    d_1 = decmon_1.e[j];
    tmp += d_1 * d_1;
  }
  tmp /= 5;
  if (tmp < 1.)
  {
    tmp = (double) 1.;
  }

  /* If n=4k+2 for any long k then update synthesis filter: */
  delta[10] = log10(tmp) * 10 - decod_1.doffst;
  if (*n % 4 == 2)
  {
    for (i = 1; i <= 105; ++i)
    {
      decmon_1.yw[(i_1 = i - 1)] = decod_1.ws[(i_2 = i - 1)] *
	ybuf[(i_3 = -10 - i + 115)];
    }
    for (k = 0; k <= 50; ++k)
    {
      tmp = (double) 0.;
      for (i = 36; i <= 55; ++i)
      {
	tmp += decmon_1.yw[i - 1] * decmon_1.yw[i + k - 1];
      }
      ryexp[k] = decod_1.rhos * ryexp[k] + tmp;
    }

    for (k = 0; k <= 50; ++k)
    {
      tmp = (double) 0.;
      for (i = 1; i <= 35; ++i)
      {
	tmp += decmon_1.yw[i - 1] * decmon_1.yw[i + k - 1];
      }
      decmon_1.ry[k] = ryexp[k] + tmp;
    }

    decmon_1.ry[0] = global_1.wnc * decmon_1.ry[0];
    if (decmon_1.ry[50] == 0.)
    {
      decmon_1.ifails = 3;
    }
    else
    {
      durbin_(decmon_1.ry, &c__50, as, rc, &decmon_1.ifails);
    }
    if (decmon_1.ifails == 0)
    {
      for (k = 1; k <= 50; ++k)
      {
	as[k] *= decod_1.lambds[k - 1];
      }
    }
  }

  /* Perform synthesis filtering: */
  for (m = 0; m <= 4; ++m)
  {
    ybuf[(i_1 = m + 115)] = decmon_1.e[m];
    for (k = 1; k <= 50; ++k)
    {
      ybuf[m + 115] = ybuf[(i_2 = m + 115)] - as[k] * ybuf[m - k + 115];
    }
  }

  /* Limit output signal when not in codebook search mode: */
  if (*mode == 1)
  {
    for (m = 0; m <= 4; ++m)
    {
      if (ybuf[m + 115] < decod_1.ymin)
      {
	ybuf[m + 115] = decod_1.ymin;
      }
      if (ybuf[m + 115] > decod_1.ymax)
      {
	ybuf[m + 115] = decod_1.ymax;
      }
    }
  }

  for (j = 0; j <= 4; ++j)
  {
    y[j] = ybuf[j + 115];
  }

  /* Update signal buffers: */
  for (m = -115; m <= -1; ++m)
  {
    ybuf[m + 115] = ybuf[m + 120];
  }

  for (m = -10; m <= -1; ++m)
  {
    delta[m + 10] = delta[m + 11];
  }

  dltbuf[34] = delta[10];
  for (m = -34; m <= -1; ++m)
  {
    dltbuf[m + 34] = dltbuf[m + 35];
  }
}
#undef dltbuf
#undef xp
#undef yp
#undef as
#undef ag
#undef ryexp
#undef rxexp
#undef q2
#undef q1
#undef ybufp
#undef rdexp
#undef delta
#undef ybuf
#undef xbuf
/* ........................... End of decode() ........................... */


/* Version: 21.11.91 / Ftj */
int             pstflt_(nvect, sx, sy)
  long           *nvect;
  double         *sx, *sy;
{
  /* System generated locals */
  long            i_1, i_2, i_3, i_4, i_5;
  double          d_1;

  /* Builtin functions */
  long            s_rnge();
  double          pow_di();
  long            i_dnnt();

  /* Local variables */
  static long     idim;
  static double   cmax;
  static long     kmax;
  static double   ptap, temp[5], ybuf[110];
  static long     itmp;
  static double   rhos, b, d[240];
  static long     i, j, k, m, n;
  static double   scale;
  static long     kpmin, kpmax;
  static double   ppfth, tapth, tiltf, stlpf[3];
  static long     kptmp;
  extern void     hybws_();
  static long     m1, m2, nfrsz;
  static double   tiltz, ryexp[51];
  static long     npwsz;
  static double   agcfac, al[3], bl[4], ap[11], gl, rc[10], sd[245];
  static long     ip, kp;
  static double   az[11], scalfl;
  static long     ifails;
  static double   ws[105], ry[51];
  static long     kpdlta;
  static double   yw[105], spfpcf;
  extern int      durbin_();
  static double   ppfzcf, cormax, spfzcf, stlpci[10], sumfil, sunfil, sppcfv[10];
  static long     icount;
  static double   stpfir[10], stpiir[10], spzcfv[10];
  static long     kp1;
  static double   dec[60], apf[11], tap, wnc, spf[5], tmp, sum, tap1;

  /* Preprocessing */
  if (*nvect == 0)
  {
    /* ...reset constant values */
    agcfac = (double) .99;
    idim = 5;
    kpdlta = 6;
    kpmin = 20;
    kpmax = 140;
    nfrsz = 20;
    npwsz = 100;
    ppfth = (double) .6;
    ppfzcf = (double) .15;
    spfpcf = (double) .75;
    spfzcf = (double) .65;
    tapth = (double) .4;
    /* ...reset initial values */
    tiltf = (double) .15;
    al[0] = (double) -2.34036589;
    al[1] = (double) 2.01190019;
    al[2] = (double) -.614109218;
    ap[0] = (double) 1.;
    az[0] = (double) 1.;
    for (i = 1; i < 11; ++i)
    {
      ap[i] = az[i] = (double) 0.;
    }

    for (i = -34; i <= 25; ++i)
    {
      dec[i + 34] = (double) 0.;
    }

    for (i = -139; i <= 100; ++i)
    {
      d[i + 139] = (double) 0.;
    }

    bl[0] = (double) .0357081667;
    bl[1] = (double) -.0069956244;
    bl[2] = (double) -.0069956244;
    bl[3] = (double) .0357081667;
    ip = 85;
    kp1 = 50;
    scalfl = (double) 1.;
    for (i = 0; i < 3; ++i)
    {
      stlpf[i] = (double) 0.;
    }
    for (i = 0; i < 10; ++i)
    {
      stlpci[i] = stpfir[i] = stpiir[i] = (double) 0.;
    }

    /* ...reset internal decoder state */
    gl = (double) 1.;
    wnc = (double) 1.00390625;
    rhos = (double) .75;
    hybws_(ws);
    for (i = 0; i <= 50; ++i)
    {
      ryexp[i] = (double) 0.;
    }
    for (i = -105; i <= 4; ++i)
    {
      ybuf[i + 105] = (double) 0.;
    }
    for (i = 1; i <= 10; ++i)
    {
      d_1 = pow(spfpcf, (double) i) * 16384;
      sppcfv[i - 1] = (double) round(d_1) / (double) 16384.;
      d_1 = pow(spfzcf, (double) i) * 16384;
      spzcfv[i - 1] = (double) round(d_1) / (double) 16384.;
    }
  }

  /* ...update 10'th order LPC coefficients at first vector in frame */
  if (*nvect % 4 == 0)
  {
    for (i = 1; i <= 105; ++i)
    {
      yw[i - 1] = ws[i - 1] * ybuf[-i + 105];
    }

    for (k = 0; k <= 50; ++k)
    {
      tmp = (double) 0.;
      for (i = 36; i <= 55; ++i)
      {
	tmp += yw[i - 1] * yw[i + k - 1];
      }
      ryexp[k] = rhos * ryexp[k] + tmp;
    }

    for (k = 0; k <= 50; ++k)
    {
      tmp = (double) 0.;
      for (i = 1; i <= 35; ++i)
      {
	tmp += yw[i - 1] * yw[i + k - 1];
      }
      ry[k] = ryexp[k] + tmp;
    }
    ry[0] = wnc * ry[0];

    if (ry[50] != 0.)
    {
      durbin_(ry, &c__10, apf, rc, &ifails);
    }
  }

  /* ...update decoded signal buffer */
  for (k = 0; k <= 4; ++k)
  {
    ybuf[k + 105] = sx[k];
  }

  for (m = -105; m <= -1; ++m)
  {
    ybuf[m + 105] = ybuf[m + 110];
  }

  /* ...transfer input vector */
  i_1 = idim;
  for (k = 1; k <= i_1; ++k)
  {
    sd[k + 239] = sx[k - 1];
  }

  /* Start processing according to postfilter description: */
  /* Compute the LPC prediction residual for the current decoded */
  /* speech vector. */
  icount = *nvect % 4 + 1;

  /* check & update IP */
  if (ip == npwsz)
  {
    ip = npwsz - nfrsz;
  }
  i_1 = idim;
  for (k = 1; k <= i_1; ++k)
  {
    itmp = ip + k;
    d[itmp + 13] = sd[k + 239];
    for (j = 10; j >= 2; --j)
    {
      /* FIR filtering. */
      d[itmp + 139] = d[itmp + 139] + stlpci[j - 1] * apf[j];
      /* Memory shift. */
      stlpci[j - 1] = stlpci[j - 2];
    }

    /* Handle last one. */
    d[itmp + 139] = d[itmp + 139] + stlpci[0] * apf[1];

    /* shift in input. */
    stlpci[0] = sd[k + 239];
  }

  /* Extract the pitch period from the LPC prediction residual */
  /* update IP. */
  ip += idim;

  /* lowpass filtering & 4:1 downsampling. */
  if (icount == 3)
  {
    i_1 = npwsz;
    for (k = npwsz - nfrsz + 1; k <= i_1; ++k)
    {
      /* IIR filter */
      tmp = d[k + 139] - stlpf[0] * al[0] - stlpf[1] * al[1] - stlpf[2] * al[2];
      if (k % 4 == 0)
      {
	/* do FIR filtering only if needed. */
	n = k / 4;
	dec[n + 34] = tmp * bl[0] + stlpf[0] * bl[1] +
	  stlpf[1] * bl[2] + stlpf[2] * bl[3];
      }

      /* shift lowpass filter memory. */
      stlpf[2] = stlpf[1];