tomwalters@0: /*************************************************************************
tomwalters@0:  sp_weights.c   Map auditory image onto spiral sector-weights contour.
tomwalters@0: --------------
tomwalters@0: 
tomwalters@0:     Read auditory image frames (gensai or genspl with output set) from a
tomwalters@0:     file or the stdin. Input must have a pipes header, (for frame dimensions).
tomwalters@0: 
tomwalters@0:     Note, you can use either gensai or genspl output, but remember that their
tomwalters@0:     default options are different, (particularly pwidth, nwidth, and dencf).
tomwalters@0: 
tomwalters@0:     For each frame, (of frameheight*framewidth binary shorts),
tomwalters@0:     write a sector-weights contour, (of m binary shorts), on the stdout.
tomwalters@0: 
tomwalters@0:     Arguments:
tomwalters@0: 	-m      The number of sectors (default SECTORS).
tomwalters@0: 	-fA-B   Process the A'th to B'th frame inclusive,
tomwalters@0: 		(or the A'th to the eof if b=='e'), (default -o1-e).
tomwalters@0: 	-fA     Process just the A'th frame.
tomwalters@0: 
tomwalters@0: *************************************************************************/
tomwalters@0: 
tomwalters@0: #include <stdio.h>
tomwalters@0: #include <math.h>
tomwalters@0: #include "header.h"
tomwalters@0: #include "options.h"
tomwalters@0: #include "units.h"
tomwalters@0: #include "strmatch.h"
tomwalters@0: 
tomwalters@0: char applic[] = "Map auditory image frames onto spiral sector-weights contours.\n           (Input and output in binary shorts)." ;
tomwalters@0: 
tomwalters@0: static char *helpstr, *debugstr, *infostr, *fstr, *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:     {   "info"      ,   "off"       ,  &infostr     ,   "print frame size information", DEBUG   },
tomwalters@0:     {   "frames"    ,   "1-max"     ,  &fstr        ,   "select frames inclusively"   , VAL     },
tomwalters@0:     {   "nsectors"  ,   "64"        ,  &mstr        ,   "number of sectors"           , SVAL    },
tomwalters@0:    ( char * ) 0 } ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: #define log2(x)     (log((double)x)/log(2.0))
tomwalters@0: #define round(x)    ((int)(x+0.5))
tomwalters@0: 
tomwalters@0: int     frameheight, framewidth ;   /* Parameters read from header */
tomwalters@0: int     frames;
tomwalters@0: double  orig_spd ;
tomwalters@0: 
tomwalters@0: main(argc, argv)
tomwalters@0: int   argc ;
tomwalters@0: char *argv[] ;
tomwalters@0: {
tomwalters@0:     FILE   *fp ;
tomwalters@0:     char   *header, *val1, *val2 ;
tomwalters@0:     char   *headeropt, *versionstr;
tomwalters@0:     int     i, a, b;
tomwalters@0:     int     m ;                 /* Number of sectors around spiral      */
tomwalters@0:     int     framebytes;         /* Number of bytes in each frame        */
tomwalters@0:     short  *frame;              /* Each frame                           */
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: 
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,"sp_weights warning: bad frame specifiers\n");
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /* Read header to find dimensions of auditory image */
tomwalters@0: 
tomwalters@0:     if ((header = ReadHeader(fp)) == (char *) 0 ) {
tomwalters@0: 	fprintf(stderr,"sp_weights: header not found\n");
tomwalters@0: 	exit(1);
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     if ( (versionstr = HeaderString( header, "Version" )) == (char *)0 ) {
tomwalters@0: 	fprintf(stderr,"sp_weights: model version-number not found in header\n");
tomwalters@0: 	exit(1);
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     if ( (headeropt = HeaderString( header, "frameheight" )) == (char *)0) {
tomwalters@0: 	fprintf(stderr,"sp_weights: frameheight not found in header\n");
tomwalters@0: 	exit(1);
tomwalters@0:     }
tomwalters@0:     else  frameheight = atoi( headeropt );
tomwalters@0: 
tomwalters@0:     if ( (headeropt = HeaderString( header, "framewidth" )) == (char *)0) {
tomwalters@0: 	fprintf(stderr,"sp_weights: framewidth not found in header\n");
tomwalters@0: 	exit(1);
tomwalters@0:     }
tomwalters@0:     else  framewidth = atoi( headeropt );
tomwalters@0: 
tomwalters@0:     if ( (headeropt = HeaderString( header, "frames" )) == (char *)0) {
tomwalters@0: 	fprintf(stderr,"sp_weights: frames not found in header\n");
tomwalters@0: 	exit(1);
tomwalters@0:     }
tomwalters@0:     else  frames = atoi( headeropt );
tomwalters@0: 
tomwalters@0:     if ( (headeropt = HeaderString( header, "orig_spd" )) == (char *)0) {
tomwalters@0: 	fprintf(stderr,"sp_weights: orig_spd not found in header\n");
tomwalters@0: 	exit(1);
tomwalters@0:     }
tomwalters@0:     else  orig_spd = atof( headeropt );
tomwalters@0: 
tomwalters@0:     if ( ison( debugstr ) ) {
tomwalters@0: 	printf("a=%d b=%d m=%d \n", a, b, m );
tomwalters@0: 	printf("frames=%d  frameheight=%d  framewidth=%d orig_spd=%.2f\n", frames,frameheight,framewidth,orig_spd ) ;
tomwalters@0: 	exit(0);
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     if ( ison( infostr ) || a==0 ) {
tomwalters@0: 	printf("sp_weights input: %s", versionstr ) ;
tomwalters@0: 	printf("                  frames=%d  frameheight=%d  framewidth=%d\n", frames,frameheight,framewidth ) ;
tomwalters@0: 	exit(0);
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     if (frames<=0) {
tomwalters@0: 	if (frames==0)  fprintf(stderr,"sp_weights: zero frames input\n");
tomwalters@0: 	if (frames<0)   fprintf(stderr,"sp_weights: garbled number of frames, (start set too big?)\n");
tomwalters@0: 	exit(1);
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     if (a>frames || b>frames)
tomwalters@0: 	fprintf(stderr,"sp_weights warning: bad frame specifiers\n");
tomwalters@0: 
tomwalters@0:     framebytes = frameheight * framewidth * sizeof(short) ;
tomwalters@0: 
tomwalters@0:     /* Allocate space for an auditory image frame */
tomwalters@0: 
tomwalters@0:     if (( frame = (short *)malloc( framebytes )) == NULL)
tomwalters@0: 	    fprintf(stderr,"sp_weights: 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:     for (i=1 ; i<a  && fread(frame,framebytes,1,fp) ; i++)
tomwalters@0: 	;
tomwalters@0:     if (b > 0)
tomwalters@0: 	for (  ; i<=b && fread(frame,framebytes,1,fp) ; i++)
tomwalters@0: 	    write_weights(frame,m);
tomwalters@0:     else
tomwalters@0: 	while ( fread(frame,framebytes,1,fp) )
tomwalters@0: 	    write_weights(frame,m);
tomwalters@0: 
tomwalters@0:     fclose(fp);
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /***************************************************************************
tomwalters@0:   Spoke sector weight routines  (originally as src/src8/spoke.c)
tomwalters@0: 
tomwalters@0:   Write m sector weights as binary shorts on the stdout.
tomwalters@0: ***************************************************************************/
tomwalters@0: 
tomwalters@0: write_weights(frame,m)
tomwalters@0: short *frame;
tomwalters@0: int    m;       /* number of sectors around the spiral plane */
tomwalters@0: {
tomwalters@0:     static  short *sector;
tomwalters@0:     static  int   *logtime;
tomwalters@0:     int     i, k, t, col;
tomwalters@0:     double  sum;
tomwalters@0:     static  int    first=1;     /* flag for first-call of this routine */
tomwalters@0: 
tomwalters@0:     /* Allocate space (first-time call only). */
tomwalters@0:     if (first) {
tomwalters@0: 	first = 0;
tomwalters@0: 	/* Allocate space and compute base-2 log time for spoke detection */
tomwalters@0: 	if ( (sector  = (short *)malloc( m*sizeof(short) ) ) == NULL)
tomwalters@0: 	    fprintf(stderr,"sp_weights: malloc out of space\n");
tomwalters@0: 	if ( (logtime = (int *)malloc( framewidth*sizeof(int) ) ) == NULL )
tomwalters@0: 	    fprintf(stderr,"sp_weights: malloc out of space\n");
tomwalters@0: 	gen_logtime(logtime,m);
tomwalters@0:     }
tomwalters@0:     /* clear sector bins */
tomwalters@0:     for (i=0 ; i<m ; i++)
tomwalters@0: 	sector[i]=0;
tomwalters@0:     /* for each time point along the SAI frame (with time origin on right) */
tomwalters@0:     /* (exclude time-origin bin to prevent overweighting by this sample)   */
tomwalters@0:     for (t=round(pow(2.0,orig_spd)) ; t<framewidth ; t++)  {
tomwalters@0: 	col = (framewidth-1)-t;
tomwalters@0: 	/* sum frame column (ie across all frequency channels) */
tomwalters@0: 	for (i=0, k=0, sum=0 ; i<frameheight ; i++, k+=framewidth)
tomwalters@0: 	    sum += frame[col+k];
tomwalters@0: 	/* Normalize (divide by frameheight) for average magnitude */
tomwalters@0: 	sum /= frameheight;
tomwalters@0: 	/* Accumulate sum in log-time sector bin */
tomwalters@0: 	sector[logtime[t]] += sum;
tomwalters@0:     }
tomwalters@0:     /* write total sum activity in spokes around the spiral from the datum */
tomwalters@0:     /* line in an anticlockwise direction. (ie in "natural order").        */
tomwalters@0:     for (i=0 ; i<m ; i++)
tomwalters@0: /*        printf("%d\n", sector[i]);    */
tomwalters@0: 	fwrite(&sector[i], sizeof(short), 1, stdout);
tomwalters@0: }
tomwalters@0: 
tomwalters@0: /****************************************************************************
tomwalters@0:   List the sector numbers, 0,1,...,m-1, which correspond to each time
tomwalters@0:   t=1,2,...,framewidth-1. The return array, "logtime", maps time index onto
tomwalters@0:   a spiral sector number.
tomwalters@0:   (Note that since the sector width increases with radius, there is an
tomwalters@0:   increasing number of adjacent frame indices within the sector bin).
tomwalters@0: ****************************************************************************/
tomwalters@0: gen_logtime(logtime, m)
tomwalters@0: int  *logtime;
tomwalters@0: int   m;
tomwalters@0: {
tomwalters@0:     int     t;
tomwalters@0:     double  theta, fptr, iptr;
tomwalters@0: 
tomwalters@0:     logtime[0] = 0; /* special case for log(0) */
tomwalters@0:     for (t=1 ; t<framewidth ; t++)  {
tomwalters@0: 	theta = log2(t) - orig_spd;
tomwalters@0: 	fptr = modf(theta,&iptr);       /* find fractional part of log time */
tomwalters@0: 	/* scale up to given resolution, and round */
tomwalters@0: 	logtime[t] = (int)(0.5+(fptr*m));
tomwalters@0:     }
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: