lib/psy.c

   1 /********************************************************************
   2  *                                                                  *
   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.       *
   7  *                                                                  *
   8  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001             *
   9  * by the XIPHOPHORUS Company http://www.xiph.org/                  *
  10  *                                                                  *
  11  ********************************************************************
  12
  13  function: psychoacoustics not including preecho
  14  last mod: $Id: psy.c,v 1.56 2001/10/02 00:14:32 segher Exp $
  15
  16  ********************************************************************/
  17
  18 #include <stdlib.h>
  19 #include <math.h>
  20 #include <string.h>
  21 #include "vorbis/codec.h"
  22 #include "codec_internal.h"
  23
  24 #include "masking.h"
  25 #include "psy.h"
  26 #include "os.h"
  27 #include "lpc.h"
  28 #include "smallft.h"
  29 #include "scales.h"
  30 #include "misc.h"
  31
  32 #define NEGINF -9999.f
  33
  34 /* Why Bark scale for encoding but not masking computation? Because
  35    masking has a strong harmonic dependency */
  36
  37 vorbis_look_psy_global *_vp_global_look(vorbis_info *vi){
  38   int i,j;
  39   codec_setup_info *ci=vi->codec_setup;
  40   vorbis_info_psy_global *gi=ci->psy_g_param;
  41   vorbis_look_psy_global *look=_ogg_calloc(1,sizeof(*look));
  42
  43   int shiftoc=rint(log(gi->eighth_octave_lines*8)/log(2))-1;
  44   look->decaylines=toOC(96000.f)*(1<<(shiftoc+1))+.5f; /* max sample
  45                                                           rate of
  46                                                           192000kHz
  47                                                           for now */
  48   look->decay=_ogg_calloc(vi->channels,sizeof(*look->decay));
  49   for(i=0;i<vi->channels;i++){
  50     look->decay[i]=_ogg_calloc(look->decaylines,sizeof(*look->decay[i]));
  51     for(j=0;j<look->decaylines;j++)
  52       look->decay[i][j]=-9999.;
  53   }
  54   look->channels=vi->channels;
  55
  56   look->ampmax=-9999.;
  57   look->gi=gi;
  58   return(look);
  59 }
  60
  61 void _vp_global_free(vorbis_look_psy_global *look){
  62   int i;
  63   if(look->decay){
  64     for(i=0;i<look->channels;i++)
  65       _ogg_free(look->decay[i]);
  66     _ogg_free(look->decay);
  67   }
  68   memset(look,0,sizeof(*look));
  69   _ogg_free(look);
  70 }
  71
  72 void _vi_psy_free(vorbis_info_psy *i){
  73   if(i){
  74     memset(i,0,sizeof(*i));
  75     _ogg_free(i);
  76   }
  77 }
  78
  79 vorbis_info_psy *_vi_psy_copy(vorbis_info_psy *i){
  80   vorbis_info_psy *ret=_ogg_malloc(sizeof(*ret));
  81   memcpy(ret,i,sizeof(*ret));
  82   return(ret);
  83 }
  84
  85 /* Set up decibel threshold slopes on a Bark frequency scale */
  86 /* ATH is the only bit left on a Bark scale.  No reason to change it
  87    right now */
  88 static void set_curve(float *ref,float *c,int n, float crate){
  89   int i,j=0;
  90
  91   for(i=0;i<MAX_BARK-1;i++){
  92     int endpos=rint(fromBARK(i+1)*2*n/crate);
  93     float base=ref[i];
  94     if(j<endpos){
  95       float delta=(ref[i+1]-base)/(endpos-j);
  96       for(;j<endpos && j<n;j++){
  97         c[j]=base;
  98         base+=delta;
  99       }
 100     }
 101   }
 102 }
 103
 104 static void min_curve(float *c,
 105                        float *c2){
 106   int i;
 107   for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i];
 108 }
 109 static void max_curve(float *c,
 110                        float *c2){
 111   int i;
 112   for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i];
 113 }
 114
 115 static void attenuate_curve(float *c,float att){
 116   int i;
 117   for(i=0;i<EHMER_MAX;i++)
 118     c[i]+=att;
 119 }
 120
 121 static void interp_curve(float *c,float *c1,float *c2,float del){
 122   int i;
 123   for(i=0;i<EHMER_MAX;i++)
 124     c[i]=c2[i]*del+c1[i]*(1.f-del);
 125 }
 126
 127 extern int analysis_noisy;
 128 static void setup_curve(float **c,
 129                         int band,
 130                         float *curveatt_dB){
 131   int i,j;
 132   float ath[EHMER_MAX];
 133   float tempc[P_LEVELS][EHMER_MAX];
 134   float *ATH=ATH_Bark_dB_lspconservative; /* just for limiting here */
 135
 136   memcpy(c[0]+2,c[4]+2,sizeof(*c[0])*EHMER_MAX);
 137   memcpy(c[2]+2,c[4]+2,sizeof(*c[2])*EHMER_MAX);
 138
 139   /* we add back in the ATH to avoid low level curves falling off to
 140      -infinity and unnecessarily cutting off high level curves in the
 141      curve limiting (last step).  But again, remember... a half-band's
 142      settings must be valid over the whole band, and it's better to
 143      mask too little than too much, so be pessimistical. */
 144
 145   for(i=0;i<EHMER_MAX;i++){
 146     float oc_min=band*.5+(i-EHMER_OFFSET)*.125;
 147     float oc_max=band*.5+(i-EHMER_OFFSET+1)*.125;
 148     float bark=toBARK(fromOC(oc_min));
 149     int ibark=floor(bark);
 150     float del=bark-ibark;
 151     float ath_min,ath_max;
 152
 153     if(ibark<26)
 154       ath_min=ATH[ibark]*(1.f-del)+ATH[ibark+1]*del;
 155     else
 156       ath_min=ATH[25];
 157
 158     bark=toBARK(fromOC(oc_max));
 159     ibark=floor(bark);
 160     del=bark-ibark;
 161
 162     if(ibark<26)
 163       ath_max=ATH[ibark]*(1.f-del)+ATH[ibark+1]*del;
 164     else
 165       ath_max=ATH[25];
 166
 167     ath[i]=min(ath_min,ath_max);
 168   }
 169
 170   /* The c array comes in as dB curves at 20 40 60 80 100 dB.
 171      interpolate intermediate dB curves */
 172   for(i=1;i<P_LEVELS;i+=2){
 173     interp_curve(c[i]+2,c[i-1]+2,c[i+1]+2,.5);
 174   }
 175
 176   /* normalize curves so the driving amplitude is 0dB */
 177   /* make temp curves with the ATH overlayed */
 178   for(i=0;i<P_LEVELS;i++){
 179     attenuate_curve(c[i]+2,curveatt_dB[i]);
 180     memcpy(tempc[i],ath,EHMER_MAX*sizeof(*tempc[i]));
 181     attenuate_curve(tempc[i],-i*10.f);
 182     max_curve(tempc[i],c[i]+2);
 183   }
 184
 185   /* Now limit the louder curves.
 186
 187      the idea is this: We don't know what the playback attenuation
 188      will be; 0dB SL moves every time the user twiddles the volume
 189      knob. So that means we have to use a single 'most pessimal' curve
 190      for all masking amplitudes, right?  Wrong.  The *loudest* sound
 191      can be in (we assume) a range of ...+100dB] SL.  However, sounds
 192      20dB down will be in a range ...+80], 40dB down is from ...+60],
 193      etc... */
 194
 195   for(j=1;j<P_LEVELS;j++){
 196     min_curve(tempc[j],tempc[j-1]);
 197     min_curve(c[j]+2,tempc[j]);
 198   }
 199
 200   /* add fenceposts */
 201   for(j=0;j<P_LEVELS;j++){
 202
 203     for(i=0;i<EHMER_OFFSET;i++)
 204       if(c[j][i+2]>-200.f)break;
 205     c[j][0]=i;
 206
 207     for(i=EHMER_MAX-1;i>EHMER_OFFSET+1;i--)
 208       if(c[j][i+2]>-200.f)
 209         break;
 210     c[j][1]=i;
 211
 212   }
 213 }
 214
 215 void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi,
 216                   vorbis_info_psy_global *gi,int n,long rate){
 217   long i,j,k,lo=-99,hi=0;
 218   long maxoc;
 219   memset(p,0,sizeof(*p));
 220
 221
 222   p->eighth_octave_lines=gi->eighth_octave_lines;
 223   p->shiftoc=rint(log(gi->eighth_octave_lines*8)/log(2))-1;
 224
 225   p->firstoc=toOC(.25f*rate/n)*(1<<(p->shiftoc+1))-gi->eighth_octave_lines;
 226   maxoc=toOC((n*.5f-.25f)*rate/n)*(1<<(p->shiftoc+1))+.5f;
 227   p->total_octave_lines=maxoc-p->firstoc+1;
 228
 229   if(vi->ath)
 230     p->ath=_ogg_malloc(n*sizeof(*p->ath));
 231   p->octave=_ogg_malloc(n*sizeof(*p->octave));
 232   p->bark=_ogg_malloc(n*sizeof(*p->bark));
 233   p->vi=vi;
 234   p->n=n;
 235   p->rate=rate;
 236
 237   /* set up the lookups for a given blocksize and sample rate */
 238   if(vi->ath)
 239     set_curve(vi->ath, p->ath,n,rate);
 240   for(i=0;i<n;i++){
 241     float bark=toBARK(rate/(2*n)*i);
 242
 243     for(;lo+vi->noisewindowlomin<i &&
 244           toBARK(rate/(2*n)*lo)<(bark-vi->noisewindowlo);lo++);
 245
 246     for(;hi<n && (hi<i+vi->noisewindowhimin ||
 247           toBARK(rate/(2*n)*hi)<(bark+vi->noisewindowhi));hi++);
 248
 249     p->bark[i]=(lo<<16)+hi;
 250
 251   }
 252
 253   for(i=0;i<n;i++)
 254     p->octave[i]=toOC((i*.5f+.25f)*rate/n)*(1<<(p->shiftoc+1))+.5f;
 255
 256   p->tonecurves=_ogg_malloc(P_BANDS*sizeof(*p->tonecurves));
 257   p->noisethresh=_ogg_malloc(n*sizeof(*p->noisethresh));
 258   p->noiseoffset=_ogg_malloc(n*sizeof(*p->noiseoffset));
 259   for(i=0;i<P_BANDS;i++)
 260     p->tonecurves[i]=_ogg_malloc(P_LEVELS*sizeof(*p->tonecurves[i]));
 261
 262   for(i=0;i<P_BANDS;i++)
 263     for(j=0;j<P_LEVELS;j++)
 264       p->tonecurves[i][j]=_ogg_malloc((EHMER_MAX+2)*sizeof(*p->tonecurves[i][j]));
 265
 266
 267   /* OK, yeah, this was a silly way to do it */
 268   memcpy(p->tonecurves[0][4]+2,tone_125_40dB_SL,sizeof(*p->tonecurves[0][4])*EHMER_MAX);
 269   memcpy(p->tonecurves[0][6]+2,tone_125_60dB_SL,sizeof(*p->tonecurves[0][6])*EHMER_MAX);
 270   memcpy(p->tonecurves[0][8]+2,tone_125_80dB_SL,sizeof(*p->tonecurves[0][8])*EHMER_MAX);
 271   memcpy(p->tonecurves[0][10]+2,tone_125_100dB_SL,sizeof(*p->tonecurves[0][10])*EHMER_MAX);
 272
 273   memcpy(p->tonecurves[2][4]+2,tone_125_40dB_SL,sizeof(*p->tonecurves[2][4])*EHMER_MAX);
 274   memcpy(p->tonecurves[2][6]+2,tone_125_60dB_SL,sizeof(*p->tonecurves[2][6])*EHMER_MAX);
 275   memcpy(p->tonecurves[2][8]+2,tone_125_80dB_SL,sizeof(*p->tonecurves[2][8])*EHMER_MAX);
 276   memcpy(p->tonecurves[2][10]+2,tone_125_100dB_SL,sizeof(*p->tonecurves[2][10])*EHMER_MAX);
 277
 278   memcpy(p->tonecurves[4][4]+2,tone_250_40dB_SL,sizeof(*p->tonecurves[4][4])*EHMER_MAX);
 279   memcpy(p->tonecurves[4][6]+2,tone_250_60dB_SL,sizeof(*p->tonecurves[4][6])*EHMER_MAX);
 280   memcpy(p->tonecurves[4][8]+2,tone_250_80dB_SL,sizeof(*p->tonecurves[4][8])*EHMER_MAX);
 281   memcpy(p->tonecurves[4][10]+2,tone_250_100dB_SL,sizeof(*p->tonecurves[4][10])*EHMER_MAX);
 282
 283   memcpy(p->tonecurves[6][4]+2,tone_500_40dB_SL,sizeof(*p->tonecurves[6][4])*EHMER_MAX);
 284   memcpy(p->tonecurves[6][6]+2,tone_500_60dB_SL,sizeof(*p->tonecurves[6][6])*EHMER_MAX);
 285   memcpy(p->tonecurves[6][8]+2,tone_500_80dB_SL,sizeof(*p->tonecurves[6][8])*EHMER_MAX);
 286   memcpy(p->tonecurves[6][10]+2,tone_500_100dB_SL,sizeof(*p->tonecurves[6][10])*EHMER_MAX);
 287
 288   memcpy(p->tonecurves[8][4]+2,tone_1000_40dB_SL,sizeof(*p->tonecurves[8][4])*EHMER_MAX);
 289   memcpy(p->tonecurves[8][6]+2,tone_1000_60dB_SL,sizeof(*p->tonecurves[8][6])*EHMER_MAX);
 290   memcpy(p->tonecurves[8][8]+2,tone_1000_80dB_SL,sizeof(*p->tonecurves[8][8])*EHMER_MAX);
 291   memcpy(p->tonecurves[8][10]+2,tone_1000_100dB_SL,sizeof(*p->tonecurves[8][10])*EHMER_MAX);
 292
 293   memcpy(p->tonecurves[10][4]+2,tone_2000_40dB_SL,sizeof(*p->tonecurves[10][4])*EHMER_MAX);
 294   memcpy(p->tonecurves[10][6]+2,tone_2000_60dB_SL,sizeof(*p->tonecurves[10][6])*EHMER_MAX);
 295   memcpy(p->tonecurves[10][8]+2,tone_2000_80dB_SL,sizeof(*p->tonecurves[10][8])*EHMER_MAX);
 296   memcpy(p->tonecurves[10][10]+2,tone_2000_100dB_SL,sizeof(*p->tonecurves[10][10])*EHMER_MAX);
 297
 298   memcpy(p->tonecurves[12][4]+2,tone_4000_40dB_SL,sizeof(*p->tonecurves[12][4])*EHMER_MAX);
 299   memcpy(p->tonecurves[12][6]+2,tone_4000_60dB_SL,sizeof(*p->tonecurves[12][6])*EHMER_MAX);
 300   memcpy(p->tonecurves[12][8]+2,tone_4000_80dB_SL,sizeof(*p->tonecurves[12][8])*EHMER_MAX);
 301   memcpy(p->tonecurves[12][10]+2,tone_4000_100dB_SL,sizeof(*p->tonecurves[12][10])*EHMER_MAX);
 302
 303   memcpy(p->tonecurves[14][4]+2,tone_8000_40dB_SL,sizeof(*p->tonecurves[14][4])*EHMER_MAX);
 304   memcpy(p->tonecurves[14][6]+2,tone_8000_60dB_SL,sizeof(*p->tonecurves[14][6])*EHMER_MAX);
 305   memcpy(p->tonecurves[14][8]+2,tone_8000_80dB_SL,sizeof(*p->tonecurves[14][8])*EHMER_MAX);
 306   memcpy(p->tonecurves[14][10]+2,tone_8000_100dB_SL,sizeof(*p->tonecurves[14][10])*EHMER_MAX);
 307
 308   memcpy(p->tonecurves[16][4]+2,tone_8000_40dB_SL,sizeof(*p->tonecurves[16][4])*EHMER_MAX);
 309   memcpy(p->tonecurves[16][6]+2,tone_8000_60dB_SL,sizeof(*p->tonecurves[16][6])*EHMER_MAX);
 310   memcpy(p->tonecurves[16][8]+2,tone_8000_80dB_SL,sizeof(*p->tonecurves[16][8])*EHMER_MAX);
 311   memcpy(p->tonecurves[16][10]+2,tone_8000_100dB_SL,sizeof(*p->tonecurves[16][10])*EHMER_MAX);
 312
 313   /* value limit the tonal masking curves; the peakatt not only
 314      optionally specifies maximum dynamic depth, but also [always]
 315      limits the masking curves to a minimum depth */
 316   for(i=0;i<P_BANDS;i+=2)
 317     for(j=4;j<P_LEVELS;j+=2)
 318       for(k=2;k<EHMER_MAX+2;k++)
 319         p->tonecurves[i][j][k]+=vi->tone_masteratt;
 320
 321   /* interpolate curves between */
 322   for(i=1;i<P_BANDS;i+=2)
 323     for(j=4;j<P_LEVELS;j+=2){
 324       memcpy(p->tonecurves[i][j]+2,p->tonecurves[i-1][j]+2,EHMER_MAX*sizeof(*p->tonecurves[i][j]));
 325       /*interp_curve(p->tonecurves[i][j],
 326                    p->tonecurves[i-1][j],
 327                    p->tonecurves[i+1][j],.5);*/
 328       min_curve(p->tonecurves[i][j]+2,p->tonecurves[i+1][j]+2);
 329     }
 330
 331   /* set up the final curves */
 332   for(i=0;i<P_BANDS;i++)
 333     setup_curve(p->tonecurves[i],i,vi->toneatt->block[i]);
 334
 335   if(vi->curvelimitp){
 336     /* value limit the tonal masking curves; the peakatt not only
 337        optionally specifies maximum dynamic depth, but also [always]
 338        limits the masking curves to a minimum depth  */
 339     for(i=0;i<P_BANDS;i++)
 340       for(j=0;j<P_LEVELS;j++){
 341         for(k=2;k<EHMER_OFFSET+2+vi->curvelimitp;k++)
 342           if(p->tonecurves[i][j][k]> vi->peakatt->block[i][j])
 343             p->tonecurves[i][j][k]=  vi->peakatt->block[i][j];
 344           else
 345             break;
 346       }
 347   }
 348
 349   if(vi->tone_guard) /* we limit depth only optionally */
 350     for(i=0;i<P_BANDS;i++)
 351       for(j=0;j<P_LEVELS;j++)
 352         if(p->tonecurves[i][j][EHMER_OFFSET+2]< vi->peakatt->block[i][j])
 353           p->tonecurves[i][j][EHMER_OFFSET+2]=  vi->peakatt->block[i][j];
 354
 355   /* but guarding is mandatory */
 356   for(i=0;i<P_BANDS;i++)
 357     for(j=0;j<P_LEVELS;j++)
 358       if(p->tonecurves[i][j][EHMER_OFFSET+2]< vi->tone_guard)
 359           p->tonecurves[i][j][EHMER_OFFSET+2]=  vi->tone_guard;
 360
 361   /* set up rolling noise median */
 362   for(i=0;i<n;i++){
 363     float halfoc=toOC((i+.5)*rate/(2.*n))*2.;
 364     int inthalfoc;
 365     float del;
 366
 367     if(halfoc<0)halfoc=0;
 368     if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1;
 369     inthalfoc=(int)halfoc;
 370     del=halfoc-inthalfoc;
 371     p->noiseoffset[i]=
 372       p->vi->noiseoff[inthalfoc]*(1.-del) +
 373       p->vi->noiseoff[inthalfoc+1]*del;
 374   }
 375
 376   analysis_noisy=1;
 377   _analysis_output("noiseoff",0,p->noiseoffset,n,1,0);
 378   _analysis_output("noisethresh",0,p->noisethresh,n,1,0);
 379
 380   for(i=0;i<P_LEVELS;i++)
 381     _analysis_output("curve_63Hz",i,p->tonecurves[0][i]+2,EHMER_MAX,0,0);
 382   for(i=0;i<P_LEVELS;i++)
 383     _analysis_output("curve_88Hz",i,p->tonecurves[1][i]+2,EHMER_MAX,0,0);
 384   for(i=0;i<P_LEVELS;i++)
 385     _analysis_output("curve_125Hz",i,p->tonecurves[2][i]+2,EHMER_MAX,0,0);
 386   for(i=0;i<P_LEVELS;i++)
 387     _analysis_output("curve_170Hz",i,p->tonecurves[3][i]+2,EHMER_MAX,0,0);
 388   for(i=0;i<P_LEVELS;i++)
 389     _analysis_output("curve_250Hz",i,p->tonecurves[4][i]+2,EHMER_MAX,0,0);
 390   for(i=0;i<P_LEVELS;i++)
 391     _analysis_output("curve_350Hz",i,p->tonecurves[5][i]+2,EHMER_MAX,0,0);
 392   for(i=0;i<P_LEVELS;i++)
 393     _analysis_output("curve_500Hz",i,p->tonecurves[6][i]+2,EHMER_MAX,0,0);
 394   for(i=0;i<P_LEVELS;i++)
 395     _analysis_output("curve_700Hz",i,p->tonecurves[7][i]+2,EHMER_MAX,0,0);
 396   for(i=0;i<P_LEVELS;i++)
 397     _analysis_output("curve_1kHz",i,p->tonecurves[8][i]+2,EHMER_MAX,0,0);
 398   for(i=0;i<P_LEVELS;i++)
 399     _analysis_output("curve_1.4Hz",i,p->tonecurves[9][i]+2,EHMER_MAX,0,0);
 400   for(i=0;i<P_LEVELS;i++)
 401     _analysis_output("curve_2kHz",i,p->tonecurves[10][i]+2,EHMER_MAX,0,0);
 402   for(i=0;i<P_LEVELS;i++)
 403     _analysis_output("curve_2.4kHz",i,p->tonecurves[11][i]+2,EHMER_MAX,0,0);
 404   for(i=0;i<P_LEVELS;i++)
 405     _analysis_output("curve_4kHz",i,p->tonecurves[12][i]+2,EHMER_MAX,0,0);
 406   for(i=0;i<P_LEVELS;i++)
 407     _analysis_output("curve_5.6kHz",i,p->tonecurves[13][i]+2,EHMER_MAX,0,0);
 408   for(i=0;i<P_LEVELS;i++)
 409     _analysis_output("curve_8kHz",i,p->tonecurves[14][i]+2,EHMER_MAX,0,0);
 410   for(i=0;i<P_LEVELS;i++)
 411     _analysis_output("curve_11.5kHz",i,p->tonecurves[15][i]+2,EHMER_MAX,0,0);
 412   for(i=0;i<P_LEVELS;i++)
 413     _analysis_output("curve_16kHz",i,p->tonecurves[16][i]+2,EHMER_MAX,0,0);
 414   analysis_noisy=1;
 415
 416 }
 417
 418 void _vp_psy_clear(vorbis_look_psy *p){
 419   int i,j;
 420   if(p){
 421     if(p->ath)_ogg_free(p->ath);
 422     if(p->octave)_ogg_free(p->octave);
 423     if(p->bark)_ogg_free(p->bark);
 424     if(p->tonecurves){
 425       for(i=0;i<P_BANDS;i++){
 426         for(j=0;j<P_LEVELS;j++){
 427           _ogg_free(p->tonecurves[i][j]);
 428         }
 429         _ogg_free(p->tonecurves[i]);
 430       }
 431       _ogg_free(p->tonecurves);
 432     }
 433     _ogg_free(p->noiseoffset);
 434     _ogg_free(p->noisethresh);
 435     memset(p,0,sizeof(*p));
 436   }
 437 }
 438
 439 /* octave/(8*eighth_octave_lines) x scale and dB y scale */
 440 static void seed_curve(float *seed,
 441                        const float **curves,
 442                        float amp,
 443                        int oc, int n,
 444                        int linesper,float dBoffset){
 445   int i,post1;
 446   int seedptr;
 447   const float *posts,*curve;
 448
 449   int choice=(int)((amp+dBoffset)*.1f);
 450   choice=max(choice,0);
 451   choice=min(choice,P_LEVELS-1);
 452   posts=curves[choice];
 453   curve=posts+2;
 454   post1=(int)posts[1];
 455   seedptr=oc+(posts[0]-16)*linesper-(linesper>>1);
 456
 457   for(i=posts[0];i<post1;i++){
 458     if(seedptr>0){
 459       float lin=amp+curve[i];
 460       if(seed[seedptr]<lin)seed[seedptr]=lin;
 461     }
 462     seedptr+=linesper;
 463     if(seedptr>=n)break;
 464   }
 465 }
 466
 467 static void seed_loop(vorbis_look_psy *p,
 468                       const float ***curves,
 469                       const float *f,
 470                       const float *flr,
 471                       float *seed,
 472                       float specmax){
 473   vorbis_info_psy *vi=p->vi;
 474   long n=p->n,i;
 475   float dBoffset=vi->max_curve_dB-specmax;
 476
 477   /* prime the working vector with peak values */
 478
 479   for(i=0;i<n;i++){
 480     float max=f[i];
 481     long oc=p->octave[i];
 482     while(i+1<n && p->octave[i+1]==oc){
 483       i++;
 484       if(f[i]>max)max=f[i];
 485     }
 486
 487     if(max+6.f>flr[i]){
 488       oc=oc>>p->shiftoc;
 489       if(oc>=P_BANDS)oc=P_BANDS-1;
 490       if(oc<0)oc=0;
 491       seed_curve(seed,
 492                  curves[oc],
 493                  max,
 494                  p->octave[i]-p->firstoc,
 495                  p->total_octave_lines,
 496                  p->eighth_octave_lines,
 497                  dBoffset);
 498     }
 499   }
 500 }
 501
 502 static void seed_chase(float *seeds, int linesper, long n){
 503   long  *posstack=alloca(n*sizeof(*posstack));
 504   float *ampstack=alloca(n*sizeof(*ampstack));
 505   long   stack=0;
 506   long   pos=0;
 507   long   i;
 508
 509   for(i=0;i<n;i++){
 510     if(stack<2){
 511       posstack[stack]=i;
 512       ampstack[stack++]=seeds[i];
 513     }else{
 514       while(1){
 515         if(seeds[i]<ampstack[stack-1]){
 516           posstack[stack]=i;
 517           ampstack[stack++]=seeds[i];
 518           break;
 519         }else{
 520           if(i<posstack[stack-1]+linesper){
 521             if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] &&
 522                i<posstack[stack-2]+linesper){
 523               /* we completely overlap, making stack-1 irrelevant.  pop it */
 524               stack--;
 525               continue;
 526             }
 527           }
 528           posstack[stack]=i;
 529           ampstack[stack++]=seeds[i];
 530           break;
 531
 532         }
 533       }
 534     }
 535   }
 536
 537   /* the stack now contains only the positions that are relevant. Scan
 538      'em straight through */
 539
 540   for(i=0;i<stack;i++){
 541     long endpos;
 542     if(i<stack-1 && ampstack[i+1]>ampstack[i]){
 543       endpos=posstack[i+1];
 544     }else{
 545       endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is
 546                                         discarded in short frames */
 547     }
 548     if(endpos>n)endpos=n;
 549     for(;pos<endpos;pos++)
 550       seeds[pos]=ampstack[i];
 551   }
 552
 553   /* there.  Linear time.  I now remember this was on a problem set I
 554      had in Grad Skool... I didn't solve it at the time ;-) */
 555
 556 }
 557
 558 /* bleaugh, this is more complicated than it needs to be */
 559 static void max_seeds(vorbis_look_psy *p,
 560                       vorbis_look_psy_global *g,
 561                       int channel,
 562                       float *seed,
 563                       float *flr){
 564   long   n=p->total_octave_lines;
 565   int    linesper=p->eighth_octave_lines;
 566   long   linpos=0;
 567   long   pos;
 568
 569   seed_chase(seed,linesper,n); /* for masking */
 570
 571   pos=p->octave[0]-p->firstoc-(linesper>>1);
 572   while(linpos+1<p->n){
 573     float minV=seed[pos];
 574     long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc;
 575     if(minV>p->vi->tone_abs_limit)minV=p->vi->tone_abs_limit;
 576     while(pos+1<=end){
 577       pos++;
 578       if((seed[pos]>NEGINF && seed[pos]<minV) || minV==NEGINF)
 579         minV=seed[pos];
 580     }
 581
 582     /* seed scale is log.  Floor is linear.  Map back to it */
 583     end=pos+p->firstoc;
 584     for(;linpos<p->n && p->octave[linpos]<=end;linpos++)
 585       if(flr[linpos]<minV)flr[linpos]=minV;
 586   }
 587
 588   {
 589     float minV=seed[p->total_octave_lines-1];
 590     for(;linpos<p->n;linpos++)
 591       if(flr[linpos]<minV)flr[linpos]=minV;
 592   }
 593
 594 }
 595
 596 static void bark_noise_hybridmp(int n,const long *b,
 597                                 const float *f,
 598                                 float *noise,
 599                                 const float offset,
 600                                 const int fixed){
 601   long i,hi=b[0]>>16,lo=b[0]>>16,hif=0,lof=0;
 602   double xa=0,xb=0;
 603   double ya=0,yb=0;
 604   double x2a=0,x2b=0;
 605   double y2a=0,y2b=0;
 606   double xya=0,xyb=0;
 607   double na=0,nb=0;
 608
 609   for(i=0;i<n;i++){
 610     if(hi<n){
 611       /* find new lo/hi */
 612       int bi=b[i]&0xffffL;
 613       for(;hi<bi;hi++){
 614         int ii=(hi<0?-hi:hi);
 615         double bin=(f[ii]<-offset?1.:f[ii]+offset);
 616         double nn= bin*bin;
 617         na  += nn;
 618         xa  += hi*nn;
 619         ya  += bin*nn;
 620         x2a += hi*hi*nn;
 621         y2a += bin*bin*nn;
 622         xya += hi*bin*nn;
 623       }
 624       bi=b[i]>>16;
 625       for(;lo<bi;lo++){
 626         int ii=(lo<0?-lo:lo);
 627         double bin=(f[ii]<-offset?1.:f[ii]+offset);
 628         double nn= bin*bin;
 629         na  -= nn;
 630         xa  -= lo*nn;
 631         ya  -= bin*nn;
 632         x2a -= lo*lo*nn;
 633         y2a -= bin*bin*nn;
 634         xya -= lo*bin*nn;
 635       }
 636     }
 637
 638     if(hif<n && fixed>0){
 639       int bi=i+fixed/2;
 640       if(bi>n)bi=n;
 641
 642       for(;hif<bi;hif++){
 643         int ii=(hif<0?-hif:hif);
 644         double bin=(f[ii]<-offset?1.:f[ii]+offset);
 645         double nn= bin*bin;
 646         nb  += nn;
 647         xb  += hif*nn;
 648         yb  += bin*nn;
 649         x2b += hif*hif*nn;
 650         y2b += bin*bin*nn;
 651         xyb += hif*bin*nn;
 652       }
 653       bi=i-(fixed+1)/2;
 654       for(;lof<bi;lof++){
 655         int ii=(lof<0?-lof:lof);
 656         double bin=(f[ii]<-offset?1.:f[ii]+offset);
 657         double nn= bin*bin;
 658         nb  -= nn;
 659         xb  -= lof*nn;
 660         yb  -= bin*nn;
 661         x2b -= lof*lof*nn;
 662         y2b -= bin*bin*nn;
 663         xyb -= lof*bin*nn;
 664       }
 665     }
 666
 667     {
 668       double va=0.f;
 669
 670       if(na>2){
 671         double denom=1./(na*x2a-xa*xa);
 672         double a=(ya*x2a-xya*xa)*denom;
 673         double b=(na*xya-xa*ya)*denom;
 674         va=a+b*i;
 675       }
 676       if(va<0.)va=0.;
 677
 678       if(fixed>0){
 679         double vb=0.f;
 680
 681         if(nb>2){
 682           double denomf=1./(nb*x2b-xb*xb);
 683           double af=(yb*x2b-xyb*xb)*denomf;
 684           double bf=(nb*xyb-xb*yb)*denomf;
 685           vb=af+bf*i;
 686         }
 687         if(vb<0.)vb=0.;
 688         if(va>vb && vb>0.)va=vb;
 689
 690       }
 691
 692       noise[i]=va-offset;
 693     }
 694   }
 695 }
 696
 697
 698 void _vp_remove_floor(vorbis_look_psy *p,
 699                       vorbis_look_psy_global *g,
 700                       float *logmdct,
 701                       float *mdct,
 702                       float *codedflr,
 703                       float *residue,
 704                       float local_specmax){
 705   int i,n=p->n;
 706
 707   for(i=0;i<n;i++)
 708     if(mdct[i]!=0.f)
 709       residue[i]=mdct[i]/codedflr[i];
 710     else
 711       residue[i]=0.f;
 712 }
 713
 714
 715 void _vp_compute_mask(vorbis_look_psy *p,
 716                        vorbis_look_psy_global *g,
 717                        int channel,
 718                        float *logfft,
 719                        float *logmdct,
 720                        float *logmask,
 721                        float global_specmax,
 722                        float local_specmax,
 723                        int lastsize){
 724   int i,n=p->n;
 725   static int seq=0;
 726
 727   float *seed=alloca(sizeof(*seed)*p->total_octave_lines);
 728   for(i=0;i<p->total_octave_lines;i++)seed[i]=NEGINF;
 729
 730   /* noise masking */
 731   if(p->vi->noisemaskp){
 732     float *work=alloca(n*sizeof(*work));
 733
 734     bark_noise_hybridmp(n,p->bark,logmdct,logmask,
 735                         140.,-1);
 736
 737     for(i=0;i<n;i++)work[i]=logmdct[i]-logmask[i];
 738
 739     _analysis_output("medianmdct",seq,work,n,1,0);
 740     bark_noise_hybridmp(n,p->bark,work,logmask,0.,
 741                         p->vi->noisewindowfixed);
 742
 743     for(i=0;i<n;i++)work[i]=logmdct[i]-work[i];
 744
 745     /* work[i] holds the median line (.5), logmask holds the upper
 746        envelope line (1.) */
 747
 748     _analysis_output("median",seq,work,n,1,0);
 749
 750     _analysis_output("medianenvelope",seq,logmask,n,1,0);
 751     for(i=0;i<n;i++)logmask[i]+=work[i];
 752     _analysis_output("envelope",seq,logmask,n,1,0);
 753     for(i=0;i<n;i++)logmask[i]-=work[i];
 754
 755     for(i=0;i<n;i++){
 756       int dB=logmask[i]+.5;
 757       if(dB>=NOISE_COMPAND_LEVELS)dB=NOISE_COMPAND_LEVELS-1;
 758       logmask[i]= work[i]+p->vi->noisecompand[dB]+p->noiseoffset[i];
 759       if(logmask[i]>p->vi->noisemaxsupp)logmask[i]=p->vi->noisemaxsupp;
 760     }
 761
 762     _analysis_output("noise",seq,logmask,n,1,0);
 763
 764   }else{
 765     for(i=0;i<n;i++)logmask[i]=NEGINF;
 766   }
 767
 768   /* set the ATH (floating below localmax, not global max by a
 769      specified att) */
 770   if(p->vi->ath){
 771     float att=local_specmax+p->vi->ath_adjatt;
 772     if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt;
 773
 774     for(i=0;i<n;i++){
 775       float av=p->ath[i]+att;
 776       if(av>logmask[i])logmask[i]=av;
 777     }
 778   }
 779
 780   /* tone/peak masking */
 781   seed_loop(p,(const float ***)p->tonecurves,logfft,logmask,seed,global_specmax);
 782   max_seeds(p,g,channel,seed,logmask);
 783
 784   /* doing this here is clean, but we need to find a faster way to do
 785      it than to just tack it on */
 786
 787   for(i=0;i<n;i++)if(logmdct[i]>=logmask[i])break;
 788   if(i==n)
 789     for(i=0;i<n;i++)logmask[i]=NEGINF;
 790   else
 791     for(i=0;i<n;i++)
 792       logfft[i]=max(logmdct[i],logfft[i]);
 793
 794   seq++;
 795
 796 }
 797
 798 float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){
 799   vorbis_info *vi=vd->vi;
 800   codec_setup_info *ci=vi->codec_setup;
 801   vorbis_info_psy_global *gi=ci->psy_g_param;
 802
 803   int n=ci->blocksizes[vd->W]/2;
 804   float secs=(float)n/vi->rate;
 805
 806   amp+=secs*gi->ampmax_att_per_sec;
 807   if(amp<-9999)amp=-9999;
 808   return(amp);
 809 }
 810
 811 static void couple_lossless(float A, float B,
 812                             float granule,float igranule,
 813                             float *mag, float *ang){
 814
 815   A=rint(A*igranule)*granule;
 816   B=rint(B*igranule)*granule;
 817
 818   if(fabs(A)>fabs(B)){
 819     *mag=A; *ang=(A>0.f?A-B:B-A);
 820   }else{
 821     *mag=B; *ang=(B>0.f?A-B:B-A);
 822   }
 823
 824   if(*ang>fabs(*mag)*1.9999f)*ang=-fabs(*mag)*2.f;
 825
 826 }
 827
 828 static void couple_8phase(float A, float B, float fA, float fB,
 829                          float granule,float igranule,
 830                          float fmag, float *mag, float *ang){
 831
 832   float origmag=FAST_HYPOT(A*fA,B*fB),corr;
 833
 834   if(fmag!=0.f){
 835     float phase=rint((A-B)/fmag);
 836
 837     if(fabs(A)>fabs(B)){
 838       *mag=A;phase=(A>0?phase:-phase);
 839     }else{
 840       *mag=B;phase=(B>0?phase:-phase);
 841     }
 842
 843     switch((int)phase){
 844     case 0:
 845       corr=origmag/FAST_HYPOT(fmag*fA,fmag*fB);
 846       *mag=rint(*mag*corr*igranule)*granule;
 847       *ang=0.f;
 848       break;
 849     case 1:
 850       *mag=rint(*mag*igranule)*granule;
 851       *ang=fabs(*mag);
 852       break;
 853     case -1:
 854       *mag=rint(*mag*igranule)*granule;
 855       *ang= -fabs(*mag);
 856       break;
 857     case -2:
 858       corr=origmag/FAST_HYPOT(fmag*fA,fmag*fB);
 859       *mag=rint(*mag*corr*igranule)*granule;
 860       *ang= -2.f*fabs(*mag);
 861       break;
 862     case 2:
 863       corr=origmag/FAST_HYPOT(fmag*fA,fmag*fB);
 864       *mag= -rint(*mag*corr*igranule)*granule;
 865       *ang= -2.f*fabs(*mag);
 866       break;
 867     }
 868   }else{
 869     *mag=0.f;
 870     *ang=0.f;
 871   }
 872 }
 873
 874 static void couple_6phase(float A, float B, float fA, float fB,
 875                          float granule,float igranule,
 876                          float fmag, float *mag, float *ang){
 877
 878   float origmag=FAST_HYPOT(A*fA,B*fB),corr;
 879
 880   if(fmag!=0.f){
 881     float phase=rint((A-B)/fmag);
 882
 883     if(fabs(A)>fabs(B)){
 884       *mag=A;phase=(A>0?phase:-phase);
 885     }else{
 886       *mag=B;phase=(B>0?phase:-phase);
 887     }
 888
 889     switch((int)phase){
 890     case 0:
 891       corr=origmag/FAST_HYPOT(fmag*fA,fmag*fB);
 892       *mag=rint(*mag*corr*igranule)*granule;
 893       *ang=0.f;
 894       break;
 895     case 1:case 2:
 896       *mag=rint(*mag*igranule)*granule;
 897       *ang= fabs(*mag);
 898       break;
 899     case -1:case -2:
 900       *mag=rint(*mag*igranule)*granule;
 901       *ang= -fabs(*mag);
 902       break;
 903     default:
 904       *mag=0.f;
 905       *ang=0.f;
 906
 907     }
 908   }else{
 909     *mag=0.f;
 910     *ang=0.f;
 911   }
 912 }
 913
 914 static void couple_4phase(float A, float B, float fA, float fB,
 915                          float granule,float igranule,
 916                          float fmag, float *mag, float *ang){
 917
 918   float origmag=FAST_HYPOT(A*fA,B*fB),corr;
 919
 920   if(fmag!=0.f){
 921     float phase=rint((A-B)*.5/fmag);
 922
 923     if(fabs(A)>fabs(B)){
 924       *mag=A;phase=(A>0?phase:-phase);
 925     }else{
 926       *mag=B;phase=(B>0?phase:-phase);
 927     }
 928
 929     corr=origmag/FAST_HYPOT(fmag*fA,fmag*fB);
 930     *mag=rint(*mag*corr*igranule)*granule;
 931     switch((int)phase){
 932     case 0:
 933       *ang=0.f;
 934       break;
 935     case -1:
 936       *ang=-2.f*fabs(*mag);
 937       break;
 938     default:
 939       *mag=-*mag;
 940       *ang=-2.f*fabs(*mag);
 941       break;
 942     }
 943   }else{
 944     *mag=0.f;
 945     *ang=0.f;
 946   }
 947 }
 948
 949 static void couple_point(float A, float B, float fA, float fB,
 950                          float granule,float igranule,
 951                          float fmag, float *mag, float *ang){
 952
 953   float origmag=FAST_HYPOT(A*fA,B*fB),corr;
 954
 955   if(fmag!=0.f){
 956     //float phase=rint((A-B)*.5/fmag);
 957
 958     if(fabs(A)>fabs(B)){
 959       *mag=A;//phase=(A>0?phase:-phase);
 960     }else{
 961       *mag=B;//phase=(B>0?phase:-phase);
 962     }
 963
 964     //switch((int)phase){
 965       //case 0:
 966       corr=origmag/FAST_HYPOT(fmag*fA,fmag*fB);
 967       *mag=rint(*mag*corr*igranule)*granule;
 968       *ang=0.f;
 969       //break;
 970       //default:
 971       //*mag=0.f;
 972       //*ang=0.f;
 973       //break;
 974       //}
 975   }else{
 976     *mag=0.f;
 977     *ang=0.f;
 978   }
 979 }
 980
 981
 982 void _vp_quantize_couple(vorbis_look_psy *p,
 983                          vorbis_info_mapping0 *vi,
 984                          float **pcm,
 985                          float **sofar,
 986                          float **quantized,
 987                          int   *nonzero,
 988                          int   passno){
 989
 990   int i,j,k,n=p->n;
 991   vorbis_info_psy *info=p->vi;
 992
 993   /* perform any requested channel coupling */
 994   for(i=0;i<vi->coupling_steps;i++){
 995     float granulem=info->couple_pass[passno].granulem;
 996     float igranulem=info->couple_pass[passno].igranulem;
 997
 998     /* make sure coupling a zero and a nonzero channel results in two
 999        nonzero channels. */
1000     if(nonzero[vi->coupling_mag[i]] ||
1001        nonzero[vi->coupling_ang[i]]){
1002
1003       float *pcmM=pcm[vi->coupling_mag[i]];
1004       float *pcmA=pcm[vi->coupling_ang[i]];
1005       float *floorM=pcm[vi->coupling_mag[i]]+n;
1006       float *floorA=pcm[vi->coupling_ang[i]]+n;
1007       float *sofarM=sofar[vi->coupling_mag[i]];
1008       float *sofarA=sofar[vi->coupling_ang[i]];
1009       float *qM=quantized[vi->coupling_mag[i]];
1010       float *qA=quantized[vi->coupling_ang[i]];
1011
1012       nonzero[vi->coupling_mag[i]]=1;
1013       nonzero[vi->coupling_ang[i]]=1;
1014
1015       for(j=0,k=0;j<n;k++){
1016         vp_couple *part=info->couple_pass[passno].couple_pass+k;
1017
1018         for(;j<part->limit && j<p->n;j++){
1019           /* partition by partition; k is our by-location partition
1020              class counter */
1021           float ang,mag,fmag=max(fabs(pcmM[j]),fabs(pcmA[j]));
1022
1023           if(fmag<part->amppost_point){
1024             couple_point(pcmM[j],pcmA[j],floorM[j],floorA[j],
1025                          granulem,igranulem,fmag,&mag,&ang);
1026           }else{
1027             if(fmag<part->amppost_6phase){
1028               couple_6phase(pcmM[j],pcmA[j],floorM[j],floorA[j],
1029                            granulem,igranulem,fmag,&mag,&ang);
1030             }else{
1031               if(fmag<part->amppost_8phase){
1032                 couple_8phase(pcmM[j],pcmA[j],floorM[j],floorA[j],
1033                               granulem,igranulem,fmag,&mag,&ang);
1034               }else{
1035                 couple_lossless(pcmM[j],pcmA[j],
1036                                 granulem,igranulem,&mag,&ang);
1037               }
1038             }
1039           }
1040
1041           qM[j]=mag-sofarM[j];
1042           qA[j]=ang-sofarA[j];
1043         }
1044       }
1045     }
1046   }
1047 }