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-2002 *
9 * by the XIPHOPHORUS Company http://www.xiph.org/ *
11 ********************************************************************
13 function: psychoacoustics not including preecho
14 last mod: $Id: psy.c,v 1.81 2002/10/21 07:00:11 xiphmont Exp $
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;
559 tN = tX = tXX = tY = tXY = 0.f;
562 if (y < 1.f) y = 1.f;
576 for (i = 1, x = 1.f; i < n; i++, x += 1.f) {
579 if (y < 1.f) y = 1.f;
596 for (i = 0, x = 0.f;; i++, x += 1.f) {
604 tXX = XX[hi] + XX[-lo];
606 tXY = XY[hi] - XY[-lo];
608 A = tY * tXX - tX * tXY;
609 B = tN * tXY - tX * tY;
610 D = tN * tXX - tX * tX;
615 noise[i] = R - offset;
618 for ( ;; i++, x += 1.f) {
626 tXX = XX[hi] - XX[lo];
628 tXY = XY[hi] - XY[lo];
630 A = tY * tXX - tX * tXY;
631 B = tN * tXY - tX * tY;
632 D = tN * tXX - tX * tX;
634 if (R < 0.f) R = 0.f;
636 noise[i] = R - offset;
638 for ( ; i < n; i++, x += 1.f) {
641 if (R < 0.f) R = 0.f;
643 noise[i] = R - offset;
646 if (fixed <= 0) return;
648 for (i = 0, x = 0.f;; i++, x += 1.f) {
655 tXX = XX[hi] + XX[-lo];
657 tXY = XY[hi] - XY[-lo];
660 A = tY * tXX - tX * tXY;
661 B = tN * tXY - tX * tY;
662 D = tN * tXX - tX * tX;
665 if (R - offset < noise[i]) noise[i] = R - offset;
667 for ( ;; i++, x += 1.f) {
675 tXX = XX[hi] - XX[lo];
677 tXY = XY[hi] - XY[lo];
679 A = tY * tXX - tX * tXY;
680 B = tN * tXY - tX * tY;
681 D = tN * tXX - tX * tX;
684 if (R - offset < noise[i]) noise[i] = R - offset;
686 for ( ; i < n; i++, x += 1.f) {
688 if (R - offset < noise[i]) noise[i] = R - offset;
692 static float FLOOR1_fromdB_INV_LOOKUP[256]={
693 0.F, 8.81683e+06F, 8.27882e+06F, 7.77365e+06F,
694 7.29930e+06F, 6.85389e+06F, 6.43567e+06F, 6.04296e+06F,
695 5.67422e+06F, 5.32798e+06F, 5.00286e+06F, 4.69759e+06F,
696 4.41094e+06F, 4.14178e+06F, 3.88905e+06F, 3.65174e+06F,
697 3.42891e+06F, 3.21968e+06F, 3.02321e+06F, 2.83873e+06F,
698 2.66551e+06F, 2.50286e+06F, 2.35014e+06F, 2.20673e+06F,
699 2.07208e+06F, 1.94564e+06F, 1.82692e+06F, 1.71544e+06F,
700 1.61076e+06F, 1.51247e+06F, 1.42018e+06F, 1.33352e+06F,
701 1.25215e+06F, 1.17574e+06F, 1.10400e+06F, 1.03663e+06F,
702 973377.F, 913981.F, 858210.F, 805842.F,
703 756669.F, 710497.F, 667142.F, 626433.F,
704 588208.F, 552316.F, 518613.F, 486967.F,
705 457252.F, 429351.F, 403152.F, 378551.F,
706 355452.F, 333762.F, 313396.F, 294273.F,
707 276316.F, 259455.F, 243623.F, 228757.F,
708 214798.F, 201691.F, 189384.F, 177828.F,
709 166977.F, 156788.F, 147221.F, 138237.F,
710 129802.F, 121881.F, 114444.F, 107461.F,
711 100903.F, 94746.3F, 88964.9F, 83536.2F,
712 78438.8F, 73652.5F, 69158.2F, 64938.1F,
713 60975.6F, 57254.9F, 53761.2F, 50480.6F,
714 47400.3F, 44507.9F, 41792.0F, 39241.9F,
715 36847.3F, 34598.9F, 32487.7F, 30505.3F,
716 28643.8F, 26896.0F, 25254.8F, 23713.7F,
717 22266.7F, 20908.0F, 19632.2F, 18434.2F,
718 17309.4F, 16253.1F, 15261.4F, 14330.1F,
719 13455.7F, 12634.6F, 11863.7F, 11139.7F,
720 10460.0F, 9821.72F, 9222.39F, 8659.64F,
721 8131.23F, 7635.06F, 7169.17F, 6731.70F,
722 6320.93F, 5935.23F, 5573.06F, 5232.99F,
723 4913.67F, 4613.84F, 4332.30F, 4067.94F,
724 3819.72F, 3586.64F, 3367.78F, 3162.28F,
725 2969.31F, 2788.13F, 2617.99F, 2458.24F,
726 2308.24F, 2167.39F, 2035.14F, 1910.95F,
727 1794.35F, 1684.85F, 1582.04F, 1485.51F,
728 1394.86F, 1309.75F, 1229.83F, 1154.78F,
729 1084.32F, 1018.15F, 956.024F, 897.687F,
730 842.910F, 791.475F, 743.179F, 697.830F,
731 655.249F, 615.265F, 577.722F, 542.469F,
732 509.367F, 478.286F, 449.101F, 421.696F,
733 395.964F, 371.803F, 349.115F, 327.812F,
734 307.809F, 289.026F, 271.390F, 254.830F,
735 239.280F, 224.679F, 210.969F, 198.096F,
736 186.008F, 174.658F, 164.000F, 153.993F,
737 144.596F, 135.773F, 127.488F, 119.708F,
738 112.404F, 105.545F, 99.1046F, 93.0572F,
739 87.3788F, 82.0469F, 77.0404F, 72.3394F,
740 67.9252F, 63.7804F, 59.8885F, 56.2341F,
741 52.8027F, 49.5807F, 46.5553F, 43.7144F,
742 41.0470F, 38.5423F, 36.1904F, 33.9821F,
743 31.9085F, 29.9614F, 28.1332F, 26.4165F,
744 24.8045F, 23.2910F, 21.8697F, 20.5352F,
745 19.2822F, 18.1056F, 17.0008F, 15.9634F,
746 14.9893F, 14.0746F, 13.2158F, 12.4094F,
747 11.6522F, 10.9411F, 10.2735F, 9.64662F,
748 9.05798F, 8.50526F, 7.98626F, 7.49894F,
749 7.04135F, 6.61169F, 6.20824F, 5.82941F,
750 5.47370F, 5.13970F, 4.82607F, 4.53158F,
751 4.25507F, 3.99542F, 3.75162F, 3.52269F,
752 3.30774F, 3.10590F, 2.91638F, 2.73842F,
753 2.57132F, 2.41442F, 2.26709F, 2.12875F,
754 1.99885F, 1.87688F, 1.76236F, 1.65482F,
755 1.55384F, 1.45902F, 1.36999F, 1.28640F,
756 1.20790F, 1.13419F, 1.06499F, 1.F
759 void _vp_remove_floor(vorbis_look_psy *p,
763 int sliding_lowpass){
767 if(sliding_lowpass>n)sliding_lowpass=n;
769 for(i=0;i<sliding_lowpass;i++){
771 mdct[i]*FLOOR1_fromdB_INV_LOOKUP[codedflr[i]];
778 void _vp_noisemask(vorbis_look_psy *p,
783 float *work=alloca(n*sizeof(*work));
785 bark_noise_hybridmp(n,p->bark,logmdct,logmask,
788 for(i=0;i<n;i++)work[i]=logmdct[i]-logmask[i];
790 bark_noise_hybridmp(n,p->bark,work,logmask,0.,
791 p->vi->noisewindowfixed);
793 for(i=0;i<n;i++)work[i]=logmdct[i]-work[i];
801 work2[i]=logmask[i]+work[i];
805 _analysis_output("median2R",seq/2,work,n,1,0,0);
807 _analysis_output("median2L",seq/2,work,n,1,0,0);
810 _analysis_output("envelope2R",seq/2,work2,n,1,0,0);
812 _analysis_output("envelope2L",seq/2,work2,n,1,0,0);
818 int dB=logmask[i]+.5;
819 if(dB>=NOISE_COMPAND_LEVELS)dB=NOISE_COMPAND_LEVELS-1;
821 logmask[i]= work[i]+p->vi->noisecompand[dB];
826 void _vp_tonemask(vorbis_look_psy *p,
829 float global_specmax,
830 float local_specmax){
834 float *seed=alloca(sizeof(*seed)*p->total_octave_lines);
835 float att=local_specmax+p->vi->ath_adjatt;
836 for(i=0;i<p->total_octave_lines;i++)seed[i]=NEGINF;
838 /* set the ATH (floating below localmax, not global max by a
840 if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt;
843 logmask[i]=p->ath[i]+att;
846 seed_loop(p,(const float ***)p->tonecurves,logfft,logmask,seed,global_specmax);
847 max_seeds(p,seed,logmask);
851 void _vp_offset_and_mix(vorbis_look_psy *p,
859 float de, coeffi, cx;/* AoTuV */
860 float toneatt=p->vi->tone_masteratt[offset_select];
865 float val= noise[i]+p->noiseoffset[offset_select][i];
866 if(val>p->vi->noisemaxsupp)val=p->vi->noisemaxsupp;
867 logmask[i]=max(val,tone[i]+toneatt);
872 The following codes improve a noise problem.
873 A fundamental idea uses the value of masking and carries out
874 the relative compensation of the MDCT.
875 However, this code is not perfect and all noise problems cannot be solved.
876 by Aoyumi @ 2004/04/18
879 if(offset_select == 1) {
880 coeffi = -17.2; /* coeffi is a -17.2dB threshold */
881 val = val - logmdct[i]; /* val == mdct line value relative to floor in dB */
884 /* mdct value is > -17.2 dB below floor */
886 de = 1.0-((val-coeffi)*0.005*cx);
887 /* pro-rated attenuation:
888 -0.00 dB boost if mdct value is -17.2dB (relative to floor)
889 -0.77 dB boost if mdct value is 0dB (relative to floor)
890 -1.64 dB boost if mdct value is +17.2dB (relative to floor)
893 if(de < 0) de = 0.0001;
895 /* mdct value is <= -17.2 dB below floor */
897 de = 1.0-((val-coeffi)*0.0003*cx);
898 /* pro-rated attenuation:
899 +0.00 dB atten if mdct value is -17.2dB (relative to floor)
900 +0.45 dB atten if mdct value is -34.4dB (relative to floor)
909 float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){
910 vorbis_info *vi=vd->vi;
911 codec_setup_info *ci=vi->codec_setup;
912 vorbis_info_psy_global *gi=&ci->psy_g_param;
914 int n=ci->blocksizes[vd->W]/2;
915 float secs=(float)n/vi->rate;
917 amp+=secs*gi->ampmax_att_per_sec;
918 if(amp<-9999)amp=-9999;
922 static void couple_lossless(float A, float B,
923 float *qA, float *qB){
924 int test1=fabs(*qA)>fabs(*qB);
925 test1-= fabs(*qA)<fabs(*qB);
927 if(!test1)test1=((fabs(A)>fabs(B))<<1)-1;
929 *qB=(*qA>0.f?*qA-*qB:*qB-*qA);
932 *qB=(*qB>0.f?*qA-*qB:*qB-*qA);
936 if(*qB>fabs(*qA)*1.9999f){
942 static float hypot_lookup[32]={
943 -0.009935, -0.011245, -0.012726, -0.014397,
944 -0.016282, -0.018407, -0.020800, -0.023494,
945 -0.026522, -0.029923, -0.033737, -0.038010,
946 -0.042787, -0.048121, -0.054064, -0.060671,
947 -0.068000, -0.076109, -0.085054, -0.094892,
948 -0.105675, -0.117451, -0.130260, -0.144134,
949 -0.159093, -0.175146, -0.192286, -0.210490,
950 -0.229718, -0.249913, -0.271001, -0.292893};
952 static void precomputed_couple_point(float premag,
953 int floorA,int floorB,
954 float *mag, float *ang){
956 int test=(floorA>floorB)-1;
957 int offset=31-abs(floorA-floorB);
958 float floormag=hypot_lookup[((offset<0)-1)&offset]+1.f;
960 floormag*=FLOOR1_fromdB_INV_LOOKUP[(floorB&test)|(floorA&(~test))];
962 *mag=premag*floormag;
966 /* just like below, this is currently set up to only do
967 single-step-depth coupling. Otherwise, we'd have to do more
968 copying (which will be inevitable later) */
970 /* doing the real circular magnitude calculation is audibly superior
972 static float dipole_hypot(float a, float b){
974 if(b>0.)return sqrt(a*a+b*b);
975 if(a>-b)return sqrt(a*a-b*b);
976 return -sqrt(b*b-a*a);
978 if(b<0.)return -sqrt(a*a+b*b);
979 if(-a>b)return -sqrt(a*a-b*b);
980 return sqrt(b*b-a*a);
982 static float round_hypot(float a, float b){
984 if(b>0.)return sqrt(a*a+b*b);
985 if(a>-b)return sqrt(a*a+b*b);
986 return -sqrt(b*b+a*a);
988 if(b<0.)return -sqrt(a*a+b*b);
989 if(-a>b)return -sqrt(a*a+b*b);
990 return sqrt(b*b+a*a);
993 /* revert to round hypot for now */
994 float **_vp_quantize_couple_memo(vorbis_block *vb,
995 vorbis_info_psy_global *g,
997 vorbis_info_mapping0 *vi,
1001 float **ret=_vorbis_block_alloc(vb,vi->coupling_steps*sizeof(*ret));
1002 int limit=g->coupling_pointlimit[p->vi->blockflag][PACKETBLOBS/2];
1004 for(i=0;i<vi->coupling_steps;i++){
1005 float *mdctM=mdct[vi->coupling_mag[i]];
1006 float *mdctA=mdct[vi->coupling_ang[i]];
1007 ret[i]=_vorbis_block_alloc(vb,n*sizeof(**ret));
1008 for(j=0;j<limit;j++)
1009 ret[i][j]=dipole_hypot(mdctM[j],mdctA[j]);
1011 ret[i][j]=round_hypot(mdctM[j],mdctA[j]);
1017 /* this is for per-channel noise normalization */
1018 static int apsort(const void *a, const void *b){
1019 float f1=fabs(**(float**)a);
1020 float f2=fabs(**(float**)b);
1021 return (f1<f2)-(f1>f2);
1024 int **_vp_quantize_couple_sort(vorbis_block *vb,
1026 vorbis_info_mapping0 *vi,
1030 if(p->vi->normal_point_p){
1032 int **ret=_vorbis_block_alloc(vb,vi->coupling_steps*sizeof(*ret));
1033 int partition=p->vi->normal_partition;
1034 float **work=alloca(sizeof(*work)*partition);
1036 for(i=0;i<vi->coupling_steps;i++){
1037 ret[i]=_vorbis_block_alloc(vb,n*sizeof(**ret));
1039 for(j=0;j<n;j+=partition){
1040 for(k=0;k<partition;k++)work[k]=mags[i]+k+j;
1041 qsort(work,partition,sizeof(*work),apsort);
1042 for(k=0;k<partition;k++)ret[i][k+j]=work[k]-mags[i];
1050 void _vp_noise_normalize_sort(vorbis_look_psy *p,
1051 float *magnitudes,int *sortedindex){
1053 vorbis_info_psy *vi=p->vi;
1054 int partition=vi->normal_partition;
1055 float **work=alloca(sizeof(*work)*partition);
1056 int start=vi->normal_start;
1058 for(j=start;j<n;j+=partition){
1059 if(j+partition>n)partition=n-j;
1060 for(i=0;i<partition;i++)work[i]=magnitudes+i+j;
1061 qsort(work,partition,sizeof(*work),apsort);
1062 for(i=0;i<partition;i++){
1063 sortedindex[i+j-start]=work[i]-magnitudes;
1068 void _vp_noise_normalize(vorbis_look_psy *p,
1069 float *in,float *out,int *sortedindex){
1070 int flag=0,i,j=0,n=p->n;
1071 vorbis_info_psy *vi=p->vi;
1072 int partition=vi->normal_partition;
1073 int start=vi->normal_start;
1077 if(vi->normal_channel_p){
1081 for(;j+partition<=n;j+=partition){
1085 for(i=j;i<j+partition;i++)
1088 for(i=0;i<partition;i++){
1089 k=sortedindex[i+j-start];
1091 if(in[k]*in[k]>=.25f){
1096 if(acc<vi->normal_thresh)break;
1097 out[k]=unitnorm(in[k]);
1102 for(;i<partition;i++){
1103 k=sortedindex[i+j-start];
1114 void _vp_couple(int blobno,
1115 vorbis_info_psy_global *g,
1117 vorbis_info_mapping0 *vi,
1123 int sliding_lowpass){
1127 /* perform any requested channel coupling */
1128 /* point stereo can only be used in a first stage (in this encoder)
1129 because of the dependency on floor lookups */
1130 for(i=0;i<vi->coupling_steps;i++){
1132 /* once we're doing multistage coupling in which a channel goes
1133 through more than one coupling step, the floor vector
1134 magnitudes will also have to be recalculated an propogated
1135 along with PCM. Right now, we're not (that will wait until 5.1
1136 most likely), so the code isn't here yet. The memory management
1137 here is all assuming single depth couplings anyway. */
1139 /* make sure coupling a zero and a nonzero channel results in two
1140 nonzero channels. */
1141 if(nonzero[vi->coupling_mag[i]] ||
1142 nonzero[vi->coupling_ang[i]]){
1145 float *rM=res[vi->coupling_mag[i]];
1146 float *rA=res[vi->coupling_ang[i]];
1149 int *floorM=ifloor[vi->coupling_mag[i]];
1150 int *floorA=ifloor[vi->coupling_ang[i]];
1151 float prepoint=stereo_threshholds[g->coupling_prepointamp[blobno]];
1152 float postpoint=stereo_threshholds[g->coupling_postpointamp[blobno]];
1153 int partition=(p->vi->normal_point_p?p->vi->normal_partition:p->n);
1154 int limit=g->coupling_pointlimit[p->vi->blockflag][blobno];
1155 int pointlimit=limit;
1157 nonzero[vi->coupling_mag[i]]=1;
1158 nonzero[vi->coupling_ang[i]]=1;
1160 /* The threshold of a stereo is changed with the size of n */
1162 postpoint=stereo_threshholds_limited[g->coupling_postpointamp[blobno]];
1164 for(j=0;j<p->n;j+=partition){
1167 for(k=0;k<partition;k++){
1170 if(l<sliding_lowpass){
1171 if((l>=limit && fabs(rM[l])<postpoint && fabs(rA[l])<postpoint) ||
1172 (fabs(rM[l])<prepoint && fabs(rA[l])<prepoint)){
1175 precomputed_couple_point(mag_memo[i][l],
1176 floorM[l],floorA[l],
1179 if(rint(qM[l])==0.f)acc+=qM[l]*qM[l];
1181 couple_lossless(rM[l],rA[l],qM+l,qA+l);
1189 if(p->vi->normal_point_p){
1190 for(k=0;k<partition && acc>=p->vi->normal_thresh;k++){
1191 int l=mag_sort[i][j+k];
1192 if(l<sliding_lowpass && l>=pointlimit && rint(qM[l])==0.f){
1193 qM[l]=unitnorm(qM[l]);
1205 The boost problem by the combination of noise normalization and point stereo is eased.
1206 However, this is a temporary patch.
1207 by Aoyumi @ 2004/04/18
1210 void hf_reduction(vorbis_info_psy_global *g,
1212 vorbis_info_mapping0 *vi,
1215 int i,j,n=p->n, de=0.3*p->m_val;
1216 int limit=g->coupling_pointlimit[p->vi->blockflag][PACKETBLOBS/2];
1217 int start=p->vi->normal_start;
1219 for(i=0; i<vi->coupling_steps; i++){
1220 /* for(j=start; j<limit; j++){} // ???*/
1221 for(j=limit; j<n; j++)
1222 mdct[i][j] *= (1.0 - de*((float)(j-limit) / (float)(n-limit)));