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-2007             *
   9  * by the Xiph.Org Foundation http://www.xiph.org/                  *
  10  *                                                                  *
  11  ********************************************************************
  12
  13  function: psychoacoustics not including preecho
  14  last mod: $Id$
  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 static double stereo_threshholds[]={0.0, .5, 1.0, 1.5, 2.5, 4.5, 8.5, 16.5, 9e10};
  34 static double stereo_threshholds_limited[]={0.0, .5, 1.0, 1.5, 2.0, 2.5, 4.5, 8.5, 9e10};
  35
  36 vorbis_look_psy_global *_vp_global_look(vorbis_info *vi){
  37   codec_setup_info *ci=vi->codec_setup;
  38   vorbis_info_psy_global *gi=&ci->psy_g_param;
  39   vorbis_look_psy_global *look=_ogg_calloc(1,sizeof(*look));
  40
  41   look->channels=vi->channels;
  42
  43   look->ampmax=-9999.;
  44   look->gi=gi;
  45   return(look);
  46 }
  47
  48 void _vp_global_free(vorbis_look_psy_global *look){
  49   if(look){
  50     memset(look,0,sizeof(*look));
  51     _ogg_free(look);
  52   }
  53 }
  54
  55 void _vi_gpsy_free(vorbis_info_psy_global *i){
  56   if(i){
  57     memset(i,0,sizeof(*i));
  58     _ogg_free(i);
  59   }
  60 }
  61
  62 void _vi_psy_free(vorbis_info_psy *i){
  63   if(i){
  64     memset(i,0,sizeof(*i));
  65     _ogg_free(i);
  66   }
  67 }
  68
  69 static void min_curve(float *c,
  70                        float *c2){
  71   int i;
  72   for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i];
  73 }
  74 static void max_curve(float *c,
  75                        float *c2){
  76   int i;
  77   for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i];
  78 }
  79
  80 static void attenuate_curve(float *c,float att){
  81   int i;
  82   for(i=0;i<EHMER_MAX;i++)
  83     c[i]+=att;
  84 }
  85
  86 static float ***setup_tone_curves(float curveatt_dB[P_BANDS],float binHz,int n,
  87                                   float center_boost, float center_decay_rate){
  88   int i,j,k,m;
  89   float ath[EHMER_MAX];
  90   float workc[P_BANDS][P_LEVELS][EHMER_MAX];
  91   float athc[P_LEVELS][EHMER_MAX];
  92   float *brute_buffer=alloca(n*sizeof(*brute_buffer));
  93
  94   float ***ret=_ogg_malloc(sizeof(*ret)*P_BANDS);
  95
  96   memset(workc,0,sizeof(workc));
  97
  98   for(i=0;i<P_BANDS;i++){
  99     /* we add back in the ATH to avoid low level curves falling off to
 100        -infinity and unnecessarily cutting off high level curves in the
 101        curve limiting (last step). */
 102
 103     /* A half-band's settings must be valid over the whole band, and
 104        it's better to mask too little than too much */
 105     int ath_offset=i*4;
 106     for(j=0;j<EHMER_MAX;j++){
 107       float min=999.;
 108       for(k=0;k<4;k++)
 109         if(j+k+ath_offset<MAX_ATH){
 110           if(min>ATH[j+k+ath_offset])min=ATH[j+k+ath_offset];
 111         }else{
 112           if(min>ATH[MAX_ATH-1])min=ATH[MAX_ATH-1];
 113         }
 114       ath[j]=min;
 115     }
 116
 117     /* copy curves into working space, replicate the 50dB curve to 30
 118        and 40, replicate the 100dB curve to 110 */
 119     for(j=0;j<6;j++)
 120       memcpy(workc[i][j+2],tonemasks[i][j],EHMER_MAX*sizeof(*tonemasks[i][j]));
 121     memcpy(workc[i][0],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0]));
 122     memcpy(workc[i][1],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0]));
 123
 124     /* apply centered curve boost/decay */
 125     for(j=0;j<P_LEVELS;j++){
 126       for(k=0;k<EHMER_MAX;k++){
 127         float adj=center_boost+abs(EHMER_OFFSET-k)*center_decay_rate;
 128         if(adj<0. && center_boost>0)adj=0.;
 129         if(adj>0. && center_boost<0)adj=0.;
 130         workc[i][j][k]+=adj;
 131       }
 132     }
 133
 134     /* normalize curves so the driving amplitude is 0dB */
 135     /* make temp curves with the ATH overlayed */
 136     for(j=0;j<P_LEVELS;j++){
 137       attenuate_curve(workc[i][j],curveatt_dB[i]+100.-(j<2?2:j)*10.-P_LEVEL_0);
 138       memcpy(athc[j],ath,EHMER_MAX*sizeof(**athc));
 139       attenuate_curve(athc[j],+100.-j*10.f-P_LEVEL_0);
 140       max_curve(athc[j],workc[i][j]);
 141     }
 142
 143     /* Now limit the louder curves.
 144
 145        the idea is this: We don't know what the playback attenuation
 146        will be; 0dB SL moves every time the user twiddles the volume
 147        knob. So that means we have to use a single 'most pessimal' curve
 148        for all masking amplitudes, right?  Wrong.  The *loudest* sound
 149        can be in (we assume) a range of ...+100dB] SL.  However, sounds
 150        20dB down will be in a range ...+80], 40dB down is from ...+60],
 151        etc... */
 152
 153     for(j=1;j<P_LEVELS;j++){
 154       min_curve(athc[j],athc[j-1]);
 155       min_curve(workc[i][j],athc[j]);
 156     }
 157   }
 158
 159   for(i=0;i<P_BANDS;i++){
 160     int hi_curve,lo_curve,bin;
 161     ret[i]=_ogg_malloc(sizeof(**ret)*P_LEVELS);
 162
 163     /* low frequency curves are measured with greater resolution than
 164        the MDCT/FFT will actually give us; we want the curve applied
 165        to the tone data to be pessimistic and thus apply the minimum
 166        masking possible for a given bin.  That means that a single bin
 167        could span more than one octave and that the curve will be a
 168        composite of multiple octaves.  It also may mean that a single
 169        bin may span > an eighth of an octave and that the eighth
 170        octave values may also be composited. */
 171
 172     /* which octave curves will we be compositing? */
 173     bin=floor(fromOC(i*.5)/binHz);
 174     lo_curve=  ceil(toOC(bin*binHz+1)*2);
 175     hi_curve=  floor(toOC((bin+1)*binHz)*2);
 176     if(lo_curve>i)lo_curve=i;
 177     if(lo_curve<0)lo_curve=0;
 178     if(hi_curve>=P_BANDS)hi_curve=P_BANDS-1;
 179
 180     for(m=0;m<P_LEVELS;m++){
 181       ret[i][m]=_ogg_malloc(sizeof(***ret)*(EHMER_MAX+2));
 182
 183       for(j=0;j<n;j++)brute_buffer[j]=999.;
 184
 185       /* render the curve into bins, then pull values back into curve.
 186          The point is that any inherent subsampling aliasing results in
 187          a safe minimum */
 188       for(k=lo_curve;k<=hi_curve;k++){
 189         int l=0;
 190
 191         for(j=0;j<EHMER_MAX;j++){
 192           int lo_bin= fromOC(j*.125+k*.5-2.0625)/binHz;
 193           int hi_bin= fromOC(j*.125+k*.5-1.9375)/binHz+1;
 194
 195           if(lo_bin<0)lo_bin=0;
 196           if(lo_bin>n)lo_bin=n;
 197           if(lo_bin<l)l=lo_bin;
 198           if(hi_bin<0)hi_bin=0;
 199           if(hi_bin>n)hi_bin=n;
 200
 201           for(;l<hi_bin && l<n;l++)
 202             if(brute_buffer[l]>workc[k][m][j])
 203               brute_buffer[l]=workc[k][m][j];
 204         }
 205
 206         for(;l<n;l++)
 207           if(brute_buffer[l]>workc[k][m][EHMER_MAX-1])
 208             brute_buffer[l]=workc[k][m][EHMER_MAX-1];
 209
 210       }
 211
 212       /* be equally paranoid about being valid up to next half ocatve */
 213       if(i+1<P_BANDS){
 214         int l=0;
 215         k=i+1;
 216         for(j=0;j<EHMER_MAX;j++){
 217           int lo_bin= fromOC(j*.125+i*.5-2.0625)/binHz;
 218           int hi_bin= fromOC(j*.125+i*.5-1.9375)/binHz+1;
 219
 220           if(lo_bin<0)lo_bin=0;
 221           if(lo_bin>n)lo_bin=n;
 222           if(lo_bin<l)l=lo_bin;
 223           if(hi_bin<0)hi_bin=0;
 224           if(hi_bin>n)hi_bin=n;
 225
 226           for(;l<hi_bin && l<n;l++)
 227             if(brute_buffer[l]>workc[k][m][j])
 228               brute_buffer[l]=workc[k][m][j];
 229         }
 230
 231         for(;l<n;l++)
 232           if(brute_buffer[l]>workc[k][m][EHMER_MAX-1])
 233             brute_buffer[l]=workc[k][m][EHMER_MAX-1];
 234
 235       }
 236
 237
 238       for(j=0;j<EHMER_MAX;j++){
 239         int bin=fromOC(j*.125+i*.5-2.)/binHz;
 240         if(bin<0){
 241           ret[i][m][j+2]=-999.;
 242         }else{
 243           if(bin>=n){
 244             ret[i][m][j+2]=-999.;
 245           }else{
 246             ret[i][m][j+2]=brute_buffer[bin];
 247           }
 248         }
 249       }
 250
 251       /* add fenceposts */
 252       for(j=0;j<EHMER_OFFSET;j++)
 253         if(ret[i][m][j+2]>-200.f)break;
 254       ret[i][m][0]=j;
 255
 256       for(j=EHMER_MAX-1;j>EHMER_OFFSET+1;j--)
 257         if(ret[i][m][j+2]>-200.f)
 258           break;
 259       ret[i][m][1]=j;
 260
 261     }
 262   }
 263
 264   return(ret);
 265 }
 266
 267 void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi,
 268                   vorbis_info_psy_global *gi,int n,long rate){
 269   long i,j,lo=-99,hi=1;
 270   long maxoc;
 271   memset(p,0,sizeof(*p));
 272
 273   p->eighth_octave_lines=gi->eighth_octave_lines;
 274   p->shiftoc=rint(log(gi->eighth_octave_lines*8.f)/log(2.f))-1;
 275
 276   p->firstoc=toOC(.25f*rate*.5/n)*(1<<(p->shiftoc+1))-gi->eighth_octave_lines;
 277   maxoc=toOC((n+.25f)*rate*.5/n)*(1<<(p->shiftoc+1))+.5f;
 278   p->total_octave_lines=maxoc-p->firstoc+1;
 279   p->ath=_ogg_malloc(n*sizeof(*p->ath));
 280
 281   p->octave=_ogg_malloc(n*sizeof(*p->octave));
 282   p->bark=_ogg_malloc(n*sizeof(*p->bark));
 283   p->vi=vi;
 284   p->n=n;
 285   p->rate=rate;
 286
 287   /* AoTuV HF weighting */
 288   p->m_val = 1.;
 289   if(rate < 26000) p->m_val = 0;
 290   else if(rate < 38000) p->m_val = .94;   /* 32kHz */
 291   else if(rate > 46000) p->m_val = 1.275; /* 48kHz */
 292
 293   /* set up the lookups for a given blocksize and sample rate */
 294
 295   for(i=0,j=0;i<MAX_ATH-1;i++){
 296     int endpos=rint(fromOC((i+1)*.125-2.)*2*n/rate);
 297     float base=ATH[i];
 298     if(j<endpos){
 299       float delta=(ATH[i+1]-base)/(endpos-j);
 300       for(;j<endpos && j<n;j++){
 301         p->ath[j]=base+100.;
 302         base+=delta;
 303       }
 304     }
 305   }
 306
 307   for(i=0;i<n;i++){
 308     float bark=toBARK(rate/(2*n)*i);
 309
 310     for(;lo+vi->noisewindowlomin<i &&
 311           toBARK(rate/(2*n)*lo)<(bark-vi->noisewindowlo);lo++);
 312
 313     for(;hi<=n && (hi<i+vi->noisewindowhimin ||
 314           toBARK(rate/(2*n)*hi)<(bark+vi->noisewindowhi));hi++);
 315
 316     p->bark[i]=((lo-1)<<16)+(hi-1);
 317
 318   }
 319
 320   for(i=0;i<n;i++)
 321     p->octave[i]=toOC((i+.25f)*.5*rate/n)*(1<<(p->shiftoc+1))+.5f;
 322
 323   p->tonecurves=setup_tone_curves(vi->toneatt,rate*.5/n,n,
 324                                   vi->tone_centerboost,vi->tone_decay);
 325
 326   /* set up rolling noise median */
 327   p->noiseoffset=_ogg_malloc(P_NOISECURVES*sizeof(*p->noiseoffset));
 328   for(i=0;i<P_NOISECURVES;i++)
 329     p->noiseoffset[i]=_ogg_malloc(n*sizeof(**p->noiseoffset));
 330
 331   for(i=0;i<n;i++){
 332     float halfoc=toOC((i+.5)*rate/(2.*n))*2.;
 333     int inthalfoc;
 334     float del;
 335
 336     if(halfoc<0)halfoc=0;
 337     if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1;
 338     inthalfoc=(int)halfoc;
 339     del=halfoc-inthalfoc;
 340
 341     for(j=0;j<P_NOISECURVES;j++)
 342       p->noiseoffset[j][i]=
 343         p->vi->noiseoff[j][inthalfoc]*(1.-del) +
 344         p->vi->noiseoff[j][inthalfoc+1]*del;
 345
 346   }
 347 #if 0
 348   {
 349     static int ls=0;
 350     _analysis_output_always("noiseoff0",ls,p->noiseoffset[0],n,1,0,0);
 351     _analysis_output_always("noiseoff1",ls,p->noiseoffset[1],n,1,0,0);
 352     _analysis_output_always("noiseoff2",ls++,p->noiseoffset[2],n,1,0,0);
 353   }
 354 #endif
 355 }
 356
 357 void _vp_psy_clear(vorbis_look_psy *p){
 358   int i,j;
 359   if(p){
 360     if(p->ath)_ogg_free(p->ath);
 361     if(p->octave)_ogg_free(p->octave);
 362     if(p->bark)_ogg_free(p->bark);
 363     if(p->tonecurves){
 364       for(i=0;i<P_BANDS;i++){
 365         for(j=0;j<P_LEVELS;j++){
 366           _ogg_free(p->tonecurves[i][j]);
 367         }
 368         _ogg_free(p->tonecurves[i]);
 369       }
 370       _ogg_free(p->tonecurves);
 371     }
 372     if(p->noiseoffset){
 373       for(i=0;i<P_NOISECURVES;i++){
 374         _ogg_free(p->noiseoffset[i]);
 375       }
 376       _ogg_free(p->noiseoffset);
 377     }
 378     memset(p,0,sizeof(*p));
 379   }
 380 }
 381
 382 /* octave/(8*eighth_octave_lines) x scale and dB y scale */
 383 static void seed_curve(float *seed,
 384                        const float **curves,
 385                        float amp,
 386                        int oc, int n,
 387                        int linesper,float dBoffset){
 388   int i,post1;
 389   int seedptr;
 390   const float *posts,*curve;
 391
 392   int choice=(int)((amp+dBoffset-P_LEVEL_0)*.1f);
 393   choice=max(choice,0);
 394   choice=min(choice,P_LEVELS-1);
 395   posts=curves[choice];
 396   curve=posts+2;
 397   post1=(int)posts[1];
 398   seedptr=oc+(posts[0]-EHMER_OFFSET)*linesper-(linesper>>1);
 399
 400   for(i=posts[0];i<post1;i++){
 401     if(seedptr>0){
 402       float lin=amp+curve[i];
 403       if(seed[seedptr]<lin)seed[seedptr]=lin;
 404     }
 405     seedptr+=linesper;
 406     if(seedptr>=n)break;
 407   }
 408 }
 409
 410 static void seed_loop(vorbis_look_psy *p,
 411                       const float ***curves,
 412                       const float *f,
 413                       const float *flr,
 414                       float *seed,
 415                       float specmax){
 416   vorbis_info_psy *vi=p->vi;
 417   long n=p->n,i;
 418   float dBoffset=vi->max_curve_dB-specmax;
 419
 420   /* prime the working vector with peak values */
 421
 422   for(i=0;i<n;i++){
 423     float max=f[i];
 424     long oc=p->octave[i];
 425     while(i+1<n && p->octave[i+1]==oc){
 426       i++;
 427       if(f[i]>max)max=f[i];
 428     }
 429
 430     if(max+6.f>flr[i]){
 431       oc=oc>>p->shiftoc;
 432
 433       if(oc>=P_BANDS)oc=P_BANDS-1;
 434       if(oc<0)oc=0;
 435
 436       seed_curve(seed,
 437                  curves[oc],
 438                  max,
 439                  p->octave[i]-p->firstoc,
 440                  p->total_octave_lines,
 441                  p->eighth_octave_lines,
 442                  dBoffset);
 443     }
 444   }
 445 }
 446
 447 static void seed_chase(float *seeds, int linesper, long n){
 448   long  *posstack=alloca(n*sizeof(*posstack));
 449   float *ampstack=alloca(n*sizeof(*ampstack));
 450   long   stack=0;
 451   long   pos=0;
 452   long   i;
 453
 454   for(i=0;i<n;i++){
 455     if(stack<2){
 456       posstack[stack]=i;
 457       ampstack[stack++]=seeds[i];
 458     }else{
 459       while(1){
 460         if(seeds[i]<ampstack[stack-1]){
 461           posstack[stack]=i;
 462           ampstack[stack++]=seeds[i];
 463           break;
 464         }else{
 465           if(i<posstack[stack-1]+linesper){
 466             if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] &&
 467                i<posstack[stack-2]+linesper){
 468               /* we completely overlap, making stack-1 irrelevant.  pop it */
 469               stack--;
 470               continue;
 471             }
 472           }
 473           posstack[stack]=i;
 474           ampstack[stack++]=seeds[i];
 475           break;
 476
 477         }
 478       }
 479     }
 480   }
 481
 482   /* the stack now contains only the positions that are relevant. Scan
 483      'em straight through */
 484
 485   for(i=0;i<stack;i++){
 486     long endpos;
 487     if(i<stack-1 && ampstack[i+1]>ampstack[i]){
 488       endpos=posstack[i+1];
 489     }else{
 490       endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is
 491                                         discarded in short frames */
 492     }
 493     if(endpos>n)endpos=n;
 494     for(;pos<endpos;pos++)
 495       seeds[pos]=ampstack[i];
 496   }
 497
 498   /* there.  Linear time.  I now remember this was on a problem set I
 499      had in Grad Skool... I didn't solve it at the time ;-) */
 500
 501 }
 502
 503 /* bleaugh, this is more complicated than it needs to be */
 504 #include<stdio.h>
 505 static void max_seeds(vorbis_look_psy *p,
 506                       float *seed,
 507                       float *flr){
 508   long   n=p->total_octave_lines;
 509   int    linesper=p->eighth_octave_lines;
 510   long   linpos=0;
 511   long   pos;
 512
 513   seed_chase(seed,linesper,n); /* for masking */
 514
 515   pos=p->octave[0]-p->firstoc-(linesper>>1);
 516
 517   while(linpos+1<p->n){
 518     float minV=seed[pos];
 519     long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc;
 520     if(minV>p->vi->tone_abs_limit)minV=p->vi->tone_abs_limit;
 521     while(pos+1<=end){
 522       pos++;
 523       if((seed[pos]>NEGINF && seed[pos]<minV) || minV==NEGINF)
 524         minV=seed[pos];
 525     }
 526
 527     end=pos+p->firstoc;
 528     for(;linpos<p->n && p->octave[linpos]<=end;linpos++)
 529       if(flr[linpos]<minV)flr[linpos]=minV;
 530   }
 531
 532   {
 533     float minV=seed[p->total_octave_lines-1];
 534     for(;linpos<p->n;linpos++)
 535       if(flr[linpos]<minV)flr[linpos]=minV;
 536   }
 537
 538 }
 539
 540 static void bark_noise_hybridmp(int n,const long *b,
 541                                 const float *f,
 542                                 float *noise,
 543                                 const float offset,
 544                                 const int fixed){
 545
 546   float *N=alloca(n*sizeof(*N));
 547   float *X=alloca(n*sizeof(*N));
 548   float *XX=alloca(n*sizeof(*N));
 549   float *Y=alloca(n*sizeof(*N));
 550   float *XY=alloca(n*sizeof(*N));
 551
 552   float tN, tX, tXX, tY, tXY;
 553   int i;
 554
 555   int lo, hi;
 556   float R=0.f;
 557   float A=0.f;
 558   float B=0.f;
 559   float D=1.f;
 560   float w, x, y;
 561
 562   tN = tX = tXX = tY = tXY = 0.f;
 563
 564   y = f[0] + offset;
 565   if (y < 1.f) y = 1.f;
 566
 567   w = y * y * .5;
 568
 569   tN += w;
 570   tX += w;
 571   tY += w * y;
 572
 573   N[0] = tN;
 574   X[0] = tX;
 575   XX[0] = tXX;
 576   Y[0] = tY;
 577   XY[0] = tXY;
 578
 579   for (i = 1, x = 1.f; i < n; i++, x += 1.f) {
 580
 581     y = f[i] + offset;
 582     if (y < 1.f) y = 1.f;
 583
 584     w = y * y;
 585
 586     tN += w;
 587     tX += w * x;
 588     tXX += w * x * x;
 589     tY += w * y;
 590     tXY += w * x * y;
 591
 592     N[i] = tN;
 593     X[i] = tX;
 594     XX[i] = tXX;
 595     Y[i] = tY;
 596     XY[i] = tXY;
 597   }
 598
 599   for (i = 0, x = 0.f;; i++, x += 1.f) {
 600
 601     lo = b[i] >> 16;
 602     if( lo>=0 ) break;
 603     hi = b[i] & 0xffff;
 604
 605     tN = N[hi] + N[-lo];
 606     tX = X[hi] - X[-lo];
 607     tXX = XX[hi] + XX[-lo];
 608     tY = Y[hi] + Y[-lo];
 609     tXY = XY[hi] - XY[-lo];
 610
 611     A = tY * tXX - tX * tXY;
 612     B = tN * tXY - tX * tY;
 613     D = tN * tXX - tX * tX;
 614     R = (A + x * B) / D;
 615     if (R < 0.f)
 616       R = 0.f;
 617
 618     noise[i] = R - offset;
 619   }
 620
 621   for ( ;; i++, x += 1.f) {
 622
 623     lo = b[i] >> 16;
 624     hi = b[i] & 0xffff;
 625     if(hi>=n)break;
 626
 627     tN = N[hi] - N[lo];
 628     tX = X[hi] - X[lo];
 629     tXX = XX[hi] - XX[lo];
 630     tY = Y[hi] - Y[lo];
 631     tXY = XY[hi] - XY[lo];
 632
 633     A = tY * tXX - tX * tXY;
 634     B = tN * tXY - tX * tY;
 635     D = tN * tXX - tX * tX;
 636     R = (A + x * B) / D;
 637     if (R < 0.f) R = 0.f;
 638
 639     noise[i] = R - offset;
 640   }
 641   for ( ; i < n; i++, x += 1.f) {
 642
 643     R = (A + x * B) / D;
 644     if (R < 0.f) R = 0.f;
 645
 646     noise[i] = R - offset;
 647   }
 648
 649   if (fixed <= 0) return;
 650
 651   for (i = 0, x = 0.f;; i++, x += 1.f) {
 652     hi = i + fixed / 2;
 653     lo = hi - fixed;
 654     if(lo>=0)break;
 655
 656     tN = N[hi] + N[-lo];
 657     tX = X[hi] - X[-lo];
 658     tXX = XX[hi] + XX[-lo];
 659     tY = Y[hi] + Y[-lo];
 660     tXY = XY[hi] - XY[-lo];
 661
 662
 663     A = tY * tXX - tX * tXY;
 664     B = tN * tXY - tX * tY;
 665     D = tN * tXX - tX * tX;
 666     R = (A + x * B) / D;
 667
 668     if (R - offset < noise[i]) noise[i] = R - offset;
 669   }
 670   for ( ;; i++, x += 1.f) {
 671
 672     hi = i + fixed / 2;
 673     lo = hi - fixed;
 674     if(hi>=n)break;
 675
 676     tN = N[hi] - N[lo];
 677     tX = X[hi] - X[lo];
 678     tXX = XX[hi] - XX[lo];
 679     tY = Y[hi] - Y[lo];
 680     tXY = XY[hi] - XY[lo];
 681
 682     A = tY * tXX - tX * tXY;
 683     B = tN * tXY - tX * tY;
 684     D = tN * tXX - tX * tX;
 685     R = (A + x * B) / D;
 686
 687     if (R - offset < noise[i]) noise[i] = R - offset;
 688   }
 689   for ( ; i < n; i++, x += 1.f) {
 690     R = (A + x * B) / D;
 691     if (R - offset < noise[i]) noise[i] = R - offset;
 692   }
 693 }
 694
 695 static float FLOOR1_fromdB_INV_LOOKUP[256]={
 696   0.F, 8.81683e+06F, 8.27882e+06F, 7.77365e+06F,
 697   7.29930e+06F, 6.85389e+06F, 6.43567e+06F, 6.04296e+06F,
 698   5.67422e+06F, 5.32798e+06F, 5.00286e+06F, 4.69759e+06F,
 699   4.41094e+06F, 4.14178e+06F, 3.88905e+06F, 3.65174e+06F,
 700   3.42891e+06F, 3.21968e+06F, 3.02321e+06F, 2.83873e+06F,
 701   2.66551e+06F, 2.50286e+06F, 2.35014e+06F, 2.20673e+06F,
 702   2.07208e+06F, 1.94564e+06F, 1.82692e+06F, 1.71544e+06F,
 703   1.61076e+06F, 1.51247e+06F, 1.42018e+06F, 1.33352e+06F,
 704   1.25215e+06F, 1.17574e+06F, 1.10400e+06F, 1.03663e+06F,
 705   973377.F, 913981.F, 858210.F, 805842.F,
 706   756669.F, 710497.F, 667142.F, 626433.F,
 707   588208.F, 552316.F, 518613.F, 486967.F,
 708   457252.F, 429351.F, 403152.F, 378551.F,
 709   355452.F, 333762.F, 313396.F, 294273.F,
 710   276316.F, 259455.F, 243623.F, 228757.F,
 711   214798.F, 201691.F, 189384.F, 177828.F,
 712   166977.F, 156788.F, 147221.F, 138237.F,
 713   129802.F, 121881.F, 114444.F, 107461.F,
 714   100903.F, 94746.3F, 88964.9F, 83536.2F,
 715   78438.8F, 73652.5F, 69158.2F, 64938.1F,
 716   60975.6F, 57254.9F, 53761.2F, 50480.6F,
 717   47400.3F, 44507.9F, 41792.0F, 39241.9F,
 718   36847.3F, 34598.9F, 32487.7F, 30505.3F,
 719   28643.8F, 26896.0F, 25254.8F, 23713.7F,
 720   22266.7F, 20908.0F, 19632.2F, 18434.2F,
 721   17309.4F, 16253.1F, 15261.4F, 14330.1F,
 722   13455.7F, 12634.6F, 11863.7F, 11139.7F,
 723   10460.0F, 9821.72F, 9222.39F, 8659.64F,
 724   8131.23F, 7635.06F, 7169.17F, 6731.70F,
 725   6320.93F, 5935.23F, 5573.06F, 5232.99F,
 726   4913.67F, 4613.84F, 4332.30F, 4067.94F,
 727   3819.72F, 3586.64F, 3367.78F, 3162.28F,
 728   2969.31F, 2788.13F, 2617.99F, 2458.24F,
 729   2308.24F, 2167.39F, 2035.14F, 1910.95F,
 730   1794.35F, 1684.85F, 1582.04F, 1485.51F,
 731   1394.86F, 1309.75F, 1229.83F, 1154.78F,
 732   1084.32F, 1018.15F, 956.024F, 897.687F,
 733   842.910F, 791.475F, 743.179F, 697.830F,
 734   655.249F, 615.265F, 577.722F, 542.469F,
 735   509.367F, 478.286F, 449.101F, 421.696F,
 736   395.964F, 371.803F, 349.115F, 327.812F,
 737   307.809F, 289.026F, 271.390F, 254.830F,
 738   239.280F, 224.679F, 210.969F, 198.096F,
 739   186.008F, 174.658F, 164.000F, 153.993F,
 740   144.596F, 135.773F, 127.488F, 119.708F,
 741   112.404F, 105.545F, 99.1046F, 93.0572F,
 742   87.3788F, 82.0469F, 77.0404F, 72.3394F,
 743   67.9252F, 63.7804F, 59.8885F, 56.2341F,
 744   52.8027F, 49.5807F, 46.5553F, 43.7144F,
 745   41.0470F, 38.5423F, 36.1904F, 33.9821F,
 746   31.9085F, 29.9614F, 28.1332F, 26.4165F,
 747   24.8045F, 23.2910F, 21.8697F, 20.5352F,
 748   19.2822F, 18.1056F, 17.0008F, 15.9634F,
 749   14.9893F, 14.0746F, 13.2158F, 12.4094F,
 750   11.6522F, 10.9411F, 10.2735F, 9.64662F,
 751   9.05798F, 8.50526F, 7.98626F, 7.49894F,
 752   7.04135F, 6.61169F, 6.20824F, 5.82941F,
 753   5.47370F, 5.13970F, 4.82607F, 4.53158F,
 754   4.25507F, 3.99542F, 3.75162F, 3.52269F,
 755   3.30774F, 3.10590F, 2.91638F, 2.73842F,
 756   2.57132F, 2.41442F, 2.26709F, 2.12875F,
 757   1.99885F, 1.87688F, 1.76236F, 1.65482F,
 758   1.55384F, 1.45902F, 1.36999F, 1.28640F,
 759   1.20790F, 1.13419F, 1.06499F, 1.F
 760 };
 761
 762 void _vp_remove_floor(vorbis_look_psy *p,
 763                       float *mdct,
 764                       int *codedflr,
 765                       float *residue,
 766                       int sliding_lowpass){
 767
 768   int i,n=p->n;
 769
 770   if(sliding_lowpass>n)sliding_lowpass=n;
 771
 772   for(i=0;i<sliding_lowpass;i++){
 773     residue[i]=
 774       mdct[i]*FLOOR1_fromdB_INV_LOOKUP[codedflr[i]];
 775   }
 776
 777   for(;i<n;i++)
 778     residue[i]=0.;
 779 }
 780
 781 void _vp_noisemask(vorbis_look_psy *p,
 782                    float *logmdct,
 783                    float *logmask){
 784
 785   int i,n=p->n;
 786   float *work=alloca(n*sizeof(*work));
 787
 788   bark_noise_hybridmp(n,p->bark,logmdct,logmask,
 789                       140.,-1);
 790
 791   for(i=0;i<n;i++)work[i]=logmdct[i]-logmask[i];
 792
 793   bark_noise_hybridmp(n,p->bark,work,logmask,0.,
 794                       p->vi->noisewindowfixed);
 795
 796   for(i=0;i<n;i++)work[i]=logmdct[i]-work[i];
 797
 798 #if 0
 799   {
 800     static int seq=0;
 801
 802     float work2[n];
 803     for(i=0;i<n;i++){
 804       work2[i]=logmask[i]+work[i];
 805     }
 806
 807     if(seq&1)
 808       _analysis_output("median2R",seq/2,work,n,1,0,0);
 809     else
 810       _analysis_output("median2L",seq/2,work,n,1,0,0);
 811
 812     if(seq&1)
 813       _analysis_output("envelope2R",seq/2,work2,n,1,0,0);
 814     else
 815       _analysis_output("envelope2L",seq/2,work2,n,1,0,0);
 816     seq++;
 817   }
 818 #endif
 819
 820   for(i=0;i<n;i++){
 821     int dB=logmask[i]+.5;
 822     if(dB>=NOISE_COMPAND_LEVELS)dB=NOISE_COMPAND_LEVELS-1;
 823     if(dB<0)dB=0;
 824     logmask[i]= work[i]+p->vi->noisecompand[dB];
 825   }
 826
 827 }
 828
 829 void _vp_tonemask(vorbis_look_psy *p,
 830                   float *logfft,
 831                   float *logmask,
 832                   float global_specmax,
 833                   float local_specmax){
 834
 835   int i,n=p->n;
 836
 837   float *seed=alloca(sizeof(*seed)*p->total_octave_lines);
 838   float att=local_specmax+p->vi->ath_adjatt;
 839   for(i=0;i<p->total_octave_lines;i++)seed[i]=NEGINF;
 840
 841   /* set the ATH (floating below localmax, not global max by a
 842      specified att) */
 843   if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt;
 844
 845   for(i=0;i<n;i++)
 846     logmask[i]=p->ath[i]+att;
 847
 848   /* tone masking */
 849   seed_loop(p,(const float ***)p->tonecurves,logfft,logmask,seed,global_specmax);
 850   max_seeds(p,seed,logmask);
 851
 852 }
 853
 854 void _vp_offset_and_mix(vorbis_look_psy *p,
 855                         float *noise,
 856                         float *tone,
 857                         int offset_select,
 858                         float *logmask,
 859                         float *mdct,
 860                         float *logmdct){
 861   int i,n=p->n;
 862   float de, coeffi, cx;/* AoTuV */
 863   float toneatt=p->vi->tone_masteratt[offset_select];
 864
 865   cx = p->m_val;
 866
 867   for(i=0;i<n;i++){
 868     float val= noise[i]+p->noiseoffset[offset_select][i];
 869     if(val>p->vi->noisemaxsupp)val=p->vi->noisemaxsupp;
 870     logmask[i]=max(val,tone[i]+toneatt);
 871
 872
 873     /* AoTuV */
 874     /** @ M1 **
 875         The following codes improve a noise problem.
 876         A fundamental idea uses the value of masking and carries out
 877         the relative compensation of the MDCT.
 878         However, this code is not perfect and all noise problems cannot be solved.
 879         by Aoyumi @ 2004/04/18
 880     */
 881
 882     if(offset_select == 1) {
 883       coeffi = -17.2;       /* coeffi is a -17.2dB threshold */
 884       val = val - logmdct[i];  /* val == mdct line value relative to floor in dB */
 885
 886       if(val > coeffi){
 887         /* mdct value is > -17.2 dB below floor */
 888
 889         de = 1.0-((val-coeffi)*0.005*cx);
 890         /* pro-rated attenuation:
 891            -0.00 dB boost if mdct value is -17.2dB (relative to floor)
 892            -0.77 dB boost if mdct value is 0dB (relative to floor)
 893            -1.64 dB boost if mdct value is +17.2dB (relative to floor)
 894            etc... */
 895
 896         if(de < 0) de = 0.0001;
 897       }else
 898         /* mdct value is <= -17.2 dB below floor */
 899
 900         de = 1.0-((val-coeffi)*0.0003*cx);
 901       /* pro-rated attenuation:
 902          +0.00 dB atten if mdct value is -17.2dB (relative to floor)
 903          +0.45 dB atten if mdct value is -34.4dB (relative to floor)
 904          etc... */
 905
 906       mdct[i] *= de;
 907
 908     }
 909   }
 910 }
 911
 912 float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){
 913   vorbis_info *vi=vd->vi;
 914   codec_setup_info *ci=vi->codec_setup;
 915   vorbis_info_psy_global *gi=&ci->psy_g_param;
 916
 917   int n=ci->blocksizes[vd->W]/2;
 918   float secs=(float)n/vi->rate;
 919
 920   amp+=secs*gi->ampmax_att_per_sec;
 921   if(amp<-9999)amp=-9999;
 922   return(amp);
 923 }
 924
 925 static void couple_lossless(float A, float B,
 926                             float *qA, float *qB){
 927   int test1=fabs(*qA)>fabs(*qB);
 928   test1-= fabs(*qA)<fabs(*qB);
 929
 930   if(!test1)test1=((fabs(A)>fabs(B))<<1)-1;
 931   if(test1==1){
 932     *qB=(*qA>0.f?*qA-*qB:*qB-*qA);
 933   }else{
 934     float temp=*qB;
 935     *qB=(*qB>0.f?*qA-*qB:*qB-*qA);
 936     *qA=temp;
 937   }
 938
 939   if(*qB>fabs(*qA)*1.9999f){
 940     *qB= -fabs(*qA)*2.f;
 941     *qA= -*qA;
 942   }
 943 }
 944
 945 static float hypot_lookup[32]={
 946   -0.009935, -0.011245, -0.012726, -0.014397,
 947   -0.016282, -0.018407, -0.020800, -0.023494,
 948   -0.026522, -0.029923, -0.033737, -0.038010,
 949   -0.042787, -0.048121, -0.054064, -0.060671,
 950   -0.068000, -0.076109, -0.085054, -0.094892,
 951   -0.105675, -0.117451, -0.130260, -0.144134,
 952   -0.159093, -0.175146, -0.192286, -0.210490,
 953   -0.229718, -0.249913, -0.271001, -0.292893};
 954
 955 static void precomputed_couple_point(float premag,
 956                                      int floorA,int floorB,
 957                                      float *mag, float *ang){
 958
 959   int test=(floorA>floorB)-1;
 960   int offset=31-abs(floorA-floorB);
 961   float floormag=hypot_lookup[((offset<0)-1)&offset]+1.f;
 962
 963   floormag*=FLOOR1_fromdB_INV_LOOKUP[(floorB&test)|(floorA&(~test))];
 964
 965   *mag=premag*floormag;
 966   *ang=0.f;
 967 }
 968
 969 /* just like below, this is currently set up to only do
 970    single-step-depth coupling.  Otherwise, we'd have to do more
 971    copying (which will be inevitable later) */
 972
 973 /* doing the real circular magnitude calculation is audibly superior
 974    to (A+B)/sqrt(2) */
 975 static float dipole_hypot(float a, float b){
 976   if(a>0.){
 977     if(b>0.)return sqrt(a*a+b*b);
 978     if(a>-b)return sqrt(a*a-b*b);
 979     return -sqrt(b*b-a*a);
 980   }
 981   if(b<0.)return -sqrt(a*a+b*b);
 982   if(-a>b)return -sqrt(a*a-b*b);
 983   return sqrt(b*b-a*a);
 984 }
 985 static float round_hypot(float a, float b){
 986   if(a>0.){
 987     if(b>0.)return sqrt(a*a+b*b);
 988     if(a>-b)return sqrt(a*a+b*b);
 989     return -sqrt(b*b+a*a);
 990   }
 991   if(b<0.)return -sqrt(a*a+b*b);
 992   if(-a>b)return -sqrt(a*a+b*b);
 993   return sqrt(b*b+a*a);
 994 }
 995
 996 /* revert to round hypot for now */
 997 float **_vp_quantize_couple_memo(vorbis_block *vb,
 998                                  vorbis_info_psy_global *g,
 999                                  vorbis_look_psy *p,
1000                                  vorbis_info_mapping0 *vi,
1001                                  float **mdct){
1002
1003   int i,j,n=p->n;
1004   float **ret=_vorbis_block_alloc(vb,vi->coupling_steps*sizeof(*ret));
1005   int limit=g->coupling_pointlimit[p->vi->blockflag][PACKETBLOBS/2];
1006
1007   for(i=0;i<vi->coupling_steps;i++){
1008     float *mdctM=mdct[vi->coupling_mag[i]];
1009     float *mdctA=mdct[vi->coupling_ang[i]];
1010     ret[i]=_vorbis_block_alloc(vb,n*sizeof(**ret));
1011     for(j=0;j<limit;j++)
1012       ret[i][j]=dipole_hypot(mdctM[j],mdctA[j]);
1013     for(;j<n;j++)
1014       ret[i][j]=round_hypot(mdctM[j],mdctA[j]);
1015   }
1016
1017   return(ret);
1018 }
1019
1020 /* this is for per-channel noise normalization */
1021 static int apsort(const void *a, const void *b){
1022   float f1=fabs(**(float**)a);
1023   float f2=fabs(**(float**)b);
1024   return (f1<f2)-(f1>f2);
1025 }
1026
1027 int **_vp_quantize_couple_sort(vorbis_block *vb,
1028                                vorbis_look_psy *p,
1029                                vorbis_info_mapping0 *vi,
1030                                float **mags){
1031
1032
1033   if(p->vi->normal_point_p){
1034     int i,j,k,n=p->n;
1035     int **ret=_vorbis_block_alloc(vb,vi->coupling_steps*sizeof(*ret));
1036     int partition=p->vi->normal_partition;
1037     float **work=alloca(sizeof(*work)*partition);
1038
1039     for(i=0;i<vi->coupling_steps;i++){
1040       ret[i]=_vorbis_block_alloc(vb,n*sizeof(**ret));
1041
1042       for(j=0;j<n;j+=partition){
1043         for(k=0;k<partition;k++)work[k]=mags[i]+k+j;
1044         qsort(work,partition,sizeof(*work),apsort);
1045         for(k=0;k<partition;k++)ret[i][k+j]=work[k]-mags[i];
1046       }
1047     }
1048     return(ret);
1049   }
1050   return(NULL);
1051 }
1052
1053 void _vp_noise_normalize_sort(vorbis_look_psy *p,
1054                               float *magnitudes,int *sortedindex){
1055   int i,j,n=p->n;
1056   vorbis_info_psy *vi=p->vi;
1057   int partition=vi->normal_partition;
1058   float **work=alloca(sizeof(*work)*partition);
1059   int start=vi->normal_start;
1060
1061   for(j=start;j<n;j+=partition){
1062     if(j+partition>n)partition=n-j;
1063     for(i=0;i<partition;i++)work[i]=magnitudes+i+j;
1064     qsort(work,partition,sizeof(*work),apsort);
1065     for(i=0;i<partition;i++){
1066       sortedindex[i+j-start]=work[i]-magnitudes;
1067     }
1068   }
1069 }
1070
1071 void _vp_noise_normalize(vorbis_look_psy *p,
1072                          float *in,float *out,int *sortedindex){
1073   int flag=0,i,j=0,n=p->n;
1074   vorbis_info_psy *vi=p->vi;
1075   int partition=vi->normal_partition;
1076   int start=vi->normal_start;
1077
1078   if(start>n)start=n;
1079
1080   if(vi->normal_channel_p){
1081     for(;j<start;j++)
1082       out[j]=rint(in[j]);
1083
1084     for(;j+partition<=n;j+=partition){
1085       float acc=0.;
1086       int k;
1087
1088       for(i=j;i<j+partition;i++)
1089         acc+=in[i]*in[i];
1090
1091       for(i=0;i<partition;i++){
1092         k=sortedindex[i+j-start];
1093
1094         if(in[k]*in[k]>=.25f){
1095           out[k]=rint(in[k]);
1096           acc-=in[k]*in[k];
1097           flag=1;
1098         }else{
1099           if(acc<vi->normal_thresh)break;
1100           out[k]=unitnorm(in[k]);
1101           acc-=1.;
1102         }
1103       }
1104
1105       for(;i<partition;i++){
1106         k=sortedindex[i+j-start];
1107         out[k]=0.;
1108       }
1109     }
1110   }
1111
1112   for(;j<n;j++)
1113     out[j]=rint(in[j]);
1114
1115 }
1116
1117 void _vp_couple(int blobno,
1118                 vorbis_info_psy_global *g,
1119                 vorbis_look_psy *p,
1120                 vorbis_info_mapping0 *vi,
1121                 float **res,
1122                 float **mag_memo,
1123                 int   **mag_sort,
1124                 int   **ifloor,
1125                 int   *nonzero,
1126                 int  sliding_lowpass){
1127
1128   int i,j,k,n=p->n;
1129
1130   /* perform any requested channel coupling */
1131   /* point stereo can only be used in a first stage (in this encoder)
1132      because of the dependency on floor lookups */
1133   for(i=0;i<vi->coupling_steps;i++){
1134
1135     /* once we're doing multistage coupling in which a channel goes
1136        through more than one coupling step, the floor vector
1137        magnitudes will also have to be recalculated an propogated
1138        along with PCM.  Right now, we're not (that will wait until 5.1
1139        most likely), so the code isn't here yet. The memory management
1140        here is all assuming single depth couplings anyway. */
1141
1142     /* make sure coupling a zero and a nonzero channel results in two
1143        nonzero channels. */
1144     if(nonzero[vi->coupling_mag[i]] ||
1145        nonzero[vi->coupling_ang[i]]){
1146
1147
1148       float *rM=res[vi->coupling_mag[i]];
1149       float *rA=res[vi->coupling_ang[i]];
1150       float *qM=rM+n;
1151       float *qA=rA+n;
1152       int *floorM=ifloor[vi->coupling_mag[i]];
1153       int *floorA=ifloor[vi->coupling_ang[i]];
1154       float prepoint=stereo_threshholds[g->coupling_prepointamp[blobno]];
1155       float postpoint=stereo_threshholds[g->coupling_postpointamp[blobno]];
1156       int partition=(p->vi->normal_point_p?p->vi->normal_partition:p->n);
1157       int limit=g->coupling_pointlimit[p->vi->blockflag][blobno];
1158       int pointlimit=limit;
1159
1160       nonzero[vi->coupling_mag[i]]=1;
1161       nonzero[vi->coupling_ang[i]]=1;
1162
1163        /* The threshold of a stereo is changed with the size of n */
1164        if(n > 1000)
1165          postpoint=stereo_threshholds_limited[g->coupling_postpointamp[blobno]];
1166
1167       for(j=0;j<p->n;j+=partition){
1168         float acc=0.f;
1169
1170         for(k=0;k<partition;k++){
1171           int l=k+j;
1172
1173           if(l<sliding_lowpass){
1174             if((l>=limit && fabs(rM[l])<postpoint && fabs(rA[l])<postpoint) ||
1175                (fabs(rM[l])<prepoint && fabs(rA[l])<prepoint)){
1176
1177
1178               precomputed_couple_point(mag_memo[i][l],
1179                                        floorM[l],floorA[l],
1180                                        qM+l,qA+l);
1181
1182               if(rint(qM[l])==0.f)acc+=qM[l]*qM[l];
1183             }else{
1184               couple_lossless(rM[l],rA[l],qM+l,qA+l);
1185             }
1186           }else{
1187             qM[l]=0.;
1188             qA[l]=0.;
1189           }
1190         }
1191
1192         if(p->vi->normal_point_p){
1193           for(k=0;k<partition && acc>=p->vi->normal_thresh;k++){
1194             int l=mag_sort[i][j+k];
1195             if(l<sliding_lowpass && l>=pointlimit && rint(qM[l])==0.f){
1196               qM[l]=unitnorm(qM[l]);
1197               acc-=1.f;
1198             }
1199           }
1200         }
1201       }
1202     }
1203   }
1204 }
1205
1206 /* AoTuV */
1207 /** @ M2 **
1208    The boost problem by the combination of noise normalization and point stereo is eased.
1209    However, this is a temporary patch.
1210    by Aoyumi @ 2004/04/18
1211 */
1212
1213 void hf_reduction(vorbis_info_psy_global *g,
1214                       vorbis_look_psy *p,
1215                       vorbis_info_mapping0 *vi,
1216                       float **mdct){
1217
1218   int i,j,n=p->n, de=0.3*p->m_val;
1219   int limit=g->coupling_pointlimit[p->vi->blockflag][PACKETBLOBS/2];
1220
1221   for(i=0; i<vi->coupling_steps; i++){
1222     /* for(j=start; j<limit; j++){} // ???*/
1223     for(j=limit; j<n; j++)
1224       mdct[i][j] *= (1.0 - de*((float)(j-limit) / (float)(n-limit)));
1225   }
1226 }