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

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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tools/pitch_strength.c	Fri May 20 15:19:45 2011 +0100
@@ -0,0 +1,449 @@
+/***************************************************************************
+  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;
+	}
+}
+