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-2007 *
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 double stereo_threshholds[]={0.0, .5, 1.0, 1.5, 2.5, 4.5, 8.5, 16.5, 9e10};
34 static 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 float bark=toBARK(rate/(2*n)*i);
310 for(;lo+vi->noisewindowlomin<i &&
311 toBARK(rate/(2*n)*lo)<(bark-vi->noisewindowlo);lo++);
313 for(;hi<=n && (hi<i+vi->noisewindowhimin ||
314 toBARK(rate/(2*n)*hi)<(bark+vi->noisewindowhi));hi++);
316 p->bark[i]=((lo-1)<<16)+(hi-1);
321 p->octave[i]=toOC((i+.25f)*.5*rate/n)*(1<<(p->shiftoc+1))+.5f;
323 p->tonecurves=setup_tone_curves(vi->toneatt,rate*.5/n,n,
324 vi->tone_centerboost,vi->tone_decay);
326 /* set up rolling noise median */
327 p->noiseoffset=_ogg_malloc(P_NOISECURVES*sizeof(*p->noiseoffset));
328 for(i=0;i<P_NOISECURVES;i++)
329 p->noiseoffset[i]=_ogg_malloc(n*sizeof(**p->noiseoffset));
332 float halfoc=toOC((i+.5)*rate/(2.*n))*2.;
336 if(halfoc<0)halfoc=0;
337 if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1;
338 inthalfoc=(int)halfoc;
339 del=halfoc-inthalfoc;
341 for(j=0;j<P_NOISECURVES;j++)
342 p->noiseoffset[j][i]=
343 p->vi->noiseoff[j][inthalfoc]*(1.-del) +
344 p->vi->noiseoff[j][inthalfoc+1]*del;
350 _analysis_output_always("noiseoff0",ls,p->noiseoffset[0],n,1,0,0);
351 _analysis_output_always("noiseoff1",ls,p->noiseoffset[1],n,1,0,0);
352 _analysis_output_always("noiseoff2",ls++,p->noiseoffset[2],n,1,0,0);
357 void _vp_psy_clear(vorbis_look_psy *p){
360 if(p->ath)_ogg_free(p->ath);
361 if(p->octave)_ogg_free(p->octave);
362 if(p->bark)_ogg_free(p->bark);
364 for(i=0;i<P_BANDS;i++){
365 for(j=0;j<P_LEVELS;j++){
366 _ogg_free(p->tonecurves[i][j]);
368 _ogg_free(p->tonecurves[i]);
370 _ogg_free(p->tonecurves);
373 for(i=0;i<P_NOISECURVES;i++){
374 _ogg_free(p->noiseoffset[i]);
376 _ogg_free(p->noiseoffset);
378 memset(p,0,sizeof(*p));
382 /* octave/(8*eighth_octave_lines) x scale and dB y scale */
383 static void seed_curve(float *seed,
384 const float **curves,
387 int linesper,float dBoffset){
390 const float *posts,*curve;
392 int choice=(int)((amp+dBoffset-P_LEVEL_0)*.1f);
393 choice=max(choice,0);
394 choice=min(choice,P_LEVELS-1);
395 posts=curves[choice];
398 seedptr=oc+(posts[0]-EHMER_OFFSET)*linesper-(linesper>>1);
400 for(i=posts[0];i<post1;i++){
402 float lin=amp+curve[i];
403 if(seed[seedptr]<lin)seed[seedptr]=lin;
410 static void seed_loop(vorbis_look_psy *p,
411 const float ***curves,
416 vorbis_info_psy *vi=p->vi;
418 float dBoffset=vi->max_curve_dB-specmax;
420 /* prime the working vector with peak values */
424 long oc=p->octave[i];
425 while(i+1<n && p->octave[i+1]==oc){
427 if(f[i]>max)max=f[i];
433 if(oc>=P_BANDS)oc=P_BANDS-1;
439 p->octave[i]-p->firstoc,
440 p->total_octave_lines,
441 p->eighth_octave_lines,
447 static void seed_chase(float *seeds, int linesper, long n){
448 long *posstack=alloca(n*sizeof(*posstack));
449 float *ampstack=alloca(n*sizeof(*ampstack));
457 ampstack[stack++]=seeds[i];
460 if(seeds[i]<ampstack[stack-1]){
462 ampstack[stack++]=seeds[i];
465 if(i<posstack[stack-1]+linesper){
466 if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] &&
467 i<posstack[stack-2]+linesper){
468 /* we completely overlap, making stack-1 irrelevant. pop it */
474 ampstack[stack++]=seeds[i];
482 /* the stack now contains only the positions that are relevant. Scan
483 'em straight through */
485 for(i=0;i<stack;i++){
487 if(i<stack-1 && ampstack[i+1]>ampstack[i]){
488 endpos=posstack[i+1];
490 endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is
491 discarded in short frames */
493 if(endpos>n)endpos=n;
494 for(;pos<endpos;pos++)
495 seeds[pos]=ampstack[i];
498 /* there. Linear time. I now remember this was on a problem set I
499 had in Grad Skool... I didn't solve it at the time ;-) */
503 /* bleaugh, this is more complicated than it needs to be */
505 static void max_seeds(vorbis_look_psy *p,
508 long n=p->total_octave_lines;
509 int linesper=p->eighth_octave_lines;
513 seed_chase(seed,linesper,n); /* for masking */
515 pos=p->octave[0]-p->firstoc-(linesper>>1);
517 while(linpos+1<p->n){
518 float minV=seed[pos];
519 long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc;
520 if(minV>p->vi->tone_abs_limit)minV=p->vi->tone_abs_limit;
523 if((seed[pos]>NEGINF && seed[pos]<minV) || minV==NEGINF)
528 for(;linpos<p->n && p->octave[linpos]<=end;linpos++)
529 if(flr[linpos]<minV)flr[linpos]=minV;
533 float minV=seed[p->total_octave_lines-1];
534 for(;linpos<p->n;linpos++)
535 if(flr[linpos]<minV)flr[linpos]=minV;
540 static void bark_noise_hybridmp(int n,const long *b,
546 float *N=alloca(n*sizeof(*N));
547 float *X=alloca(n*sizeof(*N));
548 float *XX=alloca(n*sizeof(*N));
549 float *Y=alloca(n*sizeof(*N));
550 float *XY=alloca(n*sizeof(*N));
552 float tN, tX, tXX, tY, tXY;
562 tN = tX = tXX = tY = tXY = 0.f;
565 if (y < 1.f) y = 1.f;
579 for (i = 1, x = 1.f; i < n; i++, x += 1.f) {
582 if (y < 1.f) y = 1.f;
599 for (i = 0, x = 0.f;; i++, x += 1.f) {
607 tXX = XX[hi] + XX[-lo];
609 tXY = XY[hi] - XY[-lo];
611 A = tY * tXX - tX * tXY;
612 B = tN * tXY - tX * tY;
613 D = tN * tXX - tX * tX;
618 noise[i] = R - offset;
621 for ( ;; i++, x += 1.f) {
629 tXX = XX[hi] - XX[lo];
631 tXY = XY[hi] - XY[lo];
633 A = tY * tXX - tX * tXY;
634 B = tN * tXY - tX * tY;
635 D = tN * tXX - tX * tX;
637 if (R < 0.f) R = 0.f;
639 noise[i] = R - offset;
641 for ( ; i < n; i++, x += 1.f) {
644 if (R < 0.f) R = 0.f;
646 noise[i] = R - offset;
649 if (fixed <= 0) return;
651 for (i = 0, x = 0.f;; i++, x += 1.f) {
658 tXX = XX[hi] + XX[-lo];
660 tXY = XY[hi] - XY[-lo];
663 A = tY * tXX - tX * tXY;
664 B = tN * tXY - tX * tY;
665 D = tN * tXX - tX * tX;
668 if (R - offset < noise[i]) noise[i] = R - offset;
670 for ( ;; i++, x += 1.f) {
678 tXX = XX[hi] - XX[lo];
680 tXY = XY[hi] - XY[lo];
682 A = tY * tXX - tX * tXY;
683 B = tN * tXY - tX * tY;
684 D = tN * tXX - tX * tX;
687 if (R - offset < noise[i]) noise[i] = R - offset;
689 for ( ; i < n; i++, x += 1.f) {
691 if (R - offset < noise[i]) noise[i] = R - offset;
695 static float FLOOR1_fromdB_INV_LOOKUP[256]={
696 0.F, 8.81683e+06F, 8.27882e+06F, 7.77365e+06F,
697 7.29930e+06F, 6.85389e+06F, 6.43567e+06F, 6.04296e+06F,
698 5.67422e+06F, 5.32798e+06F, 5.00286e+06F, 4.69759e+06F,
699 4.41094e+06F, 4.14178e+06F, 3.88905e+06F, 3.65174e+06F,
700 3.42891e+06F, 3.21968e+06F, 3.02321e+06F, 2.83873e+06F,
701 2.66551e+06F, 2.50286e+06F, 2.35014e+06F, 2.20673e+06F,
702 2.07208e+06F, 1.94564e+06F, 1.82692e+06F, 1.71544e+06F,
703 1.61076e+06F, 1.51247e+06F, 1.42018e+06F, 1.33352e+06F,
704 1.25215e+06F, 1.17574e+06F, 1.10400e+06F, 1.03663e+06F,
705 973377.F, 913981.F, 858210.F, 805842.F,
706 756669.F, 710497.F, 667142.F, 626433.F,
707 588208.F, 552316.F, 518613.F, 486967.F,
708 457252.F, 429351.F, 403152.F, 378551.F,
709 355452.F, 333762.F, 313396.F, 294273.F,
710 276316.F, 259455.F, 243623.F, 228757.F,
711 214798.F, 201691.F, 189384.F, 177828.F,
712 166977.F, 156788.F, 147221.F, 138237.F,
713 129802.F, 121881.F, 114444.F, 107461.F,
714 100903.F, 94746.3F, 88964.9F, 83536.2F,
715 78438.8F, 73652.5F, 69158.2F, 64938.1F,
716 60975.6F, 57254.9F, 53761.2F, 50480.6F,
717 47400.3F, 44507.9F, 41792.0F, 39241.9F,
718 36847.3F, 34598.9F, 32487.7F, 30505.3F,
719 28643.8F, 26896.0F, 25254.8F, 23713.7F,
720 22266.7F, 20908.0F, 19632.2F, 18434.2F,
721 17309.4F, 16253.1F, 15261.4F, 14330.1F,
722 13455.7F, 12634.6F, 11863.7F, 11139.7F,
723 10460.0F, 9821.72F, 9222.39F, 8659.64F,
724 8131.23F, 7635.06F, 7169.17F, 6731.70F,
725 6320.93F, 5935.23F, 5573.06F, 5232.99F,
726 4913.67F, 4613.84F, 4332.30F, 4067.94F,
727 3819.72F, 3586.64F, 3367.78F, 3162.28F,
728 2969.31F, 2788.13F, 2617.99F, 2458.24F,
729 2308.24F, 2167.39F, 2035.14F, 1910.95F,
730 1794.35F, 1684.85F, 1582.04F, 1485.51F,
731 1394.86F, 1309.75F, 1229.83F, 1154.78F,
732 1084.32F, 1018.15F, 956.024F, 897.687F,
733 842.910F, 791.475F, 743.179F, 697.830F,
734 655.249F, 615.265F, 577.722F, 542.469F,
735 509.367F, 478.286F, 449.101F, 421.696F,
736 395.964F, 371.803F, 349.115F, 327.812F,
737 307.809F, 289.026F, 271.390F, 254.830F,
738 239.280F, 224.679F, 210.969F, 198.096F,
739 186.008F, 174.658F, 164.000F, 153.993F,
740 144.596F, 135.773F, 127.488F, 119.708F,
741 112.404F, 105.545F, 99.1046F, 93.0572F,
742 87.3788F, 82.0469F, 77.0404F, 72.3394F,
743 67.9252F, 63.7804F, 59.8885F, 56.2341F,
744 52.8027F, 49.5807F, 46.5553F, 43.7144F,
745 41.0470F, 38.5423F, 36.1904F, 33.9821F,
746 31.9085F, 29.9614F, 28.1332F, 26.4165F,
747 24.8045F, 23.2910F, 21.8697F, 20.5352F,
748 19.2822F, 18.1056F, 17.0008F, 15.9634F,
749 14.9893F, 14.0746F, 13.2158F, 12.4094F,
750 11.6522F, 10.9411F, 10.2735F, 9.64662F,
751 9.05798F, 8.50526F, 7.98626F, 7.49894F,
752 7.04135F, 6.61169F, 6.20824F, 5.82941F,
753 5.47370F, 5.13970F, 4.82607F, 4.53158F,
754 4.25507F, 3.99542F, 3.75162F, 3.52269F,
755 3.30774F, 3.10590F, 2.91638F, 2.73842F,
756 2.57132F, 2.41442F, 2.26709F, 2.12875F,
757 1.99885F, 1.87688F, 1.76236F, 1.65482F,
758 1.55384F, 1.45902F, 1.36999F, 1.28640F,
759 1.20790F, 1.13419F, 1.06499F, 1.F
762 void _vp_remove_floor(vorbis_look_psy *p,
766 int sliding_lowpass){
770 if(sliding_lowpass>n)sliding_lowpass=n;
772 for(i=0;i<sliding_lowpass;i++){
774 mdct[i]*FLOOR1_fromdB_INV_LOOKUP[codedflr[i]];
781 void _vp_noisemask(vorbis_look_psy *p,
786 float *work=alloca(n*sizeof(*work));
788 bark_noise_hybridmp(n,p->bark,logmdct,logmask,
791 for(i=0;i<n;i++)work[i]=logmdct[i]-logmask[i];
793 bark_noise_hybridmp(n,p->bark,work,logmask,0.,
794 p->vi->noisewindowfixed);
796 for(i=0;i<n;i++)work[i]=logmdct[i]-work[i];
804 work2[i]=logmask[i]+work[i];
808 _analysis_output("median2R",seq/2,work,n,1,0,0);
810 _analysis_output("median2L",seq/2,work,n,1,0,0);
813 _analysis_output("envelope2R",seq/2,work2,n,1,0,0);
815 _analysis_output("envelope2L",seq/2,work2,n,1,0,0);
821 int dB=logmask[i]+.5;
822 if(dB>=NOISE_COMPAND_LEVELS)dB=NOISE_COMPAND_LEVELS-1;
824 logmask[i]= work[i]+p->vi->noisecompand[dB];
829 void _vp_tonemask(vorbis_look_psy *p,
832 float global_specmax,
833 float local_specmax){
837 float *seed=alloca(sizeof(*seed)*p->total_octave_lines);
838 float att=local_specmax+p->vi->ath_adjatt;
839 for(i=0;i<p->total_octave_lines;i++)seed[i]=NEGINF;
841 /* set the ATH (floating below localmax, not global max by a
843 if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt;
846 logmask[i]=p->ath[i]+att;
849 seed_loop(p,(const float ***)p->tonecurves,logfft,logmask,seed,global_specmax);
850 max_seeds(p,seed,logmask);
854 void _vp_offset_and_mix(vorbis_look_psy *p,
862 float de, coeffi, cx;/* AoTuV */
863 float toneatt=p->vi->tone_masteratt[offset_select];
868 float val= noise[i]+p->noiseoffset[offset_select][i];
869 if(val>p->vi->noisemaxsupp)val=p->vi->noisemaxsupp;
870 logmask[i]=max(val,tone[i]+toneatt);
875 The following codes improve a noise problem.
876 A fundamental idea uses the value of masking and carries out
877 the relative compensation of the MDCT.
878 However, this code is not perfect and all noise problems cannot be solved.
879 by Aoyumi @ 2004/04/18
882 if(offset_select == 1) {
883 coeffi = -17.2; /* coeffi is a -17.2dB threshold */
884 val = val - logmdct[i]; /* val == mdct line value relative to floor in dB */
887 /* mdct value is > -17.2 dB below floor */
889 de = 1.0-((val-coeffi)*0.005*cx);
890 /* pro-rated attenuation:
891 -0.00 dB boost if mdct value is -17.2dB (relative to floor)
892 -0.77 dB boost if mdct value is 0dB (relative to floor)
893 -1.64 dB boost if mdct value is +17.2dB (relative to floor)
896 if(de < 0) de = 0.0001;
898 /* mdct value is <= -17.2 dB below floor */
900 de = 1.0-((val-coeffi)*0.0003*cx);
901 /* pro-rated attenuation:
902 +0.00 dB atten if mdct value is -17.2dB (relative to floor)
903 +0.45 dB atten if mdct value is -34.4dB (relative to floor)
912 float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){
913 vorbis_info *vi=vd->vi;
914 codec_setup_info *ci=vi->codec_setup;
915 vorbis_info_psy_global *gi=&ci->psy_g_param;
917 int n=ci->blocksizes[vd->W]/2;
918 float secs=(float)n/vi->rate;
920 amp+=secs*gi->ampmax_att_per_sec;
921 if(amp<-9999)amp=-9999;
925 static void couple_lossless(float A, float B,
926 float *qA, float *qB){
927 int test1=fabs(*qA)>fabs(*qB);
928 test1-= fabs(*qA)<fabs(*qB);
930 if(!test1)test1=((fabs(A)>fabs(B))<<1)-1;
932 *qB=(*qA>0.f?*qA-*qB:*qB-*qA);
935 *qB=(*qB>0.f?*qA-*qB:*qB-*qA);
939 if(*qB>fabs(*qA)*1.9999f){
945 static float hypot_lookup[32]={
946 -0.009935, -0.011245, -0.012726, -0.014397,
947 -0.016282, -0.018407, -0.020800, -0.023494,
948 -0.026522, -0.029923, -0.033737, -0.038010,
949 -0.042787, -0.048121, -0.054064, -0.060671,
950 -0.068000, -0.076109, -0.085054, -0.094892,
951 -0.105675, -0.117451, -0.130260, -0.144134,
952 -0.159093, -0.175146, -0.192286, -0.210490,
953 -0.229718, -0.249913, -0.271001, -0.292893};
955 static void precomputed_couple_point(float premag,
956 int floorA,int floorB,
957 float *mag, float *ang){
959 int test=(floorA>floorB)-1;
960 int offset=31-abs(floorA-floorB);
961 float floormag=hypot_lookup[((offset<0)-1)&offset]+1.f;
963 floormag*=FLOOR1_fromdB_INV_LOOKUP[(floorB&test)|(floorA&(~test))];
965 *mag=premag*floormag;
969 /* just like below, this is currently set up to only do
970 single-step-depth coupling. Otherwise, we'd have to do more
971 copying (which will be inevitable later) */
973 /* doing the real circular magnitude calculation is audibly superior
975 static float dipole_hypot(float a, float b){
977 if(b>0.)return sqrt(a*a+b*b);
978 if(a>-b)return sqrt(a*a-b*b);
979 return -sqrt(b*b-a*a);
981 if(b<0.)return -sqrt(a*a+b*b);
982 if(-a>b)return -sqrt(a*a-b*b);
983 return sqrt(b*b-a*a);
985 static float round_hypot(float a, float b){
987 if(b>0.)return sqrt(a*a+b*b);
988 if(a>-b)return sqrt(a*a+b*b);
989 return -sqrt(b*b+a*a);
991 if(b<0.)return -sqrt(a*a+b*b);
992 if(-a>b)return -sqrt(a*a+b*b);
993 return sqrt(b*b+a*a);
996 /* revert to round hypot for now */
997 float **_vp_quantize_couple_memo(vorbis_block *vb,
998 vorbis_info_psy_global *g,
1000 vorbis_info_mapping0 *vi,
1004 float **ret=_vorbis_block_alloc(vb,vi->coupling_steps*sizeof(*ret));
1005 int limit=g->coupling_pointlimit[p->vi->blockflag][PACKETBLOBS/2];
1007 for(i=0;i<vi->coupling_steps;i++){
1008 float *mdctM=mdct[vi->coupling_mag[i]];
1009 float *mdctA=mdct[vi->coupling_ang[i]];
1010 ret[i]=_vorbis_block_alloc(vb,n*sizeof(**ret));
1011 for(j=0;j<limit;j++)
1012 ret[i][j]=dipole_hypot(mdctM[j],mdctA[j]);
1014 ret[i][j]=round_hypot(mdctM[j],mdctA[j]);
1020 /* this is for per-channel noise normalization */
1021 static int apsort(const void *a, const void *b){
1022 float f1=fabs(**(float**)a);
1023 float f2=fabs(**(float**)b);
1024 return (f1<f2)-(f1>f2);
1027 int **_vp_quantize_couple_sort(vorbis_block *vb,
1029 vorbis_info_mapping0 *vi,
1033 if(p->vi->normal_point_p){
1035 int **ret=_vorbis_block_alloc(vb,vi->coupling_steps*sizeof(*ret));
1036 int partition=p->vi->normal_partition;
1037 float **work=alloca(sizeof(*work)*partition);
1039 for(i=0;i<vi->coupling_steps;i++){
1040 ret[i]=_vorbis_block_alloc(vb,n*sizeof(**ret));
1042 for(j=0;j<n;j+=partition){
1043 for(k=0;k<partition;k++)work[k]=mags[i]+k+j;
1044 qsort(work,partition,sizeof(*work),apsort);
1045 for(k=0;k<partition;k++)ret[i][k+j]=work[k]-mags[i];
1053 void _vp_noise_normalize_sort(vorbis_look_psy *p,
1054 float *magnitudes,int *sortedindex){
1056 vorbis_info_psy *vi=p->vi;
1057 int partition=vi->normal_partition;
1058 float **work=alloca(sizeof(*work)*partition);
1059 int start=vi->normal_start;
1061 for(j=start;j<n;j+=partition){
1062 if(j+partition>n)partition=n-j;
1063 for(i=0;i<partition;i++)work[i]=magnitudes+i+j;
1064 qsort(work,partition,sizeof(*work),apsort);
1065 for(i=0;i<partition;i++){
1066 sortedindex[i+j-start]=work[i]-magnitudes;
1071 void _vp_noise_normalize(vorbis_look_psy *p,
1072 float *in,float *out,int *sortedindex){
1073 int flag=0,i,j=0,n=p->n;
1074 vorbis_info_psy *vi=p->vi;
1075 int partition=vi->normal_partition;
1076 int start=vi->normal_start;
1080 if(vi->normal_channel_p){
1084 for(;j+partition<=n;j+=partition){
1088 for(i=j;i<j+partition;i++)
1091 for(i=0;i<partition;i++){
1092 k=sortedindex[i+j-start];
1094 if(in[k]*in[k]>=.25f){
1099 if(acc<vi->normal_thresh)break;
1100 out[k]=unitnorm(in[k]);
1105 for(;i<partition;i++){
1106 k=sortedindex[i+j-start];
1117 void _vp_couple(int blobno,
1118 vorbis_info_psy_global *g,
1120 vorbis_info_mapping0 *vi,
1126 int sliding_lowpass){
1130 /* perform any requested channel coupling */
1131 /* point stereo can only be used in a first stage (in this encoder)
1132 because of the dependency on floor lookups */
1133 for(i=0;i<vi->coupling_steps;i++){
1135 /* once we're doing multistage coupling in which a channel goes
1136 through more than one coupling step, the floor vector
1137 magnitudes will also have to be recalculated an propogated
1138 along with PCM. Right now, we're not (that will wait until 5.1
1139 most likely), so the code isn't here yet. The memory management
1140 here is all assuming single depth couplings anyway. */
1142 /* make sure coupling a zero and a nonzero channel results in two
1143 nonzero channels. */
1144 if(nonzero[vi->coupling_mag[i]] ||
1145 nonzero[vi->coupling_ang[i]]){
1148 float *rM=res[vi->coupling_mag[i]];
1149 float *rA=res[vi->coupling_ang[i]];
1152 int *floorM=ifloor[vi->coupling_mag[i]];
1153 int *floorA=ifloor[vi->coupling_ang[i]];
1154 float prepoint=stereo_threshholds[g->coupling_prepointamp[blobno]];
1155 float postpoint=stereo_threshholds[g->coupling_postpointamp[blobno]];
1156 int partition=(p->vi->normal_point_p?p->vi->normal_partition:p->n);
1157 int limit=g->coupling_pointlimit[p->vi->blockflag][blobno];
1158 int pointlimit=limit;
1160 nonzero[vi->coupling_mag[i]]=1;
1161 nonzero[vi->coupling_ang[i]]=1;
1163 /* The threshold of a stereo is changed with the size of n */
1165 postpoint=stereo_threshholds_limited[g->coupling_postpointamp[blobno]];
1167 for(j=0;j<p->n;j+=partition){
1170 for(k=0;k<partition;k++){
1173 if(l<sliding_lowpass){
1174 if((l>=limit && fabs(rM[l])<postpoint && fabs(rA[l])<postpoint) ||
1175 (fabs(rM[l])<prepoint && fabs(rA[l])<prepoint)){
1178 precomputed_couple_point(mag_memo[i][l],
1179 floorM[l],floorA[l],
1182 if(rint(qM[l])==0.f)acc+=qM[l]*qM[l];
1184 couple_lossless(rM[l],rA[l],qM+l,qA+l);
1192 if(p->vi->normal_point_p){
1193 for(k=0;k<partition && acc>=p->vi->normal_thresh;k++){
1194 int l=mag_sort[i][j+k];
1195 if(l<sliding_lowpass && l>=pointlimit && rint(qM[l])==0.f){
1196 qM[l]=unitnorm(qM[l]);
1208 The boost problem by the combination of noise normalization and point stereo is eased.
1209 However, this is a temporary patch.
1210 by Aoyumi @ 2004/04/18
1213 void hf_reduction(vorbis_info_psy_global *g,
1215 vorbis_info_mapping0 *vi,
1218 int i,j,n=p->n, de=0.3*p->m_val;
1219 int limit=g->coupling_pointlimit[p->vi->blockflag][PACKETBLOBS/2];
1221 for(i=0; i<vi->coupling_steps; i++){
1222 /* for(j=start; j<limit; j++){} // ???*/
1223 for(j=limit; j<n; j++)
1224 mdct[i][j] *= (1.0 - de*((float)(j-limit) / (float)(n-limit)));