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.35 2001/01/22 01:38:25 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->peakatt);
 327     }
 328     memset(p,0,sizeof(vorbis_look_psy));
 329   }
 330 }
 331
 332 /* octave/(8*eighth_octave_lines) x scale and dB y scale */
 333 static void seed_curve(float *seed,
 334                       float **curves,
 335                       float amp,
 336                       int oc,int n,int linesper,float dBoffset){
 337   int i;
 338   long seedptr;
 339   float *posts,*curve;
 340
 341   int choice=(int)((amp+dBoffset)*.1f);
 342   choice=max(choice,0);
 343   choice=min(choice,P_LEVELS-1);
 344   posts=curves[choice];
 345   curve=posts+2;
 346   seedptr=oc+(posts[0]-16)*linesper-(linesper>>1);
 347
 348   for(i=posts[0];i<posts[1];i++){
 349     if(seedptr>0){
 350       float lin=amp+curve[i];
 351       if(seed[seedptr]<lin)seed[seedptr]=lin;
 352     }
 353     seedptr+=linesper;
 354     if(seedptr>=n)break;
 355   }
 356 }
 357
 358 static void seed_peak(float *seed,
 359                       float *att,
 360                       float amp,
 361                       int oc,
 362                       int linesper,
 363                       float dBoffset){
 364   long seedptr;
 365
 366   int choice=(int)((amp+dBoffset)*.1f);
 367   choice=max(choice,0);
 368   choice=min(choice,P_LEVELS-1);
 369   seedptr=oc-(linesper>>1);
 370
 371   amp+=att[choice];
 372   if(seed[seedptr]<amp)seed[seedptr]=amp;
 373
 374 }
 375
 376 static void seed_loop(vorbis_look_psy *p,
 377                       float ***curves,
 378                       float **att,
 379                       float *f,
 380                       float *flr,
 381                       float *minseed,
 382                       float *maxseed,
 383                       float specmax){
 384   vorbis_info_psy *vi=p->vi;
 385   long n=p->n,i;
 386   float dBoffset=vi->max_curve_dB-specmax;
 387
 388   /* prime the working vector with peak values */
 389
 390   for(i=0;i<n;i++){
 391       float max=f[i];
 392       long oc=p->octave[i];
 393       while(i+1<n && p->octave[i+1]==oc){
 394         i++;
 395         if(f[i]>max)max=f[i];
 396       }
 397
 398       if(max>flr[i]){
 399         oc=oc>>p->shiftoc;
 400         if(oc>=P_BANDS)oc=P_BANDS-1;
 401         if(oc<0)oc=0;
 402         if(vi->tonemaskp)
 403           seed_curve(minseed,
 404                      curves[oc],
 405                      max,
 406                      p->octave[i]-p->firstoc,
 407                      p->total_octave_lines,
 408                      p->eighth_octave_lines,
 409                      dBoffset);
 410         if(vi->peakattp)
 411           seed_peak(maxseed,
 412                     att[oc],
 413                     max,
 414                     p->octave[i]-p->firstoc,
 415                     p->eighth_octave_lines,
 416                     dBoffset);
 417       }
 418   }
 419 }
 420
 421 static void bound_loop(vorbis_look_psy *p,
 422                        float *f,
 423                        float *seeds,
 424                        float *flr,
 425                        float att){
 426   long n=p->n,i;
 427
 428   long off=(p->eighth_octave_lines>>1)+p->firstoc;
 429   long *ocp=p->octave;
 430
 431   for(i=0;i<n;i++){
 432     long oc=ocp[i]-off;
 433     float v=f[i]+att;
 434     if(seeds[oc]<v)seeds[oc]=v;
 435   }
 436 }
 437
 438 static void seed_chase(float *seeds, int linesper, long n){
 439   long  *posstack=alloca(n*sizeof(long));
 440   float *ampstack=alloca(n*sizeof(float));
 441   long   stack=0;
 442   long   pos=0;
 443   long   i;
 444
 445   for(i=0;i<n;i++){
 446     if(stack<2){
 447       posstack[stack]=i;
 448       ampstack[stack++]=seeds[i];
 449     }else{
 450       while(1){
 451         if(seeds[i]<ampstack[stack-1]){
 452           posstack[stack]=i;
 453           ampstack[stack++]=seeds[i];
 454           break;
 455         }else{
 456           if(i<posstack[stack-1]+linesper){
 457             if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] &&
 458                i<posstack[stack-2]+linesper){
 459               /* we completely overlap, making stack-1 irrelevant.  pop it */
 460               stack--;
 461               continue;
 462             }
 463           }
 464           posstack[stack]=i;
 465           ampstack[stack++]=seeds[i];
 466           break;
 467
 468         }
 469       }
 470     }
 471   }
 472
 473   /* the stack now contains only the positions that are relevant. Scan
 474      'em straight through */
 475
 476   for(i=0;i<stack;i++){
 477     long endpos;
 478     if(i<stack-1 && ampstack[i+1]>ampstack[i]){
 479       endpos=posstack[i+1];
 480     }else{
 481       endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is
 482                                         discarded in short frames */
 483     }
 484     if(endpos>n)endpos=n;
 485     for(;pos<endpos;pos++)
 486       seeds[pos]=ampstack[i];
 487   }
 488
 489   /* there.  Linear time.  I now remember this was on a problem set I
 490      had in Grad Skool... I didn't solve it at the time ;-) */
 491
 492 }
 493
 494 /* bleaugh, this is more complicated than it needs to be */
 495 static void max_seeds(vorbis_look_psy *p,float *minseed,float *maxseed,
 496                       float *flr){
 497   long   n=p->total_octave_lines;
 498   int    linesper=p->eighth_octave_lines;
 499   long   linpos=0;
 500   long   pos;
 501
 502   seed_chase(minseed,linesper,n); /* for masking */
 503   seed_chase(maxseed,linesper,n); /* for peak att */
 504
 505   pos=p->octave[0]-p->firstoc-(linesper>>1);
 506   while(linpos+1<p->n){
 507     float min=minseed[pos];
 508     float max=maxseed[pos];
 509     long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc;
 510     while(pos+1<=end){
 511       pos++;
 512       if((minseed[pos]>NEGINF && minseed[pos]<min) || min==NEGINF)
 513         min=minseed[pos];
 514       if(maxseed[pos]>max)max=maxseed[pos];
 515     }
 516     if(max<min)max=min;
 517
 518     /* seed scale is log.  Floor is linear.  Map back to it */
 519     end=pos+p->firstoc;
 520     for(;linpos<p->n && p->octave[linpos]<=end;linpos++)
 521       if(flr[linpos]<max)flr[linpos]=max;
 522   }
 523
 524   {
 525     float min=minseed[p->total_octave_lines-1];
 526     float max=maxseed[p->total_octave_lines-1];
 527     if(max<min)max=min;
 528     for(;linpos<p->n;linpos++)
 529       if(flr[linpos]<max)flr[linpos]=max;
 530   }
 531
 532 }
 533
 534 #define BIN(x) ((int)((x)*-4.))
 535 #define BINdB(x) ((x)*-.25)
 536 static void bark_noise_median(long n,float *b,float *f,float *noise,
 537                               float lowidth,float hiwidth,
 538                               int lomin,int himin,
 539                               float *thresh,float *off){
 540   long i=0,lo=0,hi=0;
 541   long *radix=alloca(200*4*sizeof(long)); /* quarter-dB bins */
 542
 543   long countabove=0;
 544   long median=200*4-1;
 545   long countbelow=0;
 546
 547   memset(radix,0,200*4*sizeof(long));
 548
 549   for(i=0;i<n;i++){
 550     /* find new lo/hi */
 551     for(;hi<n && (b[hi]<=b[i]+hiwidth || hi<i+himin);hi++){
 552       int bin=BIN(f[hi]);
 553       if(bin>=200*4)bin=200*4-1;
 554       radix[bin]++;
 555       if(bin<median)
 556         countabove++;
 557       else
 558         countbelow++;
 559     }
 560     for(;lo<i && b[lo]+lowidth<=b[i] && lo+lomin<i;lo++){
 561       int bin=BIN(f[lo]);
 562       if(bin>=200*4)bin=200*4-1;
 563       radix[bin]--;
 564       if(bin<median)
 565         countabove--;
 566       else
 567         countbelow--;
 568     }
 569
 570     /* move the median if needed */
 571     if(countabove+countbelow){
 572
 573       while(thresh[i]>countbelow/(float)(countabove+countbelow) && median>0){
 574         median--;
 575         countabove-=radix[median];
 576         countbelow+=radix[median];
 577       }
 578
 579       while(thresh[i]<(countbelow-radix[median])/
 580             (float)(countabove+countbelow) && median+1<200*4){
 581         countabove+=radix[median];
 582         countbelow-=radix[median];
 583         median++;
 584       }
 585     }
 586     noise[i]=BINdB(median)+off[i];
 587   }
 588
 589 }
 590
 591 float _vp_compute_mask(vorbis_look_psy *p,
 592                       float *fft,
 593                       float *mdct,
 594                       float *flr,
 595                       float *decay,
 596                       float prev_maxamp){
 597   int i,n=p->n;
 598   float specmax=NEGINF;
 599   static int seq=0;
 600
 601   float *minseed=alloca(sizeof(float)*p->total_octave_lines);
 602   float *maxseed=alloca(sizeof(float)*p->total_octave_lines);
 603   for(i=0;i<p->total_octave_lines;i++)minseed[i]=maxseed[i]=NEGINF;
 604
 605   /* go to dB scale. Also find the highest peak so we know the limits */
 606   for(i=0;i<n;i++){
 607     fft[i]=todB(fft[i]);
 608     if(fft[i]>specmax)specmax=fft[i];
 609   }
 610   if(specmax<prev_maxamp)specmax=prev_maxamp;
 611
 612
 613   for(i=0;i<n;i++){
 614     mdct[i]=todB(mdct[i]);
 615   }
 616
 617   _analysis_output("mdct",seq,mdct,n,0,0);
 618   _analysis_output("fft",seq,fft,n,0,0);
 619
 620   /* noise masking */
 621   if(p->vi->noisemaskp){
 622     bark_noise_median(n,p->bark,mdct,flr,
 623                       p->vi->noisewindowlo,
 624                       p->vi->noisewindowhi,
 625                       p->vi->noisewindowlomin,
 626                       p->vi->noisewindowhimin,
 627                       p->noisemedian,
 628                       p->noiseoffset);
 629     /* suppress any noise curve > specmax+p->vi->noisemaxsupp */
 630     for(i=0;i<n;i++)
 631       if(flr[i]>specmax+p->vi->noisemaxsupp)
 632         flr[i]=specmax+p->vi->noisemaxsupp;
 633     _analysis_output("noise",seq,flr,n,0,0);
 634   }else{
 635     for(i=0;i<n;i++)flr[i]=NEGINF;
 636   }
 637
 638   /* set the ATH (floating below specmax by a specified att) */
 639   if(p->vi->athp){
 640     float att=specmax+p->vi->ath_adjatt;
 641     if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt;
 642
 643     for(i=0;i<n;i++){
 644       float av=p->ath[i]+att;
 645       if(av>flr[i])flr[i]=av;
 646     }
 647   }
 648
 649   _analysis_output("ath",seq,flr,n,0,0);
 650
 651   /* tone/peak masking */
 652
 653   /* XXX apply decay to the fft here */
 654
 655   seed_loop(p,p->tonecurves,p->peakatt,fft,flr,minseed,maxseed,specmax);
 656   bound_loop(p,mdct,maxseed,flr,p->vi->bound_att_dB);
 657   _analysis_output("minseed",seq,minseed,p->total_octave_lines,0,0);
 658   _analysis_output("maxseed",seq,maxseed,p->total_octave_lines,0,0);
 659   max_seeds(p,minseed,maxseed,flr);
 660   _analysis_output("final",seq,flr,n,0,0);
 661
 662   /* doing this here is clean, but we need to find a faster way to do
 663      it than to just tack it on */
 664
 665   for(i=0;i<n;i++)if(mdct[i]>=flr[i])break;
 666   if(i==n)for(i=0;i<n;i++)flr[i]=NEGINF;
 667
 668
 669   seq++;
 670
 671   return(specmax);
 672 }
 673
 674
 675 /* this applies the floor and (optionally) tries to preserve noise
 676    energy in low resolution portions of the spectrum */
 677 /* f and flr are *linear* scale, not dB */
 678 void _vp_apply_floor(vorbis_look_psy *p,float *f, float *flr){
 679   float *work=alloca(p->n*sizeof(float));
 680   int j;
 681
 682   /* subtract the floor */
 683   for(j=0;j<p->n;j++){
 684     if(flr[j]<=0)
 685       work[j]=0.f;
 686     else
 687       work[j]=f[j]/flr[j];
 688   }
 689
 690   memcpy(f,work,p->n*sizeof(float));
 691 }
 692
 693 float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){
 694   vorbis_info *vi=vd->vi;
 695   codec_setup_info *ci=vi->codec_setup;
 696   int n=ci->blocksizes[vd->W]/2;
 697   float secs=(float)n/vi->rate;
 698
 699   amp+=secs*ci->ampmax_att_per_sec;
 700   if(amp<-9999)amp=-9999;
 701   return(amp);
 702 }
 703
 704
 705