Mercurial > hg > aim92

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial checkin for AIM92 aimR8.2 (last updated May 1997).
author tomwalters
date Fri, 20 May 2011 15:19:45 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:5242703e91d3
1 /***************************************************************************
2 pitch_strength.c Pitch strength measure from spiral sector-weights contour.
3 ----------------
4 Michael Allerhand, 1992.
5
6 Read m binary shorts (default m=SECTORS) for each specified frame.
7 Compute the pitch-strength of the sector-weights by projecting onto
8 the direction in the sector-weights space (given by vector `vec')
9 which best explains the variation between weak and strong pitch,
10 as exemplified by training data generated from iterated rippled noise.
11
12 The sector weights are normalized with respect to angle (pitch) and
13 magnitude (loudness), and also by extracting features using the
14 method of Fourier descriptors.
15
16 ****************************************************************************/
17
18 #include <stdio.h>
19 #include <math.h>
20 #include "options.h"
21 #include "units.h"
22 #include "strmatch.h"
23
24 char applic[] = "Pitch-strength measure from spiral sector-weights contour.\n (Input and output in binary shorts)." ;
25
26 static char *helpstr, *debugstr, *fstr, *vstr, *mstr, *shiftstr, *scalestr ;
27
28 static Options option[] = {
29 { "help" , "off" , &helpstr , "help" , DEBUG },
30 { "debug" , "off" , &debugstr , "debugging switch" , DEBUG },
31 { "frames" , "1-max" , &fstr , "select frames inclusively" , VAL },
32 { "nsectors" , "64" , &mstr , "number of sectors" , SVAL },
33 { "vector" , "off" , &vstr , "file for alternative projection vector", SVAL },
34 { "shift" , "212.0" , &shiftstr , "scale shifting by addition" , SVAL },
35 { "scale" , "0.7686" , &scalestr , "scaling by multiplication" , SVAL },
36 ( char * ) 0 } ;
37
38
39 /* projection vector for 64-sector cntour */
40
41 char *defvec[] = { "0.000", "0.150", "-0.260", "-0.190", "-0.252", "-0.372", "-0.121",
42 "-0.374", "-0.150", "-0.236", "-0.157", "-0.109", "-0.246", "-0.088",
43 "-0.045", "-0.179", "-0.164", "-0.280", "-0.028", "-0.116", "-0.013",
44 "-0.065", "-0.252", "-0.005", "-0.216", "-0.086", "0.043", "-0.061",
45 "-0.055", "0.084", "-0.045", "-0.203", "211.670" } ;
46
47
48 float shift, scale ;
49
50
51 main (argc, argv)
52 int argc;
53 char ** argv;
54 {
55 FILE *fp, *fp1 ;
56 int m, n; /* Number of sectors around spiral. */
57 short *sector; /* Sector-weights contour. */
58 float *vec;
59 int i, a, b;
60 char *val1, *val2 ;
61
62 fp = openopts( option,argc,argv ) ;
63 if ( !isoff( helpstr ) )
64 helpopts( helpstr, argv[0], applic, option ) ;
65
66 shift = atof( shiftstr ) ;
67 scale = atof( scalestr ) ;
68
69 m = atoi( mstr ) ;
70 n = m>>1 ;
71
72 /* parse bounds on number of frames */
73
74 if ( getvals( fstr, &val1, &val2, "-" ) == BADVAL ) {
75 fprintf(stderr,"loudness: bad frame selector [%s]\n", fstr ) ;
76 exit( 1 ) ;
77 }
78 a = atoi( val1 ) ;
79 if ( isempty( val2 ) ) b = a ;
80 else if ( ismax( val2 ) ) b = 0 ;
81 else b = atoi( val2 ) ;
82
83 if (b<a && b>0) fprintf(stderr,"pitch_strength: warning, bad frame specifiers\n");
84
85 if ( ison( debugstr ) ) {
86 printf("a=%d b=%d m=%d vector=%s\n",
87 a, b, m, vstr );
88 exit(0);
89 }
90
91
92 /* Allocate space for sector-weights contour and discriminant vector */
93
94 if (( sector = (short *)malloc( m*sizeof(short) )) == NULL)
95 fprintf(stderr,"pitch_strength: malloc out of space\n");
96 if (( vec = (float *)malloc( (n+1)*sizeof(float) )) == NULL)
97 fprintf(stderr,"pitch_strength: malloc out of space\n");
98
99
100 /* Assign discriminant vector */
101
102 if ( isoff( vstr) ) /* copy internal (default) vector */
103 for (i=0 ; i<=n ; i++)
104 vec[i] = atof( defvec[i] ) ;
105 else { /* read alternative vector from file */
106 if ((fp1 = fopen( vstr, "r" )) == NULL) {
107 fprintf(stderr,"pitch_strength: can't open %s\n", vstr );
108 exit(1);
109 }
110 if ( fread(vec, sizeof(float), n+1, fp1) == NULL ) {
111 fprintf(stderr,"pitch_strength: empty discriminant vector\n");
112 exit(1);
113 }
114 fclose(fp1);
115 }
116
117
118 /* Seek past a-1 frames, then read the next b-a+1 frames, or all */
119 /* remaining frames if b==0. */
120
121 if (a>0) {
122
123 for (i=1 ; i<a && fread(sector,sizeof(short),m,fp) ; i++)
124 ;
125 if (b > 0)
126 for ( ; i<=b && fread(sector,sizeof(short),m,fp) ; i++)
127 pitch_strength( sector, m, n, vec );
128 else
129 while ( fread(sector,sizeof(short),m,fp) )
130 pitch_strength( sector, m, n, vec );
131 }
132
133 fclose(fp);
134 }
135
136
137 pitch_strength( sector, m, n, vec )
138 short *sector; /* Sector-weights contour. */
139 int m, n; /* Number of sectors around spiral. */
140 float *vec;
141 {
142 float inner_product();
143 short w0;
144
145 normalize_scale(sector,m);
146 normalize_angle(sector,m);
147 Fd(sector,m,n);
148 sector[0] = 0; /* normalize for translation */
149 sector[n] = 1; /* Augment each data vector with 1 */
150 w0 = (short)( ( shift - inner_product( vec, sector, n+1 ) ) * scale ) ;
151 fwrite (&w0, sizeof(short), 1, stdout);
152 }
153
154
155 /****************************************************************************
156 Inner product (for projection).
157 ****************************************************************************/
158
159 float inner_product(W,V,m)
160 float *W;
161 short *V;
162 int m;
163 {
164 int i;
165 float w0=0;
166
167 for (i=0 ; i<m ; i++)
168 w0 += W[i] * V[i];
169 return w0;
170 }
171
172
173 /****************************************************************************
174 Normalize Scale.
175 ****************************************************************************/
176
177 /*
178 Normalize for unit variance and zero mean, and then scale for MEAN and STDDEV.
179 In general the stddev is scaled down to STDDEV, but if the stddev is already
180 less than the target STDDEV, than it is left alone.
181 */
182
183 #define MEAN 100
184 #define STDDEV 32 /* Resulting variance will be this squared */
185
186 normalize_scale(V,m)
187 short *V;
188 int m;
189 {
190 int i;
191 float mean, var;
192 float sum, sumsq, stddev, norm;
193
194 sum=0; sumsq=0;
195 for (i=0 ; i<m ; i++) {
196 sum += V[i];
197 sumsq += V[i] * V[i];
198 }
199 mean = sum/m;
200 var = sumsq/m - mean*mean;
201
202 if ((stddev = sqrt(var)) >= STDDEV)
203 for (i=0 ; i<m ; i++) {
204 norm = ( (float)V[i] - mean ) / stddev;
205 V[i] = norm * STDDEV + MEAN;
206 }
207 else
208 for (i=0 ; i<m ; i++) {
209 norm = (float)V[i]-mean;
210 V[i] = norm + MEAN;
211 }
212
213 return (int)mean;
214 }
215
216
217 /****************************************************************************
218 Normalize Angle (orientation).
219 ****************************************************************************/
220
221 /*
222 For a given spoke number k (where the primary spoke is number 1, and the
223 subsequent spokes are numbered 2,3,...,m), return the angle as a number of
224 sectors from the primary spoke.
225 In general, the angle between the primary spoke and the k'th spoke is
226 given by the fractional part of log2(2k-1). This angle is a fraction of
227 unity (the complete cycle), so in terms of m sectors, we simply multiply up
228 by m.
229 */
230
231 #define log2(x) (log((double)x)/log(2.0))
232
233 angle(k,m)
234 short k, m;
235 {
236 double theta, fptr, iptr;
237
238 if (k==0) return 0;
239 theta = log2(2*k-1);
240 fptr = modf(theta,&iptr); /* find fractional part of log time */
241 /* scale up to given resolution, and round */
242 return ((int)(0.5+(fptr*m)));
243 }
244
245
246 /*
247 Normalize orientation.
248 Find the orientation by rotating a template of four sectors set at relative
249 angles appropriate to the first four spokes, (as the angle between spokes is
250 known). The angle at which the weighted sum of these sectors is a maximum,
251 as it is rotated, gives the orientation of the sector-weights contour.
252 This angle of orientation, phi, is returned.
253 The sector-weights contour is normalized by rotating, (ie sorting), so that
254 the principal spoke, (the first of the four), is in the datum position, (at
255 sector[0]).
256 */
257
258 int normalize_angle(sector,m)
259 short *sector;
260 int m; /* number of sectors */
261 {
262 int i, k, phi;
263 float q, max;
264 short *tmp; /* temporary array for sorting */
265
266 /* Find angle of orientation phi */
267
268 phi = 0; max = 0;
269 for (k=0 ; k<m ; k++) {
270 q = (float)sector[k] + 0.8*(float)sector[(k+angle(2,m))%m] +
271 0.6*(float)sector[(k+angle(3,m))%m] + 0.4*(float)sector[(k+angle(4,m))%m] ;
272
273 if (q > max) {
274 max = q;
275 phi = k;
276 }
277 }
278
279 /* Sort sector array to rotate contour into normalized orientation */
280
281 if (( tmp = (short *)malloc(m*sizeof(short))) == NULL)
282 fprintf(stderr,"pitch_strength: malloc out of space\n");
283
284 for (i=0, k=phi ; k<m ; i++, k++)
285 tmp[i] = sector[k];
286 for (k=0 ; k<phi ; i++, k++)
287 tmp[i] = sector[k];
288 for (k=0 ; k<m ; k++)
289 sector[k] = tmp[k];
290
291 /* Return orientation, (angle increases clockwise from 0 to m-1) */
292
293 if (phi == 0) return 0;
294 else return ( m-phi );
295
296 }
297
298
299 /****************************************************************************
300 Map m sector-weights onto m/2 Fourier descriptors.
301 ****************************************************************************/
302
303 typedef struct {
304 float real;
305 float imag;
306 } complex;
307
308 #define PI 3.14159265
309 #define TWOPI 6.2831853
310 #define sq(x) ((x)*(x))
311 #define Re(C) ( (C)->real ) /* cartesian form */
312 #define Im(C) ( (C)->imag )
313 #define Mod(C) ( (C)->real ) /* polar form */
314 #define Arg(C) ( (C)->imag )
315 #define cscalar_mod(C) ( sqrt(sq(Re(C))+sq(Im(C))) )
316 #define FFT_real(V,n) ( realft((float *)(V)-1, n>>1 ,1) )
317
318
319 Fd(sector,m,n)
320 short *sector;
321 int m, n;
322 {
323 static float *V;
324 static char first=1;
325 int i, j;
326
327 if (first) { /* Allocate space first time only */
328 if (( V = (float *) malloc( m*sizeof(float) )) == NULL) {
329 fprintf(stderr,"pitch_strength: malloc out of space\n");
330 exit(1);
331 }
332 first=0;
333 }
334
335 for (i=0 ; i<m ; i++)
336 V[i] = sector[i] ;
337 FFT_real(V,m);
338 cvector_mod( sector, (complex *)V, n ) ; /* magnitude spectrum */
339
340 }
341
342 cvector_mod(V,C,m)
343 short *V;
344 complex *C;
345 int m;
346 {
347 int i;
348
349 for (i=0; i<m ; i++)
350 V[i] = (short)( cscalar_mod(&C[i]) ) ;
351 }
352
353
354 realft(data,n,isign)
355 float data[];
356 int n,isign;
357 {
358 int i,i1,i2,i3,i4,n2p3;
359 float c1=0.5,c2,h1r,h1i,h2r,h2i;
360 double wr,wi,wpr,wpi,wtemp,theta;
361
362
363 theta=3.141592653589793/(double) n;
364 if (isign == 1) {
365 c2 = -0.5;
366 four1(data,n,1);
367 } else {
368 c2=0.5;
369 theta = -theta;
370 }
371 wtemp=sin(0.5*theta);
372 wpr = -2.0*wtemp*wtemp;
373 wpi=sin(theta);
374 wr=1.0+wpr;
375 wi=wpi;
376 n2p3=2*n+3;
377 for (i=2;i<=n/2;i++) {
378 i4=1+(i3=n2p3-(i2=1+(i1=i+i-1)));
379 h1r=c1*(data[i1]+data[i3]);
380 h1i=c1*(data[i2]-data[i4]);
381 h2r = -c2*(data[i2]+data[i4]);
382 h2i=c2*(data[i1]-data[i3]);
383 data[i1]=h1r+wr*h2r-wi*h2i;
384 data[i2]=h1i+wr*h2i+wi*h2r;
385 data[i3]=h1r-wr*h2r+wi*h2i;
386 data[i4] = -h1i+wr*h2i+wi*h2r;
387 wr=(wtemp=wr)*wpr-wi*wpi+wr;
388 wi=wi*wpr+wtemp*wpi+wi;
389 }
390 if (isign == 1) {
391 data[1] = (h1r=data[1])+data[2];
392 data[2] = h1r-data[2];
393 } else {
394 data[1]=c1*((h1r=data[1])+data[2]);
395 data[2]=c1*(h1r-data[2]);
396 four1(data,n,-1);
397 }
398 }
399
400 #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
401
402 four1(data,nn,isign)
403 float data[];
404 int nn,isign;
405 {
406 int n,mmax,m,j,istep,i;
407 double wtemp,wr,wpr,wpi,wi,theta;
408 float tempr,tempi;
409
410 n=nn << 1;
411 j=1;
412 for (i=1;i<n;i+=2) {
413 if (j > i) {
414 SWAP(data[j],data[i]);
415 SWAP(data[j+1],data[i+1]);
416 }
417 m=n >> 1;
418 while (m >= 2 && j > m) {
419 j -= m;
420 m >>= 1;
421 }
422 j += m;
423 }
424 mmax=2;
425 while (n > mmax) {
426 istep=2*mmax;
427 theta=6.28318530717959/(isign*mmax);
428 wtemp=sin(0.5*theta);
429 wpr = -2.0*wtemp*wtemp;
430 wpi=sin(theta);
431 wr=1.0;
432 wi=0.0;
433 for (m=1;m<mmax;m+=2) {
434 for (i=m;i<=n;i+=istep) {
435 j=i+mmax;
436 tempr=wr*data[j]-wi*data[j+1];
437 tempi=wr*data[j+1]+wi*data[j];
438 data[j]=data[i]-tempr;
439 data[j+1]=data[i+1]-tempi;
440 data[i] += tempr;
441 data[i+1] += tempi;
442 }
443 wr=(wtemp=wr)*wpr-wi*wpi+wr;
444 wi=wi*wpr+wtemp*wpi+wi;
445 }
446 mmax=istep;
447 }
448 }
449