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