mnru.c

Reference Implementation of G.711 standard and other voice codecs

        [2] Press,W.H; Flannery,B.P; Teukolky,S.A.; Vetterling, W.T.;
           "Numerical Recipes in C: The Art of Scientific Computing";
            Cambridge University Press, Cambridge; 1990, 735 pp.
           (ISBN 0-521-35465-X)


        Prototype: MNRU.H
        ~~~~~~~~~~

        History:
        ~~~~~~~~
        27.Jan.92  1.0	Adaptation of [2]'s implementation for UGST
                        MNRU module.<[e-mail]>

=============================================================================
*/

#define MBIG 1000000000
#define MSEED 161803398
#define MZ 0
#define FAC (1.0/MBIG)
#define ITER_NO 47

float           ori_random_MNRU(mode, r, seed)
  char           *mode;
  RANDOM_state   *r;
  long           seed;
{
  long            mj, mk;
  long            i, ii, k, iter;
  float           tmp;


/*
 *   RESET OF RANDOM SEQUENCE
 */

  if (*mode == RANDOM_RESET)	/* then reset sequence */
  {
    /* Toogle mode from reset to run */
    *mode = RANDOM_RUN;

    /* Initialize ma[55] using `seed' and `MSEED' */
    mj = MSEED - (seed < 0 ? -seed : seed);
    mj %= MBIG;
    r->ma[55] = mj;
    mk = 1;

    /* Now initialize the rest of the table ma with numbers that are not
     * specially random, in a slightly random order */
    for (i = 1; i <= 54; i++)
    {
      ii = (21 * i) % 55;
      r->ma[ii] = mk;
      mk = mj - mk;
      if (mk < MZ)
	mk += MBIG;
      mj = r->ma[ii];
    }

    /* Warming-up the generator */
    for (k = 1; k <= 4; k++)
      for (i = 1; i <= 55; i++)
      {
	r->ma[i] -= r->ma[1 + (i + 30) % 55];
	if (r->ma[i] < MZ)
	  r->ma[i] += MBIG;
      }

    /* Prepar indices for 1st.generated number */
    r->inext = 0;
    r->inextp = 31;		/* The constant 31 is special; see [1] */
    r->idum = 1;
  }


/*
 *  REAL START (after initialization)
 */

  /* Accumulate samples to make an approxiamtion of the 'central limit' */
  for (tmp = 0, iter = 0; iter < ITER_NO; iter++)
  {
    /* Increment inext,inextp (mod 55) */
    if (++r->inext == 56)
      r->inext = 1;
    if (++r->inextp == 56)
      r->inextp = 1;

    /* Generate a new random number, subtractively */
    mj = r->ma[r->inext] - r->ma[r->inextp];

    /* Check range */
    if (mj < MZ)
      mj += MBIG;

    /* Save and return random number */
    r->ma[r->inext] = mj;
    tmp += (mj * FAC - 0.5);
  }
  return (tmp);
}

#undef ITER_NO
#undef MBIG
#undef MSEED
#undef MZ
#undef FAC
/*  .................... End of ori_random_MNRU() ....................... */

#endif /* *********************** STL92_RNG ****************************** */



/*
  ==========================================================================

        double *MNRU_process (char operation, MNRU_state *s,
        ~~~~~~~~~~~~~~~~~~~~  float *input, float *output,
                              long n, long seed, char mode, double Q)

        Description:
        ~~~~~~~~~~~~

        Module for addition of modulated noise to a vector of `n' samples,
        according to ITU-T Recommendation P.81, for the
        narrow-band model. Depending on the `mode', it:

        - add modulated noise to the `input' buffer at a SNR level of
          `Q' dB, saving to `output' buffer (mode==MOD_NOISE);

        - put into `output' only the noise, without adding to the original
          signal (mode==NOISE_ONLY);

        - copy to `output' the `input' samples (mode==SIGNAL_ONLY);

        There is the need of state variables, which are declared in MNRU.H.
        These are reset calling the function with the argument `operation'
        set as MNRU_START. In the last call of the function, call it with
        operation=MNRU_STOP, to release the memory allocated for the
        processing. Normal operation is followed when operation is set as
        MNRU_CONTINUE.

        Valid inputs are:
        operation:    MNRU_START, MNRU_CONTINUE, MNRU_STOP (see description
                      above; defined in MNRU.H);
        s:	      a pointer to a structure defined as MNRU_state, as in
        	      MNRU.H;
        input:        pointer to input float-data vector; must represent
                      8 kHz speech samples.
        output:	      pointer to output float-data vector; will represent
                      8 kHz speech samples.
        n:	      long with the number of samples (float) in input;
        seed:	      initial value for random number generator;
        mode:	      operation mode: MOD_NOISE, SIGNAL_ONLY, NOISE_ONLY
        	      (see description above; defined in MNRU.H);
        Q:	      double defining the desired value for the signal-to-
                      modulated-noise for the output data.

        ==================================================================
        NOTE! New values of `seed', `mode' and `Q' are considered only
              when operation==MNRU_START, because they are considered as
              INITIAL state values.
        ==================================================================

        For more details on the algorithm, see the documentation related.


        Return Value:   
        ~~~~~~~~~~~~~
        Returns a (double *)NULL if uninitialized or if initialization 
        failed; returns a (double *) to the 20 kHz data vector if reset was 
        OK and/or is in "run" (MNRU_CONTINUE) operation.


        History:
        ~~~~~~~~
        05.Feb.1992     1.10 Release of the modular version.
                             <[e-mail]>
        05.Feb.1992     2.00 Updated according to the new P.81: 
                             - no up/downsampling
			     - input signal DC removal filter
			     - output low-pass filter (instead of band-pass)
			     <[e-mail]>

  ==========================================================================
*/
/* original RPELTP: #define ALPHA 0.999 */
#define ALPHA 0.985
#define DNULL (double *)0
#ifdef STL92_RNG
#define NOISE_GAIN 0.541
#else
/* NOISE_GAIN = 0.3795 for best match with the average SNR */
/*              0.3787 for best best match with the total SNR */
/*              0.3793 for a "balanced" middle-way between both SNRs */
#define NOISE_GAIN 0.3793
#endif

double         *MNRU_process(operation, s, input, output, n, seed, mode, Q)
  char            operation, mode;
  MNRU_state     *s;
  float          *input, *output;
  long            n, seed;
  double          Q;
{
  long            count, i;
  double          noise, tmp;
  register double inp_smp, out_tmp, out_flt;


  /*
   *    ..... RESET PORTION .....
   */

  /* Check if is START of operation: reset state and allocate memory buffer */
  if (operation == MNRU_START)
  {

    /* Reset clip counter */
    s->clip = 0;

    /* Allocate memory for sample's buffer */
   if ((s->vet = (double *) calloc(n, sizeof(double))) == DNULL)
      return ((double *) DNULL);

    /* Seed for random number generation */
    s->seed = seed;

    /* Gain for signal path */
    if (mode == MOD_NOISE)
      s->signal_gain = 1.000;
    else if (mode == SIGNAL_ONLY)
      s->signal_gain = 1.000;
    else			/* (mode == NOISE_ONLY) */
      s->signal_gain = 0.000;

    /* Gain for noise path */
    if (mode == MOD_NOISE || mode == NOISE_ONLY)
      s->noise_gain = NOISE_GAIN * pow(10.0, (-0.05 * Q));
    else			/* (mode == SIGNAL_ONLY) */
      s->noise_gain = 0;

    /* Flag for random sequence initialization */
    s->rnd_mode = RANDOM_RESET;

    /* Initialization of the output low-pass filter */
    /* Cleanup memory */
    memset(s->DLY, '\0', sizeof(s->DLY));

#ifdef NBMNRU_MASK_ONLY
    /* Load numerator coefficients */
    s->A[0][0]= 0.758717518025; s->A[0][1]= 1.50771485802; s->A[0][2]= 0.758717518025;
    s->A[1][0]= 0.758717518025; s->A[1][1]= 1.46756552150; s->A[1][2]= 0.758717518025;

    /* Load denominator coefficients */
    s->B[0][0]= 1.16833932919; s->B[0][1]= 0.400250061172;
    s->B[1][0]= 1.66492368687; s->B[1][1]= 0.850653444434;
#else
    /* Load numerator coefficients */
    s->A[0][0]= 0.775841885724; s->A[0][1]= 1.54552788762; s->A[0][2]= 0.775841885724;
    s->A[1][0]= 0.775841885724; s->A[1][1]= 1.51915539326; s->A[1][2]= 0.775841885724;

    /* Load denominator coefficients */
    s->B[0][0]= 1.23307153957; s->B[0][1]= 0.430807372835;
    s->B[1][0]= 1.71128410940; s->B[1][1]= 0.859087959597;
#endif

    /* Initialization of the input DC-removal filter */
    s->last_xk = s->last_yk = 0;
  }

  /*
   *    ..... REAL MNRU WORK .....
   */

  /* Initialize memory */
  memset(s->vet, '\0', n * sizeof(double));

  for (count = 0; count < n; count++)
  {
    /* Copy sample to local variable */
    inp_smp = *input++;

#ifndef NO_DC_REMOVAL
    /* Remove DC from input sample: H(z)= (1-Z-1)/(1-a.Z-1) */
    tmp = inp_smp - s->last_xk;
    tmp += ALPHA * s->last_yk;

    /* Update for next time */
    s->last_xk = inp_smp;
    s->last_yk = tmp;

    /* Overwrite DC-removed version of the input signal */
    inp_smp = tmp;
#endif
    
    /* Random number generation */
    if (mode == SIGNAL_ONLY)
      noise = 0;
    else
    {
      noise = (double) random_MNRU(&s->rnd_mode, &s->rnd_state, s->seed);
      noise *= s->noise_gain * inp_smp;	/* noise modulated by input sample */
      if (noise>1.00 || noise <-1.00) s->clip++; /* clip counter */
    }

    /* Addition of signal and modulated noise */
    out_tmp  = noise + inp_smp * s->signal_gain;

#ifdef NO_OUT_FILTER
    out_flt = out_tmp;
#else
    /* Filter output sample by each stage of the low-pass IIR filter */
    for (i=0; i<MNRU_STAGE_OUT_FLT; i++)
    {
        out_flt = out_tmp * s->A[i][0] + s->DLY[i][1];
        s->DLY[i][1] = out_tmp * s->A[i][1] - out_flt * s->B[i][0] + 
	               s->DLY[i][0];
        s->DLY[i][0] = out_tmp * s->A[i][2] - out_flt * s->B[i][1];

        out_tmp = out_flt;  /* output becomes input for next stage */
    }
#endif

    /* Copy noise-modulated speech sample to output vector */
    *output++  = out_flt;
  }

  /* Check if is end of operation THEN release memory buffer */
  if (operation == MNRU_STOP)
  {
    free(s->rnd_state.gauss);
    free(s->vet);
    s->vet = (double *) DNULL;
  }

  /* Return address of vet: if NULL, nothing is allocated */
  return ((double *) s->vet);
}
#undef NOISE_GAIN
#undef DNULL 
#undef ALPHA

/*  .................... End of MNRU_process() ....................... */