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 http://www.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++, x += 1.f) {
 603
 604     lo = b[i] >> 16;
 605     hi = b[i] & 0xffff;
 606     if( lo>=0 ) 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)
 620       R = 0.f;
 621
 622     noise[i] = R - offset;
 623   }
 624
 625   for ( ;; i++, x += 1.f) {
 626
 627     lo = b[i] >> 16;
 628     hi = b[i] & 0xffff;
 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   for ( ; i < n; i++, x += 1.f) {
 646
 647     R = (A + x * B) / D;
 648     if (R < 0.f) R = 0.f;
 649
 650     noise[i] = R - offset;
 651   }
 652
 653   if (fixed <= 0) return;
 654
 655   for (i = 0, x = 0.f;; i++, x += 1.f) {
 656     hi = i + fixed / 2;
 657     lo = hi - fixed;
 658     if(lo>=0)break;
 659
 660     tN = N[hi] + N[-lo];
 661     tX = X[hi] - X[-lo];
 662     tXX = XX[hi] + XX[-lo];
 663     tY = Y[hi] + Y[-lo];
 664     tXY = XY[hi] - XY[-lo];
 665
 666
 667     A = tY * tXX - tX * tXY;
 668     B = tN * tXY - tX * tY;
 669     D = tN * tXX - tX * tX;
 670     R = (A + x * B) / D;
 671
 672     if (R - offset < noise[i]) noise[i] = R - offset;
 673   }
 674   for ( ;; i++, x += 1.f) {
 675
 676     hi = i + fixed / 2;
 677     lo = hi - fixed;
 678     if(hi>=n)break;
 679
 680     tN = N[hi] - N[lo];
 681     tX = X[hi] - X[lo];
 682     tXX = XX[hi] - XX[lo];
 683     tY = Y[hi] - Y[lo];
 684     tXY = XY[hi] - XY[lo];
 685
 686     A = tY * tXX - tX * tXY;
 687     B = tN * tXY - tX * tY;
 688     D = tN * tXX - tX * tX;
 689     R = (A + x * B) / D;
 690
 691     if (R - offset < noise[i]) noise[i] = R - offset;
 692   }
 693   for ( ; i < n; i++, x += 1.f) {
 694     R = (A + x * B) / D;
 695     if (R - offset < noise[i]) noise[i] = R - offset;
 696   }
 697 }
 698
 699 void _vp_noisemask(vorbis_look_psy *p,
 700                    float *logmdct,
 701                    float *logmask){
 702
 703   int i,n=p->n;
 704   float *work=alloca(n*sizeof(*work));
 705
 706   bark_noise_hybridmp(n,p->bark,logmdct,logmask,
 707                       140.,-1);
 708
 709   for(i=0;i<n;i++)work[i]=logmdct[i]-logmask[i];
 710
 711   bark_noise_hybridmp(n,p->bark,work,logmask,0.,
 712                       p->vi->noisewindowfixed);
 713
 714   for(i=0;i<n;i++)work[i]=logmdct[i]-work[i];
 715
 716 #if 0
 717   {
 718     static int seq=0;
 719
 720     float work2[n];
 721     for(i=0;i<n;i++){
 722       work2[i]=logmask[i]+work[i];
 723     }
 724
 725     if(seq&1)
 726       _analysis_output("median2R",seq/2,work,n,1,0,0);
 727     else
 728       _analysis_output("median2L",seq/2,work,n,1,0,0);
 729
 730     if(seq&1)
 731       _analysis_output("envelope2R",seq/2,work2,n,1,0,0);
 732     else
 733       _analysis_output("envelope2L",seq/2,work2,n,1,0,0);
 734     seq++;
 735   }
 736 #endif
 737
 738   for(i=0;i<n;i++){
 739     int dB=logmask[i]+.5;
 740     if(dB>=NOISE_COMPAND_LEVELS)dB=NOISE_COMPAND_LEVELS-1;
 741     if(dB<0)dB=0;
 742     logmask[i]= work[i]+p->vi->noisecompand[dB];
 743   }
 744
 745 }
 746
 747 void _vp_tonemask(vorbis_look_psy *p,
 748                   float *logfft,
 749                   float *logmask,
 750                   float global_specmax,
 751                   float local_specmax){
 752
 753   int i,n=p->n;
 754
 755   float *seed=alloca(sizeof(*seed)*p->total_octave_lines);
 756   float att=local_specmax+p->vi->ath_adjatt;
 757   for(i=0;i<p->total_octave_lines;i++)seed[i]=NEGINF;
 758
 759   /* set the ATH (floating below localmax, not global max by a
 760      specified att) */
 761   if(att<p->vi->ath_maxatt)att=p->vi->ath_maxatt;
 762
 763   for(i=0;i<n;i++)
 764     logmask[i]=p->ath[i]+att;
 765
 766   /* tone masking */
 767   seed_loop(p,(const float ***)p->tonecurves,logfft,logmask,seed,global_specmax);
 768   max_seeds(p,seed,logmask);
 769
 770 }
 771
 772 void _vp_offset_and_mix(vorbis_look_psy *p,
 773                         float *noise,
 774                         float *tone,
 775                         int offset_select,
 776                         float *logmask,
 777                         float *mdct,
 778                         float *logmdct){
 779   int i,n=p->n;
 780   float de, coeffi, cx;/* AoTuV */
 781   float toneatt=p->vi->tone_masteratt[offset_select];
 782
 783   cx = p->m_val;
 784
 785   for(i=0;i<n;i++){
 786     float val= noise[i]+p->noiseoffset[offset_select][i];
 787     if(val>p->vi->noisemaxsupp)val=p->vi->noisemaxsupp;
 788     logmask[i]=max(val,tone[i]+toneatt);
 789
 790
 791     /* AoTuV */
 792     /** @ M1 **
 793         The following codes improve a noise problem.
 794         A fundamental idea uses the value of masking and carries out
 795         the relative compensation of the MDCT.
 796         However, this code is not perfect and all noise problems cannot be solved.
 797         by Aoyumi @ 2004/04/18
 798     */
 799
 800     if(offset_select == 1) {
 801       coeffi = -17.2;       /* coeffi is a -17.2dB threshold */
 802       val = val - logmdct[i];  /* val == mdct line value relative to floor in dB */
 803
 804       if(val > coeffi){
 805         /* mdct value is > -17.2 dB below floor */
 806
 807         de = 1.0-((val-coeffi)*0.005*cx);
 808         /* pro-rated attenuation:
 809            -0.00 dB boost if mdct value is -17.2dB (relative to floor)
 810            -0.77 dB boost if mdct value is 0dB (relative to floor)
 811            -1.64 dB boost if mdct value is +17.2dB (relative to floor)
 812            etc... */
 813
 814         if(de < 0) de = 0.0001;
 815       }else
 816         /* mdct value is <= -17.2 dB below floor */
 817
 818         de = 1.0-((val-coeffi)*0.0003*cx);
 819       /* pro-rated attenuation:
 820          +0.00 dB atten if mdct value is -17.2dB (relative to floor)
 821          +0.45 dB atten if mdct value is -34.4dB (relative to floor)
 822          etc... */
 823
 824       mdct[i] *= de;
 825
 826     }
 827   }
 828 }
 829
 830 float _vp_ampmax_decay(float amp,vorbis_dsp_state *vd){
 831   vorbis_info *vi=vd->vi;
 832   codec_setup_info *ci=vi->codec_setup;
 833   vorbis_info_psy_global *gi=&ci->psy_g_param;
 834
 835   int n=ci->blocksizes[vd->W]/2;
 836   float secs=(float)n/vi->rate;
 837
 838   amp+=secs*gi->ampmax_att_per_sec;
 839   if(amp<-9999)amp=-9999;
 840   return(amp);
 841 }
 842
 843 static float FLOOR1_fromdB_LOOKUP[256]={
 844   1.0649863e-07F, 1.1341951e-07F, 1.2079015e-07F, 1.2863978e-07F,
 845   1.3699951e-07F, 1.4590251e-07F, 1.5538408e-07F, 1.6548181e-07F,
 846   1.7623575e-07F, 1.8768855e-07F, 1.9988561e-07F, 2.128753e-07F,
 847   2.2670913e-07F, 2.4144197e-07F, 2.5713223e-07F, 2.7384213e-07F,
 848   2.9163793e-07F, 3.1059021e-07F, 3.3077411e-07F, 3.5226968e-07F,
 849   3.7516214e-07F, 3.9954229e-07F, 4.2550680e-07F, 4.5315863e-07F,
 850   4.8260743e-07F, 5.1396998e-07F, 5.4737065e-07F, 5.8294187e-07F,
 851   6.2082472e-07F, 6.6116941e-07F, 7.0413592e-07F, 7.4989464e-07F,
 852   7.9862701e-07F, 8.5052630e-07F, 9.0579828e-07F, 9.6466216e-07F,
 853   1.0273513e-06F, 1.0941144e-06F, 1.1652161e-06F, 1.2409384e-06F,
 854   1.3215816e-06F, 1.4074654e-06F, 1.4989305e-06F, 1.5963394e-06F,
 855   1.7000785e-06F, 1.8105592e-06F, 1.9282195e-06F, 2.0535261e-06F,
 856   2.1869758e-06F, 2.3290978e-06F, 2.4804557e-06F, 2.6416497e-06F,
 857   2.8133190e-06F, 2.9961443e-06F, 3.1908506e-06F, 3.3982101e-06F,
 858   3.6190449e-06F, 3.8542308e-06F, 4.1047004e-06F, 4.3714470e-06F,
 859   4.6555282e-06F, 4.9580707e-06F, 5.2802740e-06F, 5.6234160e-06F,
 860   5.9888572e-06F, 6.3780469e-06F, 6.7925283e-06F, 7.2339451e-06F,
 861   7.7040476e-06F, 8.2047000e-06F, 8.7378876e-06F, 9.3057248e-06F,
 862   9.9104632e-06F, 1.0554501e-05F, 1.1240392e-05F, 1.1970856e-05F,
 863   1.2748789e-05F, 1.3577278e-05F, 1.4459606e-05F, 1.5399272e-05F,
 864   1.6400004e-05F, 1.7465768e-05F, 1.8600792e-05F, 1.9809576e-05F,
 865   2.1096914e-05F, 2.2467911e-05F, 2.3928002e-05F, 2.5482978e-05F,
 866   2.7139006e-05F, 2.8902651e-05F, 3.0780908e-05F, 3.2781225e-05F,
 867   3.4911534e-05F, 3.7180282e-05F, 3.9596466e-05F, 4.2169667e-05F,
 868   4.4910090e-05F, 4.7828601e-05F, 5.0936773e-05F, 5.4246931e-05F,
 869   5.7772202e-05F, 6.1526565e-05F, 6.5524908e-05F, 6.9783085e-05F,
 870   7.4317983e-05F, 7.9147585e-05F, 8.4291040e-05F, 8.9768747e-05F,
 871   9.5602426e-05F, 0.00010181521F, 0.00010843174F, 0.00011547824F,
 872   0.00012298267F, 0.00013097477F, 0.00013948625F, 0.00014855085F,
 873   0.00015820453F, 0.00016848555F, 0.00017943469F, 0.00019109536F,
 874   0.00020351382F, 0.00021673929F, 0.00023082423F, 0.00024582449F,
 875   0.00026179955F, 0.00027881276F, 0.00029693158F, 0.00031622787F,
 876   0.00033677814F, 0.00035866388F, 0.00038197188F, 0.00040679456F,
 877   0.00043323036F, 0.00046138411F, 0.00049136745F, 0.00052329927F,
 878   0.00055730621F, 0.00059352311F, 0.00063209358F, 0.00067317058F,
 879   0.00071691700F, 0.00076350630F, 0.00081312324F, 0.00086596457F,
 880   0.00092223983F, 0.00098217216F, 0.0010459992F, 0.0011139742F,
 881   0.0011863665F, 0.0012634633F, 0.0013455702F, 0.0014330129F,
 882   0.0015261382F, 0.0016253153F, 0.0017309374F, 0.0018434235F,
 883   0.0019632195F, 0.0020908006F, 0.0022266726F, 0.0023713743F,
 884   0.0025254795F, 0.0026895994F, 0.0028643847F, 0.0030505286F,
 885   0.0032487691F, 0.0034598925F, 0.0036847358F, 0.0039241906F,
 886   0.0041792066F, 0.0044507950F, 0.0047400328F, 0.0050480668F,
 887   0.0053761186F, 0.0057254891F, 0.0060975636F, 0.0064938176F,
 888   0.0069158225F, 0.0073652516F, 0.0078438871F, 0.0083536271F,
 889   0.0088964928F, 0.009474637F, 0.010090352F, 0.010746080F,
 890   0.011444421F, 0.012188144F, 0.012980198F, 0.013823725F,
 891   0.014722068F, 0.015678791F, 0.016697687F, 0.017782797F,
 892   0.018938423F, 0.020169149F, 0.021479854F, 0.022875735F,
 893   0.024362330F, 0.025945531F, 0.027631618F, 0.029427276F,
 894   0.031339626F, 0.033376252F, 0.035545228F, 0.037855157F,
 895   0.040315199F, 0.042935108F, 0.045725273F, 0.048696758F,
 896   0.051861348F, 0.055231591F, 0.058820850F, 0.062643361F,
 897   0.066714279F, 0.071049749F, 0.075666962F, 0.080584227F,
 898   0.085821044F, 0.091398179F, 0.097337747F, 0.10366330F,
 899   0.11039993F, 0.11757434F, 0.12521498F, 0.13335215F,
 900   0.14201813F, 0.15124727F, 0.16107617F, 0.17154380F,
 901   0.18269168F, 0.19456402F, 0.20720788F, 0.22067342F,
 902   0.23501402F, 0.25028656F, 0.26655159F, 0.28387361F,
 903   0.30232132F, 0.32196786F, 0.34289114F, 0.36517414F,
 904   0.38890521F, 0.41417847F, 0.44109412F, 0.46975890F,
 905   0.50028648F, 0.53279791F, 0.56742212F, 0.60429640F,
 906   0.64356699F, 0.68538959F, 0.72993007F, 0.77736504F,
 907   0.82788260F, 0.88168307F, 0.9389798F, 1.F,
 908 };
 909
 910 /* this is for per-channel noise normalization */
 911 static int apsort(const void *a, const void *b){
 912   float f1=**(float**)a;
 913   float f2=**(float**)b;
 914   return (f1<f2)-(f1>f2);
 915 }
 916
 917 static void flag_lossless(int limit, float prepoint, float postpoint, float *mdct,
 918                          float *floor, int *flag, int i, int jn){
 919   int j;
 920   for(j=0;j<jn;j++){
 921     float point = j>=limit-i ? postpoint : prepoint;
 922     float r = fabs(mdct[j])/floor[j];
 923     if(r<point)
 924       flag[j]=0;
 925     else
 926       flag[j]=1;
 927   }
 928 }
 929
 930 /* Overload/Side effect: On input, the *q vector holds either the
 931    quantized energy (for elements with the flag set) or the absolute
 932    values of the *r vector (for elements with flag unset).  On output,
 933    *q holds the quantized energy for all elements */
 934 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){
 935
 936   vorbis_info_psy *vi=p->vi;
 937   float **sort = alloca(n*sizeof(*sort));
 938   int j,count=0;
 939   int start = (vi->normal_p ? vi->normal_start-i : n);
 940   if(start>n)start=n;
 941
 942   /* force classic behavior where only energy in the current band is considered */
 943   acc=0.f;
 944
 945   /* still responsible for populating *out where noise norm not in
 946      effect.  There's no need to [re]populate *q in these areas */
 947   for(j=0;j<start;j++){
 948     if(!flags || !flags[j]){ /* lossless coupling already quantized.
 949                                 Don't touch; requantizing based on
 950                                 energy would be incorrect. */
 951       float ve = q[j]/f[j];
 952       if(r[j]<0)
 953         out[j] = -rint(sqrt(ve));
 954       else
 955         out[j] = rint(sqrt(ve));
 956     }
 957   }
 958
 959   /* sort magnitudes for noise norm portion of partition */
 960   for(;j<n;j++){
 961     if(!flags || !flags[j]){ /* can't noise norm elements that have
 962                                 already been loslessly coupled; we can
 963                                 only account for their energy error */
 964       float ve = q[j]/f[j];
 965       /* Despite all the new, more capable coupling code, for now we
 966          implement noise norm as it has been up to this point. Only
 967          consider promotions to unit magnitude from 0.  In addition
 968          the only energy error counted is quantizations to zero. */
 969       /* also-- the original point code only applied noise norm at > pointlimit */
 970       if(ve<.25f && (!flags || j>=limit-i)){
 971         acc += ve;
 972         sort[count++]=q+j; /* q is fabs(r) for unflagged element */
 973       }else{
 974         /* For now: no acc adjustment for nonzero quantization.  populate *out and q as this value is final. */
 975         if(r[j]<0)
 976           out[j] = -rint(sqrt(ve));
 977         else
 978           out[j] = rint(sqrt(ve));
 979         q[j] = out[j]*out[j]*f[j];
 980       }
 981     }/* else{
 982         again, no energy adjustment for error in nonzero quant-- for now
 983         }*/
 984   }
 985
 986   if(count){
 987     /* noise norm to do */
 988     qsort(sort,count,sizeof(*sort),apsort);
 989     for(j=0;j<count;j++){
 990       int k=sort[j]-q;
 991       if(acc>=vi->normal_thresh){
 992         out[k]=unitnorm(r[k]);
 993         acc-=1.f;
 994         q[k]=f[k];
 995       }else{
 996         out[k]=0;
 997         q[k]=0.f;
 998       }
 999     }
1000   }
1001
1002   return acc;
1003 }
1004
1005 /* Noise normalization, quantization and coupling are not wholly
1006    seperable processes in depth>1 coupling. */
1007 void _vp_couple_quantize_normalize(int blobno,
1008                                    vorbis_info_psy_global *g,
1009                                    vorbis_look_psy *p,
1010                                    vorbis_info_mapping0 *vi,
1011                                    float **mdct,
1012                                    int   **iwork,
1013                                    int    *nonzero,
1014                                    int     sliding_lowpass,
1015                                    int     ch){
1016
1017   int i;
1018   int n = p->n;
1019   int partition=(p->vi->normal_p ? p->vi->normal_partition : 16);
1020   int limit = g->coupling_pointlimit[p->vi->blockflag][blobno];
1021   float prepoint=stereo_threshholds[g->coupling_prepointamp[blobno]];
1022   float postpoint=stereo_threshholds[g->coupling_postpointamp[blobno]];
1023 #if 0
1024   float de=0.1*p->m_val; /* a blend of the AoTuV M2 and M3 code here and below */
1025 #endif
1026
1027   /* mdct is our raw mdct output, floor not removed. */
1028   /* inout passes in the ifloor, passes back quantized result */
1029
1030   /* unquantized energy (negative indicates amplitude has negative sign) */
1031   float **raw = alloca(ch*sizeof(*raw));
1032
1033   /* dual pupose; quantized energy (if flag set), othersize fabs(raw) */
1034   float **quant = alloca(ch*sizeof(*quant));
1035
1036   /* floor energy */
1037   float **floor = alloca(ch*sizeof(*floor));
1038
1039   /* flags indicating raw/quantized status of elements in raw vector */
1040   int   **flag  = alloca(ch*sizeof(*flag));
1041
1042   /* non-zero flag working vector */
1043   int    *nz    = alloca(ch*sizeof(*nz));
1044
1045   /* energy surplus/defecit tracking */
1046   float  *acc   = alloca((ch+vi->coupling_steps)*sizeof(*acc));
1047
1048   /* The threshold of a stereo is changed with the size of n */
1049   if(n > 1000)
1050     postpoint=stereo_threshholds_limited[g->coupling_postpointamp[blobno]];
1051
1052   raw[0]   = alloca(ch*partition*sizeof(**raw));
1053   quant[0] = alloca(ch*partition*sizeof(**quant));
1054   floor[0] = alloca(ch*partition*sizeof(**floor));
1055   flag[0]  = alloca(ch*partition*sizeof(**flag));
1056
1057   for(i=1;i<ch;i++){
1058     raw[i]   = &raw[0][partition*i];
1059     quant[i] = &quant[0][partition*i];
1060     floor[i] = &floor[0][partition*i];
1061     flag[i]  = &flag[0][partition*i];
1062   }
1063   for(i=0;i<ch+vi->coupling_steps;i++)
1064     acc[i]=0.f;
1065
1066   for(i=0;i<n;i+=partition){
1067     int k,j,jn = partition > n-i ? n-i : partition;
1068     int step,track = 0;
1069
1070     memcpy(nz,nonzero,sizeof(*nz)*ch);
1071
1072     /* prefill */
1073     memset(flag[0],0,ch*partition*sizeof(**flag));
1074     for(k=0;k<ch;k++){
1075       int *iout = &iwork[k][i];
1076       if(nz[k]){
1077
1078         for(j=0;j<jn;j++)
1079           floor[k][j] = FLOOR1_fromdB_LOOKUP[iout[j]];
1080
1081         flag_lossless(limit,prepoint,postpoint,&mdct[k][i],floor[k],flag[k],i,jn);
1082
1083         for(j=0;j<jn;j++){
1084           quant[k][j] = raw[k][j] = mdct[k][i+j]*mdct[k][i+j];
1085           if(mdct[k][i+j]<0.f) raw[k][j]*=-1.f;
1086           floor[k][j]*=floor[k][j];
1087         }
1088
1089         acc[track]=noise_normalize(p,limit,raw[k],quant[k],floor[k],NULL,acc[track],i,jn,iout);
1090
1091       }else{
1092         for(j=0;j<jn;j++){
1093           floor[k][j] = 1e-10f;
1094           raw[k][j] = 0.f;
1095           quant[k][j] = 0.f;
1096           flag[k][j] = 0;
1097           iout[j]=0;
1098         }
1099         acc[track]=0.f;
1100       }
1101       track++;
1102     }
1103
1104     /* coupling */
1105     for(step=0;step<vi->coupling_steps;step++){
1106       int Mi = vi->coupling_mag[step];
1107       int Ai = vi->coupling_ang[step];
1108       int *iM = &iwork[Mi][i];
1109       int *iA = &iwork[Ai][i];
1110       float *reM = raw[Mi];
1111       float *reA = raw[Ai];
1112       float *qeM = quant[Mi];
1113       float *qeA = quant[Ai];
1114       float *floorM = floor[Mi];
1115       float *floorA = floor[Ai];
1116       int *fM = flag[Mi];
1117       int *fA = flag[Ai];
1118
1119       if(nz[Mi] || nz[Ai]){
1120         nz[Mi] = nz[Ai] = 1;
1121
1122         for(j=0;j<jn;j++){
1123
1124           if(j<sliding_lowpass-i){
1125             if(fM[j] || fA[j]){
1126               /* lossless coupling */
1127
1128               reM[j] = fabs(reM[j])+fabs(reA[j]);
1129               qeM[j] = qeM[j]+qeA[j];
1130               fM[j]=fA[j]=1;
1131
1132               /* couple iM/iA */
1133               {
1134                 int A = iM[j];
1135                 int B = iA[j];
1136
1137                 if(abs(A)>abs(B)){
1138                   iA[j]=(A>0?A-B:B-A);
1139                 }else{
1140                   iA[j]=(B>0?A-B:B-A);
1141                   iM[j]=B;
1142                 }
1143
1144                 /* collapse two equivalent tuples to one */
1145                 if(iA[j]>=abs(iM[j])*2){
1146                   iA[j]= -iA[j];
1147                   iM[j]= -iM[j];
1148                 }
1149
1150               }
1151
1152             }else{
1153               /* lossy (point) coupling */
1154               if(j<limit-i){
1155                 /* dipole */
1156                 reM[j] += reA[j];
1157                 qeM[j] = fabs(reM[j]);
1158               }else{
1159 #if 0
1160                 /* AoTuV */
1161                 /** @ M2 **
1162                     The boost problem by the combination of noise normalization and point stereo is eased.
1163                     However, this is a temporary patch.
1164                     by Aoyumi @ 2004/04/18
1165                 */
1166                 float derate = (1.0 - de*((float)(j-limit+i) / (float)(n-limit)));
1167                 /* elliptical */
1168                 if(reM[j]+reA[j]<0){
1169                   reM[j] = - (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate);
1170                 }else{
1171                   reM[j] =   (qeM[j] = (fabs(reM[j])+fabs(reA[j]))*derate*derate);
1172                 }
1173 #else
1174                 /* elliptical */
1175                 if(reM[j]+reA[j]<0){
1176                   reM[j] = - (qeM[j] = fabs(reM[j])+fabs(reA[j]));
1177                 }else{
1178                   reM[j] =   (qeM[j] = fabs(reM[j])+fabs(reA[j]));
1179                 }
1180 #endif
1181
1182               }
1183               reA[j]=qeA[j]=0.f;
1184               fA[j]=1;
1185               iA[j]=0;
1186             }
1187           }
1188           floorM[j]=floorA[j]=floorM[j]+floorA[j];
1189         }
1190         /* normalize the resulting mag vector */
1191         acc[track]=noise_normalize(p,limit,raw[Mi],quant[Mi],floor[Mi],flag[Mi],acc[track],i,jn,iM);
1192         track++;
1193       }
1194     }
1195   }
1196
1197   for(i=0;i<vi->coupling_steps;i++){
1198     /* make sure coupling a zero and a nonzero channel results in two
1199        nonzero channels. */
1200     if(nonzero[vi->coupling_mag[i]] ||
1201        nonzero[vi->coupling_ang[i]]){
1202       nonzero[vi->coupling_mag[i]]=1;
1203       nonzero[vi->coupling_ang[i]]=1;
1204     }
1205   }
1206 }