1 /********************************************************************
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. *
8 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 *
9 * by the XIPHOPHORUS Company http://www.xiph.org/ *
11 ********************************************************************
13 function: psychoacoustics not including preecho
14 last mod: $Id: psy.c,v 1.50 2001/08/13 01:36:57 xiphmont Exp $
16 ********************************************************************/
21 #include "vorbis/codec.h"
22 #include "codec_internal.h"
32 #define NEGINF -9999.f
34 /* Why Bark scale for encoding but not masking computation? Because
35 masking has a strong harmonic dependancy */
37 vorbis_look_psy_global *_vp_global_look(vorbis_info *vi){
39 codec_setup_info *ci=vi->codec_setup;
40 vorbis_info_psy_global *gi=ci->psy_g_param;
41 vorbis_look_psy_global *look=_ogg_calloc(1,sizeof(vorbis_look_psy_global));
43 int shiftoc=rint(log(gi->eighth_octave_lines*8)/log(2))-1;
44 look->decaylines=toOC(96000.f)*(1<<(shiftoc+1))+.5f; /* max sample
48 look->decay=_ogg_calloc(vi->channels,sizeof(float *));
49 for(i=0;i<vi->channels;i++){
50 look->decay[i]=_ogg_calloc(look->decaylines,sizeof(float));
51 for(j=0;j<look->decaylines;j++)
52 look->decay[i][j]=-9999.;
54 look->channels=vi->channels;
61 void _vp_global_free(vorbis_look_psy_global *look){
64 for(i=0;i<look->channels;i++)
65 _ogg_free(look->decay[i]);
66 _ogg_free(look->decay);
68 memset(look,0,sizeof(vorbis_look_psy_global));
72 void _vi_psy_free(vorbis_info_psy *i){
74 memset(i,0,sizeof(vorbis_info_psy));
79 vorbis_info_psy *_vi_psy_copy(vorbis_info_psy *i){
80 vorbis_info_psy *ret=_ogg_malloc(sizeof(vorbis_info_psy));
81 memcpy(ret,i,sizeof(vorbis_info_psy));
85 /* Set up decibel threshold slopes on a Bark frequency scale */
86 /* ATH is the only bit left on a Bark scale. No reason to change it
88 static void set_curve(float *ref,float *c,int n, float crate){
91 for(i=0;i<MAX_BARK-1;i++){
92 int endpos=rint(fromBARK(i+1)*2*n/crate);
95 float delta=(ref[i+1]-base)/(endpos-j);
96 for(;j<endpos && j<n;j++){
104 static void min_curve(float *c,
107 for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i];
109 static void max_curve(float *c,
112 for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i];
115 static void attenuate_curve(float *c,float att){
117 for(i=0;i<EHMER_MAX;i++)
121 static void interp_curve(float *c,float *c1,float *c2,float del){
123 for(i=0;i<EHMER_MAX;i++)
124 c[i]=c2[i]*del+c1[i]*(1.f-del);
127 extern int analysis_noisy;
128 static void setup_curve(float **c,
132 float ath[EHMER_MAX];
133 float tempc[P_LEVELS][EHMER_MAX];
134 float *ATH=ATH_Bark_dB_lspconservative; /* just for limiting here */
136 memcpy(c[0]+2,c[4]+2,sizeof(float)*EHMER_MAX);
137 memcpy(c[2]+2,c[4]+2,sizeof(float)*EHMER_MAX);
139 /* we add back in the ATH to avoid low level curves falling off to
140 -infinity and unneccessarily cutting off high level curves in the
141 curve limiting (last step). But again, remember... a half-band's
142 settings must be valid over the whole band, and it's better to
143 mask too little than too much, so be pessimal. */
145 for(i=0;i<EHMER_MAX;i++){
146 float oc_min=band*.5+(i-EHMER_OFFSET)*.125;
147 float oc_max=band*.5+(i-EHMER_OFFSET+1)*.125;
148 float bark=toBARK(fromOC(oc_min));
149 int ibark=floor(bark);
150 float del=bark-ibark;
151 float ath_min,ath_max;
154 ath_min=ATH[ibark]*(1.f-del)+ATH[ibark+1]*del;
158 bark=toBARK(fromOC(oc_max));
163 ath_max=ATH[ibark]*(1.f-del)+ATH[ibark+1]*del;
167 ath[i]=min(ath_min,ath_max);
170 /* The c array is comes in as dB curves at 20 40 60 80 100 dB.
171 interpolate intermediate dB curves */
172 for(i=1;i<P_LEVELS;i+=2){
173 interp_curve(c[i]+2,c[i-1]+2,c[i+1]+2,.5);
176 /* normalize curves so the driving amplitude is 0dB */
177 /* make temp curves with the ATH overlayed */
178 for(i=0;i<P_LEVELS;i++){
179 attenuate_curve(c[i]+2,curveatt_dB[i]);
180 memcpy(tempc[i],ath,EHMER_MAX*sizeof(float));
181 attenuate_curve(tempc[i],-i*10.f);
182 max_curve(tempc[i],c[i]+2);
185 /* Now limit the louder curves.
187 the idea is this: We don't know what the playback attenuation
188 will be; 0dB SL moves every time the user twiddles the volume
189 knob. So that means we have to use a single 'most pessimal' curve
190 for all masking amplitudes, right? Wrong. The *loudest* sound
191 can be in (we assume) a range of ...+100dB] SL. However, sounds
192 20dB down will be in a range ...+80], 40dB down is from ...+60],
195 for(j=1;j<P_LEVELS;j++){
196 min_curve(tempc[j],tempc[j-1]);
197 min_curve(c[j]+2,tempc[j]);
201 for(j=0;j<P_LEVELS;j++){
203 for(i=0;i<EHMER_OFFSET;i++)
204 if(c[j][i+2]>-200.f)break;
207 for(i=EHMER_MAX-1;i>EHMER_OFFSET+1;i--)
215 void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi,
216 vorbis_info_psy_global *gi,int n,long rate){
217 long i,j,k,lo=0,hi=0;
219 memset(p,0,sizeof(vorbis_look_psy));
222 p->eighth_octave_lines=gi->eighth_octave_lines;
223 p->shiftoc=rint(log(gi->eighth_octave_lines*8)/log(2))-1;
225 p->firstoc=toOC(.25f*rate/n)*(1<<(p->shiftoc+1))-gi->eighth_octave_lines;
226 maxoc=toOC((n*.5f-.25f)*rate/n)*(1<<(p->shiftoc+1))+.5f;
227 p->total_octave_lines=maxoc-p->firstoc+1;
230 p->ath=_ogg_malloc(n*sizeof(float));
231 p->octave=_ogg_malloc(n*sizeof(long));
232 p->bark=_ogg_malloc(n*sizeof(unsigned long));
237 /* set up the lookups for a given blocksize and sample rate */
239 set_curve(vi->ath, p->ath,n,rate);
241 float bark=toBARK(rate/(2*n)*i);
243 for(;lo+vi->noisewindowlomin<i &&
244 toBARK(rate/(2*n)*lo)<(bark-vi->noisewindowlo);lo++);
246 for(;hi<n && (hi<i+vi->noisewindowhimin ||
247 toBARK(rate/(2*n)*hi)<(bark+vi->noisewindowhi));hi++);
249 p->bark[i]=(hi<<16)+lo;
254 p->octave[i]=toOC((i*.5f+.25f)*rate/n)*(1<<(p->shiftoc+1))+.5f;
256 p->tonecurves=_ogg_malloc(P_BANDS*sizeof(float **));
257 p->noisethresh=_ogg_malloc(n*sizeof(float));
258 p->noiseoffset=_ogg_malloc(n*sizeof(float));
259 for(i=0;i<P_BANDS;i++)
260 p->tonecurves[i]=_ogg_malloc(P_LEVELS*sizeof(float *));
262 for(i=0;i<P_BANDS;i++)
263 for(j=0;j<P_LEVELS;j++)
264 p->tonecurves[i][j]=_ogg_malloc((EHMER_MAX+2)*sizeof(float));
267 /* OK, yeah, this was a silly way to do it */
268 memcpy(p->tonecurves[0][4]+2,tone_125_40dB_SL,sizeof(float)*EHMER_MAX);
269 memcpy(p->tonecurves[0][6]+2,tone_125_60dB_SL,sizeof(float)*EHMER_MAX);
270 memcpy(p->tonecurves[0][8]+2,tone_125_80dB_SL,sizeof(float)*EHMER_MAX);
271 memcpy(p->tonecurves[0][10]+2,tone_125_100dB_SL,sizeof(float)*EHMER_MAX);
273 memcpy(p->tonecurves[2][4]+2,tone_125_40dB_SL,sizeof(float)*EHMER_MAX);
274 memcpy(p->tonecurves[2][6]+2,tone_125_60dB_SL,sizeof(float)*EHMER_MAX);
275 memcpy(p->tonecurves[2][8]+2,tone_125_80dB_SL,sizeof(float)*EHMER_MAX);
276 memcpy(p->tonecurves[2][10]+2,tone_125_100dB_SL,sizeof(float)*EHMER_MAX);
278 memcpy(p->tonecurves[4][4]+2,tone_250_40dB_SL,sizeof(float)*EHMER_MAX);
279 memcpy(p->tonecurves[4][6]+2,tone_250_60dB_SL,sizeof(float)*EHMER_MAX);
280 memcpy(p->tonecurves[4][8]+2,tone_250_80dB_SL,sizeof(float)*EHMER_MAX);
281 memcpy(p->tonecurves[4][10]+2,tone_250_100dB_SL,sizeof(float)*EHMER_MAX);
283 memcpy(p->tonecurves[6][4]+2,tone_500_40dB_SL,sizeof(float)*EHMER_MAX);
284 memcpy(p->tonecurves[6][6]+2,tone_500_60dB_SL,sizeof(float)*EHMER_MAX);
285 memcpy(p->tonecurves[6][8]+2,tone_500_80dB_SL,sizeof(float)*EHMER_MAX);
286 memcpy(p->tonecurves[6][10]+2,tone_500_100dB_SL,sizeof(float)*EHMER_MAX);
288 memcpy(p->tonecurves[8][4]+2,tone_1000_40dB_SL,sizeof(float)*EHMER_MAX);
289 memcpy(p->tonecurves[8][6]+2,tone_1000_60dB_SL,sizeof(float)*EHMER_MAX);
290 memcpy(p->tonecurves[8][8]+2,tone_1000_80dB_SL,sizeof(float)*EHMER_MAX);
291 memcpy(p->tonecurves[8][10]+2,tone_1000_100dB_SL,sizeof(float)*EHMER_MAX);
293 memcpy(p->tonecurves[10][4]+2,tone_2000_40dB_SL,sizeof(float)*EHMER_MAX);
294 memcpy(p->tonecurves[10][6]+2,tone_2000_60dB_SL,sizeof(float)*EHMER_MAX);
295 memcpy(p->tonecurves[10][8]+2,tone_2000_80dB_SL,sizeof(float)*EHMER_MAX);
296 memcpy(p->tonecurves[10][10]+2,tone_2000_100dB_SL,sizeof(float)*EHMER_MAX);
298 memcpy(p->tonecurves[12][4]+2,tone_4000_40dB_SL,sizeof(float)*EHMER_MAX);
299 memcpy(p->tonecurves[12][6]+2,tone_4000_60dB_SL,sizeof(float)*EHMER_MAX);
300 memcpy(p->tonecurves[12][8]+2,tone_4000_80dB_SL,sizeof(float)*EHMER_MAX);
301 memcpy(p->tonecurves[12][10]+2,tone_4000_100dB_SL,sizeof(float)*EHMER_MAX);
303 memcpy(p->tonecurves[14][4]+2,tone_8000_40dB_SL,sizeof(float)*EHMER_MAX);
304 memcpy(p->tonecurves[14][6]+2,tone_8000_60dB_SL,sizeof(float)*EHMER_MAX);
305 memcpy(p->tonecurves[14][8]+2,tone_8000_80dB_SL,sizeof(float)*EHMER_MAX);
306 memcpy(p->tonecurves[14][10]+2,tone_8000_100dB_SL,sizeof(float)*EHMER_MAX);
308 memcpy(p->tonecurves[16][4]+2,tone_8000_40dB_SL,sizeof(float)*EHMER_MAX);
309 memcpy(p->tonecurves[16][6]+2,tone_8000_60dB_SL,sizeof(float)*EHMER_MAX);
310 memcpy(p->tonecurves[16][8]+2,tone_8000_80dB_SL,sizeof(float)*EHMER_MAX);
311 memcpy(p->tonecurves[16][10]+2,tone_8000_100dB_SL,sizeof(float)*EHMER_MAX);
313 /* value limit the tonal masking curves; the peakatt not only
314 optionally specifies maximum dynamic depth, but also [always]
315 limits the masking curves to a minimum depth */
316 for(i=0;i<P_BANDS;i+=2)
317 for(j=4;j<P_LEVELS;j+=2)
318 for(k=2;k<EHMER_MAX+2;k++)
319 p->tonecurves[i][j][k]+=vi->tone_masteratt;
321 /* interpolate curves between */
322 for(i=1;i<P_BANDS;i+=2)
323 for(j=4;j<P_LEVELS;j+=2){
324 memcpy(p->tonecurves[i][j]+2,p->tonecurves[i-1][j]+2,EHMER_MAX*sizeof(float));
325 /*interp_curve(p->tonecurves[i][j],
326 p->tonecurves[i-1][j],
327 p->tonecurves[i+1][j],.5);*/
328 min_curve(p->tonecurves[i][j]+2,p->tonecurves[i+1][j]+2);
331 /* set up the final curves */
332 for(i=0;i<P_BANDS;i++)
333 setup_curve(p->tonecurves[i],i,vi->toneatt->block[i]);
336 /* value limit the tonal masking curves; the peakatt not only
337 optionally specifies maximum dynamic depth, but also [always]
338 limits the masking curves to a minimum depth */
339 for(i=0;i<P_BANDS;i++)
340 for(j=0;j<P_LEVELS;j++){
341 for(k=2;k<EHMER_OFFSET+2+vi->curvelimitp;k++)
342 if(p->tonecurves[i][j][k]> vi->peakatt->block[i][j])
343 p->tonecurves[i][j][k]= vi->peakatt->block[i][j];
349 if(vi->peakattp) /* we limit depth only optionally */
350 for(i=0;i<P_BANDS;i++)
351 for(j=0;j<P_LEVELS;j++)
352 if(p->tonecurves[i][j][EHMER_OFFSET+2]< vi->peakatt->block[i][j])
353 p->tonecurves[i][j][EHMER_OFFSET+2]= vi->peakatt->block[i][j];
355 /* but guarding is mandatory */
356 for(i=0;i<P_BANDS;i++)
357 for(j=0;j<P_LEVELS;j++)
358 if(p->tonecurves[i][j][EHMER_OFFSET+2]< vi->tone_maxatt)
359 p->tonecurves[i][j][EHMER_OFFSET+2]= vi->tone_maxatt;
361 /* set up rolling noise median */
363 float halfoc=toOC((i+.5)*rate/(2.*n))*2.;
367 if(halfoc<0)halfoc=0;
368 if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1;
369 inthalfoc=(int)halfoc;
370 del=halfoc-inthalfoc;
372 p->noisethresh[i]=((p->vi->noisethresh[inthalfoc]*(1.-del) +
373 p->vi->noisethresh[inthalfoc+1]*del))*2.f-1.f;
375 p->vi->noiseoff[inthalfoc]*(1.-del) +
376 p->vi->noiseoff[inthalfoc+1]*del;
380 _analysis_output("noiseoff",0,p->noiseoffset,n,1,0);
381 _analysis_output("noisethresh",0,p->noisethresh,n,1,0);
383 for(i=0;i<P_LEVELS;i++)
384 _analysis_output("curve_63Hz",i,p->tonecurves[0][i]+2,EHMER_MAX,0,0);
385 for(i=0;i<P_LEVELS;i++)
386 _analysis_output("curve_88Hz",i,p->tonecurves[1][i]+2,EHMER_MAX,0,0);
387 for(i=0;i<P_LEVELS;i++)
388 _analysis_output("curve_125Hz",i,p->tonecurves[2][i]+2,EHMER_MAX,0,0);
389 for(i=0;i<P_LEVELS;i++)
390 _analysis_output("curve_170Hz",i,p->tonecurves[3][i]+2,EHMER_MAX,0,0);
391 for(i=0;i<P_LEVELS;i++)
392 _analysis_output("curve_250Hz",i,p->tonecurves[4][i]+2,EHMER_MAX,0,0);
393 for(i=0;i<P_LEVELS;i++)
394 _analysis_output("curve_350Hz",i,p->tonecurves[5][i]+2,EHMER_MAX,0,0);
395 for(i=0;i<P_LEVELS;i++)
396 _analysis_output("curve_500Hz",i,p->tonecurves[6][i]+2,EHMER_MAX,0,0);
397 for(i=0;i<P_LEVELS;i++)
398 _analysis_output("curve_700Hz",i,p->tonecurves[7][i]+2,EHMER_MAX,0,0);
399 for(i=0;i<P_LEVELS;i++)
400 _analysis_output("curve_1kHz",i,p->tonecurves[8][i]+2,EHMER_MAX,0,0);
401 for(i=0;i<P_LEVELS;i++)
402 _analysis_output("curve_1.4Hz",i,p->tonecurves[9][i]+2,EHMER_MAX,0,0);
403 for(i=0;i<P_LEVELS;i++)
404 _analysis_output("curve_2kHz",i,p->tonecurves[10][i]+2,EHMER_MAX,0,0);
405 for(i=0;i<P_LEVELS;i++)
406 _analysis_output("curve_2.4kHz",i,p->tonecurves[11][i]+2,EHMER_MAX,0,0);
407 for(i=0;i<P_LEVELS;i++)
408 _analysis_output("curve_4kHz",i,p->tonecurves[12][i]+2,EHMER_MAX,0,0);
409 for(i=0;i<P_LEVELS;i++)
410 _analysis_output("curve_5.6kHz",i,p->tonecurves[13][i]+2,EHMER_MAX,0,0);
411 for(i=0;i<P_LEVELS;i++)
412 _analysis_output("curve_8kHz",i,p->tonecurves[14][i]+2,EHMER_MAX,0,0);
413 for(i=0;i<P_LEVELS;i++)
414 _analysis_output("curve_11.5kHz",i,p->tonecurves[15][i]+2,EHMER_MAX,0,0);
415 for(i=0;i<P_LEVELS;i++)
416 _analysis_output("curve_16kHz",i,p->tonecurves[16][i]+2,EHMER_MAX,0,0);
421 void _vp_psy_clear(vorbis_look_psy *p){
424 if(p->ath)_ogg_free(p->ath);
425 if(p->octave)_ogg_free(p->octave);
426 if(p->bark)_ogg_free(p->bark);
428 for(i=0;i<P_BANDS;i++){
429 for(j=0;j<P_LEVELS;j++){
430 _ogg_free(p->tonecurves[i][j]);
432 _ogg_free(p->tonecurves[i]);
434 _ogg_free(p->tonecurves);
435 _ogg_free(p->noiseoffset);
437 memset(p,0,sizeof(vorbis_look_psy));
441 /* octave/(8*eighth_octave_lines) x scale and dB y scale */
442 static void seed_curve(float *seed,
443 const float **curves,
446 int linesper,float dBoffset){
449 const float *posts,*curve;
451 int choice=(int)((amp+dBoffset)*.1f);
452 choice=max(choice,0);
453 choice=min(choice,P_LEVELS-1);
454 posts=curves[choice];
457 seedptr=oc+(posts[0]-16)*linesper-(linesper>>1);
459 for(i=posts[0];i<post1;i++){
461 float lin=amp+curve[i];
462 if(seed[seedptr]<lin)seed[seedptr]=lin;
469 static void seed_loop(vorbis_look_psy *p,
470 const float ***curves,
475 vorbis_info_psy *vi=p->vi;
477 float dBoffset=vi->max_curve_dB-specmax;
479 /* prime the working vector with peak values */
483 long oc=p->octave[i];
484 while(i+1<n && p->octave[i+1]==oc){
486 if(f[i]>max)max=f[i];
491 if(oc>=P_BANDS)oc=P_BANDS-1;
496 p->octave[i]-p->firstoc,
497 p->total_octave_lines,
498 p->eighth_octave_lines,
504 static void seed_chase(float *seeds, int linesper, long n){
505 long *posstack=alloca(n*sizeof(long));
506 float *ampstack=alloca(n*sizeof(float));
514 ampstack[stack++]=seeds[i];
517 if(seeds[i]<ampstack[stack-1]){
519 ampstack[stack++]=seeds[i];
522 if(i<posstack[stack-1]+linesper){
523 if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] &&
524 i<posstack[stack-2]+linesper){
525 /* we completely overlap, making stack-1 irrelevant. pop it */
531 ampstack[stack++]=seeds[i];
539 /* the stack now contains only the positions that are relevant. Scan
540 'em straight through */
542 for(i=0;i<stack;i++){
544 if(i<stack-1 && ampstack[i+1]>ampstack[i]){
545 endpos=posstack[i+1];
547 endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is
548 discarded in short frames */
550 if(endpos>n)endpos=n;
551 for(;pos<endpos;pos++)
552 seeds[pos]=ampstack[i];
555 /* there. Linear time. I now remember this was on a problem set I
556 had in Grad Skool... I didn't solve it at the time ;-) */
560 /* bleaugh, this is more complicated than it needs to be */
561 static void max_seeds(vorbis_look_psy *p,
562 vorbis_look_psy_global *g,
566 long n=p->total_octave_lines;
567 int linesper=p->eighth_octave_lines;
571 seed_chase(seed,linesper,n); /* for masking */
573 pos=p->octave[0]-p->firstoc-(linesper>>1);
574 while(linpos+1<p->n){
575 float minV=seed[pos];
576 long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc;
579 if((seed[pos]>NEGINF && seed[pos]<minV) || minV==NEGINF)
583 /* seed scale is log. Floor is linear. Map back to it */
585 for(;linpos<p->n && p->octave[linpos]<=end;linpos++)
586 if(flr[linpos]<minV)flr[linpos]=minV;
590 float minV=seed[p->total_octave_lines-1];
591 for(;linpos<p->n;linpos++)
592 if(flr[linpos]<minV)flr[linpos]=minV;
597 static void bark_noise_pointmp(int n,const long *b,
601 long i,hi=0,lo=0,hif=0,lof=0;
614 double bin=(f[hi]<-140.f?0.:f[hi]+140.);
625 double bin=(f[lo]<-140.f?0.:f[lo]+140.);
636 if(hif<n && fixed>0){
640 double bin=(f[hif]<-140.f?0.:f[hif]+140.);
652 double bin=(f[lof]<-140.f?0.:f[lof]+140.);
664 double denom=1./(na*x2a-xa*xa);
665 double a=(ya*x2a-xya*xa)*denom;
666 double b=(na*xya-xa*ya)*denom;
670 double denomf=1./(nb*x2b-xb*xb);
671 double af=(yb*x2b-xyb*xb)*denomf;
672 double bf=(nb*xyb-xb*yb)*denomf;
682 static void bark_noise_hybridmp(int n,const long *b,
686 long i,hi=0,lo=0,hif=0,lof=0;
693 int first=-1,firstf=-1;
714 if(first==-1)first=hi;
731 if(first<lo)first=-1;
735 for(first=lo;first<hi;first++)
736 if(f[first]>0.f)break;
737 if(first==hi)first=-1;
742 if(hif<n && fixed>0){
759 if(firstf==-1)firstf=hif;
777 if(firstf<lof)firstf=-1;
781 for(firstf=lof;firstf<hif;firstf++)
782 if(f[firstf]>0.f)break;
783 if(firstf==hif)firstf=-1;
791 if(rna>2 && (last-first)*3/2>hi-lo){
792 double denom=1./(na*x2a-xa*xa);
793 double a=(ya*x2a-xya*xa)*denom;
794 double b=(na*xya-xa*ya)*denom;
805 if(rnb>2 && (lastf-firstf)*3/2>hif-lof){
806 double denomf=1./(nb*x2b-xb*xb);
807 double af=(yb*x2b-xyb*xb)*denomf;
808 double bf=(nb*xyb-xb*yb)*denomf;
816 if(va>vb && vb>0.)va=vb;
826 void _vp_remove_floor(vorbis_look_psy *p,
827 vorbis_look_psy_global *g,
832 float local_specmax){
837 residue[i]=mdct[i]/codedflr[i];
843 void _vp_compute_mask(vorbis_look_psy *p,
844 vorbis_look_psy_global *g,
849 float global_specmax,
855 float *seed=alloca(sizeof(float)*p->total_octave_lines);
856 for(i=0;i<p->total_octave_lines;i++)seed[i]=NEGINF;
859 if(p->vi->noisemaskp){
860 float *work=alloca(n*sizeof(float));
862 bark_noise_pointmp(n,p->bark,logmdct,logmask,
865 for(i=0;i<n;i++)work[i]=logmdct[i]-logmask[i];
867 _analysis_output("medianmdct",seq,work,n,1,0);
868 bark_noise_hybridmp(n,p->bark,work,logmask,
869 p->vi->noisewindowfixed);
871 for(i=0;i<n;i++)work[i]=logmdct[i]-work[i];
873 /* work[i] holds the median line (.5), logmask holds the upper
874 envelope line (1.) */
876 _analysis_output("median",seq,work,n,1,0);
877 _analysis_output("envelope",seq,logmask,n,1,0);
879 for(i=0;i<n;i++)logmask[i]=
881 p->noisethresh[i]*logmask[i]+p->noiseoffset[i];
883 _analysis_output("noise",seq,logmask,n,1,0);
886 for(i=0;i<n;i++)logmask[i]=NEGINF;
889 /* set the ATH (floating below localmax, not global max by a
892 float att=local_specmax+p->vi->ath_adjatt;
893 if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt;
896 float av=p->ath[i]+att;
897 if(av>logmask[i])logmask[i]=av;
901 /* tone/peak masking */
902 seed_loop(p,(const float ***)p->tonecurves,logfft,logmask,seed,global_specmax);
903 max_seeds(p,g,channel,seed,logmask);
905 /* suppress any curve > p->vi->noisemaxsupp */
906 if(p->vi->noisemaxsupp<0.f)
908 if(logmask[i]>p->vi->noisemaxsupp)
909 logmask[i]=p->vi->noisemaxsupp;
912 /* doing this here is clean, but we need to find a faster way to do
913 it than to just tack it on */
915 for(i=0;i<n;i++)if(logmdct[i]>=logmask[i])break;
917 for(i=0;i<n;i++)logmask[i]=NEGINF;
920 logfft[i]=max(logmdct[i],logfft[i]);
926 float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){
927 vorbis_info *vi=vd->vi;
928 codec_setup_info *ci=vi->codec_setup;
929 vorbis_info_psy_global *gi=ci->psy_g_param;
931 int n=ci->blocksizes[vd->W]/2;
932 float secs=(float)n/vi->rate;
934 amp+=secs*gi->ampmax_att_per_sec;
935 if(amp<-9999)amp=-9999;
939 static void couple_lossless(float A, float B, float *mag, float *ang){
941 *mag=A; *ang=(A>0.f?A-B:B-A);
943 *mag=B; *ang=(B>0.f?A-B:B-A);
947 static void couple_8phase(float A, float B, float *mag, float *ang){
949 *mag=A; *ang=(A>0?A-B:B-A);
951 *mag=B; *ang=(B>0?A-B:B-A);
955 switch((int)(rint(*ang / *mag))){
971 static void couple_6phase(float A, float B, float *mag, float *ang){
973 *mag=A; *ang=(A>0?A-B:B-A);
975 *mag=B; *ang=(B>0?A-B:B-A);
979 switch((int)(rint(*ang / *mag))){
995 static void couple_point(float A, float B, float *mag, float *ang){
997 *mag=A; *ang=(A>0?A-B:B-A);
999 *mag=B; *ang=(B>0?A-B:B-A);
1003 switch((int)(rint(*ang / *mag))){
1007 case 0:case 1: case -1:
1013 void _vp_quantize_couple(vorbis_look_psy *p,
1014 vorbis_info_mapping0 *vi,
1022 vorbis_info_psy *info=p->vi;
1024 /* perform any requested channel coupling */
1025 for(i=0;i<vi->coupling_steps;i++){
1026 float granulem=info->couple_pass[passno].granulem;
1027 float igranulem=info->couple_pass[passno].igranulem;
1029 /* make sure coupling a zero and a nonzero channel results in two
1030 nonzero channels. */
1031 if(nonzero[vi->coupling_mag[i]] ||
1032 nonzero[vi->coupling_ang[i]]){
1034 float *pcmM=pcm[vi->coupling_mag[i]];
1035 float *pcmA=pcm[vi->coupling_ang[i]];
1036 float *sofarM=sofar[vi->coupling_mag[i]];
1037 float *sofarA=sofar[vi->coupling_ang[i]];
1038 float *qM=quantized[vi->coupling_mag[i]];
1039 float *qA=quantized[vi->coupling_ang[i]];
1041 nonzero[vi->coupling_mag[i]]=1;
1042 nonzero[vi->coupling_ang[i]]=1;
1044 for(j=0,k=0;j<n;k++){
1045 vp_couple *part=info->couple_pass[passno].couple_pass+k;
1047 for(;j<part->limit && j<p->n;j++){
1048 /* partition by partition; k is our by-location partition
1051 float Am=rint(pcmM[j]*igranulem)*granulem;
1052 float Bm=rint(pcmA[j]*igranulem)*granulem;
1053 float ang,mag,fmag=max(fabs(Am),fabs(Bm));
1055 if(fmag<part->amppost_point){
1056 couple_point(Am,Bm,&mag,&ang);
1058 if(fmag<part->amppost_6phase){
1059 couple_6phase(Am,Bm,&mag,&ang);
1061 if(fmag<part->amppost_8phase){
1062 couple_8phase(Am,Bm,&mag,&ang);
1064 couple_lossless(Am,Bm,&mag,&ang);
1069 if(ang>fmag*1.9999f)ang=-fmag*2.f;
1071 qM[j]=mag-sofarM[j];
1072 qA[j]=ang-sofarA[j];