stats.c

speech signal process tools

/*
  stats - compute statistics on SPS files				     
  Joseph T. Buck							    
 									   
  stats considers an SPS file to be a series of vectors. It computes mean, 
  standard deviation, maximum, and minimum values on an element-by-element 
  basis. The number of vectors used may be controlled with the -o switch.  
  The -n switch suppresses the header.					   
  									   
  From SDS version 1.3 of 10/28/85					   

  Joe Buck started converting it to SPS; Alan Parker finsished it.
  Dave Burton made ESPS modifications.

  This material contains proprietary software of Entropic Speech, Inc.   
  Any reproduction, distribution, or publication without the the prior	   
  written permission of Entropic Speech, Inc. is strictly prohibited.
  Any public distribution of copies of this work authorized in writing by
  Entropic Speech, Inc. must bear the notice			
 								
      "Copyright (c) 1987 Entropic Speech, Inc.; All rights reserved"

*/
#ifdef SCCS
	static char *sccs_id = "@(#)stats.c	1.7	6/25/87 ESI";
#endif
#include <values.h> /* maybe MASSCOMP specific */

#define MAXFLOAT ((float)3.4e+37)  /* fix around Masscomp 68010 compiler
				      bug */

#include <stdio.h>
#include <sps/sps.h>
#include <math.h>
#define MAXL 100
#define SYNTAX USAGE ("stats [-r rec_range] [-e el_range] [-m magnitude] [-n] file")
#define MAGNITUDE(x, m) (((x-m) > (m-x)) ? (x-m) : (m-x))

void lrange_switch(), range_switch();
void exit();
short get_rec_len();
main (argc, argv)
int   argc;
char  **argv;
{
    extern  optind;
    extern char *optarg;
    short l_rec;
    int     c, hflag = 1, mflag = 0;
    long s_rec, e_rec, i_rec, i;
    long    tag;
    char   *prange = NULL, *erange = NULL, *filename = NULL;
    double  data [MAXL], mins[MAXL], maxs[MAXL], avgs[MAXL], devs[MAXL], 
	    mag[MAXL];
    int     minl [MAXL], maxl[MAXL];
    int     s_ele, e_ele;
    double  scale, bias, sqrt ();
    FILE   *istrm = stdin;
    struct header *h;

    while ((c = getopt (argc, argv, "nmr:e:")) != EOF) {
	switch (c) {
	    case 'n':
		hflag = 0;
		break;
	    case 'm':
		mflag++;
		break;
	    case 'r': 
		prange = optarg;
		break;
	    case 'e':
	        erange = optarg;
		break;
	    default:
		SYNTAX;
	}
    }
    if (argc == optind + 1)
	filename = argv[optind];
    
    read_params((char *)NULL,SC_CHECK_FILE,filename);
    if (filename == NULL) {
	if(symtype("filename") == ST_UNDEF) {
	  fprintf (stderr, "stats: no input file.\n");
	  exit(1);
    	}
	filename = getsym_s("filename");
    }

    /*
     * Check for stdin
    */
    if(strcmp(filename, "-") == 0)
	filename = "<stdin>";
    else
	TRYOPEN(argv[0], filename, "r", istrm);

    if (!(h = read_header (istrm)))
	NOTSPS ("stats", filename);
    l_rec = get_rec_len(h);
    if (l_rec > MAXL) {
	fprintf (stderr, "stats: record length (%d) too long\n", l_rec);
	exit (1);
    }
    s_rec = 1;
    e_rec = MAXLONG;
    if(symtype("start") != ST_UNDEF) 
	s_rec = getsym_i("start");
    if(symtype("nan") != ST_UNDEF) 
	e_rec = s_rec + getsym_i("nan") - 1;

    lrange_switch (prange, &s_rec, &e_rec, 1);

    (void) putsym_s("filename",filename);
    (void) putsym_s("prog","stats");
    (void) putsym_i("start",(int)s_rec);

    s_ele = 1;
    e_ele = l_rec;
    range_switch (erange, &s_ele, &e_ele, 1);
    if (e_ele > l_rec) {
	fprintf (stderr, "stats: only %d elements per record\n", l_rec);
	exit (1);
    }
/* Initialize statistics */
    for (i = s_ele - 1; i < e_ele; i++) {
	avgs[i] = devs[i] = 0.0;
	mins[i] = MAXFLOAT;
	maxs[i] = -MAXFLOAT;
	mag[i] = 0.;
    }
    skiprec (istrm, s_rec - 1, size_rec (h));
    for (i_rec = s_rec; 
	 i_rec <= e_rec && get_gen_recd(data,&tag,h,istrm) != EOF; 
	 i_rec++) {

	for (i = s_ele - 1; i < e_ele; i++) {
	    if (data[i] < mins[i]) {
		mins[i] = data[i];
		minl[i] = i_rec;
	    }
	    if (data[i] > maxs[i]) {
		maxs[i] = data[i];
		maxl[i] = i_rec;
	    }
	    avgs[i] += data[i];
	    devs[i] += data[i] * data[i];
	}
    }

    e_rec = i_rec - 1; 
    (void) putsym_i("nan",(int)(e_rec-s_rec+1));
    if (hflag) printf("stats: filename: %s, start: %ld, end: %ld\n",
	filename, s_rec, e_rec);

    scale = 1.0 / (double) (e_rec-s_rec+1);
/* bias is the factor that converts the biased estimate of the variance to
   the unbiased estimate: (e_rec-s_rec+1)/((e_rec-s_rec)-1).
*/
    bias = 0.0;
    if ((e_rec-s_rec+1) > 1)
	bias = (double) (e_rec-s_rec+1) / (double) (e_rec-s_rec);
    for (i = s_ele - 1; i < e_ele; i++) {
	avgs[i] *= scale;
/*
 * for some cases this expression yields a negative number very close
 * zero that should be treated as zero.
 */
	if ((bias * (devs[i] * scale - avgs[i] * avgs[i])) <= 0)
	  devs[i] = 0;
	 else
	  devs[i] = sqrt (bias * (devs[i] * scale - avgs[i] * avgs[i]));
    }

/*
 * Now compute average adjusted magnitude, if required
*/
    if(mflag > 0){
    /* first rewind file and reread header*/
    rewind(istrm);
    h = read_header (istrm);
    skiprec (istrm, s_rec - 1, size_rec (h));
	    for (i_rec = s_rec; 
		 i_rec <= e_rec && get_gen_recd(data,&tag,h,istrm) != EOF; 
		 i_rec++) {
		for (i = s_ele - 1; i < e_ele; i++) {
		    mag[i] += MAGNITUDE(data[i], avgs[i]);
		}
	    }
    /* divide by number of records*/
    for(i=s_ele-1;i<e_ele;i++)
	mag[i] *= scale;
    }


    if (hflag)
	printf ("Element   Minimum    at     Maximum    at      Mean      Std. Dev     Magnitude\n");

        for (i = s_ele - 1; i < e_ele; i++) {
	printf ("%5d%12g%6d%12g%6d%12g%12g%12g\n", i + 1, mins[i], minl[i],
		maxs[i], maxl[i], avgs[i], devs[i], mag[i]);
	}
    
    return 0;
}