tomwalters@0: /***************************************************************************
tomwalters@0:  loudness.c     Loudness measure from spiral sector-weights contour.
tomwalters@0:  ----------
tomwalters@0: 			Michael Allerhand, 1992.
tomwalters@0: 
tomwalters@0: 	Read m binary shorts (default m=SECTORS) for each specified frame.
tomwalters@0: 	Compute the loudness as the mean value of the sector-weights.
tomwalters@0: 
tomwalters@0: ****************************************************************************/
tomwalters@0: 
tomwalters@0: #include <stdio.h>
tomwalters@0: #include <math.h>
tomwalters@0: #include "options.h"
tomwalters@0: #include "units.h"
tomwalters@0: #include "strmatch.h"
tomwalters@0: 
tomwalters@0: char applic[] = "Loudness measure from spiral sector-weights contour.\n           (Input and output in binary shorts)." ;
tomwalters@0: 
tomwalters@0: static char *helpstr, *debugstr, *fstr, *sstr, *mstr ;
tomwalters@0: 
tomwalters@0: static Options option[] = {
tomwalters@0:     {   "help"      ,   "off"       ,  &helpstr     ,   "help"                        , DEBUG   },
tomwalters@0:     {   "debug"     ,   "off"       ,  &debugstr    ,   "debugging switch"            , DEBUG   },
tomwalters@0:     {   "frames"    ,   "1-max"     ,  &fstr        ,   "select frames inclusively"   , VAL     },
tomwalters@0:     {   "scale"     ,   "1.848"     ,  &sstr        ,   "scale factor"                , SVAL    },
tomwalters@0:     {   "nsectors"  ,   "64"        ,  &mstr        ,   "number of sectors"           , SVAL    },
tomwalters@0:    ( char * ) 0 } ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: main (argc, argv)
tomwalters@0: int    argc;
tomwalters@0: char **argv;
tomwalters@0: {
tomwalters@0:     FILE   *fp;
tomwalters@0:     int     m;                  /* Number of sectors around spiral.     */
tomwalters@0:     short  *sector;             /* Sector-weights contour.              */
tomwalters@0:     float   scalar;             /* Scale factor for loudness.           */
tomwalters@0:     int     i, a, b;
tomwalters@0:     char   *val1, *val2 ;
tomwalters@0: 
tomwalters@0:     fp = openopts( option,argc,argv ) ;
tomwalters@0:     if ( !isoff( helpstr ) )
tomwalters@0: 	helpopts( helpstr, argv[0], applic, option ) ;
tomwalters@0: 
tomwalters@0:     m = atoi( mstr ) ;
tomwalters@0:     scalar = atof( sstr ) ;
tomwalters@0: 
tomwalters@0:     /* parse bounds on number of frames */
tomwalters@0: 
tomwalters@0:     if ( getvals( fstr, &val1, &val2, "-" ) == BADVAL ) {
tomwalters@0: 	fprintf(stderr,"loudness: bad frame selector [%s]\n", fstr ) ;
tomwalters@0: 	exit( 1 ) ;
tomwalters@0:     }
tomwalters@0:     a = atoi( val1 ) ;
tomwalters@0:     if ( isempty( val2 ) )    b = a ;
tomwalters@0:     else if ( ismax( val2 ) ) b = 0 ;
tomwalters@0:     else                      b = atoi( val2 ) ;
tomwalters@0: 
tomwalters@0:     if (b<a && b>0) fprintf(stderr,"loudness: warning, bad frame specifiers\n");
tomwalters@0: 
tomwalters@0:     if ( ison( debugstr ) ) {
tomwalters@0: 	printf("a=%d b=%d m=%d scalar=%.3f\n",
tomwalters@0: 		a, b, m, scalar );
tomwalters@0: 	exit(0);
tomwalters@0:     }
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /* Allocate space for sector-weights contour */
tomwalters@0: 
tomwalters@0:     if (( sector = (short *)malloc(m*sizeof(short))) == NULL)
tomwalters@0: 	    fprintf(stderr,"loudness: malloc out of space\n");
tomwalters@0: 
tomwalters@0:     /* Seek past a-1 frames, then read the next b-a+1 frames, or all    */
tomwalters@0:     /* remaining frames if b==0.                                        */
tomwalters@0: 
tomwalters@0:     if (a>0) {
tomwalters@0: 
tomwalters@0: 	for (i=1 ; i<a  && fread(sector,sizeof(short),m,fp) ; i++)
tomwalters@0: 	    ;
tomwalters@0: 	if (b > 0)
tomwalters@0: 	    for (  ; i<=b && fread(sector,sizeof(short),m,fp) ; i++)
tomwalters@0: 		loudness(sector,m,scalar);
tomwalters@0: 	else
tomwalters@0: 	    while ( fread(sector,sizeof(short),m,fp) )
tomwalters@0: 		loudness(sector,m,scalar);
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     fclose(fp);
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: loudness(sector,m,scalar)
tomwalters@0: short  *sector;         /* Sector-weights contour.                  */
tomwalters@0: int     m;              /* Number of sectors around spiral.         */
tomwalters@0: float   scalar;         /* Scale factor for loudness measure.       */
tomwalters@0: {
tomwalters@0:     short   scale;      /* Mean weights.                            */
tomwalters@0: 
tomwalters@0:     scale = (short)( normalize_scale(sector,m) * scalar );
tomwalters@0:     fwrite(&scale,sizeof(short),1,stdout);
tomwalters@0: 
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /****************************************************************************
tomwalters@0:     Normalize Scale.
tomwalters@0: ****************************************************************************/
tomwalters@0: 
tomwalters@0: /*
tomwalters@0: Normalize for unit variance and zero mean, and then scale for MEAN and STDDEV.
tomwalters@0: In general the stddev is scaled down to STDDEV, but if the stddev is already
tomwalters@0: less than the target STDDEV, than it is left alone.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: #define MEAN    100
tomwalters@0: #define STDDEV  32       /* Resulting variance will be this squared */
tomwalters@0: 
tomwalters@0: normalize_scale(V,m)
tomwalters@0: short *V;
tomwalters@0: int    m;
tomwalters@0: {
tomwalters@0:     int   i;
tomwalters@0:     float mean, var;
tomwalters@0:     float sum, sumsq, stddev, norm;
tomwalters@0: 
tomwalters@0:     sum=0; sumsq=0;
tomwalters@0:     for (i=0 ; i<m ; i++) {
tomwalters@0: 	sum   += V[i];
tomwalters@0: 	sumsq += V[i] * V[i];
tomwalters@0:     }
tomwalters@0:     mean = sum/m;
tomwalters@0:     var  = sumsq/m - mean*mean;
tomwalters@0: 
tomwalters@0:     if ((stddev = sqrt(var)) >= STDDEV)
tomwalters@0: 	for (i=0 ; i<m ; i++) {
tomwalters@0: 	    norm = ( (float)V[i] - mean ) / stddev;
tomwalters@0: 	    V[i] = norm * STDDEV + MEAN;
tomwalters@0: 	}
tomwalters@0:     else
tomwalters@0: 	for (i=0 ;  i<m ; i++) {
tomwalters@0: 	    norm = (float)V[i]-mean;
tomwalters@0: 	    V[i] = norm + MEAN;
tomwalters@0: 	}
tomwalters@0: 
tomwalters@0:     return (int)mean;
tomwalters@0: }
tomwalters@0: