vorbis_psy.c (36160B)
1 /******************************************************************** 2 * * 3 * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * 4 * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * 5 * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * 6 * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * 7 * * 8 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2010 * 9 * by the Xiph.Org Foundation https://xiph.org/ * 10 * * 11 ******************************************************************** 12 13 function: psychoacoustics not including preecho 14 15 ********************************************************************/ 16 17 #include <stdlib.h> 18 #include <math.h> 19 #include <string.h> 20 #include "vorbis/codec.h" 21 #include "codec_internal.h" 22 23 #include "masking.h" 24 #include "psy.h" 25 #include "os.h" 26 #include "lpc.h" 27 #include "smallft.h" 28 #include "scales.h" 29 #include "misc.h" 30 31 #define NEGINF -9999.f 32 static const double stereo_threshholds[]={0.0, .5, 1.0, 1.5, 2.5, 4.5, 8.5, 16.5, 9e10}; 33 static const double stereo_threshholds_limited[]={0.0, .5, 1.0, 1.5, 2.0, 2.5, 4.5, 8.5, 9e10}; 34 35 vorbis_look_psy_global *_vp_global_look(vorbis_info *vi){ 36 codec_setup_info *ci=vi->codec_setup; 37 vorbis_info_psy_global *gi=&ci->psy_g_param; 38 vorbis_look_psy_global *look=_ogg_calloc(1,sizeof(*look)); 39 40 look->channels=vi->channels; 41 42 look->ampmax=-9999.; 43 look->gi=gi; 44 return(look); 45 } 46 47 void _vp_global_free(vorbis_look_psy_global *look){ 48 if(look){ 49 memset(look,0,sizeof(*look)); 50 _ogg_free(look); 51 } 52 } 53 54 void _vi_gpsy_free(vorbis_info_psy_global *i){ 55 if(i){ 56 memset(i,0,sizeof(*i)); 57 _ogg_free(i); 58 } 59 } 60 61 void _vi_psy_free(vorbis_info_psy *i){ 62 if(i){ 63 memset(i,0,sizeof(*i)); 64 _ogg_free(i); 65 } 66 } 67 68 static void min_curve(float *c, 69 float *c2){ 70 int i; 71 for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i]; 72 } 73 static void max_curve(float *c, 74 float *c2){ 75 int i; 76 for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i]; 77 } 78 79 static void attenuate_curve(float *c,float att){ 80 int i; 81 for(i=0;i<EHMER_MAX;i++) 82 c[i]+=att; 83 } 84 85 static float ***setup_tone_curves(float curveatt_dB[P_BANDS],float binHz,int n, 86 float center_boost, float center_decay_rate){ 87 int i,j,k,m; 88 float ath[EHMER_MAX]; 89 float workc[P_BANDS][P_LEVELS][EHMER_MAX]; 90 float athc[P_LEVELS][EHMER_MAX]; 91 float *brute_buffer=alloca(n*sizeof(*brute_buffer)); 92 93 float ***ret=_ogg_malloc(sizeof(*ret)*P_BANDS); 94 95 memset(workc,0,sizeof(workc)); 96 97 for(i=0;i<P_BANDS;i++){ 98 /* we add back in the ATH to avoid low level curves falling off to 99 -infinity and unnecessarily cutting off high level curves in the 100 curve limiting (last step). */ 101 102 /* A half-band's settings must be valid over the whole band, and 103 it's better to mask too little than too much */ 104 int ath_offset=i*4; 105 for(j=0;j<EHMER_MAX;j++){ 106 float min=999.; 107 for(k=0;k<4;k++) 108 if(j+k+ath_offset<MAX_ATH){ 109 if(min>ATH[j+k+ath_offset])min=ATH[j+k+ath_offset]; 110 }else{ 111 if(min>ATH[MAX_ATH-1])min=ATH[MAX_ATH-1]; 112 } 113 ath[j]=min; 114 } 115 116 /* copy curves into working space, replicate the 50dB curve to 30 117 and 40, replicate the 100dB curve to 110 */ 118 for(j=0;j<6;j++) 119 memcpy(workc[i][j+2],tonemasks[i][j],EHMER_MAX*sizeof(*tonemasks[i][j])); 120 memcpy(workc[i][0],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0])); 121 memcpy(workc[i][1],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0])); 122 123 /* apply centered curve boost/decay */ 124 for(j=0;j<P_LEVELS;j++){ 125 for(k=0;k<EHMER_MAX;k++){ 126 float adj=center_boost+abs(EHMER_OFFSET-k)*center_decay_rate; 127 if(adj<0. && center_boost>0)adj=0.; 128 if(adj>0. && center_boost<0)adj=0.; 129 workc[i][j][k]+=adj; 130 } 131 } 132 133 /* normalize curves so the driving amplitude is 0dB */ 134 /* make temp curves with the ATH overlayed */ 135 for(j=0;j<P_LEVELS;j++){ 136 attenuate_curve(workc[i][j],curveatt_dB[i]+100.-(j<2?2:j)*10.-P_LEVEL_0); 137 memcpy(athc[j],ath,EHMER_MAX*sizeof(**athc)); 138 attenuate_curve(athc[j],+100.-j*10.f-P_LEVEL_0); 139 max_curve(athc[j],workc[i][j]); 140 } 141 142 /* Now limit the louder curves. 143 144 the idea is this: We don't know what the playback attenuation 145 will be; 0dB SL moves every time the user twiddles the volume 146 knob. So that means we have to use a single 'most pessimal' curve 147 for all masking amplitudes, right? Wrong. The *loudest* sound 148 can be in (we assume) a range of ...+100dB] SL. However, sounds 149 20dB down will be in a range ...+80], 40dB down is from ...+60], 150 etc... */ 151 152 for(j=1;j<P_LEVELS;j++){ 153 min_curve(athc[j],athc[j-1]); 154 min_curve(workc[i][j],athc[j]); 155 } 156 } 157 158 for(i=0;i<P_BANDS;i++){ 159 int hi_curve,lo_curve,bin; 160 ret[i]=_ogg_malloc(sizeof(**ret)*P_LEVELS); 161 162 /* low frequency curves are measured with greater resolution than 163 the MDCT/FFT will actually give us; we want the curve applied 164 to the tone data to be pessimistic and thus apply the minimum 165 masking possible for a given bin. That means that a single bin 166 could span more than one octave and that the curve will be a 167 composite of multiple octaves. It also may mean that a single 168 bin may span > an eighth of an octave and that the eighth 169 octave values may also be composited. */ 170 171 /* which octave curves will we be compositing? */ 172 bin=floor(fromOC(i*.5)/binHz); 173 lo_curve= ceil(toOC(bin*binHz+1)*2); 174 hi_curve= floor(toOC((bin+1)*binHz)*2); 175 if(lo_curve>i)lo_curve=i; 176 if(lo_curve<0)lo_curve=0; 177 if(hi_curve>=P_BANDS)hi_curve=P_BANDS-1; 178 179 for(m=0;m<P_LEVELS;m++){ 180 ret[i][m]=_ogg_malloc(sizeof(***ret)*(EHMER_MAX+2)); 181 182 for(j=0;j<n;j++)brute_buffer[j]=999.; 183 184 /* render the curve into bins, then pull values back into curve. 185 The point is that any inherent subsampling aliasing results in 186 a safe minimum */ 187 for(k=lo_curve;k<=hi_curve;k++){ 188 int l=0; 189 190 for(j=0;j<EHMER_MAX;j++){ 191 int lo_bin= fromOC(j*.125+k*.5-2.0625)/binHz; 192 int hi_bin= fromOC(j*.125+k*.5-1.9375)/binHz+1; 193 194 if(lo_bin<0)lo_bin=0; 195 if(lo_bin>n)lo_bin=n; 196 if(lo_bin<l)l=lo_bin; 197 if(hi_bin<0)hi_bin=0; 198 if(hi_bin>n)hi_bin=n; 199 200 for(;l<hi_bin && l<n;l++) 201 if(brute_buffer[l]>workc[k][m][j]) 202 brute_buffer[l]=workc[k][m][j]; 203 } 204 205 for(;l<n;l++) 206 if(brute_buffer[l]>workc[k][m][EHMER_MAX-1]) 207 brute_buffer[l]=workc[k][m][EHMER_MAX-1]; 208 209 } 210 211 /* be equally paranoid about being valid up to next half ocatve */ 212 if(i+1<P_BANDS){ 213 int l=0; 214 k=i+1; 215 for(j=0;j<EHMER_MAX;j++){ 216 int lo_bin= fromOC(j*.125+i*.5-2.0625)/binHz; 217 int hi_bin= fromOC(j*.125+i*.5-1.9375)/binHz+1; 218 219 if(lo_bin<0)lo_bin=0; 220 if(lo_bin>n)lo_bin=n; 221 if(lo_bin<l)l=lo_bin; 222 if(hi_bin<0)hi_bin=0; 223 if(hi_bin>n)hi_bin=n; 224 225 for(;l<hi_bin && l<n;l++) 226 if(brute_buffer[l]>workc[k][m][j]) 227 brute_buffer[l]=workc[k][m][j]; 228 } 229 230 for(;l<n;l++) 231 if(brute_buffer[l]>workc[k][m][EHMER_MAX-1]) 232 brute_buffer[l]=workc[k][m][EHMER_MAX-1]; 233 234 } 235 236 237 for(j=0;j<EHMER_MAX;j++){ 238 int bin=fromOC(j*.125+i*.5-2.)/binHz; 239 if(bin<0){ 240 ret[i][m][j+2]=-999.; 241 }else{ 242 if(bin>=n){ 243 ret[i][m][j+2]=-999.; 244 }else{ 245 ret[i][m][j+2]=brute_buffer[bin]; 246 } 247 } 248 } 249 250 /* add fenceposts */ 251 for(j=0;j<EHMER_OFFSET;j++) 252 if(ret[i][m][j+2]>-200.f)break; 253 ret[i][m][0]=j; 254 255 for(j=EHMER_MAX-1;j>EHMER_OFFSET+1;j--) 256 if(ret[i][m][j+2]>-200.f) 257 break; 258 ret[i][m][1]=j; 259 260 } 261 } 262 263 return(ret); 264 } 265 266 void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi, 267 vorbis_info_psy_global *gi,int n,long rate){ 268 long i,j,lo=-99,hi=1; 269 long maxoc; 270 memset(p,0,sizeof(*p)); 271 272 p->eighth_octave_lines=gi->eighth_octave_lines; 273 p->shiftoc=rint(log(gi->eighth_octave_lines*8.f)/log(2.f))-1; 274 275 p->firstoc=toOC(.25f*rate*.5/n)*(1<<(p->shiftoc+1))-gi->eighth_octave_lines; 276 maxoc=toOC((n+.25f)*rate*.5/n)*(1<<(p->shiftoc+1))+.5f; 277 p->total_octave_lines=maxoc-p->firstoc+1; 278 p->ath=_ogg_malloc(n*sizeof(*p->ath)); 279 280 p->octave=_ogg_malloc(n*sizeof(*p->octave)); 281 p->bark=_ogg_malloc(n*sizeof(*p->bark)); 282 p->vi=vi; 283 p->n=n; 284 p->rate=rate; 285 286 /* AoTuV HF weighting */ 287 p->m_val = 1.; 288 if(rate < 26000) p->m_val = 0; 289 else if(rate < 38000) p->m_val = .94; /* 32kHz */ 290 else if(rate > 46000) p->m_val = 1.275; /* 48kHz */ 291 292 /* set up the lookups for a given blocksize and sample rate */ 293 294 for(i=0,j=0;i<MAX_ATH-1;i++){ 295 int endpos=rint(fromOC((i+1)*.125-2.)*2*n/rate); 296 float base=ATH[i]; 297 if(j<endpos){ 298 float delta=(ATH[i+1]-base)/(endpos-j); 299 for(;j<endpos && j<n;j++){ 300 p->ath[j]=base+100.; 301 base+=delta; 302 } 303 } 304 } 305 306 for(;j<n;j++){ 307 p->ath[j]=p->ath[j-1]; 308 } 309 310 for(i=0;i<n;i++){ 311 float bark=toBARK(rate/(2*n)*i); 312 313 for(;lo+vi->noisewindowlomin<i && 314 toBARK(rate/(2*n)*lo)<(bark-vi->noisewindowlo);lo++); 315 316 for(;hi<=n && (hi<i+vi->noisewindowhimin || 317 toBARK(rate/(2*n)*hi)<(bark+vi->noisewindowhi));hi++); 318 319 p->bark[i]=((lo-1)<<16)+(hi-1); 320 321 } 322 323 for(i=0;i<n;i++) 324 p->octave[i]=toOC((i+.25f)*.5*rate/n)*(1<<(p->shiftoc+1))+.5f; 325 326 p->tonecurves=setup_tone_curves(vi->toneatt,rate*.5/n,n, 327 vi->tone_centerboost,vi->tone_decay); 328 329 /* set up rolling noise median */ 330 p->noiseoffset=_ogg_malloc(P_NOISECURVES*sizeof(*p->noiseoffset)); 331 for(i=0;i<P_NOISECURVES;i++) 332 p->noiseoffset[i]=_ogg_malloc(n*sizeof(**p->noiseoffset)); 333 334 for(i=0;i<n;i++){ 335 float halfoc=toOC((i+.5)*rate/(2.*n))*2.; 336 int inthalfoc; 337 float del; 338 339 if(halfoc<0)halfoc=0; 340 if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1; 341 inthalfoc=(int)halfoc; 342 /*If we hit the P_BANDS-1 clamp above, inthalfoc+1 will be out of bounds, 343 even though it will have an interpolation weight of 0. 344 Shift the interval so we don't read past the end of the array.*/ 345 if(inthalfoc>=P_BANDS-2)inthalfoc=P_BANDS-2; 346 del=halfoc-inthalfoc; 347 348 for(j=0;j<P_NOISECURVES;j++) 349 p->noiseoffset[j][i]= 350 p->vi->noiseoff[j][inthalfoc]*(1.-del) + 351 p->vi->noiseoff[j][inthalfoc+1]*del; 352 353 } 354 #if 0 355 { 356 static int ls=0; 357 _analysis_output_always("noiseoff0",ls,p->noiseoffset[0],n,1,0,0); 358 _analysis_output_always("noiseoff1",ls,p->noiseoffset[1],n,1,0,0); 359 _analysis_output_always("noiseoff2",ls++,p->noiseoffset[2],n,1,0,0); 360 } 361 #endif 362 } 363 364 void _vp_psy_clear(vorbis_look_psy *p){ 365 int i,j; 366 if(p){ 367 if(p->ath)_ogg_free(p->ath); 368 if(p->octave)_ogg_free(p->octave); 369 if(p->bark)_ogg_free(p->bark); 370 if(p->tonecurves){ 371 for(i=0;i<P_BANDS;i++){ 372 for(j=0;j<P_LEVELS;j++){ 373 _ogg_free(p->tonecurves[i][j]); 374 } 375 _ogg_free(p->tonecurves[i]); 376 } 377 _ogg_free(p->tonecurves); 378 } 379 if(p->noiseoffset){ 380 for(i=0;i<P_NOISECURVES;i++){ 381 _ogg_free(p->noiseoffset[i]); 382 } 383 _ogg_free(p->noiseoffset); 384 } 385 memset(p,0,sizeof(*p)); 386 } 387 } 388 389 /* octave/(8*eighth_octave_lines) x scale and dB y scale */ 390 static void seed_curve(float *seed, 391 const float **curves, 392 float amp, 393 int oc, int n, 394 int linesper,float dBoffset){ 395 int i,post1; 396 int seedptr; 397 const float *posts,*curve; 398 399 int choice=(int)((amp+dBoffset-P_LEVEL_0)*.1f); 400 choice=max(choice,0); 401 choice=min(choice,P_LEVELS-1); 402 posts=curves[choice]; 403 curve=posts+2; 404 post1=(int)posts[1]; 405 seedptr=oc+(posts[0]-EHMER_OFFSET)*linesper-(linesper>>1); 406 407 for(i=posts[0];i<post1;i++){ 408 if(seedptr>0){ 409 float lin=amp+curve[i]; 410 if(seed[seedptr]<lin)seed[seedptr]=lin; 411 } 412 seedptr+=linesper; 413 if(seedptr>=n)break; 414 } 415 } 416 417 static void seed_loop(vorbis_look_psy *p, 418 const float ***curves, 419 const float *f, 420 const float *flr, 421 float *seed, 422 float specmax){ 423 vorbis_info_psy *vi=p->vi; 424 long n=p->n,i; 425 float dBoffset=vi->max_curve_dB-specmax; 426 427 /* prime the working vector with peak values */ 428 429 for(i=0;i<n;i++){ 430 float max=f[i]; 431 long oc=p->octave[i]; 432 while(i+1<n && p->octave[i+1]==oc){ 433 i++; 434 if(f[i]>max)max=f[i]; 435 } 436 437 if(max+6.f>flr[i]){ 438 oc=oc>>p->shiftoc; 439 440 if(oc>=P_BANDS)oc=P_BANDS-1; 441 if(oc<0)oc=0; 442 443 seed_curve(seed, 444 curves[oc], 445 max, 446 p->octave[i]-p->firstoc, 447 p->total_octave_lines, 448 p->eighth_octave_lines, 449 dBoffset); 450 } 451 } 452 } 453 454 static void seed_chase(float *seeds, int linesper, long n){ 455 long *posstack=alloca(n*sizeof(*posstack)); 456 float *ampstack=alloca(n*sizeof(*ampstack)); 457 long stack=0; 458 long pos=0; 459 long i; 460 461 for(i=0;i<n;i++){ 462 if(stack<2){ 463 posstack[stack]=i; 464 ampstack[stack++]=seeds[i]; 465 }else{ 466 while(1){ 467 if(seeds[i]<ampstack[stack-1]){ 468 posstack[stack]=i; 469 ampstack[stack++]=seeds[i]; 470 break; 471 }else{ 472 if(i<posstack[stack-1]+linesper){ 473 if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] && 474 i<posstack[stack-2]+linesper){ 475 /* we completely overlap, making stack-1 irrelevant. pop it */ 476 stack--; 477 continue; 478 } 479 } 480 posstack[stack]=i; 481 ampstack[stack++]=seeds[i]; 482 break; 483 484 } 485 } 486 } 487 } 488 489 /* the stack now contains only the positions that are relevant. Scan 490 'em straight through */ 491 492 for(i=0;i<stack;i++){ 493 long endpos; 494 if(i<stack-1 && ampstack[i+1]>ampstack[i]){ 495 endpos=posstack[i+1]; 496 }else{ 497 endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is 498 discarded in short frames */ 499 } 500 if(endpos>n)endpos=n; 501 for(;pos<endpos;pos++) 502 seeds[pos]=ampstack[i]; 503 } 504 505 /* there. Linear time. I now remember this was on a problem set I 506 had in Grad Skool... I didn't solve it at the time ;-) */ 507 508 } 509 510 /* bleaugh, this is more complicated than it needs to be */ 511 #include<stdio.h> 512 static void max_seeds(vorbis_look_psy *p, 513 float *seed, 514 float *flr){ 515 long n=p->total_octave_lines; 516 int linesper=p->eighth_octave_lines; 517 long linpos=0; 518 long pos; 519 520 seed_chase(seed,linesper,n); /* for masking */ 521 522 pos=p->octave[0]-p->firstoc-(linesper>>1); 523 524 while(linpos+1<p->n){ 525 float minV=seed[pos]; 526 long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc; 527 if(minV>p->vi->tone_abs_limit)minV=p->vi->tone_abs_limit; 528 while(pos+1<=end){ 529 pos++; 530 if((seed[pos]>NEGINF && seed[pos]<minV) || minV==NEGINF) 531 minV=seed[pos]; 532 } 533 534 end=pos+p->firstoc; 535 for(;linpos<p->n && p->octave[linpos]<=end;linpos++) 536 if(flr[linpos]<minV)flr[linpos]=minV; 537 } 538 539 { 540 float minV=seed[p->total_octave_lines-1]; 541 for(;linpos<p->n;linpos++) 542 if(flr[linpos]<minV)flr[linpos]=minV; 543 } 544 545 } 546 547 static void bark_noise_hybridmp(int n,const long *b, 548 const float *f, 549 float *noise, 550 const float offset, 551 const int fixed){ 552 553 float *N=alloca(n*sizeof(*N)); 554 float *X=alloca(n*sizeof(*N)); 555 float *XX=alloca(n*sizeof(*N)); 556 float *Y=alloca(n*sizeof(*N)); 557 float *XY=alloca(n*sizeof(*N)); 558 559 float tN, tX, tXX, tY, tXY; 560 int i; 561 562 int lo, hi; 563 float R=0.f; 564 float A=0.f; 565 float B=0.f; 566 float D=1.f; 567 float w, x, y; 568 569 tN = tX = tXX = tY = tXY = 0.f; 570 571 y = f[0] + offset; 572 if (y < 1.f) y = 1.f; 573 574 w = y * y * .5; 575 576 tN += w; 577 tX += w; 578 tY += w * y; 579 580 N[0] = tN; 581 X[0] = tX; 582 XX[0] = tXX; 583 Y[0] = tY; 584 XY[0] = tXY; 585 586 for (i = 1, x = 1.f; i < n; i++, x += 1.f) { 587 588 y = f[i] + offset; 589 if (y < 1.f) y = 1.f; 590 591 w = y * y; 592 593 tN += w; 594 tX += w * x; 595 tXX += w * x * x; 596 tY += w * y; 597 tXY += w * x * y; 598 599 N[i] = tN; 600 X[i] = tX; 601 XX[i] = tXX; 602 Y[i] = tY; 603 XY[i] = tXY; 604 } 605 606 for (i = 0, x = 0.f; i < n; i++, x += 1.f) { 607 608 lo = b[i] >> 16; 609 hi = b[i] & 0xffff; 610 if( lo>=0 || -lo>=n ) break; 611 if( hi>=n ) break; 612 613 tN = N[hi] + N[-lo]; 614 tX = X[hi] - X[-lo]; 615 tXX = XX[hi] + XX[-lo]; 616 tY = Y[hi] + Y[-lo]; 617 tXY = XY[hi] - XY[-lo]; 618 619 A = tY * tXX - tX * tXY; 620 B = tN * tXY - tX * tY; 621 D = tN * tXX - tX * tX; 622 R = (A + x * B) / D; 623 if (R < 0.f) R = 0.f; 624 625 noise[i] = R - offset; 626 } 627 628 for ( ; i < n; i++, x += 1.f) { 629 630 lo = b[i] >> 16; 631 hi = b[i] & 0xffff; 632 if( lo<0 || lo>=n ) break; 633 if( hi>=n ) break; 634 635 tN = N[hi] - N[lo]; 636 tX = X[hi] - X[lo]; 637 tXX = XX[hi] - XX[lo]; 638 tY = Y[hi] - Y[lo]; 639 tXY = XY[hi] - XY[lo]; 640 641 A = tY * tXX - tX * tXY; 642 B = tN * tXY - tX * tY; 643 D = tN * tXX - tX * tX; 644 R = (A + x * B) / D; 645 if (R < 0.f) R = 0.f; 646 647 noise[i] = R - offset; 648 } 649 650 for ( ; i < n; i++, x += 1.f) { 651 652 R = (A + x * B) / D; 653 if (R < 0.f) R = 0.f; 654 655 noise[i] = R - offset; 656 } 657 658 if (fixed <= 0) return; 659 660 for (i = 0, x = 0.f; i < n; i++, x += 1.f) { 661 hi = i + fixed / 2; 662 lo = hi - fixed; 663 if ( hi>=n ) break; 664 if ( lo>=0 ) break; 665 666 tN = N[hi] + N[-lo]; 667 tX = X[hi] - X[-lo]; 668 tXX = XX[hi] + XX[-lo]; 669 tY = Y[hi] + Y[-lo]; 670 tXY = XY[hi] - XY[-lo]; 671 672 673 A = tY * tXX - tX * tXY; 674 B = tN * tXY - tX * tY; 675 D = tN * tXX - tX * tX; 676 R = (A + x * B) / D; 677 678 if (R - offset < noise[i]) noise[i] = R - offset; 679 } 680 for ( ; i < n; i++, x += 1.f) { 681 682 hi = i + fixed / 2; 683 lo = hi - fixed; 684 if ( hi>=n ) break; 685 if ( lo<0 ) break; 686 687 tN = N[hi] - N[lo]; 688 tX = X[hi] - X[lo]; 689 tXX = XX[hi] - XX[lo]; 690 tY = Y[hi] - Y[lo]; 691 tXY = XY[hi] - XY[lo]; 692 693 A = tY * tXX - tX * tXY; 694 B = tN * tXY - tX * tY; 695 D = tN * tXX - tX * tX; 696 R = (A + x * B) / D; 697 698 if (R - offset < noise[i]) noise[i] = R - offset; 699 } 700 for ( ; i < n; i++, x += 1.f) { 701 R = (A + x * B) / D; 702 if (R - offset < noise[i]) noise[i] = R - offset; 703 } 704 } 705 706 void _vp_noisemask(vorbis_look_psy *p, 707 float *logmdct, 708 float *logmask){ 709 710 int i,n=p->n; 711 float *work=alloca(n*sizeof(*work)); 712 713 bark_noise_hybridmp(n,p->bark,logmdct,logmask, 714 140.,-1); 715 716 for(i=0;i<n;i++)work[i]=logmdct[i]-logmask[i]; 717 718 bark_noise_hybridmp(n,p->bark,work,logmask,0., 719 p->vi->noisewindowfixed); 720 721 for(i=0;i<n;i++)work[i]=logmdct[i]-work[i]; 722 723 #if 0 724 { 725 static int seq=0; 726 727 float work2[n]; 728 for(i=0;i<n;i++){ 729 work2[i]=logmask[i]+work[i]; 730 } 731 732 if(seq&1) 733 _analysis_output("median2R",seq/2,work,n,1,0,0); 734 else 735 _analysis_output("median2L",seq/2,work,n,1,0,0); 736 737 if(seq&1) 738 _analysis_output("envelope2R",seq/2,work2,n,1,0,0); 739 else 740 _analysis_output("envelope2L",seq/2,work2,n,1,0,0); 741 seq++; 742 } 743 #endif 744 745 for(i=0;i<n;i++){ 746 int dB=logmask[i]+.5; 747 if(dB>=NOISE_COMPAND_LEVELS)dB=NOISE_COMPAND_LEVELS-1; 748 if(dB<0)dB=0; 749 logmask[i]= work[i]+p->vi->noisecompand[dB]; 750 } 751 752 } 753 754 void _vp_tonemask(vorbis_look_psy *p, 755 float *logfft, 756 float *logmask, 757 float global_specmax, 758 float local_specmax){ 759 760 int i,n=p->n; 761 762 float *seed=alloca(sizeof(*seed)*p->total_octave_lines); 763 float att=local_specmax+p->vi->ath_adjatt; 764 for(i=0;i<p->total_octave_lines;i++)seed[i]=NEGINF; 765 766 /* set the ATH (floating below localmax, not global max by a 767 specified att) */ 768 if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt; 769 770 for(i=0;i<n;i++) 771 logmask[i]=p->ath[i]+att; 772 773 /* tone masking */ 774 seed_loop(p,(const float ***)p->tonecurves,logfft,logmask,seed,global_specmax); 775 max_seeds(p,seed,logmask); 776 777 } 778 779 void _vp_offset_and_mix(vorbis_look_psy *p, 780 float *noise, 781 float *tone, 782 int offset_select, 783 float *logmask, 784 float *mdct, 785 float *logmdct){ 786 int i,n=p->n; 787 float de, coeffi, cx;/* AoTuV */ 788 float toneatt=p->vi->tone_masteratt[offset_select]; 789 790 cx = p->m_val; 791 792 for(i=0;i<n;i++){ 793 float val= noise[i]+p->noiseoffset[offset_select][i]; 794 if(val>p->vi->noisemaxsupp)val=p->vi->noisemaxsupp; 795 logmask[i]=max(val,tone[i]+toneatt); 796 797 798 /* AoTuV */ 799 /** @ M1 ** 800 The following codes improve a noise problem. 801 A fundamental idea uses the value of masking and carries out 802 the relative compensation of the MDCT. 803 However, this code is not perfect and all noise problems cannot be solved. 804 by Aoyumi @ 2004/04/18 805 */ 806 807 if(offset_select == 1) { 808 coeffi = -17.2; /* coeffi is a -17.2dB threshold */ 809 val = val - logmdct[i]; /* val == mdct line value relative to floor in dB */ 810 811 if(val > coeffi){ 812 /* mdct value is > -17.2 dB below floor */ 813 814 de = 1.0-((val-coeffi)*0.005*cx); 815 /* pro-rated attenuation: 816 -0.00 dB boost if mdct value is -17.2dB (relative to floor) 817 -0.77 dB boost if mdct value is 0dB (relative to floor) 818 -1.64 dB boost if mdct value is +17.2dB (relative to floor) 819 etc... */ 820 821 if(de < 0) de = 0.0001; 822 }else 823 /* mdct value is <= -17.2 dB below floor */ 824 825 de = 1.0-((val-coeffi)*0.0003*cx); 826 /* pro-rated attenuation: 827 +0.00 dB atten if mdct value is -17.2dB (relative to floor) 828 +0.45 dB atten if mdct value is -34.4dB (relative to floor) 829 etc... */ 830 831 mdct[i] *= de; 832 833 } 834 } 835 } 836 837 float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){ 838 vorbis_info *vi=vd->vi; 839 codec_setup_info *ci=vi->codec_setup; 840 vorbis_info_psy_global *gi=&ci->psy_g_param; 841 842 int n=ci->blocksizes[vd->W]/2; 843 float secs=(float)n/vi->rate; 844 845 amp+=secs*gi->ampmax_att_per_sec; 846 if(amp<-9999)amp=-9999; 847 return(amp); 848 } 849 850 static float FLOOR1_fromdB_LOOKUP[256]={ 851 1.0649863e-07F, 1.1341951e-07F, 1.2079015e-07F, 1.2863978e-07F, 852 1.3699951e-07F, 1.4590251e-07F, 1.5538408e-07F, 1.6548181e-07F, 853 1.7623575e-07F, 1.8768855e-07F, 1.9988561e-07F, 2.128753e-07F, 854 2.2670913e-07F, 2.4144197e-07F, 2.5713223e-07F, 2.7384213e-07F, 855 2.9163793e-07F, 3.1059021e-07F, 3.3077411e-07F, 3.5226968e-07F, 856 3.7516214e-07F, 3.9954229e-07F, 4.2550680e-07F, 4.5315863e-07F, 857 4.8260743e-07F, 5.1396998e-07F, 5.4737065e-07F, 5.8294187e-07F, 858 6.2082472e-07F, 6.6116941e-07F, 7.0413592e-07F, 7.4989464e-07F, 859 7.9862701e-07F, 8.5052630e-07F, 9.0579828e-07F, 9.6466216e-07F, 860 1.0273513e-06F, 1.0941144e-06F, 1.1652161e-06F, 1.2409384e-06F, 861 1.3215816e-06F, 1.4074654e-06F, 1.4989305e-06F, 1.5963394e-06F, 862 1.7000785e-06F, 1.8105592e-06F, 1.9282195e-06F, 2.0535261e-06F, 863 2.1869758e-06F, 2.3290978e-06F, 2.4804557e-06F, 2.6416497e-06F, 864 2.8133190e-06F, 2.9961443e-06F, 3.1908506e-06F, 3.3982101e-06F, 865 3.6190449e-06F, 3.8542308e-06F, 4.1047004e-06F, 4.3714470e-06F, 866 4.6555282e-06F, 4.9580707e-06F, 5.2802740e-06F, 5.6234160e-06F, 867 5.9888572e-06F, 6.3780469e-06F, 6.7925283e-06F, 7.2339451e-06F, 868 7.7040476e-06F, 8.2047000e-06F, 8.7378876e-06F, 9.3057248e-06F, 869 9.9104632e-06F, 1.0554501e-05F, 1.1240392e-05F, 1.1970856e-05F, 870 1.2748789e-05F, 1.3577278e-05F, 1.4459606e-05F, 1.5399272e-05F, 871 1.6400004e-05F, 1.7465768e-05F, 1.8600792e-05F, 1.9809576e-05F, 872 2.1096914e-05F, 2.2467911e-05F, 2.3928002e-05F, 2.5482978e-05F, 873 2.7139006e-05F, 2.8902651e-05F, 3.0780908e-05F, 3.2781225e-05F, 874 3.4911534e-05F, 3.7180282e-05F, 3.9596466e-05F, 4.2169667e-05F, 875 4.4910090e-05F, 4.7828601e-05F, 5.0936773e-05F, 5.4246931e-05F, 876 5.7772202e-05F, 6.1526565e-05F, 6.5524908e-05F, 6.9783085e-05F, 877 7.4317983e-05F, 7.9147585e-05F, 8.4291040e-05F, 8.9768747e-05F, 878 9.5602426e-05F, 0.00010181521F, 0.00010843174F, 0.00011547824F, 879 0.00012298267F, 0.00013097477F, 0.00013948625F, 0.00014855085F, 880 0.00015820453F, 0.00016848555F, 0.00017943469F, 0.00019109536F, 881 0.00020351382F, 0.00021673929F, 0.00023082423F, 0.00024582449F, 882 0.00026179955F, 0.00027881276F, 0.00029693158F, 0.00031622787F, 883 0.00033677814F, 0.00035866388F, 0.00038197188F, 0.00040679456F, 884 0.00043323036F, 0.00046138411F, 0.00049136745F, 0.00052329927F, 885 0.00055730621F, 0.00059352311F, 0.00063209358F, 0.00067317058F, 886 0.00071691700F, 0.00076350630F, 0.00081312324F, 0.00086596457F, 887 0.00092223983F, 0.00098217216F, 0.0010459992F, 0.0011139742F, 888 0.0011863665F, 0.0012634633F, 0.0013455702F, 0.0014330129F, 889 0.0015261382F, 0.0016253153F, 0.0017309374F, 0.0018434235F, 890 0.0019632195F, 0.0020908006F, 0.0022266726F, 0.0023713743F, 891 0.0025254795F, 0.0026895994F, 0.0028643847F, 0.0030505286F, 892 0.0032487691F, 0.0034598925F, 0.0036847358F, 0.0039241906F, 893 0.0041792066F, 0.0044507950F, 0.0047400328F, 0.0050480668F, 894 0.0053761186F, 0.0057254891F, 0.0060975636F, 0.0064938176F, 895 0.0069158225F, 0.0073652516F, 0.0078438871F, 0.0083536271F, 896 0.0088964928F, 0.009474637F, 0.010090352F, 0.010746080F, 897 0.011444421F, 0.012188144F, 0.012980198F, 0.013823725F, 898 0.014722068F, 0.015678791F, 0.016697687F, 0.017782797F, 899 0.018938423F, 0.020169149F, 0.021479854F, 0.022875735F, 900 0.024362330F, 0.025945531F, 0.027631618F, 0.029427276F, 901 0.031339626F, 0.033376252F, 0.035545228F, 0.037855157F, 902 0.040315199F, 0.042935108F, 0.045725273F, 0.048696758F, 903 0.051861348F, 0.055231591F, 0.058820850F, 0.062643361F, 904 0.066714279F, 0.071049749F, 0.075666962F, 0.080584227F, 905 0.085821044F, 0.091398179F, 0.097337747F, 0.10366330F, 906 0.11039993F, 0.11757434F, 0.12521498F, 0.13335215F, 907 0.14201813F, 0.15124727F, 0.16107617F, 0.17154380F, 908 0.18269168F, 0.19456402F, 0.20720788F, 0.22067342F, 909 0.23501402F, 0.25028656F, 0.26655159F, 0.28387361F, 910 0.30232132F, 0.32196786F, 0.34289114F, 0.36517414F, 911 0.38890521F, 0.41417847F, 0.44109412F, 0.46975890F, 912 0.50028648F, 0.53279791F, 0.56742212F, 0.60429640F, 913 0.64356699F, 0.68538959F, 0.72993007F, 0.77736504F, 914 0.82788260F, 0.88168307F, 0.9389798F, 1.F, 915 }; 916 917 /* this is for per-channel noise normalization */ 918 static int apsort(const void *a, const void *b){ 919 float f1=**(float**)a; 920 float f2=**(float**)b; 921 return (f1<f2)-(f1>f2); 922 } 923 924 static void flag_lossless(int limit, float prepoint, float postpoint, float *mdct, 925 float *floor, int *flag, int i, int jn){ 926 int j; 927 for(j=0;j<jn;j++){ 928 float point = j>=limit-i ? postpoint : prepoint; 929 float r = fabs(mdct[j])/floor[j]; 930 if(r<point) 931 flag[j]=0; 932 else 933 flag[j]=1; 934 } 935 } 936 937 /* Overload/Side effect: On input, the *q vector holds either the 938 quantized energy (for elements with the flag set) or the absolute 939 values of the *r vector (for elements with flag unset). On output, 940 *q holds the quantized energy for all elements */ 941 static float noise_normalize(vorbis_look_psy *p, int limit, float *r, float *q, float *f, int *flags, float acc, int i, int n, int *out){ 942 943 vorbis_info_psy *vi=p->vi; 944 float **sort = alloca(n*sizeof(*sort)); 945 int j,count=0; 946 int start = (vi->normal_p ? vi->normal_start-i : n); 947 if(start>n)start=n; 948 949 /* force classic behavior where only energy in the current band is considered */ 950 acc=0.f; 951 952 /* still responsible for populating *out where noise norm not in 953 effect. There's no need to [re]populate *q in these areas */ 954 for(j=0;j<start;j++){ 955 if(!flags || !flags[j]){ /* lossless coupling already quantized. 956 Don't touch; requantizing based on 957 energy would be incorrect. */ 958 float ve = q[j]/f[j]; 959 if(r[j]<0) 960 out[j] = -rint(sqrt(ve)); 961 else 962 out[j] = rint(sqrt(ve)); 963 } 964 } 965 966 /* sort magnitudes for noise norm portion of partition */ 967 for(;j<n;j++){ 968 if(!flags || !flags[j]){ /* can't noise norm elements that have 969 already been loslessly coupled; we can 970 only account for their energy error */ 971 float ve = q[j]/f[j]; 972 /* Despite all the new, more capable coupling code, for now we 973 implement noise norm as it has been up to this point. Only 974 consider promotions to unit magnitude from 0. In addition 975 the only energy error counted is quantizations to zero. */ 976 /* also-- the original point code only applied noise norm at > pointlimit */ 977 if(ve<.25f && (!flags || j>=limit-i)){ 978 acc += ve; 979 sort[count++]=q+j; /* q is fabs(r) for unflagged element */ 980 }else{ 981 /* For now: no acc adjustment for nonzero quantization. populate *out and q as this value is final. */ 982 if(r[j]<0) 983 out[j] = -rint(sqrt(ve)); 984 else 985 out[j] = rint(sqrt(ve)); 986 q[j] = out[j]*out[j]*f[j]; 987 } 988 }/* else{ 989 again, no energy adjustment for error in nonzero quant-- for now 990 }*/ 991 } 992 993 if(count){ 994 /* noise norm to do */ 995 qsort(sort,count,sizeof(*sort),apsort); 996 for(j=0;j<count;j++){ 997 int k=sort[j]-q; 998 if(acc>=vi->normal_thresh){ 999 out[k]=unitnorm(r[k]); 1000 acc-=1.f; 1001 q[k]=f[k]; 1002 }else{ 1003 out[k]=0; 1004 q[k]=0.f; 1005 } 1006 } 1007 } 1008 1009 return acc; 1010 } 1011 1012 /* Noise normalization, quantization and coupling are not wholly 1013 seperable processes in depth>1 coupling. */ 1014 void _vp_couple_quantize_normalize(int blobno, 1015 vorbis_info_psy_global *g, 1016 vorbis_look_psy *p, 1017 vorbis_info_mapping0 *vi, 1018 float **mdct, 1019 int **iwork, 1020 int *nonzero, 1021 int sliding_lowpass, 1022 int ch){ 1023 1024 int i; 1025 int n = p->n; 1026 int partition=(p->vi->normal_p ? p->vi->normal_partition : 16); 1027 int limit = g->coupling_pointlimit[p->vi->blockflag][blobno]; 1028 float prepoint=stereo_threshholds[g->coupling_prepointamp[blobno]]; 1029 float postpoint=stereo_threshholds[g->coupling_postpointamp[blobno]]; 1030 #if 0 1031 float de=0.1*p->m_val; /* a blend of the AoTuV M2 and M3 code here and below */ 1032 #endif 1033 1034 /* mdct is our raw mdct output, floor not removed. */ 1035 /* inout passes in the ifloor, passes back quantized result */ 1036 1037 /* unquantized energy (negative indicates amplitude has negative sign) */ 1038 float **raw = alloca(ch*sizeof(*raw)); 1039 1040 /* dual pupose; quantized energy (if flag set), othersize fabs(raw) */ 1041 float **quant = alloca(ch*sizeof(*quant)); 1042 1043 /* floor energy */ 1044 float **floor = alloca(ch*sizeof(*floor)); 1045 1046 /* flags indicating raw/quantized status of elements in raw vector */ 1047 int **flag = alloca(ch*sizeof(*flag)); 1048 1049 /* non-zero flag working vector */ 1050 int *nz = alloca(ch*sizeof(*nz)); 1051 1052 /* energy surplus/defecit tracking */ 1053 float *acc = alloca((ch+vi->coupling_steps)*sizeof(*acc)); 1054 1055 /* The threshold of a stereo is changed with the size of n */ 1056 if(n > 1000) 1057 postpoint=stereo_threshholds_limited[g->coupling_postpointamp[blobno]]; 1058 1059 raw[0] = alloca(ch*partition*sizeof(**raw)); 1060 quant[0] = alloca(ch*partition*sizeof(**quant)); 1061 floor[0] = alloca(ch*partition*sizeof(**floor)); 1062 flag[0] = alloca(ch*partition*sizeof(**flag)); 1063 1064 for(i=1;i<ch;i++){ 1065 raw[i] = &raw[0][partition*i]; 1066 quant[i] = &quant[0][partition*i]; 1067 floor[i] = &floor[0][partition*i]; 1068 flag[i] = &flag[0][partition*i]; 1069 } 1070 for(i=0;i<ch+vi->coupling_steps;i++) 1071 acc[i]=0.f; 1072 1073 for(i=0;i<n;i+=partition){ 1074 int k,j,jn = partition > n-i ? n-i : partition; 1075 int step,track = 0; 1076 1077 memcpy(nz,nonzero,sizeof(*nz)*ch); 1078 1079 /* prefill */ 1080 memset(flag[0],0,ch*partition*sizeof(**flag)); 1081 for(k=0;k<ch;k++){ 1082 int *iout = &iwork[k][i]; 1083 if(nz[k]){ 1084 1085 for(j=0;j<jn;j++) 1086 floor[k][j] = FLOOR1_fromdB_LOOKUP[iout[j]]; 1087 1088 flag_lossless(limit,prepoint,postpoint,&mdct[k][i],floor[k],flag[k],i,jn); 1089 1090 for(j=0;j<jn;j++){ 1091 quant[k][j] = raw[k][j] = mdct[k][i+j]*mdct[k][i+j]; 1092 if(mdct[k][i+j]<0.f) raw[k][j]*=-1.f; 1093 floor[k][j]*=floor[k][j]; 1094 } 1095 1096 acc[track]=noise_normalize(p,limit,raw[k],quant[k],floor[k],NULL,acc[track],i,jn,iout); 1097 1098 }else{ 1099 for(j=0;j<jn;j++){ 1100 floor[k][j] = 1e-10f; 1101 raw[k][j] = 0.f; 1102 quant[k][j] = 0.f; 1103 flag[k][j] = 0; 1104 iout[j]=0; 1105 } 1106 acc[track]=0.f; 1107 } 1108 track++; 1109 } 1110 1111 /* coupling */ 1112 for(step=0;step<vi->coupling_steps;step++){ 1113 int Mi = vi->coupling_mag[step]; 1114 int Ai = vi->coupling_ang[step]; 1115 int *iM = &iwork[Mi][i]; 1116 int *iA = &iwork[Ai][i]; 1117 float *reM = raw[Mi]; 1118 float *reA = raw[Ai]; 1119 float *qeM = quant[Mi]; 1120 float *qeA = quant[Ai]; 1121 float *floorM = floor[Mi]; 1122 float *floorA = floor[Ai]; 1123 int *fM = flag[Mi]; 1124 int *fA = flag[Ai]; 1125 1126 if(nz[Mi] || nz[Ai]){ 1127 nz[Mi] = nz[Ai] = 1; 1128 1129 for(j=0;j<jn;j++){ 1130 1131 if(j<sliding_lowpass-i){ 1132 if(fM[j] || fA[j]){ 1133 /* lossless coupling */ 1134 1135 reM[j] = fabs(reM[j])+fabs(reA[j]); 1136 qeM[j] = qeM[j]+qeA[j]; 1137 fM[j]=fA[j]=1; 1138 1139 /* couple iM/iA */ 1140 { 1141 int A = iM[j]; 1142 int B = iA[j]; 1143 1144 if(abs(A)>abs(B)){ 1145 iA[j]=(A>0?A-B:B-A); 1146 }else{ 1147 iA[j]=(B>0?A-B:B-A); 1148 iM[j]=B; 1149 } 1150 1151 /* collapse two equivalent tuples to one */ 1152 if(iA[j]>=abs(iM[j])*2){ 1153 iA[j]= -iA[j]; 1154 iM[j]= -iM[j]; 1155 } 1156 1157 } 1158 1159 }else{ 1160 /* lossy (point) coupling */ 1161 if(j<limit-i){ 1162 /* dipole */ 1163 reM[j] += reA[j]; 1164 qeM[j] = fabs(reM[j]); 1165 }else{ 1166 #if 0 1167 /* AoTuV */ 1168 /** @ M2 ** 1169 The boost problem by the combination of noise normalization and point stereo is eased. 1170 However, this is a temporary patch. 1171 by Aoyumi @ 2004/04/18 1172 */ 1173 float derate = (1.0 - de*((float)(j-limit+i) / (float)(n-limit))); 1174 /* elliptical */ 1175 if(reM[j]+reA[j]<0){ 1176 reM[j] = - (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate); 1177 }else{ 1178 reM[j] = (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate); 1179 } 1180 #else 1181 /* elliptical */ 1182 if(reM[j]+reA[j]<0){ 1183 reM[j] = - (qeM[j] = fabs(reM[j])+fabs(reA[j])); 1184 }else{ 1185 reM[j] = (qeM[j] = fabs(reM[j])+fabs(reA[j])); 1186 } 1187 #endif 1188 1189 } 1190 reA[j]=qeA[j]=0.f; 1191 fA[j]=1; 1192 iA[j]=0; 1193 } 1194 } 1195 floorM[j]=floorA[j]=floorM[j]+floorA[j]; 1196 } 1197 /* normalize the resulting mag vector */ 1198 acc[track]=noise_normalize(p,limit,raw[Mi],quant[Mi],floor[Mi],flag[Mi],acc[track],i,jn,iM); 1199 track++; 1200 } 1201 } 1202 } 1203 1204 for(i=0;i<vi->coupling_steps;i++){ 1205 /* make sure coupling a zero and a nonzero channel results in two 1206 nonzero channels. */ 1207 if(nonzero[vi->coupling_mag[i]] || 1208 nonzero[vi->coupling_ang[i]]){ 1209 nonzero[vi->coupling_mag[i]]=1; 1210 nonzero[vi->coupling_ang[i]]=1; 1211 } 1212 } 1213 }