snr.c

Reference Implementation of G.711 standard and other voice codecs

/* *********************************************************************** */
/* FORTRAN77 main program SNR */
/* Task:      Compute signal to noise ratio, SNR, between two signal files. */
/*            Global, segmental and minimum SNRs are computed */
/* ----------------------------------------------------------------------- */
/* Input:     Reference signal file */
/*            Test signal file */
/* Output:    SNR values to terminal */
/* ----------------------------------------------------------------------- */
/* File handling:      through routines FOPEN, FREAD */
/*                     files on units 101 and 102 */
/* User communication: FORTRAN terminal I/O */
/* System interaction: none */
/* ----------------------------------------------------------------------- */
/* Subroutines called:                                       from: */
/*   FOPEN, FREAD                                            FILSUB */
/*   GLBSNR,SEGSNR,MINSNR                                   this file */
/* ----------------------------------------------------------------------- */
/* Version: 1                Date: 14.10.91        By: Finn Tore Johansen */
/* *********************************************************************** */
/* snr.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"


/* Table of constant values */

static long     c__256 = 256;


int             main(argc, argv)
  int             argc;
  char           *argv[];
{
  /* Initialized data */

  static long     nlow = 0;
  static long     nseg = 0;
  static double   mn[7] = {200., 200., 200., 200., 200., 200., 200.};
  static long     imn[7] = {0, 0, 0, 0, 0, 0, 0};

  /* System generated locals */
  long            i_1;

  /* Builtin functions */
  long            s_wsfe(), do_fio(), e_wsfe(), s_rsfe(), e_rsfe(), s_wsle(), e_wsle(),
                  do_lio(), f_clos();

  /* Local variables */
  static long     i;
  static double   x[256], y[256];
  static char     infil1[80], infil2[80];
  extern int      glbsnr_();
  extern void     segsnr_();
  extern void     minsnr_();
  static long     ix1[256], ix2[256];
  static double   glb;
  static long     eof;
  static double   seg;
  static long     ifr;
  static short    buf1[256], buf2[256];
  long            samples;
  FILE           *f1, *f2;


  /* User communication */
  GET_PAR_S(1, "_Reference signal file name: ", infil1);
  GET_PAR_S(2, "_Test signal file name: ", infil2);

  /* Initialize files */
  if ((f1 = fopen(infil1, RB)) == NULL)
    KILL(infil1, 2);
  if ((f2 = fopen(infil2, RB)) == NULL)
    KILL(infil2, 3);

  ifr = 0;

  while ((i = fread(buf1, sizeof(short), (size_t) 256, f1)) > 0 &&
	 (i_1 = fread(buf2, sizeof(short), (size_t) 256, f2)) > 0)
  {
    /* Get the number of samples read in this data buffer */
    samples = (i <= i_1) ? i : i_1;

    /* Copy input signals to double buffer */
    for (i = 0; i < samples; ++i)
    {
      x[i] = (double) buf1[i];
      y[i] = (double) buf2[i] - (double) buf1[i];
    }
    /* Compute different SNRs, global, segmental (256 sample segments) and */
    /* minimum for block lengths 256,128,64,32,16,8 and 4 */
    ++ifr;
    glbsnr_(x, y, &samples, &glb);
    segsnr_(x, y, &samples, &seg, &nlow, &nseg);
    for (i = 0; i <= 6; ++i)
    {
      minsnr_(x, y, &samples, &i, &ifr, &mn[i], &imn[i]);
    }
  }

  /* Write results to terminal */
  if (glb == -1.)
  {
    printf("Zero input signal\n");
  }
  else if (glb == -2.)
  {
    printf("Files are identical\n");
  }
  else
  {
    printf("%s%s", "   SEG256  GLOB    MIN256  MIN128",
	   "  MIN64   MIN32   MIN16   MIN8    MIN4\n");
    printf("%8.2lf", seg);
    printf("%8.2lf", glb);
    for (i = 0; i <= 6; ++i)
    {
      printf("%8.2lf", mn[i]);
    }
    printf("\n   Min. segment:");
    for (i = 0; i < 7; ++i)
    {
      printf("%8ld", imn[i]);
    }
  }

  printf("\n   Active segments :%12ld", nseg);
  printf("\n   Total segments  :%12ld\n", nseg + nlow);

  /* Finalizations */
  fclose(f1);
  fclose(f2);
#ifndef VMS
  return (0);
#endif
}
/* ....................... End of main() ................................ */


/* Compute global SNR for file */
/* Return 0 if SNR, -1 if silence (no signal) and -2 if snr->oo */
int             glbsnr_(x, y, nfr, glb)
  double         *x, *y;
  long           *nfr;
  double         *glb;
{
  /* Initialized data */

  static double   xp = 0.;
  static double   yp = 0.;

  /* System generated locals */
  long            i_1;
  double          d_1;

  /* Local variables */
  static long     i;

  i_1 = *nfr - 1;
  for (i = 0; i <= i_1; ++i)
  {
    /* Computing 2nd power */
    d_1 = x[i];
    xp += d_1 * d_1 / *nfr;
    /* Computing 2nd power */
    d_1 = y[i];
    yp += d_1 * d_1 / *nfr;
    /* L10: */
  }

  /* Check limit values and calculate SNR */
  if (xp <= 0.)
  {
    /* No signal - silence */
    *glb = (double) -1.;
    return -1;
  }
  else if (yp <= 0.)
  {
    /* Signals are identical: snr -> oo */
    *glb = (double) -2.;
    return -2;
  }
  else
  {
    /* OK, they are "regular" signals */
    d_1 = xp / yp;
    *glb = log10(d_1) * 10;
  }

  /* Return value */
  return 0;
}
/* .......................... End of glbsnr() ........................... */


/* Compute average segmental SNR for segments above threshold */
void            segsnr_(x, y, nfr, seg, nlow, nseg)
  double         *x, *y;
  long           *nfr;
  double         *seg;
  long           *nlow, *nseg;
{
  /* Initialized data */
  static double   sum = 0.;

  /* System generated locals */
  long            i_1;
  double          d_1;

  /* Local variables */
  static long     i;
  static double   xp, yp, maxsnr, xth;

  xth = (double) 1e4;
  maxsnr = (double) 200.;
  xp = (double) 0.;
  yp = (double) 0.;
  i_1 = *nfr - 1;
  for (i = 0; i <= i_1; ++i)
  {
    /* Computing 2nd powers */
    d_1 = x[i];
    xp += d_1 * d_1;
    d_1 = y[i];
    yp += d_1 * d_1;
  }
  xp /= *nfr;
  yp /= *nfr;
  if (xp <= xth)
  {
    ++(*nlow);
  }
  else
  {
    ++(*nseg);
    if (yp > 0.)
    {
      d_1 = xp / yp;
      sum += log10(d_1) * 10;
    }
    else
    {
      sum += maxsnr;
    }
    *seg = sum / *nseg;
  }
}
/* ........................... End of segsnr() ........................... */


/* Compute minimum block SNR for blocks with energy above threshold */
void            minsnr_(x, y, nfr, nsub, ifr, mn, imn)
  double         *x, *y;
  long           *nfr, *nsub, *ifr;
  double         *mn;
  long           *imn;
{
  /* System generated locals */
  long            i_1, i_2;
  double          d_1;

  /* Local variables */
  static long     segl, i, n, i0;
  static double   xp, yp, maxsnr, snr, xth;

  xth = (double) 1e4;
  maxsnr = (double) 200.;
  d_1 = pow((double) 2., (double) (*nsub));
  segl = *nfr / d_1;
  i_1 = d_1 - 1;
  for (n = 0; n <= i_1; ++n)
  {
    i0 = n * segl;
    xp = (double) 0.;
    yp = (double) 0.;
    i_2 = i0 + segl - 1;
    for (i = i0; i <= i_2; ++i)
    {
      /* Computing 2nd powers */
      d_1 = x[i];
      xp += d_1 * d_1;
      d_1 = y[i];
      yp += d_1 * d_1;
    }
    xp /= segl;
    yp /= segl;

    if (xp > xth)
    {
      if (yp <= 0.)
      {
	snr = maxsnr;
      }
      else
      {
	d_1 = xp / yp;
	snr = log10(d_1) * 10;
      }
      if (snr < *mn)
      {
	*mn = snr;
	*imn = *ifr;
      }
    }
  }
}
/* ........................... End of minsnr() ........................... */