view tools/pitch_strength.c @ 0:5242703e91d3 tip

Initial checkin for AIM92 aimR8.2 (last updated May 1997).

/***************************************************************************
  pitch_strength.c      Pitch strength measure from spiral sector-weights contour.
  ----------------
			Michael Allerhand, 1992.

	Read m binary shorts (default m=SECTORS) for each specified frame.
	Compute the pitch-strength of the sector-weights by projecting onto
	the direction in the sector-weights space (given by vector `vec')
	which best explains the variation between weak and strong pitch,
	as exemplified by training data generated from iterated rippled noise.

	The sector weights are normalized with respect to angle (pitch) and
	magnitude (loudness), and also by extracting features using the
	method of Fourier descriptors.

****************************************************************************/

#include <stdio.h>
#include <math.h>
#include "options.h"
#include "units.h"
#include "strmatch.h"

char applic[] = "Pitch-strength measure from spiral sector-weights contour.\n            (Input and output in binary shorts)." ;

static char *helpstr, *debugstr, *fstr, *vstr, *mstr, *shiftstr, *scalestr ;

static Options option[] = {
    {   "help"      ,   "off"       ,  &helpstr     ,   "help"                                  , DEBUG   },
    {   "debug"     ,   "off"       ,  &debugstr    ,   "debugging switch"                      , DEBUG   },
    {   "frames"    ,   "1-max"     ,  &fstr        ,   "select frames inclusively"             , VAL     },
    {   "nsectors"  ,   "64"        ,  &mstr        ,   "number of sectors"                     , SVAL    },
    {   "vector"    ,   "off"       ,  &vstr        ,   "file for alternative projection vector", SVAL    },
    {   "shift"     ,   "212.0"     ,  &shiftstr    ,   "scale shifting by addition"            , SVAL    },
    {   "scale"     ,   "0.7686"    ,  &scalestr    ,   "scaling by multiplication"             , SVAL    },
   ( char * ) 0 } ;


/* projection vector for 64-sector cntour */

char *defvec[] = { "0.000",  "0.150", "-0.260", "-0.190", "-0.252", "-0.372", "-0.121",
		  "-0.374", "-0.150", "-0.236", "-0.157", "-0.109", "-0.246", "-0.088",
		  "-0.045", "-0.179", "-0.164", "-0.280", "-0.028", "-0.116", "-0.013",
		  "-0.065", "-0.252", "-0.005", "-0.216", "-0.086",  "0.043", "-0.061",
		  "-0.055",  "0.084", "-0.045", "-0.203", "211.670" } ;


float shift, scale ;


main (argc, argv)
int     argc;
char ** argv;
{
    FILE   *fp, *fp1 ;
    int     m, n;               /* Number of sectors around spiral.         */
    short  *sector;             /* Sector-weights contour.                  */
    float  *vec;
    int     i, a, b;
    char   *val1, *val2 ;

    fp = openopts( option,argc,argv ) ;
    if ( !isoff( helpstr ) )
	helpopts( helpstr, argv[0], applic, option ) ;

    shift = atof( shiftstr ) ;
    scale = atof( scalestr ) ;

    m = atoi( mstr ) ;
    n = m>>1 ;

    /* parse bounds on number of frames */

    if ( getvals( fstr, &val1, &val2, "-" ) == BADVAL ) {
	fprintf(stderr,"loudness: bad frame selector [%s]\n", fstr ) ;
	exit( 1 ) ;
    }
    a = atoi( val1 ) ;
    if ( isempty( val2 ) )    b = a ;
    else if ( ismax( val2 ) ) b = 0 ;
    else                      b = atoi( val2 ) ;

    if (b<a && b>0) fprintf(stderr,"pitch_strength: warning, bad frame specifiers\n");

    if ( ison( debugstr ) ) {
	printf("a=%d b=%d m=%d vector=%s\n",
		a, b, m, vstr );
	exit(0);
    }


    /* Allocate space for sector-weights contour and discriminant vector */

    if (( sector = (short *)malloc( m*sizeof(short) )) == NULL)
	    fprintf(stderr,"pitch_strength: malloc out of space\n");
    if (( vec = (float *)malloc( (n+1)*sizeof(float) )) == NULL)
	fprintf(stderr,"pitch_strength: malloc out of space\n");


    /* Assign discriminant vector */

    if ( isoff( vstr) )                 /* copy internal (default) vector */
	for (i=0 ; i<=n ; i++)
	    vec[i] = atof( defvec[i] ) ;
    else {                              /* read alternative vector from file */
	if ((fp1 = fopen( vstr, "r" )) == NULL) {
	    fprintf(stderr,"pitch_strength: can't open %s\n", vstr );
	    exit(1);
	}
	if ( fread(vec, sizeof(float), n+1, fp1) == NULL ) {
	    fprintf(stderr,"pitch_strength: empty discriminant vector\n");
	    exit(1);
	}
	fclose(fp1);
    }


    /* Seek past a-1 frames, then read the next b-a+1 frames, or all    */
    /* remaining frames if b==0.                                        */

    if (a>0) {

	for (i=1 ; i<a  && fread(sector,sizeof(short),m,fp) ; i++)
	    ;
	if (b > 0)
	    for (  ; i<=b && fread(sector,sizeof(short),m,fp) ; i++)
		pitch_strength( sector, m, n, vec );
	else
	    while ( fread(sector,sizeof(short),m,fp) )
		pitch_strength( sector, m, n, vec );
    }

    fclose(fp);
}


pitch_strength( sector, m, n, vec )
short  *sector;         /* Sector-weights contour.                  */
int     m, n;           /* Number of sectors around spiral.         */
float  *vec;
{
    float  inner_product();
    short  w0;

    normalize_scale(sector,m);
    normalize_angle(sector,m);
    Fd(sector,m,n);
    sector[0] = 0;      /* normalize for translation */
    sector[n] = 1;      /* Augment each data vector with 1 */
    w0 = (short)( ( shift - inner_product( vec, sector, n+1 ) ) * scale ) ;
    fwrite (&w0, sizeof(short), 1, stdout);
}


/****************************************************************************
    Inner product (for projection).
****************************************************************************/

float inner_product(W,V,m)
float *W;
short *V;
int    m;
{
    int   i;
    float w0=0;

    for (i=0 ; i<m ; i++)
	w0 += W[i] * V[i];
    return w0;
}


/****************************************************************************
    Normalize Scale.
****************************************************************************/

/*
Normalize for unit variance and zero mean, and then scale for MEAN and STDDEV.
In general the stddev is scaled down to STDDEV, but if the stddev is already
less than the target STDDEV, than it is left alone.
*/

#define MEAN    100
#define STDDEV  32       /* Resulting variance will be this squared */

normalize_scale(V,m)
short *V;
int    m;
{
    int   i;
    float   mean, var;
    float sum, sumsq, stddev, norm;

    sum=0; sumsq=0;
    for (i=0 ; i<m ; i++) {
	sum   += V[i];
	sumsq += V[i] * V[i];
    }
    mean = sum/m;
    var  = sumsq/m - mean*mean;

    if ((stddev = sqrt(var)) >= STDDEV)
	for (i=0 ; i<m ; i++) {
	    norm = ( (float)V[i] - mean ) / stddev;
	    V[i] = norm * STDDEV + MEAN;
	}
    else
	for (i=0 ;  i<m ; i++) {
	    norm = (float)V[i]-mean;
	    V[i] = norm + MEAN;
	}

    return (int)mean;
}


/****************************************************************************
    Normalize Angle (orientation).
****************************************************************************/

/*
  For a given spoke number k (where the primary spoke is number 1, and the
  subsequent spokes are numbered 2,3,...,m), return the angle as a number of
  sectors from the primary spoke.
  In general, the angle between the primary spoke and the k'th spoke is
  given by the fractional part of log2(2k-1). This angle is a fraction of
  unity (the complete cycle), so in terms of m sectors, we simply multiply up
  by m.
*/

#define log2(x)     (log((double)x)/log(2.0))

angle(k,m)
short k, m;
{
    double  theta, fptr, iptr;

    if (k==0) return 0;
    theta = log2(2*k-1);
    fptr = modf(theta,&iptr);       /* find fractional part of log time */
    /* scale up to given resolution, and round */
    return ((int)(0.5+(fptr*m)));
}


/*
 Normalize orientation.
 Find the orientation by rotating a template of four sectors set at relative
 angles appropriate to the first four spokes, (as the angle between spokes is
 known). The angle at which the weighted sum of these sectors is a maximum,
 as it is rotated, gives the orientation of the sector-weights contour.
 This angle of orientation, phi, is returned.
 The sector-weights contour is normalized by rotating, (ie sorting), so that
 the principal spoke, (the first of the four), is in the datum position, (at
 sector[0]).
*/

int normalize_angle(sector,m)
short *sector;
int    m;       /* number of sectors */
{
    int    i, k, phi;
    float  q, max;
    short *tmp; /* temporary array for sorting */

    /* Find angle of orientation phi */

    phi = 0; max = 0;
    for (k=0 ; k<m ; k++) {
	q = (float)sector[k] + 0.8*(float)sector[(k+angle(2,m))%m] +
	    0.6*(float)sector[(k+angle(3,m))%m] + 0.4*(float)sector[(k+angle(4,m))%m] ;

	if (q > max) {
	    max = q;
	    phi = k;
	}
    }

    /* Sort sector array to rotate contour into normalized orientation */

    if (( tmp = (short *)malloc(m*sizeof(short))) == NULL)
	fprintf(stderr,"pitch_strength: malloc out of space\n");

    for (i=0, k=phi ; k<m ; i++, k++)
	tmp[i] = sector[k];
    for (k=0 ; k<phi ; i++, k++)
	tmp[i] = sector[k];
    for (k=0 ; k<m ; k++)
	sector[k] = tmp[k];

    /* Return orientation, (angle increases clockwise from 0 to m-1) */

    if (phi == 0) return 0;
    else          return ( m-phi );

}


/****************************************************************************
 Map m sector-weights onto m/2 Fourier descriptors.
****************************************************************************/

typedef struct {
    float real;
    float imag;
} complex;

#define PI              3.14159265
#define TWOPI           6.2831853
#define sq(x)           ((x)*(x))
#define Re(C)           ( (C)->real )   /* cartesian form */
#define Im(C)           ( (C)->imag )
#define Mod(C)          ( (C)->real )   /* polar form     */
#define Arg(C)          ( (C)->imag )
#define cscalar_mod(C)  ( sqrt(sq(Re(C))+sq(Im(C))) )
#define FFT_real(V,n)   ( realft((float *)(V)-1, n>>1 ,1) )


Fd(sector,m,n)
short *sector;
int    m, n;
{
    static float   *V;
    static char first=1;
    int    i, j;

    if (first) {        /* Allocate space first time only */
	if (( V = (float *)  malloc( m*sizeof(float)   )) == NULL) {
	    fprintf(stderr,"pitch_strength: malloc out of space\n");
	    exit(1);
	}
	first=0;
    }

    for (i=0 ; i<m ; i++)
	V[i] =  sector[i] ;
    FFT_real(V,m);
    cvector_mod( sector, (complex *)V, n ) ; /* magnitude spectrum */

}

cvector_mod(V,C,m)
short   *V;
complex *C;
int      m;
{
    int   i;

    for (i=0; i<m ; i++)
	V[i] = (short)( cscalar_mod(&C[i]) ) ;
}


realft(data,n,isign)
float data[];
int n,isign;
{
	int i,i1,i2,i3,i4,n2p3;
	float c1=0.5,c2,h1r,h1i,h2r,h2i;
	double wr,wi,wpr,wpi,wtemp,theta;


	theta=3.141592653589793/(double) n;
	if (isign == 1) {
		c2 = -0.5;
		four1(data,n,1);
	} else {
		c2=0.5;
		theta = -theta;
	}
	wtemp=sin(0.5*theta);
	wpr = -2.0*wtemp*wtemp;
	wpi=sin(theta);
	wr=1.0+wpr;
	wi=wpi;
	n2p3=2*n+3;
	for (i=2;i<=n/2;i++) {
		i4=1+(i3=n2p3-(i2=1+(i1=i+i-1)));
		h1r=c1*(data[i1]+data[i3]);
		h1i=c1*(data[i2]-data[i4]);
		h2r = -c2*(data[i2]+data[i4]);
		h2i=c2*(data[i1]-data[i3]);
		data[i1]=h1r+wr*h2r-wi*h2i;
		data[i2]=h1i+wr*h2i+wi*h2r;
		data[i3]=h1r-wr*h2r+wi*h2i;
		data[i4] = -h1i+wr*h2i+wi*h2r;
		wr=(wtemp=wr)*wpr-wi*wpi+wr;
		wi=wi*wpr+wtemp*wpi+wi;
	}
	if (isign == 1) {
		data[1] = (h1r=data[1])+data[2];
		data[2] = h1r-data[2];
	} else {
		data[1]=c1*((h1r=data[1])+data[2]);
		data[2]=c1*(h1r-data[2]);
		four1(data,n,-1);
	}
}

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

four1(data,nn,isign)
float data[];
int nn,isign;
{
	int n,mmax,m,j,istep,i;
	double wtemp,wr,wpr,wpi,wi,theta;
	float tempr,tempi;

	n=nn << 1;
	j=1;
	for (i=1;i<n;i+=2) {
		if (j > i) {
			SWAP(data[j],data[i]);
			SWAP(data[j+1],data[i+1]);
		}
		m=n >> 1;
		while (m >= 2 && j > m) {
			j -= m;
			m >>= 1;
		}
		j += m;
	}
	mmax=2;
	while (n > mmax) {
		istep=2*mmax;
		theta=6.28318530717959/(isign*mmax);
		wtemp=sin(0.5*theta);
		wpr = -2.0*wtemp*wtemp;
		wpi=sin(theta);
		wr=1.0;
		wi=0.0;
		for (m=1;m<mmax;m+=2) {
			for (i=m;i<=n;i+=istep) {
				j=i+mmax;
				tempr=wr*data[j]-wi*data[j+1];
				tempi=wr*data[j+1]+wi*data[j];
				data[j]=data[i]-tempr;
				data[j+1]=data[i+1]-tempi;
				data[i] += tempr;
				data[i+1] += tempi;
			}
			wr=(wtemp=wr)*wpr-wi*wpi+wr;
			wi=wi*wpr+wtemp*wpi+wi;
		}
		mmax=istep;
	}
}