tomwalters@0: /***************************************************************************
tomwalters@0:   pitch_strength.c      Pitch strength 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 pitch-strength of the sector-weights by projecting onto
tomwalters@0: 	the direction in the sector-weights space (given by vector `vec')
tomwalters@0: 	which best explains the variation between weak and strong pitch,
tomwalters@0: 	as exemplified by training data generated from iterated rippled noise.
tomwalters@0: 
tomwalters@0: 	The sector weights are normalized with respect to angle (pitch) and
tomwalters@0: 	magnitude (loudness), and also by extracting features using the
tomwalters@0: 	method of Fourier descriptors.
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[] = "Pitch-strength measure from spiral sector-weights contour.\n            (Input and output in binary shorts)." ;
tomwalters@0: 
tomwalters@0: static char *helpstr, *debugstr, *fstr, *vstr, *mstr, *shiftstr, *scalestr ;
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:     {   "nsectors"  ,   "64"        ,  &mstr        ,   "number of sectors"                     , SVAL    },
tomwalters@0:     {   "vector"    ,   "off"       ,  &vstr        ,   "file for alternative projection vector", SVAL    },
tomwalters@0:     {   "shift"     ,   "212.0"     ,  &shiftstr    ,   "scale shifting by addition"            , SVAL    },
tomwalters@0:     {   "scale"     ,   "0.7686"    ,  &scalestr    ,   "scaling by multiplication"             , SVAL    },
tomwalters@0:    ( char * ) 0 } ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /* projection vector for 64-sector cntour */
tomwalters@0: 
tomwalters@0: char *defvec[] = { "0.000",  "0.150", "-0.260", "-0.190", "-0.252", "-0.372", "-0.121",
tomwalters@0: 		  "-0.374", "-0.150", "-0.236", "-0.157", "-0.109", "-0.246", "-0.088",
tomwalters@0: 		  "-0.045", "-0.179", "-0.164", "-0.280", "-0.028", "-0.116", "-0.013",
tomwalters@0: 		  "-0.065", "-0.252", "-0.005", "-0.216", "-0.086",  "0.043", "-0.061",
tomwalters@0: 		  "-0.055",  "0.084", "-0.045", "-0.203", "211.670" } ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: float shift, scale ;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: main (argc, argv)
tomwalters@0: int     argc;
tomwalters@0: char ** argv;
tomwalters@0: {
tomwalters@0:     FILE   *fp, *fp1 ;
tomwalters@0:     int     m, n;               /* Number of sectors around spiral.         */
tomwalters@0:     short  *sector;             /* Sector-weights contour.                  */
tomwalters@0:     float  *vec;
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:     shift = atof( shiftstr ) ;
tomwalters@0:     scale = atof( scalestr ) ;
tomwalters@0: 
tomwalters@0:     m = atoi( mstr ) ;
tomwalters@0:     n = m>>1 ;
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,"pitch_strength: warning, bad frame specifiers\n");
tomwalters@0: 
tomwalters@0:     if ( ison( debugstr ) ) {
tomwalters@0: 	printf("a=%d b=%d m=%d vector=%s\n",
tomwalters@0: 		a, b, m, vstr );
tomwalters@0: 	exit(0);
tomwalters@0:     }
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /* Allocate space for sector-weights contour and discriminant vector */
tomwalters@0: 
tomwalters@0:     if (( sector = (short *)malloc( m*sizeof(short) )) == NULL)
tomwalters@0: 	    fprintf(stderr,"pitch_strength: malloc out of space\n");
tomwalters@0:     if (( vec = (float *)malloc( (n+1)*sizeof(float) )) == NULL)
tomwalters@0: 	fprintf(stderr,"pitch_strength: malloc out of space\n");
tomwalters@0: 
tomwalters@0: 
tomwalters@0:     /* Assign discriminant vector */
tomwalters@0: 
tomwalters@0:     if ( isoff( vstr) )                 /* copy internal (default) vector */
tomwalters@0: 	for (i=0 ; i<=n ; i++)
tomwalters@0: 	    vec[i] = atof( defvec[i] ) ;
tomwalters@0:     else {                              /* read alternative vector from file */
tomwalters@0: 	if ((fp1 = fopen( vstr, "r" )) == NULL) {
tomwalters@0: 	    fprintf(stderr,"pitch_strength: can't open %s\n", vstr );
tomwalters@0: 	    exit(1);
tomwalters@0: 	}
tomwalters@0: 	if ( fread(vec, sizeof(float), n+1, fp1) == NULL ) {
tomwalters@0: 	    fprintf(stderr,"pitch_strength: empty discriminant vector\n");
tomwalters@0: 	    exit(1);
tomwalters@0: 	}
tomwalters@0: 	fclose(fp1);
tomwalters@0:     }
tomwalters@0: 
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: 		pitch_strength( sector, m, n, vec );
tomwalters@0: 	else
tomwalters@0: 	    while ( fread(sector,sizeof(short),m,fp) )
tomwalters@0: 		pitch_strength( sector, m, n, vec );
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     fclose(fp);
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: pitch_strength( sector, m, n, vec )
tomwalters@0: short  *sector;         /* Sector-weights contour.                  */
tomwalters@0: int     m, n;           /* Number of sectors around spiral.         */
tomwalters@0: float  *vec;
tomwalters@0: {
tomwalters@0:     float  inner_product();
tomwalters@0:     short  w0;
tomwalters@0: 
tomwalters@0:     normalize_scale(sector,m);
tomwalters@0:     normalize_angle(sector,m);
tomwalters@0:     Fd(sector,m,n);
tomwalters@0:     sector[0] = 0;      /* normalize for translation */
tomwalters@0:     sector[n] = 1;      /* Augment each data vector with 1 */
tomwalters@0:     w0 = (short)( ( shift - inner_product( vec, sector, n+1 ) ) * scale ) ;
tomwalters@0:     fwrite (&w0, sizeof(short), 1, stdout);
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /****************************************************************************
tomwalters@0:     Inner product (for projection).
tomwalters@0: ****************************************************************************/
tomwalters@0: 
tomwalters@0: float inner_product(W,V,m)
tomwalters@0: float *W;
tomwalters@0: short *V;
tomwalters@0: int    m;
tomwalters@0: {
tomwalters@0:     int   i;
tomwalters@0:     float w0=0;
tomwalters@0: 
tomwalters@0:     for (i=0 ; i<m ; i++)
tomwalters@0: 	w0 += W[i] * V[i];
tomwalters@0:     return w0;
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: 
tomwalters@0: 
tomwalters@0: /****************************************************************************
tomwalters@0:     Normalize Angle (orientation).
tomwalters@0: ****************************************************************************/
tomwalters@0: 
tomwalters@0: /*
tomwalters@0:   For a given spoke number k (where the primary spoke is number 1, and the
tomwalters@0:   subsequent spokes are numbered 2,3,...,m), return the angle as a number of
tomwalters@0:   sectors from the primary spoke.
tomwalters@0:   In general, the angle between the primary spoke and the k'th spoke is
tomwalters@0:   given by the fractional part of log2(2k-1). This angle is a fraction of
tomwalters@0:   unity (the complete cycle), so in terms of m sectors, we simply multiply up
tomwalters@0:   by m.
tomwalters@0: */
tomwalters@0: 
tomwalters@0: #define log2(x)     (log((double)x)/log(2.0))
tomwalters@0: 
tomwalters@0: angle(k,m)
tomwalters@0: short k, m;
tomwalters@0: {
tomwalters@0:     double  theta, fptr, iptr;
tomwalters@0: 
tomwalters@0:     if (k==0) return 0;
tomwalters@0:     theta = log2(2*k-1);
tomwalters@0:     fptr = modf(theta,&iptr);       /* find fractional part of log time */
tomwalters@0:     /* scale up to given resolution, and round */
tomwalters@0:     return ((int)(0.5+(fptr*m)));
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /*
tomwalters@0:  Normalize orientation.
tomwalters@0:  Find the orientation by rotating a template of four sectors set at relative
tomwalters@0:  angles appropriate to the first four spokes, (as the angle between spokes is
tomwalters@0:  known). The angle at which the weighted sum of these sectors is a maximum,
tomwalters@0:  as it is rotated, gives the orientation of the sector-weights contour.
tomwalters@0:  This angle of orientation, phi, is returned.
tomwalters@0:  The sector-weights contour is normalized by rotating, (ie sorting), so that
tomwalters@0:  the principal spoke, (the first of the four), is in the datum position, (at
tomwalters@0:  sector[0]).
tomwalters@0: */
tomwalters@0: 
tomwalters@0: int normalize_angle(sector,m)
tomwalters@0: short *sector;
tomwalters@0: int    m;       /* number of sectors */
tomwalters@0: {
tomwalters@0:     int    i, k, phi;
tomwalters@0:     float  q, max;
tomwalters@0:     short *tmp; /* temporary array for sorting */
tomwalters@0: 
tomwalters@0:     /* Find angle of orientation phi */
tomwalters@0: 
tomwalters@0:     phi = 0; max = 0;
tomwalters@0:     for (k=0 ; k<m ; k++) {
tomwalters@0: 	q = (float)sector[k] + 0.8*(float)sector[(k+angle(2,m))%m] +
tomwalters@0: 	    0.6*(float)sector[(k+angle(3,m))%m] + 0.4*(float)sector[(k+angle(4,m))%m] ;
tomwalters@0: 
tomwalters@0: 	if (q > max) {
tomwalters@0: 	    max = q;
tomwalters@0: 	    phi = k;
tomwalters@0: 	}
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     /* Sort sector array to rotate contour into normalized orientation */
tomwalters@0: 
tomwalters@0:     if (( tmp = (short *)malloc(m*sizeof(short))) == NULL)
tomwalters@0: 	fprintf(stderr,"pitch_strength: malloc out of space\n");
tomwalters@0: 
tomwalters@0:     for (i=0, k=phi ; k<m ; i++, k++)
tomwalters@0: 	tmp[i] = sector[k];
tomwalters@0:     for (k=0 ; k<phi ; i++, k++)
tomwalters@0: 	tmp[i] = sector[k];
tomwalters@0:     for (k=0 ; k<m ; k++)
tomwalters@0: 	sector[k] = tmp[k];
tomwalters@0: 
tomwalters@0:     /* Return orientation, (angle increases clockwise from 0 to m-1) */
tomwalters@0: 
tomwalters@0:     if (phi == 0) return 0;
tomwalters@0:     else          return ( m-phi );
tomwalters@0: 
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: /****************************************************************************
tomwalters@0:  Map m sector-weights onto m/2 Fourier descriptors.
tomwalters@0: ****************************************************************************/
tomwalters@0: 
tomwalters@0: typedef struct {
tomwalters@0:     float real;
tomwalters@0:     float imag;
tomwalters@0: } complex;
tomwalters@0: 
tomwalters@0: #define PI              3.14159265
tomwalters@0: #define TWOPI           6.2831853
tomwalters@0: #define sq(x)           ((x)*(x))
tomwalters@0: #define Re(C)           ( (C)->real )   /* cartesian form */
tomwalters@0: #define Im(C)           ( (C)->imag )
tomwalters@0: #define Mod(C)          ( (C)->real )   /* polar form     */
tomwalters@0: #define Arg(C)          ( (C)->imag )
tomwalters@0: #define cscalar_mod(C)  ( sqrt(sq(Re(C))+sq(Im(C))) )
tomwalters@0: #define FFT_real(V,n)   ( realft((float *)(V)-1, n>>1 ,1) )
tomwalters@0: 
tomwalters@0: 
tomwalters@0: Fd(sector,m,n)
tomwalters@0: short *sector;
tomwalters@0: int    m, n;
tomwalters@0: {
tomwalters@0:     static float   *V;
tomwalters@0:     static char first=1;
tomwalters@0:     int    i, j;
tomwalters@0: 
tomwalters@0:     if (first) {        /* Allocate space first time only */
tomwalters@0: 	if (( V = (float *)  malloc( m*sizeof(float)   )) == NULL) {
tomwalters@0: 	    fprintf(stderr,"pitch_strength: malloc out of space\n");
tomwalters@0: 	    exit(1);
tomwalters@0: 	}
tomwalters@0: 	first=0;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     for (i=0 ; i<m ; i++)
tomwalters@0: 	V[i] =  sector[i] ;
tomwalters@0:     FFT_real(V,m);
tomwalters@0:     cvector_mod( sector, (complex *)V, n ) ; /* magnitude spectrum */
tomwalters@0: 
tomwalters@0: }
tomwalters@0: 
tomwalters@0: cvector_mod(V,C,m)
tomwalters@0: short   *V;
tomwalters@0: complex *C;
tomwalters@0: int      m;
tomwalters@0: {
tomwalters@0:     int   i;
tomwalters@0: 
tomwalters@0:     for (i=0; i<m ; i++)
tomwalters@0: 	V[i] = (short)( cscalar_mod(&C[i]) ) ;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: 
tomwalters@0: realft(data,n,isign)
tomwalters@0: float data[];
tomwalters@0: int n,isign;
tomwalters@0: {
tomwalters@0: 	int i,i1,i2,i3,i4,n2p3;
tomwalters@0: 	float c1=0.5,c2,h1r,h1i,h2r,h2i;
tomwalters@0: 	double wr,wi,wpr,wpi,wtemp,theta;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 	theta=3.141592653589793/(double) n;
tomwalters@0: 	if (isign == 1) {
tomwalters@0: 		c2 = -0.5;
tomwalters@0: 		four1(data,n,1);
tomwalters@0: 	} else {
tomwalters@0: 		c2=0.5;
tomwalters@0: 		theta = -theta;
tomwalters@0: 	}
tomwalters@0: 	wtemp=sin(0.5*theta);
tomwalters@0: 	wpr = -2.0*wtemp*wtemp;
tomwalters@0: 	wpi=sin(theta);
tomwalters@0: 	wr=1.0+wpr;
tomwalters@0: 	wi=wpi;
tomwalters@0: 	n2p3=2*n+3;
tomwalters@0: 	for (i=2;i<=n/2;i++) {
tomwalters@0: 		i4=1+(i3=n2p3-(i2=1+(i1=i+i-1)));
tomwalters@0: 		h1r=c1*(data[i1]+data[i3]);
tomwalters@0: 		h1i=c1*(data[i2]-data[i4]);
tomwalters@0: 		h2r = -c2*(data[i2]+data[i4]);
tomwalters@0: 		h2i=c2*(data[i1]-data[i3]);
tomwalters@0: 		data[i1]=h1r+wr*h2r-wi*h2i;
tomwalters@0: 		data[i2]=h1i+wr*h2i+wi*h2r;
tomwalters@0: 		data[i3]=h1r-wr*h2r+wi*h2i;
tomwalters@0: 		data[i4] = -h1i+wr*h2i+wi*h2r;
tomwalters@0: 		wr=(wtemp=wr)*wpr-wi*wpi+wr;
tomwalters@0: 		wi=wi*wpr+wtemp*wpi+wi;
tomwalters@0: 	}
tomwalters@0: 	if (isign == 1) {
tomwalters@0: 		data[1] = (h1r=data[1])+data[2];
tomwalters@0: 		data[2] = h1r-data[2];
tomwalters@0: 	} else {
tomwalters@0: 		data[1]=c1*((h1r=data[1])+data[2]);
tomwalters@0: 		data[2]=c1*(h1r-data[2]);
tomwalters@0: 		four1(data,n,-1);
tomwalters@0: 	}
tomwalters@0: }
tomwalters@0: 
tomwalters@0: #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
tomwalters@0: 
tomwalters@0: four1(data,nn,isign)
tomwalters@0: float data[];
tomwalters@0: int nn,isign;
tomwalters@0: {
tomwalters@0: 	int n,mmax,m,j,istep,i;
tomwalters@0: 	double wtemp,wr,wpr,wpi,wi,theta;
tomwalters@0: 	float tempr,tempi;
tomwalters@0: 
tomwalters@0: 	n=nn << 1;
tomwalters@0: 	j=1;
tomwalters@0: 	for (i=1;i<n;i+=2) {
tomwalters@0: 		if (j > i) {
tomwalters@0: 			SWAP(data[j],data[i]);
tomwalters@0: 			SWAP(data[j+1],data[i+1]);
tomwalters@0: 		}
tomwalters@0: 		m=n >> 1;
tomwalters@0: 		while (m >= 2 && j > m) {
tomwalters@0: 			j -= m;
tomwalters@0: 			m >>= 1;
tomwalters@0: 		}
tomwalters@0: 		j += m;
tomwalters@0: 	}
tomwalters@0: 	mmax=2;
tomwalters@0: 	while (n > mmax) {
tomwalters@0: 		istep=2*mmax;
tomwalters@0: 		theta=6.28318530717959/(isign*mmax);
tomwalters@0: 		wtemp=sin(0.5*theta);
tomwalters@0: 		wpr = -2.0*wtemp*wtemp;
tomwalters@0: 		wpi=sin(theta);
tomwalters@0: 		wr=1.0;
tomwalters@0: 		wi=0.0;
tomwalters@0: 		for (m=1;m<mmax;m+=2) {
tomwalters@0: 			for (i=m;i<=n;i+=istep) {
tomwalters@0: 				j=i+mmax;
tomwalters@0: 				tempr=wr*data[j]-wi*data[j+1];
tomwalters@0: 				tempi=wr*data[j+1]+wi*data[j];
tomwalters@0: 				data[j]=data[i]-tempr;
tomwalters@0: 				data[j+1]=data[i+1]-tempi;
tomwalters@0: 				data[i] += tempr;
tomwalters@0: 				data[i+1] += tempi;
tomwalters@0: 			}
tomwalters@0: 			wr=(wtemp=wr)*wpr-wi*wpi+wr;
tomwalters@0: 			wi=wi*wpr+wtemp*wpi+wi;
tomwalters@0: 		}
tomwalters@0: 		mmax=istep;
tomwalters@0: 	}
tomwalters@0: }
tomwalters@0: