ldcsub.c

Reference Implementation of G.711 standard and other voice codecs

/* ldcsub.f -- translated by f2c (version of 15 October 1990  19:58:17).
   You must link the resulting object file with the libraries:
	-lF77 -lI77 -lm -lc   (in that order)
*/

#include <stdio.h>
#include <math.h>
#include "ugstdemo.h"

short           round(x)
  double          x;
{
  return (short) (x + 0.5 * ((x > 0) ? 1 : -1));
}

/* Common Block Declarations */
struct
{
  double          wnc;
}               global_;

#define global_1 global_

struct
{
  double          wg[34], rhog, lambdg[10], dltmin, dltmax, doffst, ws[105], rhos,
                  lambds[50], ymin, ymax, codebk[5][1024], dcstat[290];
}               decod_;

#define decod_1 decod_

struct
{
  long            ifailg, ifails;
  double          deltaw[34], rdelta[11], sigma, e[5], yw[105], ry[51];
}               decmon_;

#define decmon_1 decmon_

struct
{
  double          wp[60], rhop, gamma1[10], gamma2[10], dmstat[143];
}               dist_;

#define dist_1 dist_

struct
{
  long            ifailp;
  double          xw[60], rx[11], q[11];
}               dismon_;

#define dismon_1 dismon_

/* Table of constant values */

static double   c_b5 = .90625;
static double   c_b9 = .98828125;
static double   c_b43 = .9;
static double   c_b47 = .6;
static long     c__0 = 0;
static long     c__1 = 1;
static long     c__10 = 10;
static double   c_b282 = 10.;
static long     c__50 = 50;

/* File:            LDCSUB.FOR */
/* Last modified:   January 6, 1992 */
/* By: 		   Finn Tore Johansen */
/* 		   Norwegian Telecom Research */
/* 		   <[e-mail]> */
/* Content:        Fortran subroutines for the LD-CELP verification programs: */
/* 	RESET		Constants definition and state reset subroutine */
/* 	ENCODE		LD-CELP encoder subroutine (using stupid full search) */
/* 	DISTM		LD-CELP distance measure subroutine */
/*       PWFILT          Perceptual weighting filter subroutine */
/* 	DECODE		LD-CELP decoder subroutine */
/* 	PSTFLT		Postfilter subroutine */
/* 	DURBIN		Durbin recursion subroutine */
/* 	HYBWG		Hybrid window for gain predictor */
/* 	HYBWS		Hybrid window for synthesis filter */
/* 	HYBWP		Hybrid window for perceptual weighting filter */
/* 	CBOOK		Shape codebook definition subroutine */
/* Version 06.01.92 / Ftj */


/* Reset decoder and distance measure states, assign constants */
int             reset_(n)
  long           *n;
{
  /* System generated locals */
  long            i_1, i_2, i_3;
  double          d_1;

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

  /* 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 void     cbook_();
  extern void     hybwg_();
  extern void     hybwp_(), hybws_();
  static double   g[8];
  static double   shape[640];
  static long     i, j, m;


  /* Constants */
  global_1.wnc = (double) 257. / 256.;

  /* Decoder parameters */
  *n = 0;
  for (i = 1; i <= 10; ++i)
  {
    d_1 = pow(c_b5, (double) i) * (double) 16384.;
    decod_1.lambdg[i - 1] = (double) round(d_1) / (double) 16384.;
  }

  decod_1.dltmin = (double) 0.;
  decod_1.dltmax = (double) 60.;
  decod_1.doffst = (double) 32.;
  for (i = 1; i <= 50; ++i)
  {
    d_1 = pow(c_b9, (double) i) * (double) 16384.;
    decod_1.lambds[i - 1] = (double) round(d_1) / (double) 16384.;
  }

  decod_1.ymin = (double) -4095.;
  decod_1.ymax = (double) 4095.;
  decod_1.rhog = (double) .75;
  decod_1.rhos = (double) .75;

  /* ...codebook */
  hybws_(decod_1.ws);
  hybwg_(decod_1.wg);
  g[0] = (double) .515625;
  g[1] = g[0] * 7 / (double) 4.;
  g[2] = g[1] * 7 / (double) 4.;
  g[3] = g[2] * 7 / (double) 4.;
  g[4] = -g[0];
  g[5] = -g[1];
  g[6] = -g[2];
  g[7] = -g[3];
  cbook_(shape);
  for (m = 0; m <= 4; ++m)
  {
    for (i = 0; i <= 7; ++i)
    {
      for (j = 0; j <= 127; ++j)
      {
	i_1 = (8 * j) + i;
	i_2 = (5 * j) + m;
	decod_1.codebk[m][i_1] = g[i] * shape[i_2];
      }
    }

    /* State variables */
    /* ...gain predictor */
  }
  /* ...special treatment of initial values to get correct start-up */
  for (i = -34; i <= -1; ++i)
  {
    dltbuf[(i_1 = i + 34)] = (double) 0.;
  }
  dltbuf[33] = -decod_1.doffst;
  dltbuf[32] = -decod_1.doffst;
  dltbuf[31] = -decod_1.doffst;
  for (i = -10; i <= -1; ++i)
  {
    delta[(i_1 = i + 10)] = -decod_1.doffst;
  }
  for (i = 0; i <= 10; ++i)
  {
    rdexp[i] = (double) 0.;
  }
  ag[0] = (double) 1.;
  ag[1] = (double) -1.;

  /* ...synthesis filter */
  for (i = 2; i <= 10; ++i)
  {
    ag[i] = (double) 0.;
  }

  for (i = 0; i <= 50; ++i)
  {
    ryexp[i] = (double) 0.;
  }

  as[0] = (double) 1.;

  for (i = 1; i <= 50; ++i)
  {
    as[i] = (double) 0.;
  }

  /* Perceptual filter parameters */
  /* Constants */
  for (i = -115; i <= 4; ++i)
  {
    ybuf[(i_1 = i + 115)] = (double) 0.;
  }

  for (i = 1; i <= 10; ++i)
  {
    d_1 = pow(c_b43, (double) i) * (double) 16384.;
    dist_1.gamma1[(i_1 = i - 1)] = (double) i_dnnt(&d_1) / (double) 16384.;
  }

  for (i = 1; i <= 10; ++i)
  {
    d_1 = pow(c_b47, (double) i) * (double) 16384.;
    dist_1.gamma2[(i_1 = i - 1)] = (double) round(d_1) / (double) 16384.;
  }

  /* State variables */
  /* ...Perceptual filter */
  dist_1.rhop = (double) .5;
  hybwp_(dist_1.wp);
  for (i = -60; i <= 4; ++i)
  {
    xbuf[(i_1 = i + 60)] = (double) 0.;
  }
  for (i = 0; i <= 10; ++i)
  {
    rxexp[i] = (double) 0.;
  }
  for (i = -10; i <= -1; ++i)
  {
    xp[(i_1 = i + 10)] = (double) 0.;
  }
  for (i = -10; i <= -1; ++i)
  {
    ybufp[(i_1 = i + 10)] = (double) 0.;
  }
  for (i = -10; i <= -1; ++i)
  {
    yp[(i_1 = i + 10)] = (double) 0.;
  }
  q1[0] = (double) 1.;
  for (i = 1; i <= 10; ++i)
  {
    q1[i] = (double) 0.;
  }
  q2[0] = (double) 1.;
  for (i = 1; i <= 10; ++i)
  {
    q2[i] = (double) 0.;
  }
  /* Version 03.10.91 / Ftj */
}
#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 reset() ............................. */


/* LDCELP encoder */
/* Version 03.10.91 / Ftj */
int             encode_(n, x, c)
  long           *n;
  double         *x;
  long           *c;
{
  /* System generated locals */
  long            i_1, i_2;

  /* 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      decode_();
  extern int      pwfilt_();
  static double   di;
  static double   dmin_;
  static double   savdec[290], yi[5];
  static double   savdis[143];
  static long     ci;
  static long     i;

  /* Save decoder and distance measure state */
  for (i = 0; i <= 289; ++i)
  {
    savdec[i] = decod_1.dcstat[i];
  }

  for (i = 0; i <= 142; ++i)
  {
    savdis[i] = dist_1.dmstat[i];
  }

  /* Search for best codeword */
  dmin_ = (double) 1e30;
  for (ci = 0; ci <= 1023; ++ci)
  {
    decode_(n, &ci, yi, &c__0);
    pwfilt_(n, x, yi, &di);
    if (di <= dmin_)
    {
      dmin_ = di;
      *c = ci;
    }

    /* Recall saved states */
    for (i = 0; i <= 289; ++i)
    {
      decod_1.dcstat[i] = savdec[i];
    }
    for (i = 0; i <= 142; ++i)
    {
      dist_1.dmstat[i] = savdis[i];
    }
  }

  /* Pass chosen codeword through decoder and dist.measure to update state */
  decode_(n, c, yi, &c__1);
  pwfilt_(n, x, yi, &di);
}
#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 encode() ............................. */


/* Find distance measure D for a given codeword C and input vector X */
/* Version 21.11.91 / Ftj */
int             distm_(n, x, c, d)
  long           *n;
  double         *x;
  long           *c;
  double         *d;
{
  /* System generated locals */
  long            i_1, i_2;

  /* Builtin functions */
  long            s_rnge();

  /* 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      decode_();
  extern int      pwfilt_();
  static double   dum;
  static double   savdec[290];
  static double   savdis[143];
  static double   y[5];
  static long     i;

  /* Save decoder and distance measure state */
  for (i = 0; i <= 289; ++i)
  {
    savdec[i] = decod_1.dcstat[i];
  }
  for (i = 0; i <= 142; ++i)
  {
    savdis[i] = dist_1.dmstat[i];
  }

  /* Compute distance measure */
  decode_(n, c, y, &c__0);
  pwfilt_(n, x, y, d);

  /* Recall saved states */
  for (i = 0; i <= 289; ++i)
  {
    decod_1.dcstat[i] = savdec[i];
  }
  for (i = 0; i <= 142; ++i)
  {
    dist_1.dmstat[i] = savdis[i];
  }

  /* Update decoder and weighting filter state */
  decode_(n, c, y, &c__1);
  pwfilt_(n, x, y, &dum);
}
#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 distm() ........................... */


/* Compute perceptually weighted distance D between vectors X(n) and Y(n) */
/* Version 18.11.91 / Ftj */
int             pwfilt_(n, x, y, d)
  long           *n;
  double         *x, *y, *d;
{
  /* System generated locals */
  long            i_1, i_2, i_3, i_4;
  double          d_1;

  /* 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   sum;
  static long     i, j, k, m;

  for (j = 0; j <= 4; ++j)
  {
    xbuf[(i_1 = j + 60)] = x[(i_2 = j)];
  }
  for (j = 0; j <= 4; ++j)
  {
    ybufp[(i_1 = j + 10)] = y[(i_2 = j)];
  }

  /* If n=4k+2 for any long k then update perceptual filter: */
  if (*n % 4 == 2)
  {
    for (i = 1; i <= 60; ++i)
    {
      i_3 = -i + 60;
      dismon_1.xw[i - 1] = dist_1.wp[i - 1] * xbuf[i_3];
    }
    for (k = 0; k <= 10; ++k)
    {
      sum = (double) 0.;
      for (i = 31; i <= 50; ++i)
      {
	sum += dismon_1.xw[(i_1 = i - 1)] * dismon_1.xw[(i_2 = i + k - 1)];
      }
      rxexp[k] = dist_1.rhop * rxexp[k] + sum;
    }

    for (k = 0; k <= 10; ++k)
    {
      sum = (double) 0.;
      for (i = 1; i <= 30; ++i)
      {
	sum += dismon_1.xw[i - 1] * dismon_1.xw[i + k - 1];
      }
      dismon_1.rx[k] = rxexp[k] + sum;
    }
    dismon_1.rx[0] = global_1.wnc * dismon_1.rx[0];
    if (dismon_1.rx[10] == 0.)
    {
      dismon_1.ifailp = 3;
    }
    else
    {
      durbin_(dismon_1.rx, &c__10, dismon_1.q, rc, &dismon_1.ifailp);
    }
    if (dismon_1.ifailp == 0)
    {
      for (k = 1; k <= 10; ++k)
      {
	q1[k] = dismon_1.q[k] * dist_1.gamma1[k - 1];
      }
      for (k = 1; k <= 10; ++k)
      {
	q2[k] = dismon_1.q[k] * dist_1.gamma2[k - 1];
      }
    }
  }