lib/psy.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE.  *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
   5  * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE.    *
   6  * PLEASE READ THESE TERMS DISTRIBUTING.                            *
   7  *                                                                  *
   8  * THE OggSQUISH 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.26 2000/08/23 06:38:49 xiphmont Exp $
  16
  17  ********************************************************************/
  18
  19 #include <stdlib.h>
  20 #include <math.h>
  21 #include <string.h>
  22 #include "vorbis/codec.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 /* Why Bark scale for encoding but not masking computation? Because
  33    masking has a strong harmonic dependancy */
  34
  35 /* the beginnings of real psychoacoustic infrastructure.  This is
  36    still not tightly tuned */
  37 void _vi_psy_free(vorbis_info_psy *i){
  38   if(i){
  39     memset(i,0,sizeof(vorbis_info_psy));
  40     free(i);
  41   }
  42 }
  43
  44 /* Set up decibel threshhold slopes on a Bark frequency scale */
  45 /* ATH is the only bit left on a Bark scale.  No reason to change it
  46    right now */
  47 static void set_curve(double *ref,double *c,int n, double crate){
  48   int i,j=0;
  49
  50   for(i=0;i<MAX_BARK-1;i++){
  51     int endpos=rint(fromBARK(i+1)*2*n/crate);
  52     double base=ref[i];
  53     if(j<endpos){
  54       double delta=(ref[i+1]-base)/(endpos-j);
  55       for(;j<endpos && j<n;j++){
  56         c[j]=base;
  57         base+=delta;
  58       }
  59     }
  60   }
  61 }
  62
  63 static void min_curve(double *c,
  64                        double *c2){
  65   int i;
  66   for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i];
  67 }
  68 static void max_curve(double *c,
  69                        double *c2){
  70   int i;
  71   for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i];
  72 }
  73
  74 static void attenuate_curve(double *c,double att){
  75   int i;
  76   for(i=0;i<EHMER_MAX;i++)
  77     c[i]+=att;
  78 }
  79
  80 static void linear_curve(double *c){
  81   int i;
  82   for(i=0;i<EHMER_MAX;i++)
  83     if(c[i]<=-200.)
  84       c[i]=0.;
  85     else
  86       c[i]=fromdB(c[i]);
  87 }
  88
  89 static void interp_curve(double *c,double *c1,double *c2,double del){
  90   int i;
  91   for(i=0;i<EHMER_MAX;i++)
  92     c[i]=c2[i]*del+c1[i]*(1.-del);
  93 }
  94
  95 static void setup_curve(double **c,
  96                         int band,
  97                         double *curveatt_dB){
  98   int i,j;
  99   double ath[EHMER_MAX];
 100   double tempc[P_LEVELS][EHMER_MAX];
 101
 102   memcpy(c[0],c[4],sizeof(double)*EHMER_MAX);
 103   memcpy(c[2],c[4],sizeof(double)*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     double oc_min=band*.5-1+(i-EHMER_OFFSET)*.125;
 113     double oc_max=band*.5-1+(i-EHMER_OFFSET+1)*.125;
 114     double bark=toBARK(fromOC(oc_min));
 115     int ibark=floor(bark);
 116     double del=bark-ibark;
 117     double ath_min,ath_max;
 118
 119     if(ibark<26)
 120       ath_min=ATH_Bark_dB[ibark]*(1.-del)+ATH_Bark_dB[ibark+1]*del;
 121     else
 122       ath_min=200.;
 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.-del)+ATH_Bark_dB[ibark+1]*del;
 130     else
 131       ath_max=200.;
 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],c[i-1],c[i+1],.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],curveatt_dB[i]);
 146     memcpy(tempc[i],ath,EHMER_MAX*sizeof(double));
 147     attenuate_curve(tempc[i],-i*10.);
 148     max_curve(tempc[i],c[i]);
 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(i=P_LEVELS-1;i>0;i--){
 162     for(j=0;j<i;j++)
 163       min_curve(c[i],tempc[j]);
 164   }
 165
 166   /* take things out of dB domain into linear amplitude */
 167   for(i=0;i<P_LEVELS;i++)
 168     linear_curve(c[i]);
 169
 170 }
 171
 172 void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi,int n,long rate){
 173   long i,j;
 174   memset(p,0,sizeof(vorbis_look_psy));
 175   p->ath=malloc(n*sizeof(double));
 176   p->octave=malloc(n*sizeof(int));
 177   p->bark=malloc(n*sizeof(double));
 178   p->vi=vi;
 179   p->n=n;
 180
 181   /* set up the lookups for a given blocksize and sample rate */
 182   /* Vorbis max sample rate is limited by 26 Bark (54kHz) */
 183   set_curve(ATH_Bark_dB, p->ath,n,rate);
 184   for(i=0;i<n;i++)
 185     p->ath[i]=fromdB(p->ath[i]);
 186   for(i=0;i<n;i++)
 187     p->bark[i]=toBARK(rate/(2*n)*i);
 188
 189   for(i=0;i<n;i++){
 190     int oc=toOC((i+.5)*rate/(2*n))*2.+2; /* half octaves, actually */
 191     if(oc<0)oc=0;
 192     if(oc>=P_BANDS)oc=P_BANDS-1;
 193     p->octave[i]=oc;
 194   }
 195
 196   p->tonecurves=malloc(P_BANDS*sizeof(double **));
 197   p->noiseatt=malloc(P_BANDS*sizeof(double **));
 198   p->peakatt=malloc(P_BANDS*sizeof(double *));
 199   for(i=0;i<P_BANDS;i++){
 200     p->tonecurves[i]=malloc(P_LEVELS*sizeof(double *));
 201     p->noiseatt[i]=malloc(P_LEVELS*sizeof(double));
 202     p->peakatt[i]=malloc(P_LEVELS*sizeof(double));
 203   }
 204
 205   for(i=0;i<P_BANDS;i++)
 206     for(j=0;j<P_LEVELS;j++){
 207       p->tonecurves[i][j]=malloc(EHMER_MAX*sizeof(double));
 208     }
 209
 210   /* OK, yeah, this was a silly way to do it */
 211   memcpy(p->tonecurves[0][4],tone_125_40dB_SL,sizeof(double)*EHMER_MAX);
 212   memcpy(p->tonecurves[0][6],tone_125_60dB_SL,sizeof(double)*EHMER_MAX);
 213   memcpy(p->tonecurves[0][8],tone_125_80dB_SL,sizeof(double)*EHMER_MAX);
 214   memcpy(p->tonecurves[0][10],tone_125_100dB_SL,sizeof(double)*EHMER_MAX);
 215
 216   memcpy(p->tonecurves[2][4],tone_125_40dB_SL,sizeof(double)*EHMER_MAX);
 217   memcpy(p->tonecurves[2][6],tone_125_60dB_SL,sizeof(double)*EHMER_MAX);
 218   memcpy(p->tonecurves[2][8],tone_125_80dB_SL,sizeof(double)*EHMER_MAX);
 219   memcpy(p->tonecurves[2][10],tone_125_100dB_SL,sizeof(double)*EHMER_MAX);
 220
 221   memcpy(p->tonecurves[4][4],tone_250_40dB_SL,sizeof(double)*EHMER_MAX);
 222   memcpy(p->tonecurves[4][6],tone_250_60dB_SL,sizeof(double)*EHMER_MAX);
 223   memcpy(p->tonecurves[4][8],tone_250_80dB_SL,sizeof(double)*EHMER_MAX);
 224   memcpy(p->tonecurves[4][10],tone_250_100dB_SL,sizeof(double)*EHMER_MAX);
 225
 226   memcpy(p->tonecurves[6][4],tone_500_40dB_SL,sizeof(double)*EHMER_MAX);
 227   memcpy(p->tonecurves[6][6],tone_500_60dB_SL,sizeof(double)*EHMER_MAX);
 228   memcpy(p->tonecurves[6][8],tone_500_80dB_SL,sizeof(double)*EHMER_MAX);
 229   memcpy(p->tonecurves[6][10],tone_500_100dB_SL,sizeof(double)*EHMER_MAX);
 230
 231   memcpy(p->tonecurves[8][4],tone_1000_40dB_SL,sizeof(double)*EHMER_MAX);
 232   memcpy(p->tonecurves[8][6],tone_1000_60dB_SL,sizeof(double)*EHMER_MAX);
 233   memcpy(p->tonecurves[8][8],tone_1000_80dB_SL,sizeof(double)*EHMER_MAX);
 234   memcpy(p->tonecurves[8][10],tone_1000_100dB_SL,sizeof(double)*EHMER_MAX);
 235
 236   memcpy(p->tonecurves[10][4],tone_2000_40dB_SL,sizeof(double)*EHMER_MAX);
 237   memcpy(p->tonecurves[10][6],tone_2000_60dB_SL,sizeof(double)*EHMER_MAX);
 238   memcpy(p->tonecurves[10][8],tone_2000_80dB_SL,sizeof(double)*EHMER_MAX);
 239   memcpy(p->tonecurves[10][10],tone_2000_100dB_SL,sizeof(double)*EHMER_MAX);
 240
 241   memcpy(p->tonecurves[12][4],tone_4000_40dB_SL,sizeof(double)*EHMER_MAX);
 242   memcpy(p->tonecurves[12][6],tone_4000_60dB_SL,sizeof(double)*EHMER_MAX);
 243   memcpy(p->tonecurves[12][8],tone_4000_80dB_SL,sizeof(double)*EHMER_MAX);
 244   memcpy(p->tonecurves[12][10],tone_4000_100dB_SL,sizeof(double)*EHMER_MAX);
 245
 246   memcpy(p->tonecurves[14][4],tone_8000_40dB_SL,sizeof(double)*EHMER_MAX);
 247   memcpy(p->tonecurves[14][6],tone_8000_60dB_SL,sizeof(double)*EHMER_MAX);
 248   memcpy(p->tonecurves[14][8],tone_8000_80dB_SL,sizeof(double)*EHMER_MAX);
 249   memcpy(p->tonecurves[14][10],tone_8000_100dB_SL,sizeof(double)*EHMER_MAX);
 250
 251   memcpy(p->tonecurves[16][4],tone_8000_40dB_SL,sizeof(double)*EHMER_MAX);
 252   memcpy(p->tonecurves[16][6],tone_8000_60dB_SL,sizeof(double)*EHMER_MAX);
 253   memcpy(p->tonecurves[16][8],tone_8000_80dB_SL,sizeof(double)*EHMER_MAX);
 254   memcpy(p->tonecurves[16][10],tone_8000_100dB_SL,sizeof(double)*EHMER_MAX);
 255
 256   /* interpolate curves between */
 257   for(i=1;i<P_BANDS;i+=2)
 258     for(j=4;j<P_LEVELS;j+=2){
 259       memcpy(p->tonecurves[i][j],p->tonecurves[i-1][j],EHMER_MAX*sizeof(double));
 260       /*interp_curve(p->tonecurves[i][j],
 261                    p->tonecurves[i-1][j],
 262                    p->tonecurves[i+1][j],.5);*/
 263       min_curve(p->tonecurves[i][j],p->tonecurves[i+1][j]);
 264       /*min_curve(p->tonecurves[i][j],p->tonecurves[i-1][j]);*/
 265     }
 266
 267   /*for(i=0;i<P_BANDS-1;i++)
 268     for(j=4;j<P_LEVELS;j+=2)
 269     min_curve(p->tonecurves[i][j],p->tonecurves[i+1][j]);*/
 270
 271   /* set up the final curves */
 272   for(i=0;i<P_BANDS;i++)
 273     setup_curve(p->tonecurves[i],i,vi->toneatt[i]);
 274
 275   /* set up attenuation levels */
 276   for(i=0;i<P_BANDS;i++)
 277     for(j=0;j<P_LEVELS;j++){
 278       p->peakatt[i][j]=fromdB(p->vi->peakatt[i][j]);
 279       p->noiseatt[i][j]=fromdB(p->vi->noiseatt[i][j]);
 280     }
 281
 282 }
 283
 284 void _vp_psy_clear(vorbis_look_psy *p){
 285   int i,j;
 286   if(p){
 287     if(p->ath)free(p->ath);
 288     if(p->octave)free(p->octave);
 289     if(p->tonecurves){
 290       for(i=0;i<P_BANDS;i++){
 291         for(j=0;j<P_LEVELS;j++){
 292           free(p->tonecurves[i][j]);
 293         }
 294         free(p->noiseatt[i]);
 295         free(p->tonecurves[i]);
 296         free(p->peakatt[i]);
 297       }
 298       free(p->tonecurves);
 299       free(p->noiseatt);
 300       free(p->peakatt);
 301     }
 302     memset(p,0,sizeof(vorbis_look_psy));
 303   }
 304 }
 305
 306 static void compute_decay_fixed(vorbis_look_psy *p,double *f, double *decay, int n){
 307   /* handle decay */
 308   int i;
 309   double decscale=fromdB(p->vi->decay_coeff*n);
 310   double attscale=1./fromdB(p->vi->attack_coeff);
 311
 312   static int frameno=0;
 313
 314   for(i=10;i<n;i++){
 315     double pre=decay[i];
 316     if(decay[i]){
 317       double val=decay[i]*decscale;
 318       double att=fabs(f[i]/val);
 319
 320       if(att>attscale)
 321         decay[i]=fabs(f[i]/attscale);
 322       else
 323         decay[i]=val;
 324     }else{
 325       decay[i]=fabs(f[i]/attscale);
 326     }
 327     if(pre>f[i])f[i]=pre;
 328   }
 329 }
 330
 331 static long _eights[EHMER_MAX+1]={
 332   981,1069,1166,1272,
 333   1387,1512,1649,1798,
 334   1961,2139,2332,2543,
 335   2774,3025,3298,3597,
 336   3922,4277,4664,5087,
 337   5547,6049,6597,7194,
 338   7845,8555,9329,10173,
 339   11094,12098,13193,14387,
 340   15689,17109,18658,20347,
 341   22188,24196,26386,28774,
 342   31379,34219,37316,40693,
 343   44376,48393,52772,57549,
 344   62757,68437,74631,81386,
 345   88752,96785,105545,115097,
 346   125515};
 347
 348 static int seed_curve(double *flr,
 349                       double **curves,
 350                        double amp,double specmax,
 351                        int x,int n,double specatt,
 352                        int maxEH){
 353   int i;
 354   double *curve;
 355
 356   /* make this attenuation adjustable */
 357   int choice=(int)((todB(amp)-specmax+specatt)/10.+.5);
 358   choice=max(choice,0);
 359   choice=min(choice,P_LEVELS-1);
 360
 361   for(i=maxEH;i>=0;i--)
 362     if(((x*_eights[i])>>12)<n)break;
 363   maxEH=i;
 364   curve=curves[choice];
 365
 366   for(;i>=0;i--)
 367     if(curve[i]>0.)break;
 368
 369   for(;i>=0;i--){
 370     double lin=curve[i];
 371     if(lin>0.){
 372       double *fp=flr+((x*_eights[i])>>12);
 373       lin*=amp;
 374       if(*fp<lin)*fp=lin;
 375     }else break;
 376   }
 377   return(maxEH);
 378 }
 379
 380 static void seed_peak(double *flr,
 381                       double *att,
 382                       double amp,double specmax,
 383                       int x,int n,double specatt){
 384   int prevx=(x*_eights[16])>>12;
 385
 386   /* make this attenuation adjustable */
 387   int choice=rint((todB(amp)-specmax+specatt)/10.+.5);
 388   if(choice<0)choice=0;
 389   if(choice>=P_LEVELS)choice=P_LEVELS-1;
 390
 391   if(prevx<n){
 392     double lin=att[choice];
 393     if(lin){
 394       lin*=amp;
 395       if(flr[prevx]<lin)flr[prevx]=lin;
 396     }
 397   }
 398 }
 399
 400 static void seed_generic(vorbis_look_psy *p,
 401                          double ***curves,
 402                          double *f,
 403                          double *flr,
 404                          double *seeds,
 405                          double specmax){
 406   vorbis_info_psy *vi=p->vi;
 407   long n=p->n,i;
 408   int maxEH=EHMER_MAX-1;
 409
 410   /* prime the working vector with peak values */
 411   /* Use the 125 Hz curve up to 125 Hz and 8kHz curve after 8kHz. */
 412   for(i=0;i<n;i++)
 413     if(f[i]>flr[i])
 414       maxEH=seed_curve(seeds,curves[p->octave[i]],
 415                        f[i],specmax,i,n,vi->max_curve_dB,maxEH);
 416 }
 417
 418 static void seed_att(vorbis_look_psy *p,
 419                      double **att,
 420                      double *f,
 421                      double *flr,
 422                      double specmax){
 423   vorbis_info_psy *vi=p->vi;
 424   long n=p->n,i;
 425
 426   for(i=0;i<n;i++)
 427     if(f[i]>flr[i])
 428       seed_peak(flr,att[p->octave[i]],f[i],
 429                 specmax,i,n,vi->max_curve_dB);
 430 }
 431
 432 static void seed_point(vorbis_look_psy *p,
 433                      double **att,
 434                      double *f,
 435                      double *flr,
 436                      double specmax){
 437   vorbis_info_psy *vi=p->vi;
 438   long n=p->n,i;
 439
 440   for(i=0;i<n;i++){
 441     /* make this attenuation adjustable */
 442     int choice=rint((todB(f[i])-specmax+vi->max_curve_dB)/10.+.5);
 443     double lin;
 444     if(choice<0)choice=0;
 445     if(choice>=P_LEVELS)choice=P_LEVELS-1;
 446     lin=att[p->octave[i]][choice]*f[i];
 447     if(flr[i]<lin)flr[i]=lin;
 448   }
 449 }
 450
 451 /* bleaugh, this is more complicated than it needs to be */
 452 static void max_seeds(vorbis_look_psy *p,double *seeds,double *flr){
 453   long n=p->n,i,j;
 454   long *posstack=alloca(n*sizeof(long));
 455   double *ampstack=alloca(n*sizeof(double));
 456   long stack=0;
 457
 458   for(i=0;i<n;i++){
 459     if(stack<2){
 460       posstack[stack]=i;
 461       ampstack[stack++]=seeds[i];
 462     }else{
 463       while(1){
 464         if(seeds[i]<ampstack[stack-1]){
 465           posstack[stack]=i;
 466           ampstack[stack++]=seeds[i];
 467           break;
 468         }else{
 469           if(i<posstack[stack-1]*1.0905077080){
 470             if(stack>1 && ampstack[stack-1]<ampstack[stack-2] &&
 471                i<posstack[stack-2]*1.0905077080){
 472               /* we completely overlap, making stack-1 irrelevant.  pop it */
 473               stack--;
 474               continue;
 475             }
 476           }
 477           posstack[stack]=i;
 478           ampstack[stack++]=seeds[i];
 479           break;
 480
 481         }
 482       }
 483     }
 484   }
 485
 486   /* the stack now contains only the positions that are relevant. Scan
 487      'em straight through */
 488   {
 489     long pos=0;
 490     for(i=0;i<stack;i++){
 491       long endpos;
 492       if(i<stack-1 && ampstack[i+1]>ampstack[i]){
 493         endpos=posstack[i+1];
 494       }else{
 495         endpos=posstack[i]*1.0905077080+1; /* +1 is important, else bin 0 is
 496                                        discarded in short frames */
 497       }
 498       if(endpos>n)endpos=n;
 499       for(j=pos;j<endpos;j++)
 500         if(flr[j]<ampstack[i])
 501           flr[j]=ampstack[i];
 502       pos=endpos;
 503     }
 504   }
 505
 506   /* there.  Linear time.  I now remember this was on a problem set I
 507      had in Grad Skool... I didn't solve it at the time ;-) */
 508 }
 509
 510 static void bark_noise(long n,double *b,double *f,double *noise){
 511   long i=1,lo=0,hi=2;
 512   double acc=0.,val,del=0.;
 513
 514   double *norm=alloca(n*sizeof(double));
 515   double normacc=0;
 516
 517   memset(noise,0,n*sizeof(double));
 518   memset(norm,0,n*sizeof(double));
 519
 520   while(hi<n){
 521     val=todB(f[i]*f[i])+400.;
 522     del=1./(i-lo);
 523     noise[lo]+=val*del;
 524     noise[i]-=val*del;
 525     norm[lo]+=del;
 526     norm[i]-=del;
 527
 528     del=1./(hi-i);
 529     noise[i]-=val*del;
 530     noise[hi]+=val*del;
 531     norm[hi]+=del;
 532     norm[i]-=del;
 533
 534
 535     i++;
 536     for(;hi<n && b[hi]-.3<b[i];hi++);
 537     for(;lo<i-1 && b[lo]+.3<b[i];lo++);
 538     if(i==hi)hi++;
 539   }
 540
 541   {
 542     long ilo=i-lo;
 543     long hii=hi-i;
 544     long hilo=hi-lo;
 545
 546     for(;i<n;i++){
 547       val=todB(f[i]*f[i])+400.;
 548       del=1./(hii);
 549       noise[i]-=val*del;
 550       norm[i]-=del;
 551
 552       del=1./(ilo);
 553       noise[i-ilo]+=val*del;
 554       noise[i]-=val*del;
 555       norm[i-ilo]+=del;
 556       norm[i]-=del;
 557     }
 558     for(i=1,lo=n-ilo;lo<n;lo++,i++){
 559       val=todB(f[n-i]*f[n-i])+400.;
 560       del=1./ilo;
 561       noise[lo]+=val*del;
 562       norm[lo]+=del;
 563     }
 564   }
 565
 566
 567   acc=0;
 568   val=0;
 569
 570   for(i=0;i<n;i++){
 571     val+=norm[i];
 572     norm[i]=val;
 573     acc+=noise[i];
 574     noise[i]=acc;
 575   }
 576
 577   val=0;
 578   acc=0;
 579   for(i=0;i<n;i++){
 580     val+=norm[i];
 581     acc+=noise[i];
 582     if(val==0){
 583       noise[i]=0.;
 584       norm[i]=0;
 585     }else{
 586       double v=acc/val-400;
 587       noise[i]=sqrt(fromdB(v));
 588     }
 589   }
 590 }
 591
 592 /* stability doesn't matter */
 593 static int comp(const void *a,const void *b){
 594   if(fabs(**(double **)a)<fabs(**(double **)b))
 595     return(1);
 596   else
 597     return(-1);
 598 }
 599
 600 static int frameno=0;
 601 void _vp_compute_mask(vorbis_look_psy *p,double *f,
 602                       double *flr,
 603                       double *decay){
 604   double *smooth=alloca(sizeof(double)*p->n);
 605   int i,n=p->n;
 606   double specmax=0.;
 607
 608   double *seed=alloca(sizeof(double)*p->n);
 609   double *seed2=alloca(sizeof(double)*p->n);
 610
 611   memset(flr,0,n*sizeof(double));
 612
 613   /* noise masking */
 614   if(p->vi->noisemaskp){
 615     memset(seed,0,n*sizeof(double));
 616     bark_noise(n,p->bark,f,seed);
 617     seed_point(p,p->noiseatt,seed,flr,specmax);
 618
 619   }
 620
 621   /* smooth the data is that's called for ********************************/
 622   for(i=0;i<n;i++)smooth[i]=fabs(f[i]);
 623   if(p->vi->smoothp){
 624     /* compute power^.5 of three neighboring bins to smooth for peaks
 625        that get split twixt bins/peaks that nail the bin.  This evens
 626        out treatment as we're not doing additive masking any longer. */
 627     double acc=smooth[0]*smooth[0]+smooth[1]*smooth[1];
 628     double prev=smooth[0];
 629
 630     smooth[0]=sqrt(acc);
 631     for(i=1;i<n-1;i++){
 632       double this=smooth[i];
 633       acc+=smooth[i+1]*smooth[i+1];
 634       if(acc<0)acc=0; /* it can happen due to finite precision */
 635       smooth[i]=sqrt(acc);
 636       acc-=prev*prev;
 637       prev=this;
 638     }
 639     smooth[n-1]=sqrt(acc);
 640   }
 641
 642   /* find the highest peak so we know the limits *************************/
 643   for(i=0;i<n;i++){
 644     if(smooth[i]>specmax)specmax=smooth[i];
 645   }
 646   specmax=todB(specmax);
 647
 648   /* set the ATH (floating below specmax by a specified att) */
 649   if(p->vi->athp){
 650     double att=specmax+p->vi->ath_adjatt;
 651     if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt;
 652     att=fromdB(att);
 653
 654     for(i=0;i<n;i++){
 655       double av=p->ath[i]*att;
 656       if(av>flr[i])flr[i]=av;
 657     }
 658   }
 659
 660   /* peak attenuation ******/
 661   if(p->vi->peakattp){
 662     memset(seed,0,n*sizeof(double));
 663     seed_att(p,p->peakatt,smooth,seed,specmax);
 664     max_seeds(p,seed,flr);
 665   }
 666
 667   /* tone masking */
 668   if(p->vi->tonemaskp){
 669     memset(seed,0,n*sizeof(double));
 670     memset(seed2,0,n*sizeof(double));
 671
 672     seed_generic(p,p->tonecurves,smooth,flr,seed2,specmax);
 673     max_seeds(p,seed2,seed2);
 674
 675     for(i=0;i<n;i++)if(seed2[i]<flr[i])seed2[i]=flr[i];
 676     for(i=0;i<n;i++)if(seed2[i]<decay[i])seed2[i]=decay[i];
 677
 678     seed_generic(p,p->tonecurves,smooth,seed2,seed,specmax);
 679     max_seeds(p,seed,seed);
 680
 681     if(p->vi->decayp)
 682       compute_decay_fixed(p,seed,decay,n);
 683
 684     for(i=0;i<n;i++)if(flr[i]<seed[i])flr[i]=seed[i];
 685
 686   }
 687
 688   frameno++;
 689 }
 690
 691
 692 /* this applies the floor and (optionally) tries to preserve noise
 693    energy in low resolution portions of the spectrum */
 694 /* f and flr are *linear* scale, not dB */
 695 void _vp_apply_floor(vorbis_look_psy *p,double *f, double *flr){
 696   double *work=alloca(p->n*sizeof(double));
 697   int i,j,addcount=0;
 698
 699   /* subtract the floor */
 700   for(j=0;j<p->n;j++){
 701     if(flr[j]<=0)
 702       work[j]=0.;
 703     else
 704       work[j]=f[j]/flr[j];
 705   }
 706
 707   memcpy(f,work,p->n*sizeof(double));
 708 }
 709
 710