lib/psy.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
   5  * THE GNU LESSER/LIBRARY PUBLIC LICENSE, WHICH IS INCLUDED WITH    *
   6  * THIS SOURCE. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.        *
   7  *                                                                  *
   8  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2000             *
   9  * by Monty <monty@xiph.org> and the XIPHOPHORUS Company            *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: psychoacoustics not including preecho
  15  last mod: $Id: psy.c,v 1.37 2001/01/30 23:40:33 xiphmont Exp $
  16
  17  ********************************************************************/
  18
  19 #include <stdlib.h>
  20 #include <math.h>
  21 #include <string.h>
  22 #include "vorbis/codec.h"
  23 #include "codec_internal.h"
  24
  25 #include "masking.h"
  26 #include "psy.h"
  27 #include "os.h"
  28 #include "lpc.h"
  29 #include "smallft.h"
  30 #include "scales.h"
  31 #include "misc.h"
  32
  33 #define NEGINF -9999.f
  34
  35 /* Why Bark scale for encoding but not masking computation? Because
  36    masking has a strong harmonic dependancy */
  37
  38 /* the beginnings of real psychoacoustic infrastructure.  This is
  39    still not tightly tuned */
  40 void _vi_psy_free(vorbis_info_psy *i){
  41   if(i){
  42     memset(i,0,sizeof(vorbis_info_psy));
  43     _ogg_free(i);
  44   }
  45 }
  46
  47 vorbis_info_psy *_vi_psy_copy(vorbis_info_psy *i){
  48   vorbis_info_psy *ret=_ogg_malloc(sizeof(vorbis_info_psy));
  49   memcpy(ret,i,sizeof(vorbis_info_psy));
  50   return(ret);
  51 }
  52
  53 /* Set up decibel threshold slopes on a Bark frequency scale */
  54 /* ATH is the only bit left on a Bark scale.  No reason to change it
  55    right now */
  56 static void set_curve(float *ref,float *c,int n, float crate){
  57   int i,j=0;
  58
  59   for(i=0;i<MAX_BARK-1;i++){
  60     int endpos=rint(fromBARK(i+1)*2*n/crate);
  61     float base=ref[i];
  62     if(j<endpos){
  63       float delta=(ref[i+1]-base)/(endpos-j);
  64       for(;j<endpos && j<n;j++){
  65         c[j]=base;
  66         base+=delta;
  67       }
  68     }
  69   }
  70 }
  71
  72 static void min_curve(float *c,
  73                        float *c2){
  74   int i;
  75   for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i];
  76 }
  77 static void max_curve(float *c,
  78                        float *c2){
  79   int i;
  80   for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i];
  81 }
  82
  83 static void attenuate_curve(float *c,float att){
  84   int i;
  85   for(i=0;i<EHMER_MAX;i++)
  86     c[i]+=att;
  87 }
  88
  89 static void interp_curve(float *c,float *c1,float *c2,float del){
  90   int i;
  91   for(i=0;i<EHMER_MAX;i++)
  92     c[i]=c2[i]*del+c1[i]*(1.f-del);
  93 }
  94
  95 static void setup_curve(float **c,
  96                         int band,
  97                         float *curveatt_dB){
  98   int i,j;
  99   float ath[EHMER_MAX];
 100   float tempc[P_LEVELS][EHMER_MAX];
 101
 102   memcpy(c[0]+2,c[4]+2,sizeof(float)*EHMER_MAX);
 103   memcpy(c[2]+2,c[4]+2,sizeof(float)*EHMER_MAX);
 104
 105   /* we add back in the ATH to avoid low level curves falling off to
 106      -infinity and unneccessarily cutting off high level curves in the
 107      curve limiting (last step).  But again, remember... a half-band's
 108      settings must be valid over the whole band, and it's better to
 109      mask too little than too much, so be pessimal. */
 110
 111   for(i=0;i<EHMER_MAX;i++){
 112     float oc_min=band*.5+(i-EHMER_OFFSET)*.125;
 113     float oc_max=band*.5+(i-EHMER_OFFSET+1)*.125;
 114     float bark=toBARK(fromOC(oc_min));
 115     int ibark=floor(bark);
 116     float del=bark-ibark;
 117     float ath_min,ath_max;
 118
 119     if(ibark<26)
 120       ath_min=ATH_Bark_dB[ibark]*(1.f-del)+ATH_Bark_dB[ibark+1]*del;
 121     else
 122       ath_min=ATH_Bark_dB[25];
 123
 124     bark=toBARK(fromOC(oc_max));
 125     ibark=floor(bark);
 126     del=bark-ibark;
 127
 128     if(ibark<26)
 129       ath_max=ATH_Bark_dB[ibark]*(1.f-del)+ATH_Bark_dB[ibark+1]*del;
 130     else
 131       ath_max=ATH_Bark_dB[25];
 132
 133     ath[i]=min(ath_min,ath_max);
 134   }
 135
 136   /* The c array is comes in as dB curves at 20 40 60 80 100 dB.
 137      interpolate intermediate dB curves */
 138   for(i=1;i<P_LEVELS;i+=2){
 139     interp_curve(c[i]+2,c[i-1]+2,c[i+1]+2,.5);
 140   }
 141
 142   /* normalize curves so the driving amplitude is 0dB */
 143   /* make temp curves with the ATH overlayed */
 144   for(i=0;i<P_LEVELS;i++){
 145     attenuate_curve(c[i]+2,curveatt_dB[i]);
 146     memcpy(tempc[i],ath,EHMER_MAX*sizeof(float));
 147     attenuate_curve(tempc[i],-i*10.f);
 148     max_curve(tempc[i],c[i]+2);
 149   }
 150
 151   /* Now limit the louder curves.
 152
 153      the idea is this: We don't know what the playback attenuation
 154      will be; 0dB SL moves every time the user twiddles the volume
 155      knob. So that means we have to use a single 'most pessimal' curve
 156      for all masking amplitudes, right?  Wrong.  The *loudest* sound
 157      can be in (we assume) a range of ...+100dB] SL.  However, sounds
 158      20dB down will be in a range ...+80], 40dB down is from ...+60],
 159      etc... */
 160
 161   for(j=1;j<P_LEVELS;j++){
 162     min_curve(tempc[j],tempc[j-1]);
 163     min_curve(c[j]+2,tempc[j]);
 164   }
 165
 166   /* add fenceposts */
 167   for(j=0;j<P_LEVELS;j++){
 168
 169     for(i=0;i<EHMER_MAX;i++)
 170       if(c[j][i+2]>-200.f)break;
 171     c[j][0]=i;
 172
 173     for(i=EHMER_MAX-1;i>=0;i--)
 174       if(c[j][i+2]>-200.f)
 175         break;
 176     c[j][1]=i;
 177
 178   }
 179 }
 180
 181 void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi,int n,long rate){
 182   long i,j;
 183   long maxoc;
 184   memset(p,0,sizeof(vorbis_look_psy));
 185
 186
 187   p->eighth_octave_lines=vi->eighth_octave_lines;
 188   p->shiftoc=rint(log(vi->eighth_octave_lines*8)/log(2))-1;
 189
 190   p->firstoc=toOC(.25f*rate/n)*(1<<(p->shiftoc+1))-vi->eighth_octave_lines;
 191   maxoc=toOC((n*.5f-.25f)*rate/n)*(1<<(p->shiftoc+1))+.5f;
 192   p->total_octave_lines=maxoc-p->firstoc+1;
 193
 194   p->ath=_ogg_malloc(n*sizeof(float));
 195   p->octave=_ogg_malloc(n*sizeof(int));
 196   p->bark=_ogg_malloc(n*sizeof(float));
 197   p->vi=vi;
 198   p->n=n;
 199
 200   /* set up the lookups for a given blocksize and sample rate */
 201   /* Vorbis max sample rate is currently limited by 26 Bark (54kHz) */
 202   set_curve(ATH_Bark_dB, p->ath,n,rate);
 203   for(i=0;i<n;i++)
 204     p->bark[i]=toBARK(rate/(2*n)*i);
 205
 206   for(i=0;i<n;i++)
 207     p->octave[i]=toOC((i*.5f+.25f)*rate/n)*(1<<(p->shiftoc+1))+.5f;
 208
 209   p->tonecurves=_ogg_malloc(P_BANDS*sizeof(float **));
 210   p->noisemedian=_ogg_malloc(n*sizeof(float *));
 211   p->noiseoffset=_ogg_malloc(n*sizeof(float *));
 212   p->peakatt=_ogg_malloc(P_BANDS*sizeof(float *));
 213   for(i=0;i<P_BANDS;i++){
 214     p->tonecurves[i]=_ogg_malloc(P_LEVELS*sizeof(float *));
 215     p->peakatt[i]=_ogg_malloc(P_LEVELS*sizeof(float));
 216   }
 217
 218   for(i=0;i<P_BANDS;i++)
 219     for(j=0;j<P_LEVELS;j++){
 220       p->tonecurves[i][j]=_ogg_malloc((EHMER_MAX+2)*sizeof(float));
 221     }
 222
 223   /* OK, yeah, this was a silly way to do it */
 224   memcpy(p->tonecurves[0][4]+2,tone_125_40dB_SL,sizeof(float)*EHMER_MAX);
 225   memcpy(p->tonecurves[0][6]+2,tone_125_60dB_SL,sizeof(float)*EHMER_MAX);
 226   memcpy(p->tonecurves[0][8]+2,tone_125_80dB_SL,sizeof(float)*EHMER_MAX);
 227   memcpy(p->tonecurves[0][10]+2,tone_125_100dB_SL,sizeof(float)*EHMER_MAX);
 228
 229   memcpy(p->tonecurves[2][4]+2,tone_125_40dB_SL,sizeof(float)*EHMER_MAX);
 230   memcpy(p->tonecurves[2][6]+2,tone_125_60dB_SL,sizeof(float)*EHMER_MAX);
 231   memcpy(p->tonecurves[2][8]+2,tone_125_80dB_SL,sizeof(float)*EHMER_MAX);
 232   memcpy(p->tonecurves[2][10]+2,tone_125_100dB_SL,sizeof(float)*EHMER_MAX);
 233
 234   memcpy(p->tonecurves[4][4]+2,tone_250_40dB_SL,sizeof(float)*EHMER_MAX);
 235   memcpy(p->tonecurves[4][6]+2,tone_250_60dB_SL,sizeof(float)*EHMER_MAX);
 236   memcpy(p->tonecurves[4][8]+2,tone_250_80dB_SL,sizeof(float)*EHMER_MAX);
 237   memcpy(p->tonecurves[4][10]+2,tone_250_100dB_SL,sizeof(float)*EHMER_MAX);
 238
 239   memcpy(p->tonecurves[6][4]+2,tone_500_40dB_SL,sizeof(float)*EHMER_MAX);
 240   memcpy(p->tonecurves[6][6]+2,tone_500_60dB_SL,sizeof(float)*EHMER_MAX);
 241   memcpy(p->tonecurves[6][8]+2,tone_500_80dB_SL,sizeof(float)*EHMER_MAX);
 242   memcpy(p->tonecurves[6][10]+2,tone_500_100dB_SL,sizeof(float)*EHMER_MAX);
 243
 244   memcpy(p->tonecurves[8][4]+2,tone_1000_40dB_SL,sizeof(float)*EHMER_MAX);
 245   memcpy(p->tonecurves[8][6]+2,tone_1000_60dB_SL,sizeof(float)*EHMER_MAX);
 246   memcpy(p->tonecurves[8][8]+2,tone_1000_80dB_SL,sizeof(float)*EHMER_MAX);
 247   memcpy(p->tonecurves[8][10]+2,tone_1000_100dB_SL,sizeof(float)*EHMER_MAX);
 248
 249   memcpy(p->tonecurves[10][4]+2,tone_2000_40dB_SL,sizeof(float)*EHMER_MAX);
 250   memcpy(p->tonecurves[10][6]+2,tone_2000_60dB_SL,sizeof(float)*EHMER_MAX);
 251   memcpy(p->tonecurves[10][8]+2,tone_2000_80dB_SL,sizeof(float)*EHMER_MAX);
 252   memcpy(p->tonecurves[10][10]+2,tone_2000_100dB_SL,sizeof(float)*EHMER_MAX);
 253
 254   memcpy(p->tonecurves[12][4]+2,tone_4000_40dB_SL,sizeof(float)*EHMER_MAX);
 255   memcpy(p->tonecurves[12][6]+2,tone_4000_60dB_SL,sizeof(float)*EHMER_MAX);
 256   memcpy(p->tonecurves[12][8]+2,tone_4000_80dB_SL,sizeof(float)*EHMER_MAX);
 257   memcpy(p->tonecurves[12][10]+2,tone_4000_100dB_SL,sizeof(float)*EHMER_MAX);
 258
 259   memcpy(p->tonecurves[14][4]+2,tone_8000_40dB_SL,sizeof(float)*EHMER_MAX);
 260   memcpy(p->tonecurves[14][6]+2,tone_8000_60dB_SL,sizeof(float)*EHMER_MAX);
 261   memcpy(p->tonecurves[14][8]+2,tone_8000_80dB_SL,sizeof(float)*EHMER_MAX);
 262   memcpy(p->tonecurves[14][10]+2,tone_8000_100dB_SL,sizeof(float)*EHMER_MAX);
 263
 264   memcpy(p->tonecurves[16][4]+2,tone_8000_40dB_SL,sizeof(float)*EHMER_MAX);
 265   memcpy(p->tonecurves[16][6]+2,tone_8000_60dB_SL,sizeof(float)*EHMER_MAX);
 266   memcpy(p->tonecurves[16][8]+2,tone_8000_80dB_SL,sizeof(float)*EHMER_MAX);
 267   memcpy(p->tonecurves[16][10]+2,tone_8000_100dB_SL,sizeof(float)*EHMER_MAX);
 268
 269   /* interpolate curves between */
 270   for(i=1;i<P_BANDS;i+=2)
 271     for(j=4;j<P_LEVELS;j+=2){
 272       memcpy(p->tonecurves[i][j]+2,p->tonecurves[i-1][j]+2,EHMER_MAX*sizeof(float));
 273       /*interp_curve(p->tonecurves[i][j],
 274                    p->tonecurves[i-1][j],
 275                    p->tonecurves[i+1][j],.5);*/
 276       min_curve(p->tonecurves[i][j]+2,p->tonecurves[i+1][j]+2);
 277     }
 278
 279   /* set up the final curves */
 280   for(i=0;i<P_BANDS;i++)
 281     setup_curve(p->tonecurves[i],i,vi->toneatt[i]);
 282
 283   /* set up attenuation levels */
 284   for(i=0;i<P_BANDS;i++)
 285     for(j=0;j<P_LEVELS;j++){
 286       p->peakatt[i][j]=p->vi->peakatt[i][j];
 287     }
 288
 289   /* set up rolling noise median */
 290   for(i=0;i<n;i++){
 291     float halfoc=toOC((i+.5)*rate/(2.*n))*2.+2.;
 292     int inthalfoc;
 293     float del;
 294
 295     if(halfoc<0)halfoc=0;
 296     if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1;
 297     inthalfoc=(int)halfoc;
 298     del=halfoc-inthalfoc;
 299
 300     p->noisemedian[i]=
 301       p->vi->noisemedian[inthalfoc*2]*(1.-del) +
 302       p->vi->noisemedian[inthalfoc*2+2]*del;
 303     p->noiseoffset[i]=
 304       p->vi->noisemedian[inthalfoc*2+1]*(1.-del) +
 305       p->vi->noisemedian[inthalfoc*2+3]*del;
 306   }
 307   /*_analysis_output("mediancurve",0,p->noisemedian,n,0,0);*/
 308 }
 309
 310 void _vp_psy_clear(vorbis_look_psy *p){
 311   int i,j;
 312   if(p){
 313     if(p->ath)_ogg_free(p->ath);
 314     if(p->octave)_ogg_free(p->octave);
 315     if(p->bark)_ogg_free(p->bark);
 316     if(p->tonecurves){
 317       for(i=0;i<P_BANDS;i++){
 318         for(j=0;j<P_LEVELS;j++){
 319           _ogg_free(p->tonecurves[i][j]);
 320         }
 321         _ogg_free(p->tonecurves[i]);
 322         _ogg_free(p->peakatt[i]);
 323       }
 324       _ogg_free(p->tonecurves);
 325       _ogg_free(p->noisemedian);
 326       _ogg_free(p->noiseoffset);
 327       _ogg_free(p->peakatt);
 328     }
 329     memset(p,0,sizeof(vorbis_look_psy));
 330   }
 331 }
 332
 333 /* octave/(8*eighth_octave_lines) x scale and dB y scale */
 334 static void seed_curve(float *seed,
 335                       float **curves,
 336                       float amp,
 337                       int oc,int n,int linesper,float dBoffset){
 338   int i;
 339   long seedptr;
 340   float *posts,*curve;
 341
 342   int choice=(int)((amp+dBoffset)*.1f);
 343   choice=max(choice,0);
 344   choice=min(choice,P_LEVELS-1);
 345   posts=curves[choice];
 346   curve=posts+2;
 347   seedptr=oc+(posts[0]-16)*linesper-(linesper>>1);
 348
 349   for(i=posts[0];i<posts[1];i++){
 350     if(seedptr>0){
 351       float lin=amp+curve[i];
 352       if(seed[seedptr]<lin)seed[seedptr]=lin;
 353     }
 354     seedptr+=linesper;
 355     if(seedptr>=n)break;
 356   }
 357 }
 358
 359 static void seed_peak(float *seed,
 360                       float *att,
 361                       float amp,
 362                       int oc,
 363                       int linesper,
 364                       float dBoffset){
 365   long seedptr;
 366
 367   int choice=(int)((amp+dBoffset)*.1f);
 368   choice=max(choice,0);
 369   choice=min(choice,P_LEVELS-1);
 370   seedptr=oc-(linesper>>1);
 371
 372   amp+=att[choice];
 373   if(seed[seedptr]<amp)seed[seedptr]=amp;
 374
 375 }
 376
 377 static void seed_loop(vorbis_look_psy *p,
 378                       float ***curves,
 379                       float **att,
 380                       float *f,
 381                       float *flr,
 382                       float *minseed,
 383                       float *maxseed,
 384                       float specmax){
 385   vorbis_info_psy *vi=p->vi;
 386   long n=p->n,i;
 387   float dBoffset=vi->max_curve_dB-specmax;
 388
 389   /* prime the working vector with peak values */
 390
 391   for(i=0;i<n;i++){
 392       float max=f[i];
 393       long oc=p->octave[i];
 394       while(i+1<n && p->octave[i+1]==oc){
 395         i++;
 396         if(f[i]>max)max=f[i];
 397       }
 398
 399       if(max>flr[i]){
 400         oc=oc>>p->shiftoc;
 401         if(oc>=P_BANDS)oc=P_BANDS-1;
 402         if(oc<0)oc=0;
 403         if(vi->tonemaskp)
 404           seed_curve(minseed,
 405                      curves[oc],
 406                      max,
 407                      p->octave[i]-p->firstoc,
 408                      p->total_octave_lines,
 409                      p->eighth_octave_lines,
 410                      dBoffset);
 411         if(vi->peakattp)
 412           seed_peak(maxseed,
 413                     att[oc],
 414                     max,
 415                     p->octave[i]-p->firstoc,
 416                     p->eighth_octave_lines,
 417                     dBoffset);
 418       }
 419   }
 420 }
 421
 422 static void bound_loop(vorbis_look_psy *p,
 423                        float *f,
 424                        float *seeds,
 425                        float *flr,
 426                        float att){
 427   long n=p->n,i;
 428
 429   long off=(p->eighth_octave_lines>>1)+p->firstoc;
 430   long *ocp=p->octave;
 431
 432   for(i=0;i<n;i++){
 433     long oc=ocp[i]-off;
 434     float v=f[i]+att;
 435     if(seeds[oc]<v)seeds[oc]=v;
 436   }
 437 }
 438
 439 static void seed_chase(float *seeds, int linesper, long n){
 440   long  *posstack=alloca(n*sizeof(long));
 441   float *ampstack=alloca(n*sizeof(float));
 442   long   stack=0;
 443   long   pos=0;
 444   long   i;
 445
 446   for(i=0;i<n;i++){
 447     if(stack<2){
 448       posstack[stack]=i;
 449       ampstack[stack++]=seeds[i];
 450     }else{
 451       while(1){
 452         if(seeds[i]<ampstack[stack-1]){
 453           posstack[stack]=i;
 454           ampstack[stack++]=seeds[i];
 455           break;
 456         }else{
 457           if(i<posstack[stack-1]+linesper){
 458             if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] &&
 459                i<posstack[stack-2]+linesper){
 460               /* we completely overlap, making stack-1 irrelevant.  pop it */
 461               stack--;
 462               continue;
 463             }
 464           }
 465           posstack[stack]=i;
 466           ampstack[stack++]=seeds[i];
 467           break;
 468
 469         }
 470       }
 471     }
 472   }
 473
 474   /* the stack now contains only the positions that are relevant. Scan
 475      'em straight through */
 476
 477   for(i=0;i<stack;i++){
 478     long endpos;
 479     if(i<stack-1 && ampstack[i+1]>ampstack[i]){
 480       endpos=posstack[i+1];
 481     }else{
 482       endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is
 483                                         discarded in short frames */
 484     }
 485     if(endpos>n)endpos=n;
 486     for(;pos<endpos;pos++)
 487       seeds[pos]=ampstack[i];
 488   }
 489
 490   /* there.  Linear time.  I now remember this was on a problem set I
 491      had in Grad Skool... I didn't solve it at the time ;-) */
 492
 493 }
 494
 495 /* bleaugh, this is more complicated than it needs to be */
 496 static void max_seeds(vorbis_look_psy *p,float *minseed,float *maxseed,
 497                       float *flr){
 498   long   n=p->total_octave_lines;
 499   int    linesper=p->eighth_octave_lines;
 500   long   linpos=0;
 501   long   pos;
 502
 503   seed_chase(minseed,linesper,n); /* for masking */
 504   seed_chase(maxseed,linesper,n); /* for peak att */
 505
 506   pos=p->octave[0]-p->firstoc-(linesper>>1);
 507   while(linpos+1<p->n){
 508     float min=minseed[pos];
 509     float max=maxseed[pos];
 510     long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc;
 511     while(pos+1<=end){
 512       pos++;
 513       if((minseed[pos]>NEGINF && minseed[pos]<min) || min==NEGINF)
 514         min=minseed[pos];
 515       if(maxseed[pos]>max)max=maxseed[pos];
 516     }
 517     if(max<min)max=min;
 518
 519     /* seed scale is log.  Floor is linear.  Map back to it */
 520     end=pos+p->firstoc;
 521     for(;linpos<p->n && p->octave[linpos]<=end;linpos++)
 522       if(flr[linpos]<max)flr[linpos]=max;
 523   }
 524
 525   {
 526     float min=minseed[p->total_octave_lines-1];
 527     float max=maxseed[p->total_octave_lines-1];
 528     if(max<min)max=min;
 529     for(;linpos<p->n;linpos++)
 530       if(flr[linpos]<max)flr[linpos]=max;
 531   }
 532
 533 }
 534
 535 /* quarter-dB bins */
 536 #define BIN(x)   ((int)((x)*negFour))
 537 #define BINdB(x) ((x)*negQuarter)
 538 #define BINCOUNT (200*4)
 539 #define LASTBIN  (BINCOUNT-1)
 540
 541 static void bark_noise_median(long n,float *b,float *f,float *noise,
 542                               float lowidth,float hiwidth,
 543                               int lomin,int himin,
 544                               float *thresh,float *off){
 545   long i=0,lo=0,hi=0;
 546   float bi,threshi;
 547   long median=LASTBIN;
 548   float negFour = -4.0f;
 549   float negQuarter = -0.25f;
 550
 551    /* these are really integral values, but we store them in floats to
 552       avoid excessive float/int conversions, which GCC and MSVC are
 553       farily poor at optimizing. */
 554
 555   float radix[BINCOUNT];
 556   float countabove=0;
 557   float countbelow=0;
 558
 559   memset(radix,0,sizeof(radix));
 560
 561   for(i=0;i<n;i++){
 562     /* find new lo/hi */
 563     bi=b[i];
 564     for(;hi<n && (hi<i+himin || b[hi]<=bi+hiwidth);hi++){
 565       int bin=BIN(f[hi]);
 566       if(bin>LASTBIN)bin=LASTBIN;
 567       radix[bin]++;
 568       if(bin<median)
 569         countabove++;
 570       else
 571         countbelow++;
 572     }
 573     for(;lo<i && lo+lomin<i && b[lo]+lowidth<=bi;lo++){
 574       int bin=BIN(f[lo]);
 575       if(bin>LASTBIN)bin=LASTBIN;
 576       radix[bin]--;
 577       if(bin<median)
 578         countabove--;
 579       else
 580         countbelow--;
 581     }
 582
 583     /* move the median if needed */
 584     if(countabove+countbelow){
 585       threshi = thresh[i];
 586
 587       while((countabove+countbelow)*threshi>countbelow && median>0){
 588         median--;
 589         countabove-=radix[median];
 590         countbelow+=radix[median];
 591       }
 592
 593       while((countabove+countbelow)*thresh[i]<(countbelow-radix[median])
 594             && median+1<BINCOUNT){
 595         countabove+=radix[median];
 596         countbelow-=radix[median];
 597         median++;
 598       }
 599     }
 600     noise[i]=BINdB(median)+off[i];
 601   }
 602
 603 }
 604
 605 float _vp_compute_mask(vorbis_look_psy *p,
 606                       float *fft,
 607                       float *mdct,
 608                       float *flr,
 609                       float *decay,
 610                       float specmax){
 611   int i,n=p->n;
 612   float localmax=NEGINF;
 613   static int seq=0;
 614
 615   float *minseed=alloca(sizeof(float)*p->total_octave_lines);
 616   float *maxseed=alloca(sizeof(float)*p->total_octave_lines);
 617   for(i=0;i<p->total_octave_lines;i++)minseed[i]=maxseed[i]=NEGINF;
 618
 619   /* go to dB scale. Also find the highest peak so we know the limits */
 620   for(i=0;i<n;i++){
 621     fft[i]=todB(fft[i]);
 622     if(fft[i]>localmax)localmax=fft[i];
 623   }
 624   if(specmax<localmax)specmax=localmax;
 625
 626
 627   for(i=0;i<n;i++){
 628     mdct[i]=todB(mdct[i]);
 629   }
 630
 631   _analysis_output("mdct",seq,mdct,n,0,0);
 632   _analysis_output("fft",seq,fft,n,0,0);
 633
 634   /* noise masking */
 635   if(p->vi->noisemaskp){
 636     bark_noise_median(n,p->bark,mdct,flr,
 637                       p->vi->noisewindowlo,
 638                       p->vi->noisewindowhi,
 639                       p->vi->noisewindowlomin,
 640                       p->vi->noisewindowhimin,
 641                       p->noisemedian,
 642                       p->noiseoffset);
 643     /* suppress any noise curve > specmax+p->vi->noisemaxsupp */
 644     for(i=0;i<n;i++)
 645       if(flr[i]>specmax+p->vi->noisemaxsupp)
 646         flr[i]=specmax+p->vi->noisemaxsupp;
 647     _analysis_output("noise",seq,flr,n,0,0);
 648   }else{
 649     for(i=0;i<n;i++)flr[i]=NEGINF;
 650   }
 651
 652   /* set the ATH (floating below localmax, not global max by a
 653      specified att) */
 654   if(p->vi->athp){
 655     float att=localmax+p->vi->ath_adjatt;
 656     if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt;
 657
 658     for(i=0;i<n;i++){
 659       float av=p->ath[i]+att;
 660       if(av>flr[i])flr[i]=av;
 661     }
 662   }
 663
 664   _analysis_output("ath",seq,flr,n,0,0);
 665
 666   /* tone/peak masking */
 667
 668   /* XXX apply decay to the fft here */
 669
 670   seed_loop(p,p->tonecurves,p->peakatt,fft,flr,minseed,maxseed,specmax);
 671   bound_loop(p,mdct,maxseed,flr,p->vi->bound_att_dB);
 672   _analysis_output("minseed",seq,minseed,p->total_octave_lines,0,0);
 673   _analysis_output("maxseed",seq,maxseed,p->total_octave_lines,0,0);
 674   max_seeds(p,minseed,maxseed,flr);
 675   _analysis_output("final",seq,flr,n,0,0);
 676
 677   /* doing this here is clean, but we need to find a faster way to do
 678      it than to just tack it on */
 679
 680   for(i=0;i<n;i++)if(mdct[i]>=flr[i])break;
 681   if(i==n)for(i=0;i<n;i++)flr[i]=NEGINF;
 682
 683
 684   seq++;
 685
 686   return(specmax);
 687 }
 688
 689
 690 /* this applies the floor and (optionally) tries to preserve noise
 691    energy in low resolution portions of the spectrum */
 692 /* f and flr are *linear* scale, not dB */
 693 void _vp_apply_floor(vorbis_look_psy *p,float *f, float *flr){
 694   float *work=alloca(p->n*sizeof(float));
 695   int j;
 696
 697   /* subtract the floor */
 698   for(j=0;j<p->n;j++){
 699     if(flr[j]<=0)
 700       work[j]=0.f;
 701     else
 702       work[j]=f[j]/flr[j];
 703   }
 704
 705   memcpy(f,work,p->n*sizeof(float));
 706 }
 707
 708 float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){
 709   vorbis_info *vi=vd->vi;
 710   codec_setup_info *ci=vi->codec_setup;
 711   int n=ci->blocksizes[vd->W]/2;
 712   float secs=(float)n/vi->rate;
 713
 714   amp+=secs*ci->ampmax_att_per_sec;
 715   if(amp<-9999)amp=-9999;
 716   return(amp);
 717 }
 718
 719
 720