mnru.c

Reference Implementation of G.711 standard and other voice codecs

/*                                                        31.JULY.1995 v.2.00
=============================================================================

                          U    U   GGG    SSSS  TTTTT
                          U    U  G       S       T
                          U    U  G  GG   SSSS    T
                          U    U  G   G       S   T
                           UUU     GG     SSS     T

                   ========================================
                    ITU-T - USER'S GROUP ON SOFTWARE TOOLS
                   ========================================


       =============================================================
       COPYRIGHT NOTE: This source code, and all of its derivations,
       is subject to the "ITU-T General Public License". Please have
       it  read  in    the  distribution  disk,   or  in  the  ITU-T
       Recommendation G.191 on "SOFTWARE TOOLS FOR SPEECH AND  AUDIO 
       CODING STANDARDS".
       =============================================================


MODULE:         NEW VERSION OF THE MNRU.C, MODULATED NOISE REFERENCE 
                UNIT'S MODULE, MNRU (ACCORDING P.81, 1995).

ORIGINAL BY:    Simao Ferraz de Campos Neto <[e-mail]>,
                based on ITU-T STL92 MNRU, which was based on 
                CSELT's MNRU.FOR algorithm.

FUNCTIONS:

MNRU_process: ......	Processes the `input' buffer of `n' samples (in
                        float format), adding to it modulated noise at a
                        `Q' dB SNR level, if `mode' is modulated noise, and
                        saving to the float buffer `output'. Otherwise, if
                        mode is noise-only, then saves only noise to
                        `output' buffer, or the filtered signal, if
                        signal-only `mode'. Depending on the `operation'
                        chosen, state variables in `s' are reset, as well
                        as memory allocated (start), kept as are (continue)
                        or memory is released (stop). Written for narrow-band
                        model (input data is sampled at 8 kHz). Its prototype
                        is in mnru.h.

random_MNRU: .......... Generates gaussian-like noise samples for use by the
                        MNRU_process function. Depends on a seed when `*mode'
                        is 1 (RANDOM_RESET), causing the initialization of
                        the generator. Then `*mode' is changed to 0
                        (RANDOM_RUN), and on a state variable that should
                        not be changed by the user. Its prototype is found
                        in mnru.h.

HISTORY:

  25.Set.91  v1.0F      Fortran version released to UGST by CSELT/Italy.
  23.Jan.92  v1.0C      C version, bit-exact with Fortran impl. in VAX
                        <[e-mail]>
  27.Jan.92  v1.1       Modular version with portable RNG
                        <[e-mail]>
  18.May.92  v1.2       Removal of of addition of +-0.5, at the exit of 
                        MNRU_process, needed to make it work with data in the 
                        normalized range. <[e-mail]>
  31.Jul.95  v2.0       Redifinition of the module due to the revision of P.81:
                        no more 1:5 up/down-sampling, inclusion of a DC filter
                        and a low-pass (instead of a band-pass) output filter.
                        To increase speed, a new random number generator
                        has been included. Works for both narrow-band and 
                        wideband speech.
=============================================================================
*/

/* Definitions for the algorithm itself */
#include "mnru.h"

/* General includes */
#include <math.h>
#include <stdlib.h> /* for calloc(), free() */
#include <string.h> /* for memset() */

#ifndef STL92_RNG /* Uses the new Random Number Generator */
#define random_MNRU new_random_MNRU

/* Local function prototypes */
float new_random_MNRU ARGS((char *mode, new_RANDOM_state *r, long seed));
float ran_vax ARGS((void));
unsigned long ran16_32c ARGS((void));

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

	new_random_MNRU (char *mode, RANDOM_state *r, long seed)
        ~~~~~~~~~~~~~~~

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

        Random number generator based on a gaussian random 
	number table accessed by a uniform random number generator.

	The gaussian random sample table is generated at start-up time
        by a Monte-Carlo algorithm using a linear congruential generator
        (LCG). During run-time the algorithm accesses randomly this table
	using indeces generated by another LCG.

	To (re)initialize the sequence, use mode=RANDOM_RESET (the routine
	will change mode to RANDOM_RUN).

	Functions used:
	~~~~~~~~~~~~~~~
        ran_vax():    generate uniformly distributed samples in the range
                      0..1. Return a float number.
	ran16_32C():  generate uniformly distributed samples in the range
	              0..65535 (2^16-1)

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

        History:
        ~~~~~~~~
        01.Jul.95  1.0	Created based on the random number generator 
                        implemented by Aachen University and used by  
                        the ITU-T 8kbit/s speech codec host laboratory  
                        (in hardware). <[e-mail]>

=============================================================================
*/
#define S1    -8.0      /* s1 = my - 4 * sigma (=-8.0 for gaussian noise) */
#define S2     8.0      /* s2 = my + 4 * sigma (= 8.0 for gaussian noise) */
#define DIF    16.0     /* s2 - s1                                        */
#define MO     8.0      /* mo = 2 * (sigma)^2  (= 8.0 for gaussian noise) */
#define BIT15  32767.0
#define TABLE_SIZE 8192 /* 2^13 */
#define ITER_NO 8
#define FACTOR 8  /* = 65536(max.no returned by ran16_32c) div.by TABLE_SIZE */
float           new_random_MNRU(mode, r, seed)
  char           *mode;
  RANDOM_state   *r;
  long           seed;
{

  long    i;
  double  z1;            /* white random number   -8...8          */
                         /* weighted with a gaussian distribution */
  double  z2;            /* white random number    0...1          */
  double  phi;           /* gauss curve                           */

  extern float ran_vax();
  extern unsigned long ran16_32c();

  long index;

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

    /* Allocate memory for gaussian table */
    r->gauss = (float *)calloc(TABLE_SIZE, sizeof(float));

    /* Generate gaussian random number table */
    for(i=0L; i<TABLE_SIZE; i++)
    {
      /* Interact until find gaussian sample */
      do
      {
	z1 = S1 + DIF*(double)ran_vax();
	phi= exp( -(z1)*(z1)/MO);
	z2 = (double)ran_vax();
      } while(z2 > phi);

      /* Save gaussian-distributed sample in table */
      r->gauss[i] = (float)z1;
    }
  }

  /* ***  REAL GENERATOR (after initialization) ***/
  for (z1=0, i=0;i<ITER_NO;i++) 
  {
    index = ran16_32c()/FACTOR;
    z1 += r->gauss[index];
  }
  z1 /= 2; /* provisional */

  /* Return gaussian sample */
  return ((float)z1);
}
#undef TABLE_SIZE
#undef BIT15
#undef MO
#undef DIF
#undef S2
#undef S1
/*  .................... End of new_random_MNRU() ....................... */


/*
  ===========================================================================
  float ran_vax(void);
  ~~~~~~~~~~~~~

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

  Function that simulates the VAX Fortran function RAN(x), that returns 
  a number uniformly distributed between 0.0 and 1.0. This implementation 
  is based on Aachen University's randm() function of the narrow-band 
  MNRU table-generation program montrand.c by CA (6.3.90).

  Parameters: none.
  ~~~~~~~~~~~

  Return value:
  ~~~~~~~~~~~~~
  An float number uniformly distributed in the range 0.0 and 1.0.

  Author:
  ~~~~~~~
  Simao Ferraz de Campos Neto          
  Comsat Laboratories                  Tel:    [phone]
  22300 Comsat Drive                   Fax:    [phone]
  Clarksburg MD 20871 - USA            E-mail: [e-mail]
  
  History:
  ~~~~~~~~
  01.Jul.95  v1.00  Created, adapted from montrand.c

  ===========================================================================
*/
#define CONST         69069
#define INIT          314159265L
#define BIT32         4294967296.0
float ran_vax()
{
  static unsigned long  seed, buffer;
  static float          ran;
  static short          firsttime=0;

  if(firsttime == 0)
  {
  	firsttime = 1;
  	seed      = INIT;
  }

  seed    = seed * CONST + 1;         /* includes the mod 2**32 operation  */
  buffer  = seed & 0xFFFFFF00;        /* mask the first 24 bit             */
  ran     = (float)buffer / BIT32;    /* and divide by 2**32 to get random */

  return(ran);
}
/*  ......................... End of ran_vax() ............................ */


/*
  ===========================================================================
  unsigned long ran16_32c(void);
  ~~~~~~~~~~~~~~~~~~~~~~~

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

  Function that simulates the DSP32C function RAN24(), modified to return 
  a number between 0 and 2^16-1. This is based on Aachen University's 
  randm() of the narrow-band MNRU program mnrusim.c by PB (08.04.1991).

  Parameters: none.
  ~~~~~~~~~~~

  Return value:
  ~~~~~~~~~~~~~
  An unsigned long number in the range 0 and 2^16-1.

  Author:
  ~~~~~~~
  Simao Ferraz de Campos Neto          
  Comsat Laboratories                  Tel:    [phone]
  22300 Comsat Drive                   Fax:    [phone]
  Clarksburg MD 20871 - USA            E-mail: [e-mail]
  
  History:
  ~~~~~~~~
  01.Jul.95  v1.00  Created, adapted from mnrusim.c

  ===========================================================================
*/
#define BIT24	16777216.0
#define BIT8    256.0
unsigned long ran16_32c()
{
  static float   seed = 12345.0;
  double         buffer1, buffer2;
  long           seedl;
  unsigned long  result;

  buffer1 = ((253.0 * seed) + 1.0);
  buffer2 = (buffer1/BIT24) ;
  seedl   = ((long)buffer2) & 0x00FFFFFFL;
  seed    = buffer1 = buffer1 - (float)seedl * BIT24;
  result  = buffer1 / BIT8;

  return result;
}
#undef BIT8
#undef BIT24
/*  .................... End of ran16_32c() ....................... */

#else /* Use the original MNRU noise generator */

#define random_MNRU ori_random_MNRU

/* Local function prototypes */
float ori_random_MNRU ARGS((char *mode, ori_RANDOM_state *r, long seed));

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

	ori_random_MNRU (char *mode, RANDOM_state *r, long seed)
        ~~~~~~~~~~~~~~~

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

        Random number generator based on Donald Knuth subtractive
        method [1] and in Press [2] C implementation (ran3()).
        The core of the routine generates a uniform deviate between -0.5 
        and 0.5. By calling this core several times (after the initialization
        phase), the cumulation of the uniform noise samples leads to a 
        gaussian noise (central limit theorem), thus generating a noise 
        suitable for the MNRU module. Tests showed that the SNR of the 
        signal processed by the MNRU is very close to the desired value 
        of Q when the number of cumulations (defined by ITER_NO) is 47.

	To (re)initialize the sequence, use mode=RANDOM_RESET (the routine
	will change mode to RANDOM_RUN).

        According to [1], any large MBIG and any smaller (but still
        large) MSEED can be used; this keeps the values chose by [2],
        pg.212. The dimension of 56 to ma MUST be kept as is, as well
        as inextp=31 [1].


        [1] Knuth, D.; "Seminumerical Algorithms", 2nd. ed., Vol.2 of
           "The Art of Computer Programming"; Addison-Wesley, Mass;
            1981, Parts 3.2-3.3.