Changeset eb7f743 for src/mathutils.c
 Timestamp:
 Oct 2, 2009, 6:11:07 AM (10 years ago)
 Branches:
 feature/autosink, feature/constantq, feature/pitchshift, feature/pydocstrings, feature/timestretch, master, pitchshift, sampler, timestretch, yinfft+
 Children:
 352fd5f
 Parents:
 1498ced
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

src/mathutils.c
r1498ced reb7f743 26 26 #include "config.h" 27 27 28 fvec_t * new_aubio_window(uint_t size, aubio_window_type wintype) { 28 fvec_t * 29 new_aubio_window (uint_t size, aubio_window_type wintype) 30 { 29 31 // create fvec of size x 1 channel 30 32 fvec_t * win = new_fvec( size, 1); … … 79 81 } 80 82 81 smpl_t aubio_unwrap2pi(smpl_t phase) { 83 smpl_t 84 aubio_unwrap2pi (smpl_t phase) 85 { 82 86 /* mod(phase+pi,2pi)+pi */ 83 return phase + TWO_PI * (1. + FLOOR((phase+PI)/TWO_PI)); 84 } 85 86 smpl_t fvec_mean(fvec_t *s) { 87 uint_t i,j; 87 return phase + TWO_PI * (1. + FLOOR ((phase + PI) / TWO_PI)); 88 } 89 90 smpl_t 91 fvec_mean (fvec_t * s) 92 { 93 uint_t i, j; 88 94 smpl_t tmp = 0.0; 89 for (i =0; i < s>channels; i++)90 for (j =0; j < s>length; j++)95 for (i = 0; i < s>channels; i++) 96 for (j = 0; j < s>length; j++) 91 97 tmp += s>data[i][j]; 92 return tmp/(smpl_t)(s>length); 93 } 94 95 smpl_t fvec_sum(fvec_t *s) { 96 uint_t i,j; 98 return tmp / (smpl_t) (s>length); 99 } 100 101 smpl_t 102 fvec_sum (fvec_t * s) 103 { 104 uint_t i, j; 97 105 smpl_t tmp = 0.0; 98 for (i =0; i < s>channels; i++)99 for (j =0; j < s>length; j++)106 for (i = 0; i < s>channels; i++) { 107 for (j = 0; j < s>length; j++) { 100 108 tmp += s>data[i][j]; 109 } 110 } 101 111 return tmp; 102 112 } 103 113 104 smpl_t fvec_max(fvec_t *s) { 105 uint_t i,j; 114 smpl_t 115 fvec_max (fvec_t * s) 116 { 117 uint_t i, j; 106 118 smpl_t tmp = 0.0; 107 for (i=0; i < s>channels; i++) 108 for (j=0; j < s>length; j++) 109 tmp = (tmp > s>data[i][j])? tmp : s>data[i][j]; 119 for (i = 0; i < s>channels; i++) { 120 for (j = 0; j < s>length; j++) { 121 tmp = (tmp > s>data[i][j]) ? tmp : s>data[i][j]; 122 } 123 } 110 124 return tmp; 111 125 } 112 126 113 smpl_t fvec_min(fvec_t *s) { 114 uint_t i,j; 127 smpl_t 128 fvec_min (fvec_t * s) 129 { 130 uint_t i, j; 115 131 smpl_t tmp = s>data[0][0]; 116 for (i=0; i < s>channels; i++) 117 for (j=0; j < s>length; j++) 118 tmp = (tmp < s>data[i][j])? tmp : s>data[i][j] ; 132 for (i = 0; i < s>channels; i++) { 133 for (j = 0; j < s>length; j++) { 134 tmp = (tmp < s>data[i][j]) ? tmp : s>data[i][j]; 135 } 136 } 119 137 return tmp; 120 138 } 121 139 122 uint_t fvec_min_elem(fvec_t *s) { 123 uint_t i,j=0, pos=0.; 140 uint_t 141 fvec_min_elem (fvec_t * s) 142 { 143 uint_t i, j = 0, pos = 0.; 124 144 smpl_t tmp = s>data[0][0]; 125 for (i=0; i < s>channels; i++) 126 for (j=0; j < s>length; j++) { 127 pos = (tmp < s>data[i][j])? pos : j; 128 tmp = (tmp < s>data[i][j])? tmp : s>data[i][j] ; 129 } 145 for (i = 0; i < s>channels; i++) { 146 for (j = 0; j < s>length; j++) { 147 pos = (tmp < s>data[i][j]) ? pos : j; 148 tmp = (tmp < s>data[i][j]) ? tmp : s>data[i][j]; 149 } 150 } 130 151 return pos; 131 152 } 132 153 133 uint_t fvec_max_elem(fvec_t *s) { 134 uint_t i,j=0, pos=0.; 154 uint_t 155 fvec_max_elem (fvec_t * s) 156 { 157 uint_t i, j, pos; 135 158 smpl_t tmp = 0.0; 136 for (i=0; i < s>channels; i++) 137 for (j=0; j < s>length; j++) { 138 pos = (tmp > s>data[i][j])? pos : j; 139 tmp = (tmp > s>data[i][j])? tmp : s>data[i][j] ; 140 } 159 for (i = 0; i < s>channels; i++) { 160 for (j = 0; j < s>length; j++) { 161 pos = (tmp > s>data[i][j]) ? pos : j; 162 tmp = (tmp > s>data[i][j]) ? tmp : s>data[i][j]; 163 } 164 } 141 165 return pos; 142 166 } 143 167 144 void fvec_shift(fvec_t *s) { 145 uint_t i,j; 146 //smpl_t tmp = 0.0; 147 for (i=0; i < s>channels; i++) 148 for (j=0; j < s>length / 2 ; j++) { 149 //tmp = s>data[i][j]; 150 //s>data[i][j] = s>data[i][j+s>length/2]; 151 //s>data[i][j+s>length/2] = tmp; 152 ELEM_SWAP(s>data[i][j],s>data[i][j+s>length/2]); 153 } 154 } 155 156 smpl_t fvec_local_energy(fvec_t * f) { 157 smpl_t locE = 0.; 158 uint_t i,j; 159 for (i=0;i<f>channels;i++) 160 for (j=0;j<f>length;j++) 161 locE+=SQR(f>data[i][j]); 162 return locE; 163 } 164 165 smpl_t fvec_local_hfc(fvec_t * f) { 166 smpl_t locE = 0.; 167 uint_t i,j; 168 for (i=0;i<f>channels;i++) 169 for (j=0;j<f>length;j++) 170 locE+=(i+1)*f>data[i][j]; 171 return locE; 172 } 173 174 smpl_t fvec_alpha_norm(fvec_t * DF, smpl_t alpha) { 168 void 169 fvec_shift (fvec_t * s) 170 { 171 uint_t i, j; 172 for (i = 0; i < s>channels; i++) { 173 for (j = 0; j < s>length / 2; j++) { 174 ELEM_SWAP (s>data[i][j], s>data[i][j + s>length / 2]); 175 } 176 } 177 } 178 179 smpl_t 180 fvec_local_energy (fvec_t * f) 181 { 182 smpl_t energy = 0.; 183 uint_t i, j; 184 for (i = 0; i < f>channels; i++) { 185 for (j = 0; j < f>length; j++) { 186 energy += SQR (f>data[i][j]); 187 } 188 } 189 return energy; 190 } 191 192 smpl_t 193 fvec_local_hfc (fvec_t * v) 194 { 195 smpl_t hfc = 0.; 196 uint_t i, j; 197 for (i = 0; i < v>channels; i++) { 198 for (j = 0; j < v>length; j++) { 199 hfc += (i + 1) * v>data[i][j]; 200 } 201 } 202 return hfc; 203 } 204 205 void 206 fvec_min_removal (fvec_t * v) 207 { 208 smpl_t v_min = fvec_min (v); 209 fvec_add (v,  v_min ); 210 } 211 212 smpl_t 213 fvec_alpha_norm (fvec_t * o, smpl_t alpha) 214 { 215 uint_t i, j; 175 216 smpl_t tmp = 0.; 176 uint_t i,j;177 for (i=0;i<DF>channels;i++)178 for (j=0;j<DF>length;j++)179 tmp += POW(ABS(DF>data[i][j]),alpha);180 return POW(tmp/DF>length,1./alpha);181 }182 183 void184 fvec_min_removal (fvec_t * o)185 {186 uint_t i, j;187 smpl_t mini = fvec_min (o);188 217 for (i = 0; i < o>channels; i++) { 189 218 for (j = 0; j < o>length; j++) { 190 o>data[i][j] = mini; 191 } 192 } 193 } 194 195 void fvec_alpha_normalise(fvec_t * mag, uint_t alpha) { 196 smpl_t alphan = 1.; 197 uint_t length = mag>length, i=0, j; 198 alphan = fvec_alpha_norm(mag,alpha); 199 for (j=0;j<length;j++){ 200 mag>data[i][j] /= alphan; 201 } 202 } 203 204 void fvec_add(fvec_t * mag, smpl_t threshold) { 205 uint_t length = mag>length, i=0, j; 206 for (j=0;j<length;j++) { 207 mag>data[i][j] += threshold; 219 tmp += POW (ABS (o>data[i][j]), alpha); 220 } 221 } 222 return POW (tmp / o>length, 1. / alpha); 223 } 224 225 void 226 fvec_alpha_normalise (fvec_t * o, smpl_t alpha) 227 { 228 uint_t i, j; 229 smpl_t norm = fvec_alpha_norm (o, alpha); 230 for (i = 0; i < o>channels; i++) { 231 for (j = 0; j < o>length; j++) { 232 o>data[i][j] /= norm; 233 } 234 } 235 } 236 237 void 238 fvec_add (fvec_t * o, smpl_t val) 239 { 240 uint_t i, j; 241 for (i = 0; i < o>channels; i++) { 242 for (j = 0; j < o>length; j++) { 243 o>data[i][j] += val; 244 } 208 245 } 209 246 } … … 217 254 } 218 255 219 smpl_t fvec_moving_thres(fvec_t * vec, fvec_t * tmpvec, 220 uint_t post, uint_t pre, uint_t pos) { 221 smpl_t * medar = (smpl_t *)tmpvec>data[0]; 256 smpl_t 257 fvec_moving_thres (fvec_t * vec, fvec_t * tmpvec, 258 uint_t post, uint_t pre, uint_t pos) 259 { 260 smpl_t *medar = (smpl_t *) tmpvec>data[0]; 222 261 uint_t k; 223 uint_t win_length = post+pre+1;224 uint_t length = 262 uint_t win_length = post + pre + 1; 263 uint_t length = vec>length; 225 264 /* post part of the buffer does not exist */ 226 if (pos <post+1) {227 for (k =0;k<post+1pos;k++)228 medar[k] = 0.; /* 0padding at the beginning */229 for (k =post+1pos;k<win_length;k++)230 medar[k] = vec>data[0][k +pospost];231 /* the buffer is fully defined */232 } else if (pos +pre<length) {233 for (k =0;k<win_length;k++)234 medar[k] = vec>data[0][k +pospost];235 /* pre part of the buffer does not exist */265 if (pos < post + 1) { 266 for (k = 0; k < post + 1  pos; k++) 267 medar[k] = 0.; /* 0padding at the beginning */ 268 for (k = post + 1  pos; k < win_length; k++) 269 medar[k] = vec>data[0][k + pos  post]; 270 /* the buffer is fully defined */ 271 } else if (pos + pre < length) { 272 for (k = 0; k < win_length; k++) 273 medar[k] = vec>data[0][k + pos  post]; 274 /* pre part of the buffer does not exist */ 236 275 } else { 237 for (k =0;k<lengthpos+post;k++)238 medar[k] = vec>data[0][k +pospost];239 for (k =lengthpos+post;k<win_length;k++)240 medar[k] = 0.; /* 0padding at the end */241 } 242 return fvec_median (tmpvec);276 for (k = 0; k < length  pos + post; k++) 277 medar[k] = vec>data[0][k + pos  post]; 278 for (k = length  pos + post; k < win_length; k++) 279 medar[k] = 0.; /* 0padding at the end */ 280 } 281 return fvec_median (tmpvec); 243 282 } 244 283 … … 301 340 if (x2 == pos) return (x>data[0][pos] <= x>data[0][x0]) ? pos : x0; 302 341 s0 = x>data[0][x0]; 303 s1 = x>data[0][pos] 342 s1 = x>data[0][pos]; 304 343 s2 = x>data[0][x2]; 305 344 return pos + 0.5 * (s2  s0 ) / (s2  2.* s1 + s0); 306 }307 308 smpl_t aubio_quadfrac(smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf) {309 smpl_t tmp = s0 + (pf/2.) * (pf * ( s0  2.*s1 + s2 )  3.*s0 + 4.*s1  s2);310 return tmp;311 345 } 312 346 … … 320 354 } 321 355 322 smpl_t aubio_freqtomidi(smpl_t freq) { 356 smpl_t 357 aubio_quadfrac (smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf) 358 { 359 smpl_t tmp = 360 s0 + (pf / 2.) * (pf * (s0  2. * s1 + s2)  3. * s0 + 4. * s1  s2); 361 return tmp; 362 } 363 364 smpl_t 365 aubio_freqtomidi (smpl_t freq) 366 { 323 367 /* log(freq/A2)/log(2) */ 324 smpl_t midi = freq /6.875;325 midi = LOG (midi)/0.69314718055995;368 smpl_t midi = freq / 6.875; 369 midi = LOG (midi) / 0.69314718055995; 326 370 midi *= 12; 327 371 midi = 3; … … 329 373 } 330 374 331 smpl_t aubio_miditofreq(smpl_t midi) { 332 smpl_t freq = (midi+3.)/12.; 333 freq = EXP(freq*0.69314718055995); 375 smpl_t 376 aubio_miditofreq (smpl_t midi) 377 { 378 smpl_t freq = (midi + 3.) / 12.; 379 freq = EXP (freq * 0.69314718055995); 334 380 freq *= 6.875; 335 381 return freq; 336 382 } 337 383 338 smpl_t aubio_bintofreq(smpl_t bin, smpl_t samplerate, smpl_t fftsize) { 339 smpl_t freq = samplerate/fftsize; 340 return freq*bin; 341 } 342 343 smpl_t aubio_bintomidi(smpl_t bin, smpl_t samplerate, smpl_t fftsize) { 344 smpl_t midi = aubio_bintofreq(bin,samplerate,fftsize); 345 return aubio_freqtomidi(midi); 346 } 347 348 smpl_t aubio_freqtobin(smpl_t freq, smpl_t samplerate, smpl_t fftsize) { 349 smpl_t bin = fftsize/samplerate; 350 return freq*bin; 351 } 352 353 smpl_t aubio_miditobin(smpl_t midi, smpl_t samplerate, smpl_t fftsize) { 354 smpl_t freq = aubio_miditofreq(midi); 355 return aubio_freqtobin(freq,samplerate,fftsize); 356 } 357 358 /** returns 1 if wassilence is 0 and RMS(ibuf)<threshold 359 * \bug mono 360 */ 361 uint_t aubio_silence_detection(fvec_t * ibuf, smpl_t threshold) { 362 smpl_t loudness = 0; 363 uint_t i=0,j; 364 for (j=0;j<ibuf>length;j++) { 365 loudness += SQR(ibuf>data[i][j]); 366 } 367 loudness = SQRT(loudness); 368 loudness /= (smpl_t)ibuf>length; 369 loudness = LIN2DB(loudness); 370 371 return (loudness < threshold); 372 } 373 374 /** returns level log(RMS(ibuf)) if < threshold, 1 otherwise 375 * \bug mono 376 */ 377 smpl_t aubio_level_detection(fvec_t * ibuf, smpl_t threshold) { 378 smpl_t loudness = 0; 379 uint_t i=0,j; 380 for (j=0;j<ibuf>length;j++) { 381 loudness += SQR(ibuf>data[i][j]); 382 } 383 loudness = SQRT(loudness); 384 loudness /= (smpl_t)ibuf>length; 385 loudness = LIN2DB(loudness); 386 387 if (loudness < threshold) 384 smpl_t 385 aubio_bintofreq (smpl_t bin, smpl_t samplerate, smpl_t fftsize) 386 { 387 smpl_t freq = samplerate / fftsize; 388 return freq * bin; 389 } 390 391 smpl_t 392 aubio_bintomidi (smpl_t bin, smpl_t samplerate, smpl_t fftsize) 393 { 394 smpl_t midi = aubio_bintofreq (bin, samplerate, fftsize); 395 return aubio_freqtomidi (midi); 396 } 397 398 smpl_t 399 aubio_freqtobin (smpl_t freq, smpl_t samplerate, smpl_t fftsize) 400 { 401 smpl_t bin = fftsize / samplerate; 402 return freq * bin; 403 } 404 405 smpl_t 406 aubio_miditobin (smpl_t midi, smpl_t samplerate, smpl_t fftsize) 407 { 408 smpl_t freq = aubio_miditofreq (midi); 409 return aubio_freqtobin (freq, samplerate, fftsize); 410 } 411 412 smpl_t 413 aubio_db_spl (fvec_t * o) 414 { 415 smpl_t val = SQRT (fvec_local_energy (o)); 416 val /= (smpl_t) o>length; 417 return LIN2DB (val); 418 } 419 420 uint_t 421 aubio_silence_detection (fvec_t * o, smpl_t threshold) 422 { 423 return (aubio_db_spl (o) < threshold); 424 } 425 426 smpl_t 427 aubio_level_detection (fvec_t * o, smpl_t threshold) 428 { 429 smpl_t db_spl = aubio_db_spl (o); 430 if (db_spl < threshold) { 388 431 return 1.; 389 else 390 return loudness; 391 } 392 393 smpl_t aubio_zero_crossing_rate(fvec_t * input) { 394 uint_t i=0,j; 432 } else { 433 return db_spl; 434 } 435 } 436 437 smpl_t 438 aubio_zero_crossing_rate (fvec_t * input) 439 { 440 uint_t i = 0, j; 395 441 uint_t zcr = 0; 396 for ( j = 1; j < input>length; j++) {442 for (j = 1; j < input>length; j++) { 397 443 // previous was strictly negative 398 if ( input>data[i][j1] < 0.) {444 if (input>data[i][j  1] < 0.) { 399 445 // current is positive or null 400 if ( input>data[i][j] >= 0.) {446 if (input>data[i][j] >= 0.) { 401 447 zcr += 1; 402 448 } 403 // previous was positive or null449 // previous was positive or null 404 450 } else { 405 451 // current is strictly negative 406 if ( input>data[i][j] < 0.) {452 if (input>data[i][j] < 0.) { 407 453 zcr += 1; 408 454 } 409 455 } 410 456 } 411 return zcr/(smpl_t)input>length; 412 } 413 414 void aubio_autocorr(fvec_t * input, fvec_t * output) { 415 uint_t i = 0, j = 0, length = input>length; 416 smpl_t * data = input>data[0]; 417 smpl_t * acf = output>data[0]; 418 smpl_t tmp =0.; 419 for(i=0;i<length;i++){ 420 for(j=i;j<length;j++){ 421 tmp += data[ji]*data[j]; 422 } 423 acf[i] = tmp /(smpl_t)(lengthi); 424 tmp = 0.0; 425 } 426 } 427 428 void aubio_cleanup(void) { 457 return zcr / (smpl_t) input>length; 458 } 459 460 void 461 aubio_autocorr (fvec_t * input, fvec_t * output) 462 { 463 uint_t i, j, k, length = input>length; 464 smpl_t *data, *acf; 465 smpl_t tmp = 0; 466 for (k = 0; k < input>channels; k++) { 467 data = input>data[k]; 468 acf = output>data[k]; 469 for (i = 0; i < length; i++) { 470 tmp = 0.; 471 for (j = i; j < length; j++) { 472 tmp += data[j  i] * data[j]; 473 } 474 acf[i] = tmp / (smpl_t) (length  i); 475 } 476 } 477 } 478 479 void 480 aubio_cleanup (void) 481 { 429 482 #if HAVE_FFTW3 430 fftw_cleanup ();483 fftw_cleanup (); 431 484 #else 432 485 #if HAVE_FFTW3F 433 fftwf_cleanup ();486 fftwf_cleanup (); 434 487 #endif 435 488 #endif
Note: See TracChangeset
for help on using the changeset viewer.