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 http://www.xiph.org/ *
11 ********************************************************************
13 function: psychoacoustics not including preecho
16 ********************************************************************/
21 #include "vorbis/codec.h"
22 #include "codec_internal.h"
32 #define NEGINF -9999.f
33 static const double stereo_threshholds[]={0.0, .5, 1.0, 1.5, 2.5, 4.5, 8.5, 16.5, 9e10};
34 static const double stereo_threshholds_limited[]={0.0, .5, 1.0, 1.5, 2.0, 2.5, 4.5, 8.5, 9e10};
36 vorbis_look_psy_global *_vp_global_look(vorbis_info *vi){
37 codec_setup_info *ci=vi->codec_setup;
38 vorbis_info_psy_global *gi=&ci->psy_g_param;
39 vorbis_look_psy_global *look=_ogg_calloc(1,sizeof(*look));
41 look->channels=vi->channels;
48 void _vp_global_free(vorbis_look_psy_global *look){
50 memset(look,0,sizeof(*look));
55 void _vi_gpsy_free(vorbis_info_psy_global *i){
57 memset(i,0,sizeof(*i));
62 void _vi_psy_free(vorbis_info_psy *i){
64 memset(i,0,sizeof(*i));
69 static void min_curve(float *c,
72 for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i];
74 static void max_curve(float *c,
77 for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i];
80 static void attenuate_curve(float *c,float att){
82 for(i=0;i<EHMER_MAX;i++)
86 static float ***setup_tone_curves(float curveatt_dB[P_BANDS],float binHz,int n,
87 float center_boost, float center_decay_rate){
90 float workc[P_BANDS][P_LEVELS][EHMER_MAX];
91 float athc[P_LEVELS][EHMER_MAX];
92 float *brute_buffer=alloca(n*sizeof(*brute_buffer));
94 float ***ret=_ogg_malloc(sizeof(*ret)*P_BANDS);
96 memset(workc,0,sizeof(workc));
98 for(i=0;i<P_BANDS;i++){
99 /* we add back in the ATH to avoid low level curves falling off to
100 -infinity and unnecessarily cutting off high level curves in the
101 curve limiting (last step). */
103 /* A half-band's settings must be valid over the whole band, and
104 it's better to mask too little than too much */
106 for(j=0;j<EHMER_MAX;j++){
109 if(j+k+ath_offset<MAX_ATH){
110 if(min>ATH[j+k+ath_offset])min=ATH[j+k+ath_offset];
112 if(min>ATH[MAX_ATH-1])min=ATH[MAX_ATH-1];
117 /* copy curves into working space, replicate the 50dB curve to 30
118 and 40, replicate the 100dB curve to 110 */
120 memcpy(workc[i][j+2],tonemasks[i][j],EHMER_MAX*sizeof(*tonemasks[i][j]));
121 memcpy(workc[i][0],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0]));
122 memcpy(workc[i][1],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0]));
124 /* apply centered curve boost/decay */
125 for(j=0;j<P_LEVELS;j++){
126 for(k=0;k<EHMER_MAX;k++){
127 float adj=center_boost+abs(EHMER_OFFSET-k)*center_decay_rate;
128 if(adj<0. && center_boost>0)adj=0.;
129 if(adj>0. && center_boost<0)adj=0.;
134 /* normalize curves so the driving amplitude is 0dB */
135 /* make temp curves with the ATH overlayed */
136 for(j=0;j<P_LEVELS;j++){
137 attenuate_curve(workc[i][j],curveatt_dB[i]+100.-(j<2?2:j)*10.-P_LEVEL_0);
138 memcpy(athc[j],ath,EHMER_MAX*sizeof(**athc));
139 attenuate_curve(athc[j],+100.-j*10.f-P_LEVEL_0);
140 max_curve(athc[j],workc[i][j]);
143 /* Now limit the louder curves.
145 the idea is this: We don't know what the playback attenuation
146 will be; 0dB SL moves every time the user twiddles the volume
147 knob. So that means we have to use a single 'most pessimal' curve
148 for all masking amplitudes, right? Wrong. The *loudest* sound
149 can be in (we assume) a range of ...+100dB] SL. However, sounds
150 20dB down will be in a range ...+80], 40dB down is from ...+60],
153 for(j=1;j<P_LEVELS;j++){
154 min_curve(athc[j],athc[j-1]);
155 min_curve(workc[i][j],athc[j]);
159 for(i=0;i<P_BANDS;i++){
160 int hi_curve,lo_curve,bin;
161 ret[i]=_ogg_malloc(sizeof(**ret)*P_LEVELS);
163 /* low frequency curves are measured with greater resolution than
164 the MDCT/FFT will actually give us; we want the curve applied
165 to the tone data to be pessimistic and thus apply the minimum
166 masking possible for a given bin. That means that a single bin
167 could span more than one octave and that the curve will be a
168 composite of multiple octaves. It also may mean that a single
169 bin may span > an eighth of an octave and that the eighth
170 octave values may also be composited. */
172 /* which octave curves will we be compositing? */
173 bin=floor(fromOC(i*.5)/binHz);
174 lo_curve= ceil(toOC(bin*binHz+1)*2);
175 hi_curve= floor(toOC((bin+1)*binHz)*2);
176 if(lo_curve>i)lo_curve=i;
177 if(lo_curve<0)lo_curve=0;
178 if(hi_curve>=P_BANDS)hi_curve=P_BANDS-1;
180 for(m=0;m<P_LEVELS;m++){
181 ret[i][m]=_ogg_malloc(sizeof(***ret)*(EHMER_MAX+2));
183 for(j=0;j<n;j++)brute_buffer[j]=999.;
185 /* render the curve into bins, then pull values back into curve.
186 The point is that any inherent subsampling aliasing results in
188 for(k=lo_curve;k<=hi_curve;k++){
191 for(j=0;j<EHMER_MAX;j++){
192 int lo_bin= fromOC(j*.125+k*.5-2.0625)/binHz;
193 int hi_bin= fromOC(j*.125+k*.5-1.9375)/binHz+1;
195 if(lo_bin<0)lo_bin=0;
196 if(lo_bin>n)lo_bin=n;
197 if(lo_bin<l)l=lo_bin;
198 if(hi_bin<0)hi_bin=0;
199 if(hi_bin>n)hi_bin=n;
201 for(;l<hi_bin && l<n;l++)
202 if(brute_buffer[l]>workc[k][m][j])
203 brute_buffer[l]=workc[k][m][j];
207 if(brute_buffer[l]>workc[k][m][EHMER_MAX-1])
208 brute_buffer[l]=workc[k][m][EHMER_MAX-1];
212 /* be equally paranoid about being valid up to next half ocatve */
216 for(j=0;j<EHMER_MAX;j++){
217 int lo_bin= fromOC(j*.125+i*.5-2.0625)/binHz;
218 int hi_bin= fromOC(j*.125+i*.5-1.9375)/binHz+1;
220 if(lo_bin<0)lo_bin=0;
221 if(lo_bin>n)lo_bin=n;
222 if(lo_bin<l)l=lo_bin;
223 if(hi_bin<0)hi_bin=0;
224 if(hi_bin>n)hi_bin=n;
226 for(;l<hi_bin && l<n;l++)
227 if(brute_buffer[l]>workc[k][m][j])
228 brute_buffer[l]=workc[k][m][j];
232 if(brute_buffer[l]>workc[k][m][EHMER_MAX-1])
233 brute_buffer[l]=workc[k][m][EHMER_MAX-1];
238 for(j=0;j<EHMER_MAX;j++){
239 int bin=fromOC(j*.125+i*.5-2.)/binHz;
241 ret[i][m][j+2]=-999.;
244 ret[i][m][j+2]=-999.;
246 ret[i][m][j+2]=brute_buffer[bin];
252 for(j=0;j<EHMER_OFFSET;j++)
253 if(ret[i][m][j+2]>-200.f)break;
256 for(j=EHMER_MAX-1;j>EHMER_OFFSET+1;j--)
257 if(ret[i][m][j+2]>-200.f)
267 void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi,
268 vorbis_info_psy_global *gi,int n,long rate){
269 long i,j,lo=-99,hi=1;
271 memset(p,0,sizeof(*p));
273 p->eighth_octave_lines=gi->eighth_octave_lines;
274 p->shiftoc=rint(log(gi->eighth_octave_lines*8.f)/log(2.f))-1;
276 p->firstoc=toOC(.25f*rate*.5/n)*(1<<(p->shiftoc+1))-gi->eighth_octave_lines;
277 maxoc=toOC((n+.25f)*rate*.5/n)*(1<<(p->shiftoc+1))+.5f;
278 p->total_octave_lines=maxoc-p->firstoc+1;
279 p->ath=_ogg_malloc(n*sizeof(*p->ath));
281 p->octave=_ogg_malloc(n*sizeof(*p->octave));
282 p->bark=_ogg_malloc(n*sizeof(*p->bark));
287 /* AoTuV HF weighting */
289 if(rate < 26000) p->m_val = 0;
290 else if(rate < 38000) p->m_val = .94; /* 32kHz */
291 else if(rate > 46000) p->m_val = 1.275; /* 48kHz */
293 /* set up the lookups for a given blocksize and sample rate */
295 for(i=0,j=0;i<MAX_ATH-1;i++){
296 int endpos=rint(fromOC((i+1)*.125-2.)*2*n/rate);
299 float delta=(ATH[i+1]-base)/(endpos-j);
300 for(;j<endpos && j<n;j++){
308 p->ath[j]=p->ath[j-1];
312 float bark=toBARK(rate/(2*n)*i);
314 for(;lo+vi->noisewindowlomin<i &&
315 toBARK(rate/(2*n)*lo)<(bark-vi->noisewindowlo);lo++);
317 for(;hi<=n && (hi<i+vi->noisewindowhimin ||
318 toBARK(rate/(2*n)*hi)<(bark+vi->noisewindowhi));hi++);
320 p->bark[i]=((lo-1)<<16)+(hi-1);
325 p->octave[i]=toOC((i+.25f)*.5*rate/n)*(1<<(p->shiftoc+1))+.5f;
327 p->tonecurves=setup_tone_curves(vi->toneatt,rate*.5/n,n,
328 vi->tone_centerboost,vi->tone_decay);
330 /* set up rolling noise median */
331 p->noiseoffset=_ogg_malloc(P_NOISECURVES*sizeof(*p->noiseoffset));
332 for(i=0;i<P_NOISECURVES;i++)
333 p->noiseoffset[i]=_ogg_malloc(n*sizeof(**p->noiseoffset));
336 float halfoc=toOC((i+.5)*rate/(2.*n))*2.;
340 if(halfoc<0)halfoc=0;
341 if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1;
342 inthalfoc=(int)halfoc;
343 del=halfoc-inthalfoc;
345 for(j=0;j<P_NOISECURVES;j++)
346 p->noiseoffset[j][i]=
347 p->vi->noiseoff[j][inthalfoc]*(1.-del) +
348 p->vi->noiseoff[j][inthalfoc+1]*del;
354 _analysis_output_always("noiseoff0",ls,p->noiseoffset[0],n,1,0,0);
355 _analysis_output_always("noiseoff1",ls,p->noiseoffset[1],n,1,0,0);
356 _analysis_output_always("noiseoff2",ls++,p->noiseoffset[2],n,1,0,0);
361 void _vp_psy_clear(vorbis_look_psy *p){
364 if(p->ath)_ogg_free(p->ath);
365 if(p->octave)_ogg_free(p->octave);
366 if(p->bark)_ogg_free(p->bark);
368 for(i=0;i<P_BANDS;i++){
369 for(j=0;j<P_LEVELS;j++){
370 _ogg_free(p->tonecurves[i][j]);
372 _ogg_free(p->tonecurves[i]);
374 _ogg_free(p->tonecurves);
377 for(i=0;i<P_NOISECURVES;i++){
378 _ogg_free(p->noiseoffset[i]);
380 _ogg_free(p->noiseoffset);
382 memset(p,0,sizeof(*p));
386 /* octave/(8*eighth_octave_lines) x scale and dB y scale */
387 static void seed_curve(float *seed,
388 const float **curves,
391 int linesper,float dBoffset){
394 const float *posts,*curve;
396 int choice=(int)((amp+dBoffset-P_LEVEL_0)*.1f);
397 choice=max(choice,0);
398 choice=min(choice,P_LEVELS-1);
399 posts=curves[choice];
402 seedptr=oc+(posts[0]-EHMER_OFFSET)*linesper-(linesper>>1);
404 for(i=posts[0];i<post1;i++){
406 float lin=amp+curve[i];
407 if(seed[seedptr]<lin)seed[seedptr]=lin;
414 static void seed_loop(vorbis_look_psy *p,
415 const float ***curves,
420 vorbis_info_psy *vi=p->vi;
422 float dBoffset=vi->max_curve_dB-specmax;
424 /* prime the working vector with peak values */
428 long oc=p->octave[i];
429 while(i+1<n && p->octave[i+1]==oc){
431 if(f[i]>max)max=f[i];
437 if(oc>=P_BANDS)oc=P_BANDS-1;
443 p->octave[i]-p->firstoc,
444 p->total_octave_lines,
445 p->eighth_octave_lines,
451 static void seed_chase(float *seeds, int linesper, long n){
452 long *posstack=alloca(n*sizeof(*posstack));
453 float *ampstack=alloca(n*sizeof(*ampstack));
461 ampstack[stack++]=seeds[i];
464 if(seeds[i]<ampstack[stack-1]){
466 ampstack[stack++]=seeds[i];
469 if(i<posstack[stack-1]+linesper){
470 if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] &&
471 i<posstack[stack-2]+linesper){
472 /* we completely overlap, making stack-1 irrelevant. pop it */
478 ampstack[stack++]=seeds[i];
486 /* the stack now contains only the positions that are relevant. Scan
487 'em straight through */
489 for(i=0;i<stack;i++){
491 if(i<stack-1 && ampstack[i+1]>ampstack[i]){
492 endpos=posstack[i+1];
494 endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is
495 discarded in short frames */
497 if(endpos>n)endpos=n;
498 for(;pos<endpos;pos++)
499 seeds[pos]=ampstack[i];
502 /* there. Linear time. I now remember this was on a problem set I
503 had in Grad Skool... I didn't solve it at the time ;-) */
507 /* bleaugh, this is more complicated than it needs to be */
509 static void max_seeds(vorbis_look_psy *p,
512 long n=p->total_octave_lines;
513 int linesper=p->eighth_octave_lines;
517 seed_chase(seed,linesper,n); /* for masking */
519 pos=p->octave[0]-p->firstoc-(linesper>>1);
521 while(linpos+1<p->n){
522 float minV=seed[pos];
523 long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc;
524 if(minV>p->vi->tone_abs_limit)minV=p->vi->tone_abs_limit;
527 if((seed[pos]>NEGINF && seed[pos]<minV) || minV==NEGINF)
532 for(;linpos<p->n && p->octave[linpos]<=end;linpos++)
533 if(flr[linpos]<minV)flr[linpos]=minV;
537 float minV=seed[p->total_octave_lines-1];
538 for(;linpos<p->n;linpos++)
539 if(flr[linpos]<minV)flr[linpos]=minV;
544 static void bark_noise_hybridmp(int n,const long *b,
550 float *N=alloca(n*sizeof(*N));
551 float *X=alloca(n*sizeof(*N));
552 float *XX=alloca(n*sizeof(*N));
553 float *Y=alloca(n*sizeof(*N));
554 float *XY=alloca(n*sizeof(*N));
556 float tN, tX, tXX, tY, tXY;
566 tN = tX = tXX = tY = tXY = 0.f;
569 if (y < 1.f) y = 1.f;
583 for (i = 1, x = 1.f; i < n; i++, x += 1.f) {
586 if (y < 1.f) y = 1.f;
603 for (i = 0, x = 0.f;; i++, x += 1.f) {
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;
622 noise[i] = R - offset;
625 for ( ;; i++, x += 1.f) {
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;
645 for ( ; i < n; i++, x += 1.f) {
648 if (R < 0.f) R = 0.f;
650 noise[i] = R - offset;
653 if (fixed <= 0) return;
655 for (i = 0, x = 0.f;; i++, x += 1.f) {
662 tXX = XX[hi] + XX[-lo];
664 tXY = XY[hi] - XY[-lo];
667 A = tY * tXX - tX * tXY;
668 B = tN * tXY - tX * tY;
669 D = tN * tXX - tX * tX;
672 if (R - offset < noise[i]) noise[i] = R - offset;
674 for ( ;; i++, x += 1.f) {
682 tXX = XX[hi] - XX[lo];
684 tXY = XY[hi] - XY[lo];
686 A = tY * tXX - tX * tXY;
687 B = tN * tXY - tX * tY;
688 D = tN * tXX - tX * tX;
691 if (R - offset < noise[i]) noise[i] = R - offset;
693 for ( ; i < n; i++, x += 1.f) {
695 if (R - offset < noise[i]) noise[i] = R - offset;
699 void _vp_noisemask(vorbis_look_psy *p,
704 float *work=alloca(n*sizeof(*work));
706 bark_noise_hybridmp(n,p->bark,logmdct,logmask,
709 for(i=0;i<n;i++)work[i]=logmdct[i]-logmask[i];
711 bark_noise_hybridmp(n,p->bark,work,logmask,0.,
712 p->vi->noisewindowfixed);
714 for(i=0;i<n;i++)work[i]=logmdct[i]-work[i];
722 work2[i]=logmask[i]+work[i];
726 _analysis_output("median2R",seq/2,work,n,1,0,0);
728 _analysis_output("median2L",seq/2,work,n,1,0,0);
731 _analysis_output("envelope2R",seq/2,work2,n,1,0,0);
733 _analysis_output("envelope2L",seq/2,work2,n,1,0,0);
739 int dB=logmask[i]+.5;
740 if(dB>=NOISE_COMPAND_LEVELS)dB=NOISE_COMPAND_LEVELS-1;
742 logmask[i]= work[i]+p->vi->noisecompand[dB];
747 void _vp_tonemask(vorbis_look_psy *p,
750 float global_specmax,
751 float local_specmax){
755 float *seed=alloca(sizeof(*seed)*p->total_octave_lines);
756 float att=local_specmax+p->vi->ath_adjatt;
757 for(i=0;i<p->total_octave_lines;i++)seed[i]=NEGINF;
759 /* set the ATH (floating below localmax, not global max by a
761 if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt;
764 logmask[i]=p->ath[i]+att;
767 seed_loop(p,(const float ***)p->tonecurves,logfft,logmask,seed,global_specmax);
768 max_seeds(p,seed,logmask);
772 void _vp_offset_and_mix(vorbis_look_psy *p,
780 float de, coeffi, cx;/* AoTuV */
781 float toneatt=p->vi->tone_masteratt[offset_select];
786 float val= noise[i]+p->noiseoffset[offset_select][i];
787 if(val>p->vi->noisemaxsupp)val=p->vi->noisemaxsupp;
788 logmask[i]=max(val,tone[i]+toneatt);
793 The following codes improve a noise problem.
794 A fundamental idea uses the value of masking and carries out
795 the relative compensation of the MDCT.
796 However, this code is not perfect and all noise problems cannot be solved.
797 by Aoyumi @ 2004/04/18
800 if(offset_select == 1) {
801 coeffi = -17.2; /* coeffi is a -17.2dB threshold */
802 val = val - logmdct[i]; /* val == mdct line value relative to floor in dB */
805 /* mdct value is > -17.2 dB below floor */
807 de = 1.0-((val-coeffi)*0.005*cx);
808 /* pro-rated attenuation:
809 -0.00 dB boost if mdct value is -17.2dB (relative to floor)
810 -0.77 dB boost if mdct value is 0dB (relative to floor)
811 -1.64 dB boost if mdct value is +17.2dB (relative to floor)
814 if(de < 0) de = 0.0001;
816 /* mdct value is <= -17.2 dB below floor */
818 de = 1.0-((val-coeffi)*0.0003*cx);
819 /* pro-rated attenuation:
820 +0.00 dB atten if mdct value is -17.2dB (relative to floor)
821 +0.45 dB atten if mdct value is -34.4dB (relative to floor)
830 float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){
831 vorbis_info *vi=vd->vi;
832 codec_setup_info *ci=vi->codec_setup;
833 vorbis_info_psy_global *gi=&ci->psy_g_param;
835 int n=ci->blocksizes[vd->W]/2;
836 float secs=(float)n/vi->rate;
838 amp+=secs*gi->ampmax_att_per_sec;
839 if(amp<-9999)amp=-9999;
843 static float FLOOR1_fromdB_LOOKUP[256]={
844 1.0649863e-07F, 1.1341951e-07F, 1.2079015e-07F, 1.2863978e-07F,
845 1.3699951e-07F, 1.4590251e-07F, 1.5538408e-07F, 1.6548181e-07F,
846 1.7623575e-07F, 1.8768855e-07F, 1.9988561e-07F, 2.128753e-07F,
847 2.2670913e-07F, 2.4144197e-07F, 2.5713223e-07F, 2.7384213e-07F,
848 2.9163793e-07F, 3.1059021e-07F, 3.3077411e-07F, 3.5226968e-07F,
849 3.7516214e-07F, 3.9954229e-07F, 4.2550680e-07F, 4.5315863e-07F,
850 4.8260743e-07F, 5.1396998e-07F, 5.4737065e-07F, 5.8294187e-07F,
851 6.2082472e-07F, 6.6116941e-07F, 7.0413592e-07F, 7.4989464e-07F,
852 7.9862701e-07F, 8.5052630e-07F, 9.0579828e-07F, 9.6466216e-07F,
853 1.0273513e-06F, 1.0941144e-06F, 1.1652161e-06F, 1.2409384e-06F,
854 1.3215816e-06F, 1.4074654e-06F, 1.4989305e-06F, 1.5963394e-06F,
855 1.7000785e-06F, 1.8105592e-06F, 1.9282195e-06F, 2.0535261e-06F,
856 2.1869758e-06F, 2.3290978e-06F, 2.4804557e-06F, 2.6416497e-06F,
857 2.8133190e-06F, 2.9961443e-06F, 3.1908506e-06F, 3.3982101e-06F,
858 3.6190449e-06F, 3.8542308e-06F, 4.1047004e-06F, 4.3714470e-06F,
859 4.6555282e-06F, 4.9580707e-06F, 5.2802740e-06F, 5.6234160e-06F,
860 5.9888572e-06F, 6.3780469e-06F, 6.7925283e-06F, 7.2339451e-06F,
861 7.7040476e-06F, 8.2047000e-06F, 8.7378876e-06F, 9.3057248e-06F,
862 9.9104632e-06F, 1.0554501e-05F, 1.1240392e-05F, 1.1970856e-05F,
863 1.2748789e-05F, 1.3577278e-05F, 1.4459606e-05F, 1.5399272e-05F,
864 1.6400004e-05F, 1.7465768e-05F, 1.8600792e-05F, 1.9809576e-05F,
865 2.1096914e-05F, 2.2467911e-05F, 2.3928002e-05F, 2.5482978e-05F,
866 2.7139006e-05F, 2.8902651e-05F, 3.0780908e-05F, 3.2781225e-05F,
867 3.4911534e-05F, 3.7180282e-05F, 3.9596466e-05F, 4.2169667e-05F,
868 4.4910090e-05F, 4.7828601e-05F, 5.0936773e-05F, 5.4246931e-05F,
869 5.7772202e-05F, 6.1526565e-05F, 6.5524908e-05F, 6.9783085e-05F,
870 7.4317983e-05F, 7.9147585e-05F, 8.4291040e-05F, 8.9768747e-05F,
871 9.5602426e-05F, 0.00010181521F, 0.00010843174F, 0.00011547824F,
872 0.00012298267F, 0.00013097477F, 0.00013948625F, 0.00014855085F,
873 0.00015820453F, 0.00016848555F, 0.00017943469F, 0.00019109536F,
874 0.00020351382F, 0.00021673929F, 0.00023082423F, 0.00024582449F,
875 0.00026179955F, 0.00027881276F, 0.00029693158F, 0.00031622787F,
876 0.00033677814F, 0.00035866388F, 0.00038197188F, 0.00040679456F,
877 0.00043323036F, 0.00046138411F, 0.00049136745F, 0.00052329927F,
878 0.00055730621F, 0.00059352311F, 0.00063209358F, 0.00067317058F,
879 0.00071691700F, 0.00076350630F, 0.00081312324F, 0.00086596457F,
880 0.00092223983F, 0.00098217216F, 0.0010459992F, 0.0011139742F,
881 0.0011863665F, 0.0012634633F, 0.0013455702F, 0.0014330129F,
882 0.0015261382F, 0.0016253153F, 0.0017309374F, 0.0018434235F,
883 0.0019632195F, 0.0020908006F, 0.0022266726F, 0.0023713743F,
884 0.0025254795F, 0.0026895994F, 0.0028643847F, 0.0030505286F,
885 0.0032487691F, 0.0034598925F, 0.0036847358F, 0.0039241906F,
886 0.0041792066F, 0.0044507950F, 0.0047400328F, 0.0050480668F,
887 0.0053761186F, 0.0057254891F, 0.0060975636F, 0.0064938176F,
888 0.0069158225F, 0.0073652516F, 0.0078438871F, 0.0083536271F,
889 0.0088964928F, 0.009474637F, 0.010090352F, 0.010746080F,
890 0.011444421F, 0.012188144F, 0.012980198F, 0.013823725F,
891 0.014722068F, 0.015678791F, 0.016697687F, 0.017782797F,
892 0.018938423F, 0.020169149F, 0.021479854F, 0.022875735F,
893 0.024362330F, 0.025945531F, 0.027631618F, 0.029427276F,
894 0.031339626F, 0.033376252F, 0.035545228F, 0.037855157F,
895 0.040315199F, 0.042935108F, 0.045725273F, 0.048696758F,
896 0.051861348F, 0.055231591F, 0.058820850F, 0.062643361F,
897 0.066714279F, 0.071049749F, 0.075666962F, 0.080584227F,
898 0.085821044F, 0.091398179F, 0.097337747F, 0.10366330F,
899 0.11039993F, 0.11757434F, 0.12521498F, 0.13335215F,
900 0.14201813F, 0.15124727F, 0.16107617F, 0.17154380F,
901 0.18269168F, 0.19456402F, 0.20720788F, 0.22067342F,
902 0.23501402F, 0.25028656F, 0.26655159F, 0.28387361F,
903 0.30232132F, 0.32196786F, 0.34289114F, 0.36517414F,
904 0.38890521F, 0.41417847F, 0.44109412F, 0.46975890F,
905 0.50028648F, 0.53279791F, 0.56742212F, 0.60429640F,
906 0.64356699F, 0.68538959F, 0.72993007F, 0.77736504F,
907 0.82788260F, 0.88168307F, 0.9389798F, 1.F,
910 /* this is for per-channel noise normalization */
911 static int apsort(const void *a, const void *b){
912 float f1=**(float**)a;
913 float f2=**(float**)b;
914 return (f1<f2)-(f1>f2);
917 static void flag_lossless(int limit, float prepoint, float postpoint, float *mdct,
918 float *floor, int *flag, int i, int jn){
921 float point = j>=limit-i ? postpoint : prepoint;
922 float r = fabs(mdct[j])/floor[j];
930 /* Overload/Side effect: On input, the *q vector holds either the
931 quantized energy (for elements with the flag set) or the absolute
932 values of the *r vector (for elements with flag unset). On output,
933 *q holds the quantized energy for all elements */
934 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){
936 vorbis_info_psy *vi=p->vi;
937 float **sort = alloca(n*sizeof(*sort));
939 int start = (vi->normal_p ? vi->normal_start-i : n);
942 /* still responsible for populating *out where noise norm not in
943 effect. There's no need to [re]populate *q in these areas */
944 for(j=0;j<start;j++){
945 if(!flags || !flags[j]){ /* lossless coupling already quantized.
946 Don't touch; requantizing based on
947 energy would be incorrect. */
948 float ve = q[j]/f[j];
950 out[j] = -rint(sqrt(ve));
952 out[j] = rint(sqrt(ve));
956 /* sort magnitudes for noise norm portion of partition */
958 if(!flags || !flags[j]){ /* can't noise norm elements that have
959 already been loslessly coupled; we can
960 only account for their energy error */
961 float ve = q[j]/f[j];
962 /* Despite all the new, more capable coupling code, for now we
963 implement noise norm as it has been up to this point. Only
964 consider promotions to unit magnitude from 0. In addition
965 the only energy error counted is quantizations to zero. */
966 /* also-- the original point code only applied noise norm at > pointlimit */
967 if(ve<.25f && (!flags || j>=limit-i)){
969 sort[count++]=q+j; /* q is fabs(r) for unflagged element */
971 /* For now: no acc adjustment for nonzero quantization. populate *out and q as this value is final. */
973 out[j] = -rint(sqrt(ve));
975 out[j] = rint(sqrt(ve));
976 q[j] = out[j]*out[j]*f[j];
979 again, no energy adjustment for error in nonzero quant-- for now
984 /* noise norm to do */
985 qsort(sort,count,sizeof(*sort),apsort);
986 for(j=0;j<count;j++){
988 if(acc>=vi->normal_thresh){
989 out[k]=unitnorm(r[k]);
1002 /* Noise normalization, quantization and coupling are not wholly
1003 seperable processes in depth>1 coupling. */
1004 void _vp_couple_quantize_normalize(int blobno,
1005 vorbis_info_psy_global *g,
1007 vorbis_info_mapping0 *vi,
1011 int sliding_lowpass,
1016 int partition=(p->vi->normal_p ? p->vi->normal_partition : 16);
1017 int limit = g->coupling_pointlimit[p->vi->blockflag][blobno];
1018 float prepoint=stereo_threshholds[g->coupling_prepointamp[blobno]];
1019 float postpoint=stereo_threshholds[g->coupling_postpointamp[blobno]];
1020 float de=0.1*p->m_val; /* a blend of the AoTuV M2 and M3 code here and below */
1022 /* mdct is our raw mdct output, floor not removed. */
1023 /* inout passes in the ifloor, passes back quantized result */
1025 /* unquantized energy (negative indicates amplitude has negative sign) */
1026 float **raw = alloca(ch*sizeof(*raw));
1028 /* dual pupose; quantized energy (if flag set), othersize fabs(raw) */
1029 float **quant = alloca(ch*sizeof(*quant));
1032 float **floor = alloca(ch*sizeof(*floor));
1034 /* flags indicating raw/quantized status of elements in raw vector */
1035 int **flag = alloca(ch*sizeof(*flag));
1037 /* non-zero flag working vector */
1038 int *nz = alloca(ch*sizeof(*nz));
1040 /* energy surplus/defecit tracking */
1041 float *acc = alloca((ch+vi->coupling_steps)*sizeof(*acc));
1043 /* The threshold of a stereo is changed with the size of n */
1045 postpoint=stereo_threshholds_limited[g->coupling_postpointamp[blobno]];
1047 raw[0] = alloca(ch*partition*sizeof(**raw));
1048 quant[0] = alloca(ch*partition*sizeof(**quant));
1049 floor[0] = alloca(ch*partition*sizeof(**floor));
1050 flag[0] = alloca(ch*partition*sizeof(**flag));
1053 raw[i] = &raw[0][partition*i];
1054 quant[i] = &quant[0][partition*i];
1055 floor[i] = &floor[0][partition*i];
1056 flag[i] = &flag[0][partition*i];
1058 for(i=0;i<ch+vi->coupling_steps;i++)
1061 for(i=0;i<n;i+=partition){
1062 int k,j,jn = partition > n-i ? n-i : partition;
1065 memcpy(nz,nonzero,sizeof(*nz)*ch);
1068 memset(flag[0],0,ch*partition*sizeof(**flag));
1070 int *iout = &iwork[k][i];
1074 floor[k][j] = FLOOR1_fromdB_LOOKUP[iout[j]];
1076 flag_lossless(limit,prepoint,postpoint,&mdct[k][i],floor[k],flag[k],i,jn);
1079 quant[k][j] = raw[k][j] = mdct[k][i+j]*mdct[k][i+j];
1080 if(mdct[k][i+j]<0.f) raw[k][j]*=-1.f;
1081 floor[k][j]*=floor[k][j];
1084 acc[track]=noise_normalize(p,limit,raw[k],quant[k],floor[k],NULL,acc[track],i,jn,iout);
1088 floor[k][j] = 1e-10f;
1100 for(step=0;step<vi->coupling_steps;step++){
1101 int Mi = vi->coupling_mag[step];
1102 int Ai = vi->coupling_ang[step];
1103 int *iM = &iwork[Mi][i];
1104 int *iA = &iwork[Ai][i];
1105 float *reM = raw[Mi];
1106 float *reA = raw[Ai];
1107 float *qeM = quant[Mi];
1108 float *qeA = quant[Ai];
1109 float *floorM = floor[Mi];
1110 float *floorA = floor[Ai];
1114 if(nz[Mi] || nz[Ai]){
1115 nz[Mi] = nz[Ai] = 1;
1119 if(j<sliding_lowpass-i){
1121 /* lossless coupling */
1123 reM[j] = fabs(reM[j])+fabs(reA[j]);
1124 qeM[j] = qeM[j]+qeA[j];
1133 iA[j]=(A>0?A-B:B-A);
1135 iA[j]=(B>0?A-B:B-A);
1139 /* collapse two equivalent tuples to one */
1140 if(iA[j]>=abs(iM[j])*2){
1148 /* lossy (point) coupling */
1152 qeM[j] = fabs(reM[j]);
1156 The boost problem by the combination of noise normalization and point stereo is eased.
1157 However, this is a temporary patch.
1158 by Aoyumi @ 2004/04/18
1160 float derate = (1.0 - de*((float)(j-limit+i) / (float)(n-limit)));
1163 if(reM[j]+reA[j]<0){
1164 reM[j] = - (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate);
1166 reM[j] = (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate);
1174 floorM[j]=floorA[j]=floorM[j]+floorA[j];
1176 /* normalize the resulting mag vector */
1177 acc[track]=noise_normalize(p,limit,raw[Mi],quant[Mi],floor[Mi],flag[Mi],acc[track],i,jn,iM);
1183 for(i=0;i<vi->coupling_steps;i++){
1184 /* make sure coupling a zero and a nonzero channel results in two
1185 nonzero channels. */
1186 if(nonzero[vi->coupling_mag[i]] ||
1187 nonzero[vi->coupling_ang[i]]){
1188 nonzero[vi->coupling_mag[i]]=1;
1189 nonzero[vi->coupling_ang[i]]=1;