sys/src/cmd/audio/libvorbis/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=malloc(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   free(brute_buffer);
 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(;j<n;j++){
 308     p->ath[j]=p->ath[j-1];
 309   }
 310
 311   for(i=0;i<n;i++){
 312     float bark=toBARK(rate/(2*n)*i);
 313
 314     for(;lo+vi->noisewindowlomin<i &&
 315           toBARK(rate/(2*n)*lo)<(bark-vi->noisewindowlo);lo++);
 316
 317     for(;hi<=n && (hi<i+vi->noisewindowhimin ||
 318           toBARK(rate/(2*n)*hi)<(bark+vi->noisewindowhi));hi++);
 319
 320     p->bark[i]=((lo-1)<<16)+(hi-1);
 321
 322   }
 323
 324   for(i=0;i<n;i++)
 325     p->octave[i]=toOC((i+.25f)*.5*rate/n)*(1<<(p->shiftoc+1))+.5f;
 326
 327   p->tonecurves=setup_tone_curves(vi->toneatt,rate*.5/n,n,
 328                                   vi->tone_centerboost,vi->tone_decay);
 329
 330   /* set up rolling noise median */
 331   p->noiseoffset=_ogg_malloc(P_NOISECURVES*sizeof(*p->noiseoffset));
 332   for(i=0;i<P_NOISECURVES;i++)
 333     p->noiseoffset[i]=_ogg_malloc(n*sizeof(**p->noiseoffset));
 334
 335   for(i=0;i<n;i++){
 336     float halfoc=toOC((i+.5)*rate/(2.*n))*2.;
 337     int inthalfoc;
 338     float del;
 339
 340     if(halfoc<0)halfoc=0;
 341     if(halfoc>=P_BANDS-1)halfoc=P_BANDS-1;
 342     inthalfoc=(int)halfoc;
 343     del=halfoc-inthalfoc;
 344
 345     for(j=0;j<P_NOISECURVES;j++)
 346       p->noiseoffset[j][i]=
 347         p->vi->noiseoff[j][inthalfoc]*(1.-del) +
 348         p->vi->noiseoff[j][inthalfoc+1]*del;
 349
 350   }
 351 #if 0
 352   {
 353     static int ls=0;
 354     _analysis_output_always("noiseoff0",ls,p->noiseoffset[0],n,1,0,0);
 355     _analysis_output_always("noiseoff1",ls,p->noiseoffset[1],n,1,0,0);
 356     _analysis_output_always("noiseoff2",ls++,p->noiseoffset[2],n,1,0,0);
 357   }
 358 #endif
 359 }
 360
 361 void _vp_psy_clear(vorbis_look_psy *p){
 362   int i,j;
 363   if(p){
 364     if(p->ath)_ogg_free(p->ath);
 365     if(p->octave)_ogg_free(p->octave);
 366     if(p->bark)_ogg_free(p->bark);
 367     if(p->tonecurves){
 368       for(i=0;i<P_BANDS;i++){
 369         for(j=0;j<P_LEVELS;j++){
 370           _ogg_free(p->tonecurves[i][j]);
 371         }
 372         _ogg_free(p->tonecurves[i]);
 373       }
 374       _ogg_free(p->tonecurves);
 375     }
 376     if(p->noiseoffset){
 377       for(i=0;i<P_NOISECURVES;i++){
 378         _ogg_free(p->noiseoffset[i]);
 379       }
 380       _ogg_free(p->noiseoffset);
 381     }
 382     memset(p,0,sizeof(*p));
 383   }
 384 }
 385
 386 /* octave/(8*eighth_octave_lines) x scale and dB y scale */
 387 static void seed_curve(float *seed,
 388                        const float **curves,
 389                        float amp,
 390                        int oc, int n,
 391                        int linesper,float dBoffset){
 392   int i,post1;
 393   int seedptr;
 394   const float *posts,*curve;
 395
 396   int choice=(int)((amp+dBoffset-P_LEVEL_0)*.1f);
 397   choice=max(choice,0);
 398   choice=min(choice,P_LEVELS-1);
 399   posts=curves[choice];
 400   curve=posts+2;
 401   post1=(int)posts[1];
 402   seedptr=oc+(posts[0]-EHMER_OFFSET)*linesper-(linesper>>1);
 403
 404   for(i=posts[0];i<post1;i++){
 405     if(seedptr>0){
 406       float lin=amp+curve[i];
 407       if(seed[seedptr]<lin)seed[seedptr]=lin;
 408     }
 409     seedptr+=linesper;
 410     if(seedptr>=n)break;
 411   }
 412 }
 413
 414 static void seed_loop(vorbis_look_psy *p,
 415                       const float ***curves,
 416                       const float *f,
 417                       const float *flr,
 418                       float *seed,
 419                       float specmax){
 420   vorbis_info_psy *vi=p->vi;
 421   long n=p->n,i;
 422   float dBoffset=vi->max_curve_dB-specmax;
 423
 424   /* prime the working vector with peak values */
 425
 426   for(i=0;i<n;i++){
 427     float max=f[i];
 428     long oc=p->octave[i];
 429     while(i+1<n && p->octave[i+1]==oc){
 430       i++;
 431       if(f[i]>max)max=f[i];
 432     }
 433
 434     if(max+6.f>flr[i]){
 435       oc=oc>>p->shiftoc;
 436
 437       if(oc>=P_BANDS)oc=P_BANDS-1;
 438       if(oc<0)oc=0;
 439
 440       seed_curve(seed,
 441                  curves[oc],
 442                  max,
 443                  p->octave[i]-p->firstoc,
 444                  p->total_octave_lines,
 445                  p->eighth_octave_lines,
 446                  dBoffset);
 447     }
 448   }
 449 }
 450
 451 static void seed_chase(float *seeds, int linesper, long n){
 452   long  *posstack=malloc(n*sizeof(*posstack));
 453   float *ampstack=malloc(n*sizeof(*ampstack));
 454   long   stack=0;
 455   long   pos=0;
 456   long   i;
 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]+linesper){
 470             if(stack>1 && ampstack[stack-1]<=ampstack[stack-2] &&
 471                i<posstack[stack-2]+linesper){
 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   for(i=0;i<stack;i++){
 490     long endpos;
 491     if(i<stack-1 && ampstack[i+1]>ampstack[i]){
 492       endpos=posstack[i+1];
 493     }else{
 494       endpos=posstack[i]+linesper+1; /* +1 is important, else bin 0 is
 495                                         discarded in short frames */
 496     }
 497     if(endpos>n)endpos=n;
 498     for(;pos<endpos;pos++)
 499       seeds[pos]=ampstack[i];
 500   }
 501
 502   /* there.  Linear time.  I now remember this was on a problem set I
 503      had in Grad Skool... I didn't solve it at the time ;-) */
 504   free(posstack);
 505   free(ampstack);
 506 }
 507
 508 /* bleaugh, this is more complicated than it needs to be */
 509 #include<stdio.h>
 510 static void max_seeds(vorbis_look_psy *p,
 511                       float *seed,
 512                       float *flr){
 513   long   n=p->total_octave_lines;
 514   int    linesper=p->eighth_octave_lines;
 515   long   linpos=0;
 516   long   pos;
 517
 518   seed_chase(seed,linesper,n); /* for masking */
 519
 520   pos=p->octave[0]-p->firstoc-(linesper>>1);
 521
 522   while(linpos+1<p->n){
 523     float minV=seed[pos];
 524     long end=((p->octave[linpos]+p->octave[linpos+1])>>1)-p->firstoc;
 525     if(minV>p->vi->tone_abs_limit)minV=p->vi->tone_abs_limit;
 526     while(pos+1<=end){
 527       pos++;
 528       if((seed[pos]>NEGINF && seed[pos]<minV) || minV==NEGINF)
 529         minV=seed[pos];
 530     }
 531
 532     end=pos+p->firstoc;
 533     for(;linpos<p->n && p->octave[linpos]<=end;linpos++)
 534       if(flr[linpos]<minV)flr[linpos]=minV;
 535   }
 536
 537   {
 538     float minV=seed[p->total_octave_lines-1];
 539     for(;linpos<p->n;linpos++)
 540       if(flr[linpos]<minV)flr[linpos]=minV;
 541   }
 542
 543 }
 544
 545 static void bark_noise_hybridmp(int n,const long *b,
 546                                 const float *f,
 547                                 float *noise,
 548                                 const float offset,
 549                                 const int fixed){
 550
 551   float *N=malloc(5*n*sizeof(*N));
 552   float *X=&N[n];
 553   float *XX=&N[2*n];
 554   float *Y=&N[3*n];
 555   float *XY=&N[4*n];
 556
 557   float tN, tX, tXX, tY, tXY;
 558   int i;
 559
 560   int lo, hi;
 561   float R=0.f;
 562   float A=0.f;
 563   float B=0.f;
 564   float D=1.f;
 565   float w, x, y;
 566
 567   tN = tX = tXX = tY = tXY = 0.f;
 568
 569   y = f[0] + offset;
 570   if (y < 1.f) y = 1.f;
 571
 572   w = y * y * .5;
 573
 574   tN += w;
 575   tX += w;
 576   tY += w * y;
 577
 578   N[0] = tN;
 579   X[0] = tX;
 580   XX[0] = tXX;
 581   Y[0] = tY;
 582   XY[0] = tXY;
 583
 584   for (i = 1, x = 1.f; i < n; i++, x += 1.f) {
 585
 586     y = f[i] + offset;
 587     if (y < 1.f) y = 1.f;
 588
 589     w = y * y;
 590
 591     tN += w;
 592     tX += w * x;
 593     tXX += w * x * x;
 594     tY += w * y;
 595     tXY += w * x * y;
 596
 597     N[i] = tN;
 598     X[i] = tX;
 599     XX[i] = tXX;
 600     Y[i] = tY;
 601     XY[i] = tXY;
 602   }
 603
 604   for (i = 0, x = 0.f; i < n; i++, x += 1.f) {
 605
 606     lo = b[i] >> 16;
 607     hi = b[i] & 0xffff;
 608     if( lo>=0 || -lo>=n ) break;
 609     if( hi>=n ) break;
 610
 611     tN = N[hi] + N[-lo];
 612     tX = X[hi] - X[-lo];
 613     tXX = XX[hi] + XX[-lo];
 614     tY = Y[hi] + Y[-lo];
 615     tXY = XY[hi] - XY[-lo];
 616
 617     A = tY * tXX - tX * tXY;
 618     B = tN * tXY - tX * tY;
 619     D = tN * tXX - tX * tX;
 620     R = (A + x * B) / D;
 621     if (R < 0.f) R = 0.f;
 622
 623     noise[i] = R - offset;
 624   }
 625
 626   for ( ; i < n; i++, x += 1.f) {
 627
 628     lo = b[i] >> 16;
 629     hi = b[i] & 0xffff;
 630     if( lo<0 || lo>=n ) break;
 631     if( hi>=n ) break;
 632
 633     tN = N[hi] - N[lo];
 634     tX = X[hi] - X[lo];
 635     tXX = XX[hi] - XX[lo];
 636     tY = Y[hi] - Y[lo];
 637     tXY = XY[hi] - XY[lo];
 638
 639     A = tY * tXX - tX * tXY;
 640     B = tN * tXY - tX * tY;
 641     D = tN * tXX - tX * tX;
 642     R = (A + x * B) / D;
 643     if (R < 0.f) R = 0.f;
 644
 645     noise[i] = R - offset;
 646   }
 647
 648   for ( ; i < n; i++, x += 1.f) {
 649
 650     R = (A + x * B) / D;
 651     if (R < 0.f) R = 0.f;
 652
 653     noise[i] = R - offset;
 654   }
 655
 656   if (fixed <= 0) {
 657     free(N);
 658     return;
 659   }
 660
 661   for (i = 0, x = 0.f; i < n; i++, x += 1.f) {
 662     hi = i + fixed / 2;
 663     lo = hi - fixed;
 664     if ( hi>=n ) break;
 665     if ( lo>=0 ) break;
 666
 667     tN = N[hi] + N[-lo];
 668     tX = X[hi] - X[-lo];
 669     tXX = XX[hi] + XX[-lo];
 670     tY = Y[hi] + Y[-lo];
 671     tXY = XY[hi] - XY[-lo];
 672
 673
 674     A = tY * tXX - tX * tXY;
 675     B = tN * tXY - tX * tY;
 676     D = tN * tXX - tX * tX;
 677     R = (A + x * B) / D;
 678
 679     if (R - offset < noise[i]) noise[i] = R - offset;
 680   }
 681   for ( ; i < n; i++, x += 1.f) {
 682
 683     hi = i + fixed / 2;
 684     lo = hi - fixed;
 685     if ( hi>=n ) break;
 686     if ( lo<0 ) break;
 687
 688     tN = N[hi] - N[lo];
 689     tX = X[hi] - X[lo];
 690     tXX = XX[hi] - XX[lo];
 691     tY = Y[hi] - Y[lo];
 692     tXY = XY[hi] - XY[lo];
 693
 694     A = tY * tXX - tX * tXY;
 695     B = tN * tXY - tX * tY;
 696     D = tN * tXX - tX * tX;
 697     R = (A + x * B) / D;
 698
 699     if (R - offset < noise[i]) noise[i] = R - offset;
 700   }
 701   for ( ; i < n; i++, x += 1.f) {
 702     R = (A + x * B) / D;
 703     if (R - offset < noise[i]) noise[i] = R - offset;
 704   }
 705   free(N);
 706 }
 707
 708 void _vp_noisemask(vorbis_look_psy *p,
 709                    float *logmdct,
 710                    float *logmask){
 711
 712   int i,n=p->n;
 713   float *work=malloc(n*sizeof(*work));
 714
 715   bark_noise_hybridmp(n,p->bark,logmdct,logmask,
 716                       140.,-1);
 717
 718   for(i=0;i<n;i++)work[i]=logmdct[i]-logmask[i];
 719
 720   bark_noise_hybridmp(n,p->bark,work,logmask,0.,
 721                       p->vi->noisewindowfixed);
 722
 723   for(i=0;i<n;i++)work[i]=logmdct[i]-work[i];
 724
 725 #if 0
 726   {
 727     static int seq=0;
 728
 729     float work2[n];
 730     for(i=0;i<n;i++){
 731       work2[i]=logmask[i]+work[i];
 732     }
 733
 734     if(seq&1)
 735       _analysis_output("median2R",seq/2,work,n,1,0,0);
 736     else
 737       _analysis_output("median2L",seq/2,work,n,1,0,0);
 738
 739     if(seq&1)
 740       _analysis_output("envelope2R",seq/2,work2,n,1,0,0);
 741     else
 742       _analysis_output("envelope2L",seq/2,work2,n,1,0,0);
 743     seq++;
 744   }
 745 #endif
 746
 747   for(i=0;i<n;i++){
 748     int dB=logmask[i]+.5;
 749     if(dB>=NOISE_COMPAND_LEVELS)dB=NOISE_COMPAND_LEVELS-1;
 750     if(dB<0)dB=0;
 751     logmask[i]= work[i]+p->vi->noisecompand[dB];
 752   }
 753   free(work);
 754 }
 755
 756 void _vp_tonemask(vorbis_look_psy *p,
 757                   float *logfft,
 758                   float *logmask,
 759                   float global_specmax,
 760                   float local_specmax){
 761
 762   int i,n=p->n;
 763
 764   float *seed=malloc(sizeof(*seed)*p->total_octave_lines);
 765   float att=local_specmax+p->vi->ath_adjatt;
 766   for(i=0;i<p->total_octave_lines;i++)seed[i]=NEGINF;
 767
 768   /* set the ATH (floating below localmax, not global max by a
 769      specified att) */
 770   if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt;
 771
 772   for(i=0;i<n;i++)
 773     logmask[i]=p->ath[i]+att;
 774
 775   /* tone masking */
 776   seed_loop(p,(const float ***)p->tonecurves,logfft,logmask,seed,global_specmax);
 777   max_seeds(p,seed,logmask);
 778   free(seed);
 779 }
 780
 781 void _vp_offset_and_mix(vorbis_look_psy *p,
 782                         float *noise,
 783                         float *tone,
 784                         int offset_select,
 785                         float *logmask,
 786                         float *mdct,
 787                         float *logmdct){
 788   int i,n=p->n;
 789   float de, coeffi, cx;/* AoTuV */
 790   float toneatt=p->vi->tone_masteratt[offset_select];
 791
 792   cx = p->m_val;
 793
 794   for(i=0;i<n;i++){
 795     float val= noise[i]+p->noiseoffset[offset_select][i];
 796     if(val>p->vi->noisemaxsupp)val=p->vi->noisemaxsupp;
 797     logmask[i]=max(val,tone[i]+toneatt);
 798
 799
 800     /* AoTuV */
 801     /** @ M1 **
 802         The following codes improve a noise problem.
 803         A fundamental idea uses the value of masking and carries out
 804         the relative compensation of the MDCT.
 805         However, this code is not perfect and all noise problems cannot be solved.
 806         by Aoyumi @ 2004/04/18
 807     */
 808
 809     if(offset_select == 1) {
 810       coeffi = -17.2;       /* coeffi is a -17.2dB threshold */
 811       val = val - logmdct[i];  /* val == mdct line value relative to floor in dB */
 812
 813       if(val > coeffi){
 814         /* mdct value is > -17.2 dB below floor */
 815
 816         de = 1.0-((val-coeffi)*0.005*cx);
 817         /* pro-rated attenuation:
 818            -0.00 dB boost if mdct value is -17.2dB (relative to floor)
 819            -0.77 dB boost if mdct value is 0dB (relative to floor)
 820            -1.64 dB boost if mdct value is +17.2dB (relative to floor)
 821            etc... */
 822
 823         if(de < 0) de = 0.0001;
 824       }else
 825         /* mdct value is <= -17.2 dB below floor */
 826
 827         de = 1.0-((val-coeffi)*0.0003*cx);
 828       /* pro-rated attenuation:
 829          +0.00 dB atten if mdct value is -17.2dB (relative to floor)
 830          +0.45 dB atten if mdct value is -34.4dB (relative to floor)
 831          etc... */
 832
 833       mdct[i] *= de;
 834
 835     }
 836   }
 837 }
 838
 839 float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){
 840   vorbis_info *vi=vd->vi;
 841   codec_setup_info *ci=vi->codec_setup;
 842   vorbis_info_psy_global *gi=&ci->psy_g_param;
 843
 844   int n=ci->blocksizes[vd->W]/2;
 845   float secs=(float)n/vi->rate;
 846
 847   amp+=secs*gi->ampmax_att_per_sec;
 848   if(amp<-9999)amp=-9999;
 849   return(amp);
 850 }
 851
 852 static float FLOOR1_fromdB_LOOKUP[256]={
 853   1.0649863e-07F, 1.1341951e-07F, 1.2079015e-07F, 1.2863978e-07F,
 854   1.3699951e-07F, 1.4590251e-07F, 1.5538408e-07F, 1.6548181e-07F,
 855   1.7623575e-07F, 1.8768855e-07F, 1.9988561e-07F, 2.128753e-07F,
 856   2.2670913e-07F, 2.4144197e-07F, 2.5713223e-07F, 2.7384213e-07F,
 857   2.9163793e-07F, 3.1059021e-07F, 3.3077411e-07F, 3.5226968e-07F,
 858   3.7516214e-07F, 3.9954229e-07F, 4.2550680e-07F, 4.5315863e-07F,
 859   4.8260743e-07F, 5.1396998e-07F, 5.4737065e-07F, 5.8294187e-07F,
 860   6.2082472e-07F, 6.6116941e-07F, 7.0413592e-07F, 7.4989464e-07F,
 861   7.9862701e-07F, 8.5052630e-07F, 9.0579828e-07F, 9.6466216e-07F,
 862   1.0273513e-06F, 1.0941144e-06F, 1.1652161e-06F, 1.2409384e-06F,
 863   1.3215816e-06F, 1.4074654e-06F, 1.4989305e-06F, 1.5963394e-06F,
 864   1.7000785e-06F, 1.8105592e-06F, 1.9282195e-06F, 2.0535261e-06F,
 865   2.1869758e-06F, 2.3290978e-06F, 2.4804557e-06F, 2.6416497e-06F,
 866   2.8133190e-06F, 2.9961443e-06F, 3.1908506e-06F, 3.3982101e-06F,
 867   3.6190449e-06F, 3.8542308e-06F, 4.1047004e-06F, 4.3714470e-06F,
 868   4.6555282e-06F, 4.9580707e-06F, 5.2802740e-06F, 5.6234160e-06F,
 869   5.9888572e-06F, 6.3780469e-06F, 6.7925283e-06F, 7.2339451e-06F,
 870   7.7040476e-06F, 8.2047000e-06F, 8.7378876e-06F, 9.3057248e-06F,
 871   9.9104632e-06F, 1.0554501e-05F, 1.1240392e-05F, 1.1970856e-05F,
 872   1.2748789e-05F, 1.3577278e-05F, 1.4459606e-05F, 1.5399272e-05F,
 873   1.6400004e-05F, 1.7465768e-05F, 1.8600792e-05F, 1.9809576e-05F,
 874   2.1096914e-05F, 2.2467911e-05F, 2.3928002e-05F, 2.5482978e-05F,
 875   2.7139006e-05F, 2.8902651e-05F, 3.0780908e-05F, 3.2781225e-05F,
 876   3.4911534e-05F, 3.7180282e-05F, 3.9596466e-05F, 4.2169667e-05F,
 877   4.4910090e-05F, 4.7828601e-05F, 5.0936773e-05F, 5.4246931e-05F,
 878   5.7772202e-05F, 6.1526565e-05F, 6.5524908e-05F, 6.9783085e-05F,
 879   7.4317983e-05F, 7.9147585e-05F, 8.4291040e-05F, 8.9768747e-05F,
 880   9.5602426e-05F, 0.00010181521F, 0.00010843174F, 0.00011547824F,
 881   0.00012298267F, 0.00013097477F, 0.00013948625F, 0.00014855085F,
 882   0.00015820453F, 0.00016848555F, 0.00017943469F, 0.00019109536F,
 883   0.00020351382F, 0.00021673929F, 0.00023082423F, 0.00024582449F,
 884   0.00026179955F, 0.00027881276F, 0.00029693158F, 0.00031622787F,
 885   0.00033677814F, 0.00035866388F, 0.00038197188F, 0.00040679456F,
 886   0.00043323036F, 0.00046138411F, 0.00049136745F, 0.00052329927F,
 887   0.00055730621F, 0.00059352311F, 0.00063209358F, 0.00067317058F,
 888   0.00071691700F, 0.00076350630F, 0.00081312324F, 0.00086596457F,
 889   0.00092223983F, 0.00098217216F, 0.0010459992F, 0.0011139742F,
 890   0.0011863665F, 0.0012634633F, 0.0013455702F, 0.0014330129F,
 891   0.0015261382F, 0.0016253153F, 0.0017309374F, 0.0018434235F,
 892   0.0019632195F, 0.0020908006F, 0.0022266726F, 0.0023713743F,
 893   0.0025254795F, 0.0026895994F, 0.0028643847F, 0.0030505286F,
 894   0.0032487691F, 0.0034598925F, 0.0036847358F, 0.0039241906F,
 895   0.0041792066F, 0.0044507950F, 0.0047400328F, 0.0050480668F,
 896   0.0053761186F, 0.0057254891F, 0.0060975636F, 0.0064938176F,
 897   0.0069158225F, 0.0073652516F, 0.0078438871F, 0.0083536271F,
 898   0.0088964928F, 0.009474637F, 0.010090352F, 0.010746080F,
 899   0.011444421F, 0.012188144F, 0.012980198F, 0.013823725F,
 900   0.014722068F, 0.015678791F, 0.016697687F, 0.017782797F,
 901   0.018938423F, 0.020169149F, 0.021479854F, 0.022875735F,
 902   0.024362330F, 0.025945531F, 0.027631618F, 0.029427276F,
 903   0.031339626F, 0.033376252F, 0.035545228F, 0.037855157F,
 904   0.040315199F, 0.042935108F, 0.045725273F, 0.048696758F,
 905   0.051861348F, 0.055231591F, 0.058820850F, 0.062643361F,
 906   0.066714279F, 0.071049749F, 0.075666962F, 0.080584227F,
 907   0.085821044F, 0.091398179F, 0.097337747F, 0.10366330F,
 908   0.11039993F, 0.11757434F, 0.12521498F, 0.13335215F,
 909   0.14201813F, 0.15124727F, 0.16107617F, 0.17154380F,
 910   0.18269168F, 0.19456402F, 0.20720788F, 0.22067342F,
 911   0.23501402F, 0.25028656F, 0.26655159F, 0.28387361F,
 912   0.30232132F, 0.32196786F, 0.34289114F, 0.36517414F,
 913   0.38890521F, 0.41417847F, 0.44109412F, 0.46975890F,
 914   0.50028648F, 0.53279791F, 0.56742212F, 0.60429640F,
 915   0.64356699F, 0.68538959F, 0.72993007F, 0.77736504F,
 916   0.82788260F, 0.88168307F, 0.9389798F, 1.F,
 917 };
 918
 919 /* this is for per-channel noise normalization */
 920 static int apsort(const void *a, const void *b){
 921   float f1=**(float**)a;
 922   float f2=**(float**)b;
 923   return (f1<f2)-(f1>f2);
 924 }
 925
 926 static void flag_lossless(int limit, float prepoint, float postpoint, float *mdct,
 927                          float *floor, int *flag, int i, int jn){
 928   int j;
 929   for(j=0;j<jn;j++){
 930     float point = j>=limit-i ? postpoint : prepoint;
 931     float r = fabs(mdct[j])/floor[j];
 932     if(r<point)
 933       flag[j]=0;
 934     else
 935       flag[j]=1;
 936   }
 937 }
 938
 939 /* Overload/Side effect: On input, the *q vector holds either the
 940    quantized energy (for elements with the flag set) or the absolute
 941    values of the *r vector (for elements with flag unset).  On output,
 942    *q holds the quantized energy for all elements */
 943 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){
 944
 945   vorbis_info_psy *vi=p->vi;
 946   float **sort = malloc(n*sizeof(*sort));
 947   int j,count=0;
 948   int start = (vi->normal_p ? vi->normal_start-i : n);
 949   if(start>n)start=n;
 950
 951   /* force classic behavior where only energy in the current band is considered */
 952   acc=0.f;
 953
 954   /* still responsible for populating *out where noise norm not in
 955      effect.  There's no need to [re]populate *q in these areas */
 956   for(j=0;j<start;j++){
 957     if(!flags || !flags[j]){ /* lossless coupling already quantized.
 958                                 Don't touch; requantizing based on
 959                                 energy would be incorrect. */
 960       float ve = q[j]/f[j];
 961       if(r[j]<0)
 962         out[j] = -rint(sqrt(ve));
 963       else
 964         out[j] = rint(sqrt(ve));
 965     }
 966   }
 967
 968   /* sort magnitudes for noise norm portion of partition */
 969   for(;j<n;j++){
 970     if(!flags || !flags[j]){ /* can't noise norm elements that have
 971                                 already been loslessly coupled; we can
 972                                 only account for their energy error */
 973       float ve = q[j]/f[j];
 974       /* Despite all the new, more capable coupling code, for now we
 975          implement noise norm as it has been up to this point. Only
 976          consider promotions to unit magnitude from 0.  In addition
 977          the only energy error counted is quantizations to zero. */
 978       /* also-- the original point code only applied noise norm at > pointlimit */
 979       if(ve<.25f && (!flags || j>=limit-i)){
 980         acc += ve;
 981         sort[count++]=q+j; /* q is fabs(r) for unflagged element */
 982       }else{
 983         /* For now: no acc adjustment for nonzero quantization.  populate *out and q as this value is final. */
 984         if(r[j]<0)
 985           out[j] = -rint(sqrt(ve));
 986         else
 987           out[j] = rint(sqrt(ve));
 988         q[j] = out[j]*out[j]*f[j];
 989       }
 990     }/* else{
 991         again, no energy adjustment for error in nonzero quant-- for now
 992         }*/
 993   }
 994
 995   if(count){
 996     /* noise norm to do */
 997     qsort(sort,count,sizeof(*sort),apsort);
 998     for(j=0;j<count;j++){
 999       int k=sort[j]-q;
1000       if(acc>=vi->normal_thresh){
1001         out[k]=unitnorm(r[k]);
1002         acc-=1.f;
1003         q[k]=f[k];
1004       }else{
1005         out[k]=0;
1006         q[k]=0.f;
1007       }
1008     }
1009   }
1010
1011   free(sort);
1012   return acc;
1013 }
1014
1015 /* Noise normalization, quantization and coupling are not wholly
1016    seperable processes in depth>1 coupling. */
1017 void _vp_couple_quantize_normalize(int blobno,
1018                                    vorbis_info_psy_global *g,
1019                                    vorbis_look_psy *p,
1020                                    vorbis_info_mapping0 *vi,
1021                                    float **mdct,
1022                                    int   **iwork,
1023                                    int    *nonzero,
1024                                    int     sliding_lowpass,
1025                                    int     ch){
1026
1027   int i;
1028   int n = p->n;
1029   int partition=(p->vi->normal_p ? p->vi->normal_partition : 16);
1030   int limit = g->coupling_pointlimit[p->vi->blockflag][blobno];
1031   float prepoint=stereo_threshholds[g->coupling_prepointamp[blobno]];
1032   float postpoint=stereo_threshholds[g->coupling_postpointamp[blobno]];
1033 #if 0
1034   float de=0.1*p->m_val; /* a blend of the AoTuV M2 and M3 code here and below */
1035 #endif
1036
1037   /* mdct is our raw mdct output, floor not removed. */
1038   /* inout passes in the ifloor, passes back quantized result */
1039
1040   /* unquantized energy (negative indicates amplitude has negative sign) */
1041   float **raw = malloc(ch*sizeof(*raw));
1042
1043   /* dual pupose; quantized energy (if flag set), othersize fabs(raw) */
1044   float **quant = malloc(ch*sizeof(*quant));
1045
1046   /* floor energy */
1047   float **floor = malloc(ch*sizeof(*floor));
1048
1049   /* flags indicating raw/quantized status of elements in raw vector */
1050   int   **flag  = malloc(ch*sizeof(*flag));
1051
1052   /* non-zero flag working vector */
1053   int    *nz    = malloc(ch*sizeof(*nz));
1054
1055   /* energy surplus/defecit tracking */
1056   float  *acc   = malloc((ch+vi->coupling_steps)*sizeof(*acc));
1057
1058   /* The threshold of a stereo is changed with the size of n */
1059   if(n > 1000)
1060     postpoint=stereo_threshholds_limited[g->coupling_postpointamp[blobno]];
1061
1062   raw[0]   = malloc(ch*partition*sizeof(**raw));
1063   quant[0] = malloc(ch*partition*sizeof(**quant));
1064   floor[0] = malloc(ch*partition*sizeof(**floor));
1065   flag[0]  = malloc(ch*partition*sizeof(**flag));
1066
1067   for(i=1;i<ch;i++){
1068     raw[i]   = &raw[0][partition*i];
1069     quant[i] = &quant[0][partition*i];
1070     floor[i] = &floor[0][partition*i];
1071     flag[i]  = &flag[0][partition*i];
1072   }
1073   for(i=0;i<ch+vi->coupling_steps;i++)
1074     acc[i]=0.f;
1075
1076   for(i=0;i<n;i+=partition){
1077     int k,j,jn = partition > n-i ? n-i : partition;
1078     int step,track = 0;
1079
1080     memcpy(nz,nonzero,sizeof(*nz)*ch);
1081
1082     /* prefill */
1083     memset(flag[0],0,ch*partition*sizeof(**flag));
1084     for(k=0;k<ch;k++){
1085       int *iout = &iwork[k][i];
1086       if(nz[k]){
1087
1088         for(j=0;j<jn;j++)
1089           floor[k][j] = FLOOR1_fromdB_LOOKUP[iout[j]];
1090
1091         flag_lossless(limit,prepoint,postpoint,&mdct[k][i],floor[k],flag[k],i,jn);
1092
1093         for(j=0;j<jn;j++){
1094           quant[k][j] = raw[k][j] = mdct[k][i+j]*mdct[k][i+j];
1095           if(mdct[k][i+j]<0.f) raw[k][j]*=-1.f;
1096           floor[k][j]*=floor[k][j];
1097         }
1098
1099         acc[track]=noise_normalize(p,limit,raw[k],quant[k],floor[k],NULL,acc[track],i,jn,iout);
1100
1101       }else{
1102         for(j=0;j<jn;j++){
1103           floor[k][j] = 1e-10f;
1104           raw[k][j] = 0.f;
1105           quant[k][j] = 0.f;
1106           flag[k][j] = 0;
1107           iout[j]=0;
1108         }
1109         acc[track]=0.f;
1110       }
1111       track++;
1112     }
1113
1114     /* coupling */
1115     for(step=0;step<vi->coupling_steps;step++){
1116       int Mi = vi->coupling_mag[step];
1117       int Ai = vi->coupling_ang[step];
1118       int *iM = &iwork[Mi][i];
1119       int *iA = &iwork[Ai][i];
1120       float *reM = raw[Mi];
1121       float *reA = raw[Ai];
1122       float *qeM = quant[Mi];
1123       float *qeA = quant[Ai];
1124       float *floorM = floor[Mi];
1125       float *floorA = floor[Ai];
1126       int *fM = flag[Mi];
1127       int *fA = flag[Ai];
1128
1129       if(nz[Mi] || nz[Ai]){
1130         nz[Mi] = nz[Ai] = 1;
1131
1132         for(j=0;j<jn;j++){
1133
1134           if(j<sliding_lowpass-i){
1135             if(fM[j] || fA[j]){
1136               /* lossless coupling */
1137
1138               reM[j] = fabs(reM[j])+fabs(reA[j]);
1139               qeM[j] = qeM[j]+qeA[j];
1140               fM[j]=fA[j]=1;
1141
1142               /* couple iM/iA */
1143               {
1144                 int A = iM[j];
1145                 int B = iA[j];
1146
1147                 if(abs(A)>abs(B)){
1148                   iA[j]=(A>0?A-B:B-A);
1149                 }else{
1150                   iA[j]=(B>0?A-B:B-A);
1151                   iM[j]=B;
1152                 }
1153
1154                 /* collapse two equivalent tuples to one */
1155                 if(iA[j]>=abs(iM[j])*2){
1156                   iA[j]= -iA[j];
1157                   iM[j]= -iM[j];
1158                 }
1159
1160               }
1161
1162             }else{
1163               /* lossy (point) coupling */
1164               if(j<limit-i){
1165                 /* dipole */
1166                 reM[j] += reA[j];
1167                 qeM[j] = fabs(reM[j]);
1168               }else{
1169 #if 0
1170                 /* AoTuV */
1171                 /** @ M2 **
1172                     The boost problem by the combination of noise normalization and point stereo is eased.
1173                     However, this is a temporary patch.
1174                     by Aoyumi @ 2004/04/18
1175                 */
1176                 float derate = (1.0 - de*((float)(j-limit+i) / (float)(n-limit)));
1177                 /* elliptical */
1178                 if(reM[j]+reA[j]<0){
1179                   reM[j] = - (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate);
1180                 }else{
1181                   reM[j] =   (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate);
1182                 }
1183 #else
1184                 /* elliptical */
1185                 if(reM[j]+reA[j]<0){
1186                   reM[j] = - (qeM[j] = fabs(reM[j])+fabs(reA[j]));
1187                 }else{
1188                   reM[j] =   (qeM[j] = fabs(reM[j])+fabs(reA[j]));
1189                 }
1190 #endif
1191
1192               }
1193               reA[j]=qeA[j]=0.f;
1194               fA[j]=1;
1195               iA[j]=0;
1196             }
1197           }
1198           floorM[j]=floorA[j]=floorM[j]+floorA[j];
1199         }
1200         /* normalize the resulting mag vector */
1201         acc[track]=noise_normalize(p,limit,raw[Mi],quant[Mi],floor[Mi],flag[Mi],acc[track],i,jn,iM);
1202         track++;
1203       }
1204     }
1205   }
1206
1207   for(i=0;i<vi->coupling_steps;i++){
1208     /* make sure coupling a zero and a nonzero channel results in two
1209        nonzero channels. */
1210     if(nonzero[vi->coupling_mag[i]] ||
1211        nonzero[vi->coupling_ang[i]]){
1212       nonzero[vi->coupling_mag[i]]=1;
1213       nonzero[vi->coupling_ang[i]]=1;
1214     }
1215   }
1216   free(raw[0]);
1217   free(quant[0]);
1218   free(floor[0]);
1219   free(flag[0]);
1220   free(raw);
1221   free(quant);
1222   free(floor);
1223   free(flag);
1224   free(nz);
1225   free(acc);
1226 }