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-2010 *
9 * by the Xiph.Org Foundation https://xiph.org/ *
11 ********************************************************************
13 function: psychoacoustics not including preecho
15 ********************************************************************/
20 #include "vorbis/codec.h"
21 #include "codec_internal.h"
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};
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));
40 look->channels=vi->channels;
47 void _vp_global_free(vorbis_look_psy_global *look){
49 memset(look,0,sizeof(*look));
54 void _vi_gpsy_free(vorbis_info_psy_global *i){
56 memset(i,0,sizeof(*i));
61 void _vi_psy_free(vorbis_info_psy *i){
63 memset(i,0,sizeof(*i));
68 static void min_curve(float *c,
71 for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i];
73 static void max_curve(float *c,
76 for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i];
79 static void attenuate_curve(float *c,float att){
81 for(i=0;i<EHMER_MAX;i++)
85 static float ***setup_tone_curves(float curveatt_dB[P_BANDS],float binHz,int n,
86 float center_boost, float center_decay_rate){
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));
93 float ***ret=_ogg_malloc(sizeof(*ret)*P_BANDS);
95 memset(workc,0,sizeof(workc));
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). */
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 */
105 for(j=0;j<EHMER_MAX;j++){
108 if(j+k+ath_offset<MAX_ATH){
109 if(min>ATH[j+k+ath_offset])min=ATH[j+k+ath_offset];
111 if(min>ATH[MAX_ATH-1])min=ATH[MAX_ATH-1];
116 /* copy curves into working space, replicate the 50dB curve to 30
117 and 40, replicate the 100dB curve to 110 */
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]));
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.;
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]);
142 /* Now limit the louder curves.
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],
152 for(j=1;j<P_LEVELS;j++){
153 min_curve(athc[j],athc[j-1]);
154 min_curve(workc[i][j],athc[j]);
158 for(i=0;i<P_BANDS;i++){
159 int hi_curve,lo_curve,bin;
160 ret[i]=_ogg_malloc(sizeof(**ret)*P_LEVELS);
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. */
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;
179 for(m=0;m<P_LEVELS;m++){
180 ret[i][m]=_ogg_malloc(sizeof(***ret)*(EHMER_MAX+2));
182 for(j=0;j<n;j++)brute_buffer[j]=999.;
184 /* render the curve into bins, then pull values back into curve.
185 The point is that any inherent subsampling aliasing results in
187 for(k=lo_curve;k<=hi_curve;k++){
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;
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;
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];
206 if(brute_buffer[l]>workc[k][m][EHMER_MAX-1])
207 brute_buffer[l]=workc[k][m][EHMER_MAX-1];
211 /* be equally paranoid about being valid up to next half ocatve */
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;
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;
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];
231 if(brute_buffer[l]>workc[k][m][EHMER_MAX-1])
232 brute_buffer[l]=workc[k][m][EHMER_MAX-1];
237 for(j=0;j<EHMER_MAX;j++){
238 int bin=fromOC(j*.125+i*.5-2.)/binHz;
240 ret[i][m][j+2]=-999.;
243 ret[i][m][j+2]=-999.;
245 ret[i][m][j+2]=brute_buffer[bin];
251 for(j=0;j<EHMER_OFFSET;j++)
252 if(ret[i][m][j+2]>-200.f)break;
255 for(j=EHMER_MAX-1;j>EHMER_OFFSET+1;j--)
256 if(ret[i][m][j+2]>-200.f)
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;
270 memset(p,0,sizeof(*p));
272 p->eighth_octave_lines=gi->eighth_octave_lines;
273 p->shiftoc=rint(log(gi->eighth_octave_lines*8.f)/log(2.f))-1;
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));
280 p->octave=_ogg_malloc(n*sizeof(*p->octave));
281 p->bark=_ogg_malloc(n*sizeof(*p->bark));
286 /* AoTuV HF weighting */
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 */
292 /* set up the lookups for a given blocksize and sample rate */
294 for(i=0,j=0;i<MAX_ATH-1;i++){
295 int endpos=rint(fromOC((i+1)*.125-2.)*2*n/rate);
298 float delta=(ATH[i+1]-base)/(endpos-j);
299 for(;j<endpos && j<n;j++){
307 p->ath[j]=p->ath[j-1];
311 float bark=toBARK(rate/(2*n)*i);
313 for(;lo+vi->noisewindowlomin<i &&
314 toBARK(rate/(2*n)*lo)<(bark-vi->noisewindowlo);lo++);
316 for(;hi<=n && (hi<i+vi->noisewindowhimin ||
317 toBARK(rate/(2*n)*hi)<(bark+vi->noisewindowhi));hi++);
319 p->bark[i]=((lo-1)<<16)+(hi-1);
324 p->octave[i]=toOC((i+.25f)*.5*rate/n)*(1<<(p->shiftoc+1))+.5f;
326 p->tonecurves=setup_tone_curves(vi->toneatt,rate*.5/n,n,
327 vi->tone_centerboost,vi->tone_decay);
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));
335 float halfoc=toOC((i+.5)*rate/(2.*n))*2.;
339 if(halfoc<0)halfoc=0;
340 if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1;
341 inthalfoc=(int)halfoc;
342 del=halfoc-inthalfoc;
344 for(j=0;j<P_NOISECURVES;j++)
345 p->noiseoffset[j][i]=
346 p->vi->noiseoff[j][inthalfoc]*(1.-del) +
347 p->vi->noiseoff[j][inthalfoc+1]*del;
353 _analysis_output_always("noiseoff0",ls,p->noiseoffset[0],n,1,0,0);
354 _analysis_output_always("noiseoff1",ls,p->noiseoffset[1],n,1,0,0);
355 _analysis_output_always("noiseoff2",ls++,p->noiseoffset[2],n,1,0,0);
360 void _vp_psy_clear(vorbis_look_psy *p){
363 if(p->ath)_ogg_free(p->ath);
364 if(p->octave)_ogg_free(p->octave);
365 if(p->bark)_ogg_free(p->bark);
367 for(i=0;i<P_BANDS;i++){
368 for(j=0;j<P_LEVELS;j++){
369 _ogg_free(p->tonecurves[i][j]);
371 _ogg_free(p->tonecurves[i]);
373 _ogg_free(p->tonecurves);
376 for(i=0;i<P_NOISECURVES;i++){
377 _ogg_free(p->noiseoffset[i]);
379 _ogg_free(p->noiseoffset);
381 memset(p,0,sizeof(*p));
385 /* octave/(8*eighth_octave_lines) x scale and dB y scale */
386 static void seed_curve(float *seed,
387 const float **curves,
390 int linesper,float dBoffset){
393 const float *posts,*curve;
395 int choice=(int)((amp+dBoffset-P_LEVEL_0)*.1f);
396 choice=max(choice,0);
397 choice=min(choice,P_LEVELS-1);
398 posts=curves[choice];
401 seedptr=oc+(posts[0]-EHMER_OFFSET)*linesper-(linesper>>1);
403 for(i=posts[0];i<post1;i++){
405 float lin=amp+curve[i];
406 if(seed[seedptr]<lin)seed[seedptr]=lin;
413 static void seed_loop(vorbis_look_psy *p,
414 const float ***curves,
419 vorbis_info_psy *vi=p->vi;
421 float dBoffset=vi->max_curve_dB-specmax;
423 /* prime the working vector with peak values */
427 long oc=p->octave[i];
428 while(i+1<n && p->octave[i+1]==oc){
430 if(f[i]>max)max=f[i];
436 if(oc>=P_BANDS)oc=P_BANDS-1;
442 p->octave[i]-p->firstoc,
443 p->total_octave_lines,
444 p->eighth_octave_lines,
450 static void seed_chase(float *seeds, int linesper, long n){
451 long *posstack=alloca(n*sizeof(*posstack));
452 float *ampstack=alloca(n*sizeof(*ampstack));
460 ampstack[stack++]=seeds[i];
463 if(seeds[i]<ampstack[stack-1]){
465 ampstack[stack++]=seeds[i];
468 if(i<posstack[stack-1]+linesper){
469 if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] &&
470 i<posstack[stack-2]+linesper){
471 /* we completely overlap, making stack-1 irrelevant. pop it */
477 ampstack[stack++]=seeds[i];
485 /* the stack now contains only the positions that are relevant. Scan
486 'em straight through */
488 for(i=0;i<stack;i++){
490 if(i<stack-1 && ampstack[i+1]>ampstack[i]){
491 endpos=posstack[i+1];
493 endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is
494 discarded in short frames */
496 if(endpos>n)endpos=n;
497 for(;pos<endpos;pos++)
498 seeds[pos]=ampstack[i];
501 /* there. Linear time. I now remember this was on a problem set I
502 had in Grad Skool... I didn't solve it at the time ;-) */
506 /* bleaugh, this is more complicated than it needs to be */
508 static void max_seeds(vorbis_look_psy *p,
511 long n=p->total_octave_lines;
512 int linesper=p->eighth_octave_lines;
516 seed_chase(seed,linesper,n); /* for masking */
518 pos=p->octave[0]-p->firstoc-(linesper>>1);
520 while(linpos+1<p->n){
521 float minV=seed[pos];
522 long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc;
523 if(minV>p->vi->tone_abs_limit)minV=p->vi->tone_abs_limit;
526 if((seed[pos]>NEGINF && seed[pos]<minV) || minV==NEGINF)
531 for(;linpos<p->n && p->octave[linpos]<=end;linpos++)
532 if(flr[linpos]<minV)flr[linpos]=minV;
536 float minV=seed[p->total_octave_lines-1];
537 for(;linpos<p->n;linpos++)
538 if(flr[linpos]<minV)flr[linpos]=minV;
543 static void bark_noise_hybridmp(int n,const long *b,
549 float *N=alloca(n*sizeof(*N));
550 float *X=alloca(n*sizeof(*N));
551 float *XX=alloca(n*sizeof(*N));
552 float *Y=alloca(n*sizeof(*N));
553 float *XY=alloca(n*sizeof(*N));
555 float tN, tX, tXX, tY, tXY;
565 tN = tX = tXX = tY = tXY = 0.f;
568 if (y < 1.f) y = 1.f;
582 for (i = 1, x = 1.f; i < n; i++, x += 1.f) {
585 if (y < 1.f) y = 1.f;
602 for (i = 0, x = 0.f; i < n; i++, x += 1.f) {
606 if( lo>=0 || -lo>=n ) break;
611 tXX = XX[hi] + XX[-lo];
613 tXY = XY[hi] - XY[-lo];
615 A = tY * tXX - tX * tXY;
616 B = tN * tXY - tX * tY;
617 D = tN * tXX - tX * tX;
619 if (R < 0.f) R = 0.f;
621 noise[i] = R - offset;
624 for ( ; i < n; i++, x += 1.f) {
628 if( lo<0 || lo>=n ) break;
633 tXX = XX[hi] - XX[lo];
635 tXY = XY[hi] - XY[lo];
637 A = tY * tXX - tX * tXY;
638 B = tN * tXY - tX * tY;
639 D = tN * tXX - tX * tX;
641 if (R < 0.f) R = 0.f;
643 noise[i] = R - offset;
646 for ( ; i < n; i++, x += 1.f) {
649 if (R < 0.f) R = 0.f;
651 noise[i] = R - offset;
654 if (fixed <= 0) return;
656 for (i = 0, x = 0.f; i < n; i++, x += 1.f) {
664 tXX = XX[hi] + XX[-lo];
666 tXY = XY[hi] - XY[-lo];
669 A = tY * tXX - tX * tXY;
670 B = tN * tXY - tX * tY;
671 D = tN * tXX - tX * tX;
674 if (R - offset < noise[i]) noise[i] = R - offset;
676 for ( ; i < n; i++, x += 1.f) {
685 tXX = XX[hi] - XX[lo];
687 tXY = XY[hi] - XY[lo];
689 A = tY * tXX - tX * tXY;
690 B = tN * tXY - tX * tY;
691 D = tN * tXX - tX * tX;
694 if (R - offset < noise[i]) noise[i] = R - offset;
696 for ( ; i < n; i++, x += 1.f) {
698 if (R - offset < noise[i]) noise[i] = R - offset;
702 void _vp_noisemask(vorbis_look_psy *p,
707 float *work=alloca(n*sizeof(*work));
709 bark_noise_hybridmp(n,p->bark,logmdct,logmask,
712 for(i=0;i<n;i++)work[i]=logmdct[i]-logmask[i];
714 bark_noise_hybridmp(n,p->bark,work,logmask,0.,
715 p->vi->noisewindowfixed);
717 for(i=0;i<n;i++)work[i]=logmdct[i]-work[i];
725 work2[i]=logmask[i]+work[i];
729 _analysis_output("median2R",seq/2,work,n,1,0,0);
731 _analysis_output("median2L",seq/2,work,n,1,0,0);
734 _analysis_output("envelope2R",seq/2,work2,n,1,0,0);
736 _analysis_output("envelope2L",seq/2,work2,n,1,0,0);
742 int dB=logmask[i]+.5;
743 if(dB>=NOISE_COMPAND_LEVELS)dB=NOISE_COMPAND_LEVELS-1;
745 logmask[i]= work[i]+p->vi->noisecompand[dB];
750 void _vp_tonemask(vorbis_look_psy *p,
753 float global_specmax,
754 float local_specmax){
758 float *seed=alloca(sizeof(*seed)*p->total_octave_lines);
759 float att=local_specmax+p->vi->ath_adjatt;
760 for(i=0;i<p->total_octave_lines;i++)seed[i]=NEGINF;
762 /* set the ATH (floating below localmax, not global max by a
764 if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt;
767 logmask[i]=p->ath[i]+att;
770 seed_loop(p,(const float ***)p->tonecurves,logfft,logmask,seed,global_specmax);
771 max_seeds(p,seed,logmask);
775 void _vp_offset_and_mix(vorbis_look_psy *p,
783 float de, coeffi, cx;/* AoTuV */
784 float toneatt=p->vi->tone_masteratt[offset_select];
789 float val= noise[i]+p->noiseoffset[offset_select][i];
790 if(val>p->vi->noisemaxsupp)val=p->vi->noisemaxsupp;
791 logmask[i]=max(val,tone[i]+toneatt);
796 The following codes improve a noise problem.
797 A fundamental idea uses the value of masking and carries out
798 the relative compensation of the MDCT.
799 However, this code is not perfect and all noise problems cannot be solved.
800 by Aoyumi @ 2004/04/18
803 if(offset_select == 1) {
804 coeffi = -17.2; /* coeffi is a -17.2dB threshold */
805 val = val - logmdct[i]; /* val == mdct line value relative to floor in dB */
808 /* mdct value is > -17.2 dB below floor */
810 de = 1.0-((val-coeffi)*0.005*cx);
811 /* pro-rated attenuation:
812 -0.00 dB boost if mdct value is -17.2dB (relative to floor)
813 -0.77 dB boost if mdct value is 0dB (relative to floor)
814 -1.64 dB boost if mdct value is +17.2dB (relative to floor)
817 if(de < 0) de = 0.0001;
819 /* mdct value is <= -17.2 dB below floor */
821 de = 1.0-((val-coeffi)*0.0003*cx);
822 /* pro-rated attenuation:
823 +0.00 dB atten if mdct value is -17.2dB (relative to floor)
824 +0.45 dB atten if mdct value is -34.4dB (relative to floor)
833 float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){
834 vorbis_info *vi=vd->vi;
835 codec_setup_info *ci=vi->codec_setup;
836 vorbis_info_psy_global *gi=&ci->psy_g_param;
838 int n=ci->blocksizes[vd->W]/2;
839 float secs=(float)n/vi->rate;
841 amp+=secs*gi->ampmax_att_per_sec;
842 if(amp<-9999)amp=-9999;
846 static float FLOOR1_fromdB_LOOKUP[256]={
847 1.0649863e-07F, 1.1341951e-07F, 1.2079015e-07F, 1.2863978e-07F,
848 1.3699951e-07F, 1.4590251e-07F, 1.5538408e-07F, 1.6548181e-07F,
849 1.7623575e-07F, 1.8768855e-07F, 1.9988561e-07F, 2.128753e-07F,
850 2.2670913e-07F, 2.4144197e-07F, 2.5713223e-07F, 2.7384213e-07F,
851 2.9163793e-07F, 3.1059021e-07F, 3.3077411e-07F, 3.5226968e-07F,
852 3.7516214e-07F, 3.9954229e-07F, 4.2550680e-07F, 4.5315863e-07F,
853 4.8260743e-07F, 5.1396998e-07F, 5.4737065e-07F, 5.8294187e-07F,
854 6.2082472e-07F, 6.6116941e-07F, 7.0413592e-07F, 7.4989464e-07F,
855 7.9862701e-07F, 8.5052630e-07F, 9.0579828e-07F, 9.6466216e-07F,
856 1.0273513e-06F, 1.0941144e-06F, 1.1652161e-06F, 1.2409384e-06F,
857 1.3215816e-06F, 1.4074654e-06F, 1.4989305e-06F, 1.5963394e-06F,
858 1.7000785e-06F, 1.8105592e-06F, 1.9282195e-06F, 2.0535261e-06F,
859 2.1869758e-06F, 2.3290978e-06F, 2.4804557e-06F, 2.6416497e-06F,
860 2.8133190e-06F, 2.9961443e-06F, 3.1908506e-06F, 3.3982101e-06F,
861 3.6190449e-06F, 3.8542308e-06F, 4.1047004e-06F, 4.3714470e-06F,
862 4.6555282e-06F, 4.9580707e-06F, 5.2802740e-06F, 5.6234160e-06F,
863 5.9888572e-06F, 6.3780469e-06F, 6.7925283e-06F, 7.2339451e-06F,
864 7.7040476e-06F, 8.2047000e-06F, 8.7378876e-06F, 9.3057248e-06F,
865 9.9104632e-06F, 1.0554501e-05F, 1.1240392e-05F, 1.1970856e-05F,
866 1.2748789e-05F, 1.3577278e-05F, 1.4459606e-05F, 1.5399272e-05F,
867 1.6400004e-05F, 1.7465768e-05F, 1.8600792e-05F, 1.9809576e-05F,
868 2.1096914e-05F, 2.2467911e-05F, 2.3928002e-05F, 2.5482978e-05F,
869 2.7139006e-05F, 2.8902651e-05F, 3.0780908e-05F, 3.2781225e-05F,
870 3.4911534e-05F, 3.7180282e-05F, 3.9596466e-05F, 4.2169667e-05F,
871 4.4910090e-05F, 4.7828601e-05F, 5.0936773e-05F, 5.4246931e-05F,
872 5.7772202e-05F, 6.1526565e-05F, 6.5524908e-05F, 6.9783085e-05F,
873 7.4317983e-05F, 7.9147585e-05F, 8.4291040e-05F, 8.9768747e-05F,
874 9.5602426e-05F, 0.00010181521F, 0.00010843174F, 0.00011547824F,
875 0.00012298267F, 0.00013097477F, 0.00013948625F, 0.00014855085F,
876 0.00015820453F, 0.00016848555F, 0.00017943469F, 0.00019109536F,
877 0.00020351382F, 0.00021673929F, 0.00023082423F, 0.00024582449F,
878 0.00026179955F, 0.00027881276F, 0.00029693158F, 0.00031622787F,
879 0.00033677814F, 0.00035866388F, 0.00038197188F, 0.00040679456F,
880 0.00043323036F, 0.00046138411F, 0.00049136745F, 0.00052329927F,
881 0.00055730621F, 0.00059352311F, 0.00063209358F, 0.00067317058F,
882 0.00071691700F, 0.00076350630F, 0.00081312324F, 0.00086596457F,
883 0.00092223983F, 0.00098217216F, 0.0010459992F, 0.0011139742F,
884 0.0011863665F, 0.0012634633F, 0.0013455702F, 0.0014330129F,
885 0.0015261382F, 0.0016253153F, 0.0017309374F, 0.0018434235F,
886 0.0019632195F, 0.0020908006F, 0.0022266726F, 0.0023713743F,
887 0.0025254795F, 0.0026895994F, 0.0028643847F, 0.0030505286F,
888 0.0032487691F, 0.0034598925F, 0.0036847358F, 0.0039241906F,
889 0.0041792066F, 0.0044507950F, 0.0047400328F, 0.0050480668F,
890 0.0053761186F, 0.0057254891F, 0.0060975636F, 0.0064938176F,
891 0.0069158225F, 0.0073652516F, 0.0078438871F, 0.0083536271F,
892 0.0088964928F, 0.009474637F, 0.010090352F, 0.010746080F,
893 0.011444421F, 0.012188144F, 0.012980198F, 0.013823725F,
894 0.014722068F, 0.015678791F, 0.016697687F, 0.017782797F,
895 0.018938423F, 0.020169149F, 0.021479854F, 0.022875735F,
896 0.024362330F, 0.025945531F, 0.027631618F, 0.029427276F,
897 0.031339626F, 0.033376252F, 0.035545228F, 0.037855157F,
898 0.040315199F, 0.042935108F, 0.045725273F, 0.048696758F,
899 0.051861348F, 0.055231591F, 0.058820850F, 0.062643361F,
900 0.066714279F, 0.071049749F, 0.075666962F, 0.080584227F,
901 0.085821044F, 0.091398179F, 0.097337747F, 0.10366330F,
902 0.11039993F, 0.11757434F, 0.12521498F, 0.13335215F,
903 0.14201813F, 0.15124727F, 0.16107617F, 0.17154380F,
904 0.18269168F, 0.19456402F, 0.20720788F, 0.22067342F,
905 0.23501402F, 0.25028656F, 0.26655159F, 0.28387361F,
906 0.30232132F, 0.32196786F, 0.34289114F, 0.36517414F,
907 0.38890521F, 0.41417847F, 0.44109412F, 0.46975890F,
908 0.50028648F, 0.53279791F, 0.56742212F, 0.60429640F,
909 0.64356699F, 0.68538959F, 0.72993007F, 0.77736504F,
910 0.82788260F, 0.88168307F, 0.9389798F, 1.F,
913 /* this is for per-channel noise normalization */
914 static int apsort(const void *a, const void *b){
915 float f1=**(float**)a;
916 float f2=**(float**)b;
917 return (f1<f2)-(f1>f2);
920 static void flag_lossless(int limit, float prepoint, float postpoint, float *mdct,
921 float *floor, int *flag, int i, int jn){
924 float point = j>=limit-i ? postpoint : prepoint;
925 float r = fabs(mdct[j])/floor[j];
933 /* Overload/Side effect: On input, the *q vector holds either the
934 quantized energy (for elements with the flag set) or the absolute
935 values of the *r vector (for elements with flag unset). On output,
936 *q holds the quantized energy for all elements */
937 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){
939 vorbis_info_psy *vi=p->vi;
940 float **sort = alloca(n*sizeof(*sort));
942 int start = (vi->normal_p ? vi->normal_start-i : n);
945 /* force classic behavior where only energy in the current band is considered */
948 /* still responsible for populating *out where noise norm not in
949 effect. There's no need to [re]populate *q in these areas */
950 for(j=0;j<start;j++){
951 if(!flags || !flags[j]){ /* lossless coupling already quantized.
952 Don't touch; requantizing based on
953 energy would be incorrect. */
954 float ve = q[j]/f[j];
956 out[j] = -rint(sqrt(ve));
958 out[j] = rint(sqrt(ve));
962 /* sort magnitudes for noise norm portion of partition */
964 if(!flags || !flags[j]){ /* can't noise norm elements that have
965 already been loslessly coupled; we can
966 only account for their energy error */
967 float ve = q[j]/f[j];
968 /* Despite all the new, more capable coupling code, for now we
969 implement noise norm as it has been up to this point. Only
970 consider promotions to unit magnitude from 0. In addition
971 the only energy error counted is quantizations to zero. */
972 /* also-- the original point code only applied noise norm at > pointlimit */
973 if(ve<.25f && (!flags || j>=limit-i)){
975 sort[count++]=q+j; /* q is fabs(r) for unflagged element */
977 /* For now: no acc adjustment for nonzero quantization. populate *out and q as this value is final. */
979 out[j] = -rint(sqrt(ve));
981 out[j] = rint(sqrt(ve));
982 q[j] = out[j]*out[j]*f[j];
985 again, no energy adjustment for error in nonzero quant-- for now
990 /* noise norm to do */
991 qsort(sort,count,sizeof(*sort),apsort);
992 for(j=0;j<count;j++){
994 if(acc>=vi->normal_thresh){
995 out[k]=unitnorm(r[k]);
1008 /* Noise normalization, quantization and coupling are not wholly
1009 seperable processes in depth>1 coupling. */
1010 void _vp_couple_quantize_normalize(int blobno,
1011 vorbis_info_psy_global *g,
1013 vorbis_info_mapping0 *vi,
1017 int sliding_lowpass,
1022 int partition=(p->vi->normal_p ? p->vi->normal_partition : 16);
1023 int limit = g->coupling_pointlimit[p->vi->blockflag][blobno];
1024 float prepoint=stereo_threshholds[g->coupling_prepointamp[blobno]];
1025 float postpoint=stereo_threshholds[g->coupling_postpointamp[blobno]];
1027 float de=0.1*p->m_val; /* a blend of the AoTuV M2 and M3 code here and below */
1030 /* mdct is our raw mdct output, floor not removed. */
1031 /* inout passes in the ifloor, passes back quantized result */
1033 /* unquantized energy (negative indicates amplitude has negative sign) */
1034 float **raw = alloca(ch*sizeof(*raw));
1036 /* dual pupose; quantized energy (if flag set), othersize fabs(raw) */
1037 float **quant = alloca(ch*sizeof(*quant));
1040 float **floor = alloca(ch*sizeof(*floor));
1042 /* flags indicating raw/quantized status of elements in raw vector */
1043 int **flag = alloca(ch*sizeof(*flag));
1045 /* non-zero flag working vector */
1046 int *nz = alloca(ch*sizeof(*nz));
1048 /* energy surplus/defecit tracking */
1049 float *acc = alloca((ch+vi->coupling_steps)*sizeof(*acc));
1051 /* The threshold of a stereo is changed with the size of n */
1053 postpoint=stereo_threshholds_limited[g->coupling_postpointamp[blobno]];
1055 raw[0] = alloca(ch*partition*sizeof(**raw));
1056 quant[0] = alloca(ch*partition*sizeof(**quant));
1057 floor[0] = alloca(ch*partition*sizeof(**floor));
1058 flag[0] = alloca(ch*partition*sizeof(**flag));
1061 raw[i] = &raw[0][partition*i];
1062 quant[i] = &quant[0][partition*i];
1063 floor[i] = &floor[0][partition*i];
1064 flag[i] = &flag[0][partition*i];
1066 for(i=0;i<ch+vi->coupling_steps;i++)
1069 for(i=0;i<n;i+=partition){
1070 int k,j,jn = partition > n-i ? n-i : partition;
1073 memcpy(nz,nonzero,sizeof(*nz)*ch);
1076 memset(flag[0],0,ch*partition*sizeof(**flag));
1078 int *iout = &iwork[k][i];
1082 floor[k][j] = FLOOR1_fromdB_LOOKUP[iout[j]];
1084 flag_lossless(limit,prepoint,postpoint,&mdct[k][i],floor[k],flag[k],i,jn);
1087 quant[k][j] = raw[k][j] = mdct[k][i+j]*mdct[k][i+j];
1088 if(mdct[k][i+j]<0.f) raw[k][j]*=-1.f;
1089 floor[k][j]*=floor[k][j];
1092 acc[track]=noise_normalize(p,limit,raw[k],quant[k],floor[k],NULL,acc[track],i,jn,iout);
1096 floor[k][j] = 1e-10f;
1108 for(step=0;step<vi->coupling_steps;step++){
1109 int Mi = vi->coupling_mag[step];
1110 int Ai = vi->coupling_ang[step];
1111 int *iM = &iwork[Mi][i];
1112 int *iA = &iwork[Ai][i];
1113 float *reM = raw[Mi];
1114 float *reA = raw[Ai];
1115 float *qeM = quant[Mi];
1116 float *qeA = quant[Ai];
1117 float *floorM = floor[Mi];
1118 float *floorA = floor[Ai];
1122 if(nz[Mi] || nz[Ai]){
1123 nz[Mi] = nz[Ai] = 1;
1127 if(j<sliding_lowpass-i){
1129 /* lossless coupling */
1131 reM[j] = fabs(reM[j])+fabs(reA[j]);
1132 qeM[j] = qeM[j]+qeA[j];
1141 iA[j]=(A>0?A-B:B-A);
1143 iA[j]=(B>0?A-B:B-A);
1147 /* collapse two equivalent tuples to one */
1148 if(iA[j]>=abs(iM[j])*2){
1156 /* lossy (point) coupling */
1160 qeM[j] = fabs(reM[j]);
1165 The boost problem by the combination of noise normalization and point stereo is eased.
1166 However, this is a temporary patch.
1167 by Aoyumi @ 2004/04/18
1169 float derate = (1.0 - de*((float)(j-limit+i) / (float)(n-limit)));
1171 if(reM[j]+reA[j]<0){
1172 reM[j] = - (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate);
1174 reM[j] = (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate);
1178 if(reM[j]+reA[j]<0){
1179 reM[j] = - (qeM[j] = fabs(reM[j])+fabs(reA[j]));
1181 reM[j] = (qeM[j] = fabs(reM[j])+fabs(reA[j]));
1191 floorM[j]=floorA[j]=floorM[j]+floorA[j];
1193 /* normalize the resulting mag vector */
1194 acc[track]=noise_normalize(p,limit,raw[Mi],quant[Mi],floor[Mi],flag[Mi],acc[track],i,jn,iM);
1200 for(i=0;i<vi->coupling_steps;i++){
1201 /* make sure coupling a zero and a nonzero channel results in two
1202 nonzero channels. */
1203 if(nonzero[vi->coupling_mag[i]] ||
1204 nonzero[vi->coupling_ang[i]]){
1205 nonzero[vi->coupling_mag[i]]=1;
1206 nonzero[vi->coupling_ang[i]]=1;