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-2010             *
   9  * by the Xiph.Org Foundation https://xiph.org/                     *
  10  *                                                                  *
  11  ********************************************************************
  12
  13  function: psychoacoustics not including preecho
  14
  15  ********************************************************************/
  16
  17 #include <stdlib.h>
  18 #include <math.h>
  19 #include <string.h>
  20 #include "vorbis/codec.h"
  21 #include "codec_internal.h"
  22
  23 #include "masking.h"
  24 #include "psy.h"
  25 #include "os.h"
  26 #include "lpc.h"
  27 #include "smallft.h"
  28 #include "scales.h"
  29 #include "misc.h"
  30
  31 #define NEGINF -9999.f
  32 static const double stereo_threshholds[]={0.0, .5, 1.0, 1.5, 2.5, 4.5, 8.5, 16.5, 9e10};
  33 static const double stereo_threshholds_limited[]={0.0, .5, 1.0, 1.5, 2.0, 2.5, 4.5, 8.5, 9e10};
  34
  35 vorbis_look_psy_global *_vp_global_look(vorbis_info *vi){
  36   codec_setup_info *ci=vi->codec_setup;
  37   vorbis_info_psy_global *gi=&ci->psy_g_param;
  38   vorbis_look_psy_global *look=_ogg_calloc(1,sizeof(*look));
  39
  40   look->channels=vi->channels;
  41
  42   look->ampmax=-9999.;
  43   look->gi=gi;
  44   return(look);
  45 }
  46
  47 void _vp_global_free(vorbis_look_psy_global *look){
  48   if(look){
  49     memset(look,0,sizeof(*look));
  50     _ogg_free(look);
  51   }
  52 }
  53
  54 void _vi_gpsy_free(vorbis_info_psy_global *i){
  55   if(i){
  56     memset(i,0,sizeof(*i));
  57     _ogg_free(i);
  58   }
  59 }
  60
  61 void _vi_psy_free(vorbis_info_psy *i){
  62   if(i){
  63     memset(i,0,sizeof(*i));
  64     _ogg_free(i);
  65   }
  66 }
  67
  68 static void min_curve(float *c,
  69                        float *c2){
  70   int i;
  71   for(i=0;i<EHMER_MAX;i++)if(c2[i]<c[i])c[i]=c2[i];
  72 }
  73 static void max_curve(float *c,
  74                        float *c2){
  75   int i;
  76   for(i=0;i<EHMER_MAX;i++)if(c2[i]>c[i])c[i]=c2[i];
  77 }
  78
  79 static void attenuate_curve(float *c,float att){
  80   int i;
  81   for(i=0;i<EHMER_MAX;i++)
  82     c[i]+=att;
  83 }
  84
  85 static float ***setup_tone_curves(float curveatt_dB[P_BANDS],float binHz,int n,
  86                                   float center_boost, float center_decay_rate){
  87   int i,j,k,m;
  88   float ath[EHMER_MAX];
  89   float workc[P_BANDS][P_LEVELS][EHMER_MAX];
  90   float athc[P_LEVELS][EHMER_MAX];
  91   float *brute_buffer=alloca(n*sizeof(*brute_buffer));
  92
  93   float ***ret=_ogg_malloc(sizeof(*ret)*P_BANDS);
  94
  95   memset(workc,0,sizeof(workc));
  96
  97   for(i=0;i<P_BANDS;i++){
  98     /* we add back in the ATH to avoid low level curves falling off to
  99        -infinity and unnecessarily cutting off high level curves in the
 100        curve limiting (last step). */
 101
 102     /* A half-band's settings must be valid over the whole band, and
 103        it's better to mask too little than too much */
 104     int ath_offset=i*4;
 105     for(j=0;j<EHMER_MAX;j++){
 106       float min=999.;
 107       for(k=0;k<4;k++)
 108         if(j+k+ath_offset<MAX_ATH){
 109           if(min>ATH[j+k+ath_offset])min=ATH[j+k+ath_offset];
 110         }else{
 111           if(min>ATH[MAX_ATH-1])min=ATH[MAX_ATH-1];
 112         }
 113       ath[j]=min;
 114     }
 115
 116     /* copy curves into working space, replicate the 50dB curve to 30
 117        and 40, replicate the 100dB curve to 110 */
 118     for(j=0;j<6;j++)
 119       memcpy(workc[i][j+2],tonemasks[i][j],EHMER_MAX*sizeof(*tonemasks[i][j]));
 120     memcpy(workc[i][0],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0]));
 121     memcpy(workc[i][1],tonemasks[i][0],EHMER_MAX*sizeof(*tonemasks[i][0]));
 122
 123     /* apply centered curve boost/decay */
 124     for(j=0;j<P_LEVELS;j++){
 125       for(k=0;k<EHMER_MAX;k++){
 126         float adj=center_boost+abs(EHMER_OFFSET-k)*center_decay_rate;
 127         if(adj<0. && center_boost>0)adj=0.;
 128         if(adj>0. && center_boost<0)adj=0.;
 129         workc[i][j][k]+=adj;
 130       }
 131     }
 132
 133     /* normalize curves so the driving amplitude is 0dB */
 134     /* make temp curves with the ATH overlayed */
 135     for(j=0;j<P_LEVELS;j++){
 136       attenuate_curve(workc[i][j],curveatt_dB[i]+100.-(j<2?2:j)*10.-P_LEVEL_0);
 137       memcpy(athc[j],ath,EHMER_MAX*sizeof(**athc));
 138       attenuate_curve(athc[j],+100.-j*10.f-P_LEVEL_0);
 139       max_curve(athc[j],workc[i][j]);
 140     }
 141
 142     /* Now limit the louder curves.
 143
 144        the idea is this: We don't know what the playback attenuation
 145        will be; 0dB SL moves every time the user twiddles the volume
 146        knob. So that means we have to use a single 'most pessimal' curve
 147        for all masking amplitudes, right?  Wrong.  The *loudest* sound
 148        can be in (we assume) a range of ...+100dB] SL.  However, sounds
 149        20dB down will be in a range ...+80], 40dB down is from ...+60],
 150        etc... */
 151
 152     for(j=1;j<P_LEVELS;j++){
 153       min_curve(athc[j],athc[j-1]);
 154       min_curve(workc[i][j],athc[j]);
 155     }
 156   }
 157
 158   for(i=0;i<P_BANDS;i++){
 159     int hi_curve,lo_curve,bin;
 160     ret[i]=_ogg_malloc(sizeof(**ret)*P_LEVELS);
 161
 162     /* low frequency curves are measured with greater resolution than
 163        the MDCT/FFT will actually give us; we want the curve applied
 164        to the tone data to be pessimistic and thus apply the minimum
 165        masking possible for a given bin.  That means that a single bin
 166        could span more than one octave and that the curve will be a
 167        composite of multiple octaves.  It also may mean that a single
 168        bin may span > an eighth of an octave and that the eighth
 169        octave values may also be composited. */
 170
 171     /* which octave curves will we be compositing? */
 172     bin=floor(fromOC(i*.5)/binHz);
 173     lo_curve=  ceil(toOC(bin*binHz+1)*2);
 174     hi_curve=  floor(toOC((bin+1)*binHz)*2);
 175     if(lo_curve>i)lo_curve=i;
 176     if(lo_curve<0)lo_curve=0;
 177     if(hi_curve>=P_BANDS)hi_curve=P_BANDS-1;
 178
 179     for(m=0;m<P_LEVELS;m++){
 180       ret[i][m]=_ogg_malloc(sizeof(***ret)*(EHMER_MAX+2));
 181
 182       for(j=0;j<n;j++)brute_buffer[j]=999.;
 183
 184       /* render the curve into bins, then pull values back into curve.
 185          The point is that any inherent subsampling aliasing results in
 186          a safe minimum */
 187       for(k=lo_curve;k<=hi_curve;k++){
 188         int l=0;
 189
 190         for(j=0;j<EHMER_MAX;j++){
 191           int lo_bin= fromOC(j*.125+k*.5-2.0625)/binHz;
 192           int hi_bin= fromOC(j*.125+k*.5-1.9375)/binHz+1;
 193
 194           if(lo_bin<0)lo_bin=0;
 195           if(lo_bin>n)lo_bin=n;
 196           if(lo_bin<l)l=lo_bin;
 197           if(hi_bin<0)hi_bin=0;
 198           if(hi_bin>n)hi_bin=n;
 199
 200           for(;l<hi_bin && l<n;l++)
 201             if(brute_buffer[l]>workc[k][m][j])
 202               brute_buffer[l]=workc[k][m][j];
 203         }
 204
 205         for(;l<n;l++)
 206           if(brute_buffer[l]>workc[k][m][EHMER_MAX-1])
 207             brute_buffer[l]=workc[k][m][EHMER_MAX-1];
 208
 209       }
 210
 211       /* be equally paranoid about being valid up to next half ocatve */
 212       if(i+1<P_BANDS){
 213         int l=0;
 214         k=i+1;
 215         for(j=0;j<EHMER_MAX;j++){
 216           int lo_bin= fromOC(j*.125+i*.5-2.0625)/binHz;
 217           int hi_bin= fromOC(j*.125+i*.5-1.9375)/binHz+1;
 218
 219           if(lo_bin<0)lo_bin=0;
 220           if(lo_bin>n)lo_bin=n;
 221           if(lo_bin<l)l=lo_bin;
 222           if(hi_bin<0)hi_bin=0;
 223           if(hi_bin>n)hi_bin=n;
 224
 225           for(;l<hi_bin && l<n;l++)
 226             if(brute_buffer[l]>workc[k][m][j])
 227               brute_buffer[l]=workc[k][m][j];
 228         }
 229
 230         for(;l<n;l++)
 231           if(brute_buffer[l]>workc[k][m][EHMER_MAX-1])
 232             brute_buffer[l]=workc[k][m][EHMER_MAX-1];
 233
 234       }
 235
 236
 237       for(j=0;j<EHMER_MAX;j++){
 238         int bin=fromOC(j*.125+i*.5-2.)/binHz;
 239         if(bin<0){
 240           ret[i][m][j+2]=-999.;
 241         }else{
 242           if(bin>=n){
 243             ret[i][m][j+2]=-999.;
 244           }else{
 245             ret[i][m][j+2]=brute_buffer[bin];
 246           }
 247         }
 248       }
 249
 250       /* add fenceposts */
 251       for(j=0;j<EHMER_OFFSET;j++)
 252         if(ret[i][m][j+2]>-200.f)break;
 253       ret[i][m][0]=j;
 254
 255       for(j=EHMER_MAX-1;j>EHMER_OFFSET+1;j--)
 256         if(ret[i][m][j+2]>-200.f)
 257           break;
 258       ret[i][m][1]=j;
 259
 260     }
 261   }
 262
 263   return(ret);
 264 }
 265
 266 void _vp_psy_init(vorbis_look_psy *p,vorbis_info_psy *vi,
 267                   vorbis_info_psy_global *gi,int n,long rate){
 268   long i,j,lo=-99,hi=1;
 269   long maxoc;
 270   memset(p,0,sizeof(*p));
 271
 272   p->eighth_octave_lines=gi->eighth_octave_lines;
 273   p->shiftoc=rint(log(gi->eighth_octave_lines*8.f)/log(2.f))-1;
 274
 275   p->firstoc=toOC(.25f*rate*.5/n)*(1<<(p->shiftoc+1))-gi->eighth_octave_lines;
 276   maxoc=toOC((n+.25f)*rate*.5/n)*(1<<(p->shiftoc+1))+.5f;
 277   p->total_octave_lines=maxoc-p->firstoc+1;
 278   p->ath=_ogg_malloc(n*sizeof(*p->ath));
 279
 280   p->octave=_ogg_malloc(n*sizeof(*p->octave));
 281   p->bark=_ogg_malloc(n*sizeof(*p->bark));
 282   p->vi=vi;
 283   p->n=n;
 284   p->rate=rate;
 285
 286   /* AoTuV HF weighting */
 287   p->m_val = 1.;
 288   if(rate < 26000) p->m_val = 0;
 289   else if(rate < 38000) p->m_val = .94;   /* 32kHz */
 290   else if(rate > 46000) p->m_val = 1.275; /* 48kHz */
 291
 292   /* set up the lookups for a given blocksize and sample rate */
 293
 294   for(i=0,j=0;i<MAX_ATH-1;i++){
 295     int endpos=rint(fromOC((i+1)*.125-2.)*2*n/rate);
 296     float base=ATH[i];
 297     if(j<endpos){
 298       float delta=(ATH[i+1]-base)/(endpos-j);
 299       for(;j<endpos && j<n;j++){
 300         p->ath[j]=base+100.;
 301         base+=delta;
 302       }
 303     }
 304   }
 305
 306   for(;j<n;j++){
 307     p->ath[j]=p->ath[j-1];
 308   }
 309
 310   for(i=0;i<n;i++){
 311     float bark=toBARK(rate/(2*n)*i);
 312
 313     for(;lo+vi->noisewindowlomin<i &&
 314           toBARK(rate/(2*n)*lo)<(bark-vi->noisewindowlo);lo++);
 315
 316     for(;hi<=n && (hi<i+vi->noisewindowhimin ||
 317           toBARK(rate/(2*n)*hi)<(bark+vi->noisewindowhi));hi++);
 318
 319     p->bark[i]=((lo-1)<<16)+(hi-1);
 320
 321   }
 322
 323   for(i=0;i<n;i++)
 324     p->octave[i]=toOC((i+.25f)*.5*rate/n)*(1<<(p->shiftoc+1))+.5f;
 325
 326   p->tonecurves=setup_tone_curves(vi->toneatt,rate*.5/n,n,
 327                                   vi->tone_centerboost,vi->tone_decay);
 328
 329   /* set up rolling noise median */
 330   p->noiseoffset=_ogg_malloc(P_NOISECURVES*sizeof(*p->noiseoffset));
 331   for(i=0;i<P_NOISECURVES;i++)
 332     p->noiseoffset[i]=_ogg_malloc(n*sizeof(**p->noiseoffset));
 333
 334   for(i=0;i<n;i++){
 335     float halfoc=toOC((i+.5)*rate/(2.*n))*2.;
 336     int inthalfoc;
 337     float del;
 338
 339     if(halfoc<0)halfoc=0;
 340     if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1;
 341     inthalfoc=(int)halfoc;
 342     del=halfoc-inthalfoc;
 343
 344     for(j=0;j<P_NOISECURVES;j++)
 345       p->noiseoffset[j][i]=
 346         p->vi->noiseoff[j][inthalfoc]*(1.-del) +
 347         p->vi->noiseoff[j][inthalfoc+1]*del;
 348
 349   }
 350 #if 0
 351   {
 352     static int ls=0;
 353     _analysis_output_always("noiseoff0",ls,p->noiseoffset[0],n,1,0,0);
 354     _analysis_output_always("noiseoff1",ls,p->noiseoffset[1],n,1,0,0);
 355     _analysis_output_always("noiseoff2",ls++,p->noiseoffset[2],n,1,0,0);
 356   }
 357 #endif
 358 }
 359
 360 void _vp_psy_clear(vorbis_look_psy *p){
 361   int i,j;
 362   if(p){
 363     if(p->ath)_ogg_free(p->ath);
 364     if(p->octave)_ogg_free(p->octave);
 365     if(p->bark)_ogg_free(p->bark);
 366     if(p->tonecurves){
 367       for(i=0;i<P_BANDS;i++){
 368         for(j=0;j<P_LEVELS;j++){
 369           _ogg_free(p->tonecurves[i][j]);
 370         }
 371         _ogg_free(p->tonecurves[i]);
 372       }
 373       _ogg_free(p->tonecurves);
 374     }
 375     if(p->noiseoffset){
 376       for(i=0;i<P_NOISECURVES;i++){
 377         _ogg_free(p->noiseoffset[i]);
 378       }
 379       _ogg_free(p->noiseoffset);
 380     }
 381     memset(p,0,sizeof(*p));
 382   }
 383 }
 384
 385 /* octave/(8*eighth_octave_lines) x scale and dB y scale */
 386 static void seed_curve(float *seed,
 387                        const float **curves,
 388                        float amp,
 389                        int oc, int n,
 390                        int linesper,float dBoffset){
 391   int i,post1;
 392   int seedptr;
 393   const float *posts,*curve;
 394
 395   int choice=(int)((amp+dBoffset-P_LEVEL_0)*.1f);
 396   choice=max(choice,0);
 397   choice=min(choice,P_LEVELS-1);
 398   posts=curves[choice];
 399   curve=posts+2;
 400   post1=(int)posts[1];
 401   seedptr=oc+(posts[0]-EHMER_OFFSET)*linesper-(linesper>>1);
 402
 403   for(i=posts[0];i<post1;i++){
 404     if(seedptr>0){
 405       float lin=amp+curve[i];
 406       if(seed[seedptr]<lin)seed[seedptr]=lin;
 407     }
 408     seedptr+=linesper;
 409     if(seedptr>=n)break;
 410   }
 411 }
 412
 413 static void seed_loop(vorbis_look_psy *p,
 414                       const float ***curves,
 415                       const float *f,
 416                       const float *flr,
 417                       float *seed,
 418                       float specmax){
 419   vorbis_info_psy *vi=p->vi;
 420   long n=p->n,i;
 421   float dBoffset=vi->max_curve_dB-specmax;
 422
 423   /* prime the working vector with peak values */
 424
 425   for(i=0;i<n;i++){
 426     float max=f[i];
 427     long oc=p->octave[i];
 428     while(i+1<n && p->octave[i+1]==oc){
 429       i++;
 430       if(f[i]>max)max=f[i];
 431     }
 432
 433     if(max+6.f>flr[i]){
 434       oc=oc>>p->shiftoc;
 435
 436       if(oc>=P_BANDS)oc=P_BANDS-1;
 437       if(oc<0)oc=0;
 438
 439       seed_curve(seed,
 440                  curves[oc],
 441                  max,
 442                  p->octave[i]-p->firstoc,
 443                  p->total_octave_lines,
 444                  p->eighth_octave_lines,
 445                  dBoffset);
 446     }
 447   }
 448 }
 449
 450 static void seed_chase(float *seeds, int linesper, long n){
 451   long  *posstack=alloca(n*sizeof(*posstack));
 452   float *ampstack=alloca(n*sizeof(*ampstack));
 453   long   stack=0;
 454   long   pos=0;
 455   long   i;
 456
 457   for(i=0;i<n;i++){
 458     if(stack<2){
 459       posstack[stack]=i;
 460       ampstack[stack++]=seeds[i];
 461     }else{
 462       while(1){
 463         if(seeds[i]<ampstack[stack-1]){
 464           posstack[stack]=i;
 465           ampstack[stack++]=seeds[i];
 466           break;
 467         }else{
 468           if(i<posstack[stack-1]+linesper){
 469             if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] &&
 470                i<posstack[stack-2]+linesper){
 471               /* we completely overlap, making stack-1 irrelevant.  pop it */
 472               stack--;
 473               continue;
 474             }
 475           }
 476           posstack[stack]=i;
 477           ampstack[stack++]=seeds[i];
 478           break;
 479
 480         }
 481       }
 482     }
 483   }
 484
 485   /* the stack now contains only the positions that are relevant. Scan
 486      'em straight through */
 487
 488   for(i=0;i<stack;i++){
 489     long endpos;
 490     if(i<stack-1 && ampstack[i+1]>ampstack[i]){
 491       endpos=posstack[i+1];
 492     }else{
 493       endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is
 494                                         discarded in short frames */
 495     }
 496     if(endpos>n)endpos=n;
 497     for(;pos<endpos;pos++)
 498       seeds[pos]=ampstack[i];
 499   }
 500
 501   /* there.  Linear time.  I now remember this was on a problem set I
 502      had in Grad Skool... I didn't solve it at the time ;-) */
 503
 504 }
 505
 506 /* bleaugh, this is more complicated than it needs to be */
 507 #include<stdio.h>
 508 static void max_seeds(vorbis_look_psy *p,
 509                       float *seed,
 510                       float *flr){
 511   long   n=p->total_octave_lines;
 512   int    linesper=p->eighth_octave_lines;
 513   long   linpos=0;
 514   long   pos;
 515
 516   seed_chase(seed,linesper,n); /* for masking */
 517
 518   pos=p->octave[0]-p->firstoc-(linesper>>1);
 519
 520   while(linpos+1<p->n){
 521     float minV=seed[pos];
 522     long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc;
 523     if(minV>p->vi->tone_abs_limit)minV=p->vi->tone_abs_limit;
 524     while(pos+1<=end){
 525       pos++;
 526       if((seed[pos]>NEGINF && seed[pos]<minV) || minV==NEGINF)
 527         minV=seed[pos];
 528     }
 529
 530     end=pos+p->firstoc;
 531     for(;linpos<p->n && p->octave[linpos]<=end;linpos++)
 532       if(flr[linpos]<minV)flr[linpos]=minV;
 533   }
 534
 535   {
 536     float minV=seed[p->total_octave_lines-1];
 537     for(;linpos<p->n;linpos++)
 538       if(flr[linpos]<minV)flr[linpos]=minV;
 539   }
 540
 541 }
 542
 543 static void bark_noise_hybridmp(int n,const long *b,
 544                                 const float *f,
 545                                 float *noise,
 546                                 const float offset,
 547                                 const int fixed){
 548
 549   float *N=alloca(n*sizeof(*N));
 550   float *X=alloca(n*sizeof(*N));
 551   float *XX=alloca(n*sizeof(*N));
 552   float *Y=alloca(n*sizeof(*N));
 553   float *XY=alloca(n*sizeof(*N));
 554
 555   float tN, tX, tXX, tY, tXY;
 556   int i;
 557
 558   int lo, hi;
 559   float R=0.f;
 560   float A=0.f;
 561   float B=0.f;
 562   float D=1.f;
 563   float w, x, y;
 564
 565   tN = tX = tXX = tY = tXY = 0.f;
 566
 567   y = f[0] + offset;
 568   if (y < 1.f) y = 1.f;
 569
 570   w = y * y * .5;
 571
 572   tN += w;
 573   tX += w;
 574   tY += w * y;
 575
 576   N[0] = tN;
 577   X[0] = tX;
 578   XX[0] = tXX;
 579   Y[0] = tY;
 580   XY[0] = tXY;
 581
 582   for (i = 1, x = 1.f; i < n; i++, x += 1.f) {
 583
 584     y = f[i] + offset;
 585     if (y < 1.f) y = 1.f;
 586
 587     w = y * y;
 588
 589     tN += w;
 590     tX += w * x;
 591     tXX += w * x * x;
 592     tY += w * y;
 593     tXY += w * x * y;
 594
 595     N[i] = tN;
 596     X[i] = tX;
 597     XX[i] = tXX;
 598     Y[i] = tY;
 599     XY[i] = tXY;
 600   }
 601
 602   for (i = 0, x = 0.f; i < n; i++, x += 1.f) {
 603
 604     lo = b[i] >> 16;
 605     hi = b[i] & 0xffff;
 606     if( lo>=0 || -lo>=n ) break;
 607     if( hi>=n ) break;
 608
 609     tN = N[hi] + N[-lo];
 610     tX = X[hi] - X[-lo];
 611     tXX = XX[hi] + XX[-lo];
 612     tY = Y[hi] + Y[-lo];
 613     tXY = XY[hi] - XY[-lo];
 614
 615     A = tY * tXX - tX * tXY;
 616     B = tN * tXY - tX * tY;
 617     D = tN * tXX - tX * tX;
 618     R = (A + x * B) / D;
 619     if (R < 0.f) R = 0.f;
 620
 621     noise[i] = R - offset;
 622   }
 623
 624   for ( ; i < n; i++, x += 1.f) {
 625
 626     lo = b[i] >> 16;
 627     hi = b[i] & 0xffff;
 628     if( lo<0 || lo>=n ) break;
 629     if( hi>=n ) break;
 630
 631     tN = N[hi] - N[lo];
 632     tX = X[hi] - X[lo];
 633     tXX = XX[hi] - XX[lo];
 634     tY = Y[hi] - Y[lo];
 635     tXY = XY[hi] - XY[lo];
 636
 637     A = tY * tXX - tX * tXY;
 638     B = tN * tXY - tX * tY;
 639     D = tN * tXX - tX * tX;
 640     R = (A + x * B) / D;
 641     if (R < 0.f) R = 0.f;
 642
 643     noise[i] = R - offset;
 644   }
 645
 646   for ( ; i < n; i++, x += 1.f) {
 647
 648     R = (A + x * B) / D;
 649     if (R < 0.f) R = 0.f;
 650
 651     noise[i] = R - offset;
 652   }
 653
 654   if (fixed <= 0) return;
 655
 656   for (i = 0, x = 0.f; i < n; i++, x += 1.f) {
 657     hi = i + fixed / 2;
 658     lo = hi - fixed;
 659     if ( hi>=n ) break;
 660     if ( lo>=0 ) break;
 661
 662     tN = N[hi] + N[-lo];
 663     tX = X[hi] - X[-lo];
 664     tXX = XX[hi] + XX[-lo];
 665     tY = Y[hi] + Y[-lo];
 666     tXY = XY[hi] - XY[-lo];
 667
 668
 669     A = tY * tXX - tX * tXY;
 670     B = tN * tXY - tX * tY;
 671     D = tN * tXX - tX * tX;
 672     R = (A + x * B) / D;
 673
 674     if (R - offset < noise[i]) noise[i] = R - offset;
 675   }
 676   for ( ; i < n; i++, x += 1.f) {
 677
 678     hi = i + fixed / 2;
 679     lo = hi - fixed;
 680     if ( hi>=n ) break;
 681     if ( lo<0 ) break;
 682
 683     tN = N[hi] - N[lo];
 684     tX = X[hi] - X[lo];
 685     tXX = XX[hi] - XX[lo];
 686     tY = Y[hi] - Y[lo];
 687     tXY = XY[hi] - XY[lo];
 688
 689     A = tY * tXX - tX * tXY;
 690     B = tN * tXY - tX * tY;
 691     D = tN * tXX - tX * tX;
 692     R = (A + x * B) / D;
 693
 694     if (R - offset < noise[i]) noise[i] = R - offset;
 695   }
 696   for ( ; i < n; i++, x += 1.f) {
 697     R = (A + x * B) / D;
 698     if (R - offset < noise[i]) noise[i] = R - offset;
 699   }
 700 }
 701
 702 void _vp_noisemask(vorbis_look_psy *p,
 703                    float *logmdct,
 704                    float *logmask){
 705
 706   int i,n=p->n;
 707   float *work=alloca(n*sizeof(*work));
 708
 709   bark_noise_hybridmp(n,p->bark,logmdct,logmask,
 710                       140.,-1);
 711
 712   for(i=0;i<n;i++)work[i]=logmdct[i]-logmask[i];
 713
 714   bark_noise_hybridmp(n,p->bark,work,logmask,0.,
 715                       p->vi->noisewindowfixed);
 716
 717   for(i=0;i<n;i++)work[i]=logmdct[i]-work[i];
 718
 719 #if 0
 720   {
 721     static int seq=0;
 722
 723     float work2[n];
 724     for(i=0;i<n;i++){
 725       work2[i]=logmask[i]+work[i];
 726     }
 727
 728     if(seq&1)
 729       _analysis_output("median2R",seq/2,work,n,1,0,0);
 730     else
 731       _analysis_output("median2L",seq/2,work,n,1,0,0);
 732
 733     if(seq&1)
 734       _analysis_output("envelope2R",seq/2,work2,n,1,0,0);
 735     else
 736       _analysis_output("envelope2L",seq/2,work2,n,1,0,0);
 737     seq++;
 738   }
 739 #endif
 740
 741   for(i=0;i<n;i++){
 742     int dB=logmask[i]+.5;
 743     if(dB>=NOISE_COMPAND_LEVELS)dB=NOISE_COMPAND_LEVELS-1;
 744     if(dB<0)dB=0;
 745     logmask[i]= work[i]+p->vi->noisecompand[dB];
 746   }
 747
 748 }
 749
 750 void _vp_tonemask(vorbis_look_psy *p,
 751                   float *logfft,
 752                   float *logmask,
 753                   float global_specmax,
 754                   float local_specmax){
 755
 756   int i,n=p->n;
 757
 758   float *seed=alloca(sizeof(*seed)*p->total_octave_lines);
 759   float att=local_specmax+p->vi->ath_adjatt;
 760   for(i=0;i<p->total_octave_lines;i++)seed[i]=NEGINF;
 761
 762   /* set the ATH (floating below localmax, not global max by a
 763      specified att) */
 764   if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt;
 765
 766   for(i=0;i<n;i++)
 767     logmask[i]=p->ath[i]+att;
 768
 769   /* tone masking */
 770   seed_loop(p,(const float ***)p->tonecurves,logfft,logmask,seed,global_specmax);
 771   max_seeds(p,seed,logmask);
 772
 773 }
 774
 775 void _vp_offset_and_mix(vorbis_look_psy *p,
 776                         float *noise,
 777                         float *tone,
 778                         int offset_select,
 779                         float *logmask,
 780                         float *mdct,
 781                         float *logmdct){
 782   int i,n=p->n;
 783   float de, coeffi, cx;/* AoTuV */
 784   float toneatt=p->vi->tone_masteratt[offset_select];
 785
 786   cx = p->m_val;
 787
 788   for(i=0;i<n;i++){
 789     float val= noise[i]+p->noiseoffset[offset_select][i];
 790     if(val>p->vi->noisemaxsupp)val=p->vi->noisemaxsupp;
 791     logmask[i]=max(val,tone[i]+toneatt);
 792
 793
 794     /* AoTuV */
 795     /** @ M1 **
 796         The following codes improve a noise problem.
 797         A fundamental idea uses the value of masking and carries out
 798         the relative compensation of the MDCT.
 799         However, this code is not perfect and all noise problems cannot be solved.
 800         by Aoyumi @ 2004/04/18
 801     */
 802
 803     if(offset_select == 1) {
 804       coeffi = -17.2;       /* coeffi is a -17.2dB threshold */
 805       val = val - logmdct[i];  /* val == mdct line value relative to floor in dB */
 806
 807       if(val > coeffi){
 808         /* mdct value is > -17.2 dB below floor */
 809
 810         de = 1.0-((val-coeffi)*0.005*cx);
 811         /* pro-rated attenuation:
 812            -0.00 dB boost if mdct value is -17.2dB (relative to floor)
 813            -0.77 dB boost if mdct value is 0dB (relative to floor)
 814            -1.64 dB boost if mdct value is +17.2dB (relative to floor)
 815            etc... */
 816
 817         if(de < 0) de = 0.0001;
 818       }else
 819         /* mdct value is <= -17.2 dB below floor */
 820
 821         de = 1.0-((val-coeffi)*0.0003*cx);
 822       /* pro-rated attenuation:
 823          +0.00 dB atten if mdct value is -17.2dB (relative to floor)
 824          +0.45 dB atten if mdct value is -34.4dB (relative to floor)
 825          etc... */
 826
 827       mdct[i] *= de;
 828
 829     }
 830   }
 831 }
 832
 833 float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){
 834   vorbis_info *vi=vd->vi;
 835   codec_setup_info *ci=vi->codec_setup;
 836   vorbis_info_psy_global *gi=&ci->psy_g_param;
 837
 838   int n=ci->blocksizes[vd->W]/2;
 839   float secs=(float)n/vi->rate;
 840
 841   amp+=secs*gi->ampmax_att_per_sec;
 842   if(amp<-9999)amp=-9999;
 843   return(amp);
 844 }
 845
 846 static float FLOOR1_fromdB_LOOKUP[256]={
 847   1.0649863e-07F, 1.1341951e-07F, 1.2079015e-07F, 1.2863978e-07F,
 848   1.3699951e-07F, 1.4590251e-07F, 1.5538408e-07F, 1.6548181e-07F,
 849   1.7623575e-07F, 1.8768855e-07F, 1.9988561e-07F, 2.128753e-07F,
 850   2.2670913e-07F, 2.4144197e-07F, 2.5713223e-07F, 2.7384213e-07F,
 851   2.9163793e-07F, 3.1059021e-07F, 3.3077411e-07F, 3.5226968e-07F,
 852   3.7516214e-07F, 3.9954229e-07F, 4.2550680e-07F, 4.5315863e-07F,
 853   4.8260743e-07F, 5.1396998e-07F, 5.4737065e-07F, 5.8294187e-07F,
 854   6.2082472e-07F, 6.6116941e-07F, 7.0413592e-07F, 7.4989464e-07F,
 855   7.9862701e-07F, 8.5052630e-07F, 9.0579828e-07F, 9.6466216e-07F,
 856   1.0273513e-06F, 1.0941144e-06F, 1.1652161e-06F, 1.2409384e-06F,
 857   1.3215816e-06F, 1.4074654e-06F, 1.4989305e-06F, 1.5963394e-06F,
 858   1.7000785e-06F, 1.8105592e-06F, 1.9282195e-06F, 2.0535261e-06F,
 859   2.1869758e-06F, 2.3290978e-06F, 2.4804557e-06F, 2.6416497e-06F,
 860   2.8133190e-06F, 2.9961443e-06F, 3.1908506e-06F, 3.3982101e-06F,
 861   3.6190449e-06F, 3.8542308e-06F, 4.1047004e-06F, 4.3714470e-06F,
 862   4.6555282e-06F, 4.9580707e-06F, 5.2802740e-06F, 5.6234160e-06F,
 863   5.9888572e-06F, 6.3780469e-06F, 6.7925283e-06F, 7.2339451e-06F,
 864   7.7040476e-06F, 8.2047000e-06F, 8.7378876e-06F, 9.3057248e-06F,
 865   9.9104632e-06F, 1.0554501e-05F, 1.1240392e-05F, 1.1970856e-05F,
 866   1.2748789e-05F, 1.3577278e-05F, 1.4459606e-05F, 1.5399272e-05F,
 867   1.6400004e-05F, 1.7465768e-05F, 1.8600792e-05F, 1.9809576e-05F,
 868   2.1096914e-05F, 2.2467911e-05F, 2.3928002e-05F, 2.5482978e-05F,
 869   2.7139006e-05F, 2.8902651e-05F, 3.0780908e-05F, 3.2781225e-05F,
 870   3.4911534e-05F, 3.7180282e-05F, 3.9596466e-05F, 4.2169667e-05F,
 871   4.4910090e-05F, 4.7828601e-05F, 5.0936773e-05F, 5.4246931e-05F,
 872   5.7772202e-05F, 6.1526565e-05F, 6.5524908e-05F, 6.9783085e-05F,
 873   7.4317983e-05F, 7.9147585e-05F, 8.4291040e-05F, 8.9768747e-05F,
 874   9.5602426e-05F, 0.00010181521F, 0.00010843174F, 0.00011547824F,
 875   0.00012298267F, 0.00013097477F, 0.00013948625F, 0.00014855085F,
 876   0.00015820453F, 0.00016848555F, 0.00017943469F, 0.00019109536F,
 877   0.00020351382F, 0.00021673929F, 0.00023082423F, 0.00024582449F,
 878   0.00026179955F, 0.00027881276F, 0.00029693158F, 0.00031622787F,
 879   0.00033677814F, 0.00035866388F, 0.00038197188F, 0.00040679456F,
 880   0.00043323036F, 0.00046138411F, 0.00049136745F, 0.00052329927F,
 881   0.00055730621F, 0.00059352311F, 0.00063209358F, 0.00067317058F,
 882   0.00071691700F, 0.00076350630F, 0.00081312324F, 0.00086596457F,
 883   0.00092223983F, 0.00098217216F, 0.0010459992F, 0.0011139742F,
 884   0.0011863665F, 0.0012634633F, 0.0013455702F, 0.0014330129F,
 885   0.0015261382F, 0.0016253153F, 0.0017309374F, 0.0018434235F,
 886   0.0019632195F, 0.0020908006F, 0.0022266726F, 0.0023713743F,
 887   0.0025254795F, 0.0026895994F, 0.0028643847F, 0.0030505286F,
 888   0.0032487691F, 0.0034598925F, 0.0036847358F, 0.0039241906F,
 889   0.0041792066F, 0.0044507950F, 0.0047400328F, 0.0050480668F,
 890   0.0053761186F, 0.0057254891F, 0.0060975636F, 0.0064938176F,
 891   0.0069158225F, 0.0073652516F, 0.0078438871F, 0.0083536271F,
 892   0.0088964928F, 0.009474637F, 0.010090352F, 0.010746080F,
 893   0.011444421F, 0.012188144F, 0.012980198F, 0.013823725F,
 894   0.014722068F, 0.015678791F, 0.016697687F, 0.017782797F,
 895   0.018938423F, 0.020169149F, 0.021479854F, 0.022875735F,
 896   0.024362330F, 0.025945531F, 0.027631618F, 0.029427276F,
 897   0.031339626F, 0.033376252F, 0.035545228F, 0.037855157F,
 898   0.040315199F, 0.042935108F, 0.045725273F, 0.048696758F,
 899   0.051861348F, 0.055231591F, 0.058820850F, 0.062643361F,
 900   0.066714279F, 0.071049749F, 0.075666962F, 0.080584227F,
 901   0.085821044F, 0.091398179F, 0.097337747F, 0.10366330F,
 902   0.11039993F, 0.11757434F, 0.12521498F, 0.13335215F,
 903   0.14201813F, 0.15124727F, 0.16107617F, 0.17154380F,
 904   0.18269168F, 0.19456402F, 0.20720788F, 0.22067342F,
 905   0.23501402F, 0.25028656F, 0.26655159F, 0.28387361F,
 906   0.30232132F, 0.32196786F, 0.34289114F, 0.36517414F,
 907   0.38890521F, 0.41417847F, 0.44109412F, 0.46975890F,
 908   0.50028648F, 0.53279791F, 0.56742212F, 0.60429640F,
 909   0.64356699F, 0.68538959F, 0.72993007F, 0.77736504F,
 910   0.82788260F, 0.88168307F, 0.9389798F, 1.F,
 911 };
 912
 913 /* this is for per-channel noise normalization */
 914 static int apsort(const void *a, const void *b){
 915   float f1=**(float**)a;
 916   float f2=**(float**)b;
 917   return (f1<f2)-(f1>f2);
 918 }
 919
 920 static void flag_lossless(int limit, float prepoint, float postpoint, float *mdct,
 921                          float *floor, int *flag, int i, int jn){
 922   int j;
 923   for(j=0;j<jn;j++){
 924     float point = j>=limit-i ? postpoint : prepoint;
 925     float r = fabs(mdct[j])/floor[j];
 926     if(r<point)
 927       flag[j]=0;
 928     else
 929       flag[j]=1;
 930   }
 931 }
 932
 933 /* Overload/Side effect: On input, the *q vector holds either the
 934    quantized energy (for elements with the flag set) or the absolute
 935    values of the *r vector (for elements with flag unset).  On output,
 936    *q holds the quantized energy for all elements */
 937 static float noise_normalize(vorbis_look_psy *p, int limit, float *r, float *q, float *f, int *flags, float acc, int i, int n, int *out){
 938
 939   vorbis_info_psy *vi=p->vi;
 940   float **sort = alloca(n*sizeof(*sort));
 941   int j,count=0;
 942   int start = (vi->normal_p ? vi->normal_start-i : n);
 943   if(start>n)start=n;
 944
 945   /* force classic behavior where only energy in the current band is considered */
 946   acc=0.f;
 947
 948   /* still responsible for populating *out where noise norm not in
 949      effect.  There's no need to [re]populate *q in these areas */
 950   for(j=0;j<start;j++){
 951     if(!flags || !flags[j]){ /* lossless coupling already quantized.
 952                                 Don't touch; requantizing based on
 953                                 energy would be incorrect. */
 954       float ve = q[j]/f[j];
 955       if(r[j]<0)
 956         out[j] = -rint(sqrt(ve));
 957       else
 958         out[j] = rint(sqrt(ve));
 959     }
 960   }
 961
 962   /* sort magnitudes for noise norm portion of partition */
 963   for(;j<n;j++){
 964     if(!flags || !flags[j]){ /* can't noise norm elements that have
 965                                 already been loslessly coupled; we can
 966                                 only account for their energy error */
 967       float ve = q[j]/f[j];
 968       /* Despite all the new, more capable coupling code, for now we
 969          implement noise norm as it has been up to this point. Only
 970          consider promotions to unit magnitude from 0.  In addition
 971          the only energy error counted is quantizations to zero. */
 972       /* also-- the original point code only applied noise norm at > pointlimit */
 973       if(ve<.25f && (!flags || j>=limit-i)){
 974         acc += ve;
 975         sort[count++]=q+j; /* q is fabs(r) for unflagged element */
 976       }else{
 977         /* For now: no acc adjustment for nonzero quantization.  populate *out and q as this value is final. */
 978         if(r[j]<0)
 979           out[j] = -rint(sqrt(ve));
 980         else
 981           out[j] = rint(sqrt(ve));
 982         q[j] = out[j]*out[j]*f[j];
 983       }
 984     }/* else{
 985         again, no energy adjustment for error in nonzero quant-- for now
 986         }*/
 987   }
 988
 989   if(count){
 990     /* noise norm to do */
 991     qsort(sort,count,sizeof(*sort),apsort);
 992     for(j=0;j<count;j++){
 993       int k=sort[j]-q;
 994       if(acc>=vi->normal_thresh){
 995         out[k]=unitnorm(r[k]);
 996         acc-=1.f;
 997         q[k]=f[k];
 998       }else{
 999         out[k]=0;
1000         q[k]=0.f;
1001       }
1002     }
1003   }
1004
1005   return acc;
1006 }
1007
1008 /* Noise normalization, quantization and coupling are not wholly
1009    seperable processes in depth>1 coupling. */
1010 void _vp_couple_quantize_normalize(int blobno,
1011                                    vorbis_info_psy_global *g,
1012                                    vorbis_look_psy *p,
1013                                    vorbis_info_mapping0 *vi,
1014                                    float **mdct,
1015                                    int   **iwork,
1016                                    int    *nonzero,
1017                                    int     sliding_lowpass,
1018                                    int     ch){
1019
1020   int i;
1021   int n = p->n;
1022   int partition=(p->vi->normal_p ? p->vi->normal_partition : 16);
1023   int limit = g->coupling_pointlimit[p->vi->blockflag][blobno];
1024   float prepoint=stereo_threshholds[g->coupling_prepointamp[blobno]];
1025   float postpoint=stereo_threshholds[g->coupling_postpointamp[blobno]];
1026 #if 0
1027   float de=0.1*p->m_val; /* a blend of the AoTuV M2 and M3 code here and below */
1028 #endif
1029
1030   /* mdct is our raw mdct output, floor not removed. */
1031   /* inout passes in the ifloor, passes back quantized result */
1032
1033   /* unquantized energy (negative indicates amplitude has negative sign) */
1034   float **raw = alloca(ch*sizeof(*raw));
1035
1036   /* dual pupose; quantized energy (if flag set), othersize fabs(raw) */
1037   float **quant = alloca(ch*sizeof(*quant));
1038
1039   /* floor energy */
1040   float **floor = alloca(ch*sizeof(*floor));
1041
1042   /* flags indicating raw/quantized status of elements in raw vector */
1043   int   **flag  = alloca(ch*sizeof(*flag));
1044
1045   /* non-zero flag working vector */
1046   int    *nz    = alloca(ch*sizeof(*nz));
1047
1048   /* energy surplus/defecit tracking */
1049   float  *acc   = alloca((ch+vi->coupling_steps)*sizeof(*acc));
1050
1051   /* The threshold of a stereo is changed with the size of n */
1052   if(n > 1000)
1053     postpoint=stereo_threshholds_limited[g->coupling_postpointamp[blobno]];
1054
1055   raw[0]   = alloca(ch*partition*sizeof(**raw));
1056   quant[0] = alloca(ch*partition*sizeof(**quant));
1057   floor[0] = alloca(ch*partition*sizeof(**floor));
1058   flag[0]  = alloca(ch*partition*sizeof(**flag));
1059
1060   for(i=1;i<ch;i++){
1061     raw[i]   = &raw[0][partition*i];
1062     quant[i] = &quant[0][partition*i];
1063     floor[i] = &floor[0][partition*i];
1064     flag[i]  = &flag[0][partition*i];
1065   }
1066   for(i=0;i<ch+vi->coupling_steps;i++)
1067     acc[i]=0.f;
1068
1069   for(i=0;i<n;i+=partition){
1070     int k,j,jn = partition > n-i ? n-i : partition;
1071     int step,track = 0;
1072
1073     memcpy(nz,nonzero,sizeof(*nz)*ch);
1074
1075     /* prefill */
1076     memset(flag[0],0,ch*partition*sizeof(**flag));
1077     for(k=0;k<ch;k++){
1078       int *iout = &iwork[k][i];
1079       if(nz[k]){
1080
1081         for(j=0;j<jn;j++)
1082           floor[k][j] = FLOOR1_fromdB_LOOKUP[iout[j]];
1083
1084         flag_lossless(limit,prepoint,postpoint,&mdct[k][i],floor[k],flag[k],i,jn);
1085
1086         for(j=0;j<jn;j++){
1087           quant[k][j] = raw[k][j] = mdct[k][i+j]*mdct[k][i+j];
1088           if(mdct[k][i+j]<0.f) raw[k][j]*=-1.f;
1089           floor[k][j]*=floor[k][j];
1090         }
1091
1092         acc[track]=noise_normalize(p,limit,raw[k],quant[k],floor[k],NULL,acc[track],i,jn,iout);
1093
1094       }else{
1095         for(j=0;j<jn;j++){
1096           floor[k][j] = 1e-10f;
1097           raw[k][j] = 0.f;
1098           quant[k][j] = 0.f;
1099           flag[k][j] = 0;
1100           iout[j]=0;
1101         }
1102         acc[track]=0.f;
1103       }
1104       track++;
1105     }
1106
1107     /* coupling */
1108     for(step=0;step<vi->coupling_steps;step++){
1109       int Mi = vi->coupling_mag[step];
1110       int Ai = vi->coupling_ang[step];
1111       int *iM = &iwork[Mi][i];
1112       int *iA = &iwork[Ai][i];
1113       float *reM = raw[Mi];
1114       float *reA = raw[Ai];
1115       float *qeM = quant[Mi];
1116       float *qeA = quant[Ai];
1117       float *floorM = floor[Mi];
1118       float *floorA = floor[Ai];
1119       int *fM = flag[Mi];
1120       int *fA = flag[Ai];
1121
1122       if(nz[Mi] || nz[Ai]){
1123         nz[Mi] = nz[Ai] = 1;
1124
1125         for(j=0;j<jn;j++){
1126
1127           if(j<sliding_lowpass-i){
1128             if(fM[j] || fA[j]){
1129               /* lossless coupling */
1130
1131               reM[j] = fabs(reM[j])+fabs(reA[j]);
1132               qeM[j] = qeM[j]+qeA[j];
1133               fM[j]=fA[j]=1;
1134
1135               /* couple iM/iA */
1136               {
1137                 int A = iM[j];
1138                 int B = iA[j];
1139
1140                 if(abs(A)>abs(B)){
1141                   iA[j]=(A>0?A-B:B-A);
1142                 }else{
1143                   iA[j]=(B>0?A-B:B-A);
1144                   iM[j]=B;
1145                 }
1146
1147                 /* collapse two equivalent tuples to one */
1148                 if(iA[j]>=abs(iM[j])*2){
1149                   iA[j]= -iA[j];
1150                   iM[j]= -iM[j];
1151                 }
1152
1153               }
1154
1155             }else{
1156               /* lossy (point) coupling */
1157               if(j<limit-i){
1158                 /* dipole */
1159                 reM[j] += reA[j];
1160                 qeM[j] = fabs(reM[j]);
1161               }else{
1162 #if 0
1163                 /* AoTuV */
1164                 /** @ M2 **
1165                     The boost problem by the combination of noise normalization and point stereo is eased.
1166                     However, this is a temporary patch.
1167                     by Aoyumi @ 2004/04/18
1168                 */
1169                 float derate = (1.0 - de*((float)(j-limit+i) / (float)(n-limit)));
1170                 /* elliptical */
1171                 if(reM[j]+reA[j]<0){
1172                   reM[j] = - (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate);
1173                 }else{
1174                   reM[j] =   (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate);
1175                 }
1176 #else
1177                 /* elliptical */
1178                 if(reM[j]+reA[j]<0){
1179                   reM[j] = - (qeM[j] = fabs(reM[j])+fabs(reA[j]));
1180                 }else{
1181                   reM[j] =   (qeM[j] = fabs(reM[j])+fabs(reA[j]));
1182                 }
1183 #endif
1184
1185               }
1186               reA[j]=qeA[j]=0.f;
1187               fA[j]=1;
1188               iA[j]=0;
1189             }
1190           }
1191           floorM[j]=floorA[j]=floorM[j]+floorA[j];
1192         }
1193         /* normalize the resulting mag vector */
1194         acc[track]=noise_normalize(p,limit,raw[Mi],quant[Mi],floor[Mi],flag[Mi],acc[track],i,jn,iM);
1195         track++;
1196       }
1197     }
1198   }
1199
1200   for(i=0;i<vi->coupling_steps;i++){
1201     /* make sure coupling a zero and a nonzero channel results in two
1202        nonzero channels. */
1203     if(nonzero[vi->coupling_mag[i]] ||
1204        nonzero[vi->coupling_ang[i]]){
1205       nonzero[vi->coupling_mag[i]]=1;
1206       nonzero[vi->coupling_ang[i]]=1;
1207     }
1208   }
1209 }