lib/lame/psymodel.c

   1 /*
   2  *      psymodel.c
   3  *
   4  *      Copyright (c) 1999 Mark Taylor
   5  *
   6  * This library is free software; you can redistribute it and/or
   7  * modify it under the terms of the GNU Library General Public
   8  * License as published by the Free Software Foundation; either
   9  * version 2 of the License, or (at your option) any later version.
  10  *
  11  * This library is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * Library General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Library General Public
  17  * License along with this library; if not, write to the
  18  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  19  * Boston, MA 02111-1307, USA.
  20  */
  21
  22 /* $Id: psymodel.c,v 1.1 2002/04/28 17:30:23 kramm Exp $ */
  23
  24
  25 /*
  26 PSYCHO ACOUSTICS
  27
  28
  29 This routine computes the psycho acoustics, delayed by
  30 one granule.
  31
  32 Input: buffer of PCM data (1024 samples).
  33
  34 This window should be centered over the 576 sample granule window.
  35 The routine will compute the psycho acoustics for
  36 this granule, but return the psycho acoustics computed
  37 for the *previous* granule.  This is because the block
  38 type of the previous granule can only be determined
  39 after we have computed the psycho acoustics for the following
  40 granule.
  41
  42 Output:  maskings and energies for each scalefactor band.
  43 block type, PE, and some correlation measures.
  44 The PE is used by CBR modes to determine if extra bits
  45 from the bit reservoir should be used.  The correlation
  46 measures are used to determine mid/side or regular stereo.
  47
  48
  49 Notation:
  50
  51 barks:  a non-linear frequency scale.  Mapping from frequency to
  52         barks is given by freq2bark()
  53
  54 scalefactor bands: The spectrum (frequencies) are broken into
  55                    SBMAX "scalefactor bands".  Thes bands
  56                    are determined by the MPEG ISO spec.  In
  57                    the noise shaping/quantization code, we allocate
  58                    bits among the partition bands to achieve the
  59                    best possible quality
  60
  61 partition bands:   The spectrum is also broken into about
  62                    64 "partition bands".  Each partition
  63                    band is about .34 barks wide.  There are about 2-5
  64                    partition bands for each scalefactor band.
  65
  66 LAME computes all psycho acoustic information for each partition
  67 band.  Then at the end of the computations, this information
  68 is mapped to scalefactor bands.  The energy in each scalefactor
  69 band is taken as the sum of the energy in all partition bands
  70 which overlap the scalefactor band.  The maskings can be computed
  71 in the same way (and thus represent the average masking in that band)
  72 or by taking the minmum value multiplied by the number of
  73 partition bands used (which represents a minimum masking in that band).
  74
  75
  76 The general outline is as follows:
  77
  78
  79 1. compute the energy in each partition band
  80 2. compute the tonality in each partition band
  81 3. compute the strength of each partion band "masker"
  82 4. compute the masking (via the spreading function applied to each masker)
  83 5. Modifications for mid/side masking.
  84
  85 Each partition band is considiered a "masker".  The strength
  86 of the i'th masker in band j is given by:
  87
  88     s3(bark(i)-bark(j))*strength(i)
  89
  90 The strength of the masker is a function of the energy and tonality.
  91 The more tonal, the less masking.  LAME uses a simple linear formula
  92 (controlled by NMT and TMN) which says the strength is given by the
  93 energy divided by a linear function of the tonality.
  94
  95
  96 s3() is the "spreading function".  It is given by a formula
  97 determined via listening tests.
  98
  99 The total masking in the j'th partition band is the sum over
 100 all maskings i.  It is thus given by the convolution of
 101 the strength with s3(), the "spreading function."
 102
 103 masking(j) = sum_over_i  s3(i-j)*strength(i)  = s3 o strength
 104
 105 where "o" = convolution operator.  s3 is given by a formula determined
 106 via listening tests.  It is normalized so that s3 o 1 = 1.
 107
 108 Note: instead of a simple convolution, LAME also has the
 109 option of using "additive masking"
 110
 111 The most critical part is step 2, computing the tonality of each
 112 partition band.  LAME has two tonality estimators.  The first
 113 is based on the ISO spec, and measures how predictiable the
 114 signal is over time.  The more predictable, the more tonal.
 115 The second measure is based on looking at the spectrum of
 116 a single granule.  The more peaky the spectrum, the more
 117 tonal.  By most indications, the latter approach is better.
 118
 119 Finally, in step 5, the maskings for the mid and side
 120 channel are possibly increased.  Under certain circumstances,
 121 noise in the mid & side channels is assumed to also
 122 be masked by strong maskers in the L or R channels.
 123
 124
 125 Other data computed by the psy-model:
 126
 127 ms_ratio        side-channel / mid-channel masking ratio (for previous granule)
 128 ms_ratio_next   side-channel / mid-channel masking ratio for this granule
 129
 130 percep_entropy[2]     L and R values (prev granule) of PE - A measure of how
 131                       much pre-echo is in the previous granule
 132 percep_entropy_MS[2]  mid and side channel values (prev granule) of percep_entropy
 133 energy[4]             L,R,M,S energy in each channel, prev granule
 134 blocktype_d[2]        block type to use for previous granule
 135
 136
 137 */
 138
 139
 140
 141
 142 #include "config_static.h"
 143
 144 #include "util.h"
 145 #include "encoder.h"
 146 #include "psymodel.h"
 147 #include "l3side.h"
 148 #include <assert.h>
 149 #include "tables.h"
 150 #include "fft.h"
 151
 152 #ifdef WITH_DMALLOC
 153 #include <dmalloc.h>
 154 #endif
 155
 156 #define NSFIRLEN 21
 157 #define rpelev 2
 158 #define rpelev2 16
 159 #define rpelev_s 2
 160 #define rpelev2_s 16
 161
 162 /* size of each partition band, in barks: */
 163 #define DELBARK .34
 164
 165
 166 #if 1
 167     /* AAC values, results in more masking over MP3 values */
 168 # define TMN 18
 169 # define NMT 6
 170 #else
 171     /* MP3 values */
 172 # define TMN 29
 173 # define NMT 6
 174 #endif
 175
 176 #define NBPSY_l  (SBMAX_l)
 177 #define NBPSY_s  (SBMAX_s)
 178
 179
 180 #ifdef M_LN10
 181 #define         LN_TO_LOG10             (M_LN10/10)
 182 #else
 183 #define         LN_TO_LOG10             0.2302585093
 184 #endif
 185
 186 FLOAT
 187 psycho_loudness_approx( FLOAT *energy, lame_global_flags *gfp );
 188
 189
 190
 191
 192
 193 /*
 194    L3psycho_anal.  Compute psycho acoustics.
 195
 196    Data returned to the calling program must be delayed by one
 197    granule.
 198
 199    This is done in two places.
 200    If we do not need to know the blocktype, the copying
 201    can be done here at the top of the program: we copy the data for
 202    the last granule (computed during the last call) before it is
 203    overwritten with the new data.  It looks like this:
 204
 205    0. static psymodel_data
 206    1. calling_program_data = psymodel_data
 207    2. compute psymodel_data
 208
 209    For data which needs to know the blocktype, the copying must be
 210    done at the end of this loop, and the old values must be saved:
 211
 212    0. static psymodel_data_old
 213    1. compute psymodel_data
 214    2. compute possible block type of this granule
 215    3. compute final block type of previous granule based on #2.
 216    4. calling_program_data = psymodel_data_old
 217    5. psymodel_data_old = psymodel_data
 218 */
 219 int L3psycho_anal( lame_global_flags * gfp,
 220                     const sample_t *buffer[2], int gr_out,
 221                     FLOAT8 *ms_ratio,
 222                     FLOAT8 *ms_ratio_next,
 223                     III_psy_ratio masking_ratio[2][2],
 224                     III_psy_ratio masking_MS_ratio[2][2],
 225                     FLOAT8 percep_entropy[2],FLOAT8 percep_MS_entropy[2],
 226                     FLOAT8 energy[4],
 227                     int blocktype_d[2])
 228 {
 229 /* to get a good cache performance, one has to think about
 230  * the sequence, in which the variables are used.
 231  * (Note: these static variables have been moved to the gfc-> struct,
 232  * and their order in memory is layed out in util.h)
 233  */
 234   lame_internal_flags *gfc=gfp->internal_flags;
 235
 236
 237   /* fft and energy calculation   */
 238   FLOAT (*wsamp_l)[BLKSIZE];
 239   FLOAT (*wsamp_s)[3][BLKSIZE_s];
 240
 241   /* convolution   */
 242   FLOAT8 eb[CBANDS];
 243   FLOAT8 cb[CBANDS];
 244   FLOAT8 thr[CBANDS];
 245
 246   /* ratios    */
 247   FLOAT8 ms_ratio_l=0,ms_ratio_s=0;
 248
 249   /* block type  */
 250   int blocktype[2],uselongblock[2];
 251
 252   /* usual variables like loop indices, etc..    */
 253   int numchn, chn;
 254   int   b, i, j, k;
 255   int   sb,sblock;
 256
 257
 258   if(gfc->psymodel_init==0) {
 259       psymodel_init(gfp);
 260       init_fft(gfc);
 261       gfc->psymodel_init=1;
 262
 263       for (chn = 0; chn < 4; ++chn )
 264       for (b = 0; b < CBANDS; ++b )
 265         gfc->nb_s1[chn][b] = gfc->nb_s2[chn][b] = 1.0;
 266
 267   }
 268
 269
 270
 271
 272
 273   numchn = gfc->channels_out;
 274   /* chn=2 and 3 = Mid and Side channels */
 275   if (gfp->mode == JOINT_STEREO) numchn=4;
 276
 277   for (chn=0; chn<numchn; chn++) {
 278       for (i=0; i<numchn; ++i) {
 279           energy[chn]=gfc->tot_ener[chn];
 280       }
 281   }
 282
 283   for (chn=0; chn<numchn; chn++) {
 284
 285       /* there is a one granule delay.  Copy maskings computed last call
 286        * into masking_ratio to return to calling program.
 287        */
 288       if (chn < 2) {
 289           /* LR maskings  */
 290           percep_entropy            [chn]       = gfc -> pe  [chn];
 291           masking_ratio    [gr_out] [chn]  .en  = gfc -> en  [chn];
 292           masking_ratio    [gr_out] [chn]  .thm = gfc -> thm [chn];
 293       } else {
 294           /* MS maskings  */
 295           percep_MS_entropy         [chn-2]     = gfc -> pe  [chn];
 296           masking_MS_ratio [gr_out] [chn-2].en  = gfc -> en  [chn];
 297           masking_MS_ratio [gr_out] [chn-2].thm = gfc -> thm [chn];
 298       }
 299
 300
 301
 302     /**********************************************************************
 303      *  compute FFTs
 304      **********************************************************************/
 305     wsamp_s = gfc->wsamp_S+(chn & 1);
 306     wsamp_l = gfc->wsamp_L+(chn & 1);
 307     if (chn<2) {
 308       fft_long ( gfc, *wsamp_l, chn, buffer);
 309       fft_short( gfc, *wsamp_s, chn, buffer);
 310     }
 311     /* FFT data for mid and side channel is derived from L & R */
 312     if (chn == 2)
 313       {
 314         for (j = BLKSIZE-1; j >=0 ; --j)
 315         {
 316           FLOAT l = gfc->wsamp_L[0][j];
 317           FLOAT r = gfc->wsamp_L[1][j];
 318           gfc->wsamp_L[0][j] = (l+r)*(FLOAT)(SQRT2*0.5);
 319           gfc->wsamp_L[1][j] = (l-r)*(FLOAT)(SQRT2*0.5);
 320         }
 321         for (b = 2; b >= 0; --b)
 322         {
 323           for (j = BLKSIZE_s-1; j >= 0 ; --j)
 324           {
 325             FLOAT l = gfc->wsamp_S[0][b][j];
 326             FLOAT r = gfc->wsamp_S[1][b][j];
 327             gfc->wsamp_S[0][b][j] = (l+r)*(FLOAT)(SQRT2*0.5);
 328             gfc->wsamp_S[1][b][j] = (l-r)*(FLOAT)(SQRT2*0.5);
 329           }
 330         }
 331       }
 332
 333
 334     /**********************************************************************
 335      *  compute energies
 336      **********************************************************************/
 337
 338
 339
 340     gfc->energy[0]  = (*wsamp_l)[0];
 341     gfc->energy[0] *= gfc->energy[0];
 342
 343     gfc->tot_ener[chn] = gfc->energy[0]; /* sum total energy at nearly no extra cost */
 344
 345     for (j=BLKSIZE/2-1; j >= 0; --j)
 346     {
 347       FLOAT re = (*wsamp_l)[BLKSIZE/2-j];
 348       FLOAT im = (*wsamp_l)[BLKSIZE/2+j];
 349       gfc->energy[BLKSIZE/2-j] = (re * re + im * im) * (FLOAT)0.5;
 350
 351       if (BLKSIZE/2-j > 10)
 352         gfc->tot_ener[chn] += gfc->energy[BLKSIZE/2-j];
 353     }
 354     for (b = 2; b >= 0; --b)
 355     {
 356       gfc->energy_s[b][0]  = (*wsamp_s)[b][0];
 357       gfc->energy_s[b][0] *=  gfc->energy_s [b][0];
 358       for (j=BLKSIZE_s/2-1; j >= 0; --j)
 359       {
 360         FLOAT re = (*wsamp_s)[b][BLKSIZE_s/2-j];
 361         FLOAT im = (*wsamp_s)[b][BLKSIZE_s/2+j];
 362         gfc->energy_s[b][BLKSIZE_s/2-j] = (re * re + im * im) * (FLOAT)0.5;
 363       }
 364     }
 365
 366 #if defined(HAVE_GTK)
 367     if (gfp->analysis) {
 368       for (j=0; j<HBLKSIZE ; j++) {
 369         gfc->pinfo->energy[gr_out][chn][j]=gfc->energy_save[chn][j];
 370         gfc->energy_save[chn][j]=gfc->energy[j];
 371       }
 372     }
 373 #endif
 374
 375
 376
 377     /**********************************************************************
 378      * compute loudness approximation (used for ATH auto-level adjustment)
 379      **********************************************************************/
 380     if( gfp->athaa_loudapprox == 2 ) {
 381       if( chn < 2 ) {           /* no loudness for mid and side channels */
 382         gfc->loudness_sq[gr_out][chn] = gfc->loudness_sq_save[chn];
 383         gfc->loudness_sq_save[chn]
 384           = psycho_loudness_approx( gfc->energy, gfp);
 385       }
 386     }
 387
 388
 389
 390     /**********************************************************************
 391      *    compute unpredicatability of first six spectral lines            *
 392      **********************************************************************/
 393     for ( j = 0; j < gfc->cw_lower_index; j++ )
 394       {  /* calculate unpredictability measure cw */
 395         FLOAT an, a1, a2;
 396         FLOAT bn, b1, b2;
 397         FLOAT rn, r1, r2;
 398         FLOAT numre, numim, den;
 399
 400         a2 = gfc-> ax_sav[chn][1][j];
 401         b2 = gfc-> bx_sav[chn][1][j];
 402         r2 = gfc-> rx_sav[chn][1][j];
 403         a1 = gfc-> ax_sav[chn][1][j] = gfc-> ax_sav[chn][0][j];
 404         b1 = gfc-> bx_sav[chn][1][j] = gfc-> bx_sav[chn][0][j];
 405         r1 = gfc-> rx_sav[chn][1][j] = gfc-> rx_sav[chn][0][j];
 406         an = gfc-> ax_sav[chn][0][j] = (*wsamp_l)[j];
 407         bn = gfc-> bx_sav[chn][0][j] = j==0 ? (*wsamp_l)[0] : (*wsamp_l)[BLKSIZE-j];
 408         rn = gfc-> rx_sav[chn][0][j] = sqrt(gfc->energy[j]);
 409
 410         { /* square (x1,y1) */
 411           if( r1 != 0 ) {
 412             numre = (a1*b1);
 413             numim = (a1*a1-b1*b1)*(FLOAT)0.5;
 414             den = r1*r1;
 415           } else {
 416             numre = 1;
 417             numim = 0;
 418             den = 1;
 419           }
 420         }
 421
 422         { /* multiply by (x2,-y2) */
 423           if( r2 != 0 ) {
 424             FLOAT tmp2 = (numim+numre)*(a2+b2)*(FLOAT)0.5;
 425             FLOAT tmp1 = -a2*numre+tmp2;
 426             numre =       -b2*numim+tmp2;
 427             numim = tmp1;
 428             den *= r2;
 429           } else {
 430             /* do nothing */
 431           }
 432         }
 433
 434         { /* r-prime factor */
 435           FLOAT tmp = (2*r1-r2)/den;
 436           numre *= tmp;
 437           numim *= tmp;
 438         }
 439         den=rn+fabs(2*r1-r2);
 440         if( den != 0 ) {
 441           numre = (an+bn)*(FLOAT)0.5-numre;
 442           numim = (an-bn)*(FLOAT)0.5-numim;
 443           den = sqrt(numre*numre+numim*numim)/den;
 444         }
 445         gfc->cw[j] = den;
 446       }
 447
 448
 449
 450     /**********************************************************************
 451      *     compute unpredicatibility of next 200 spectral lines            *
 452      **********************************************************************/
 453     for ( j = gfc->cw_lower_index; j < gfc->cw_upper_index; j += 4 )
 454       {/* calculate unpredictability measure cw */
 455         FLOAT rn, r1, r2;
 456         FLOAT numre, numim, den;
 457
 458         k = (j+2) / 4;
 459
 460         { /* square (x1,y1) */
 461           r1 = gfc->energy_s[0][k];
 462           if( r1 != 0 ) {
 463             FLOAT a1 = (*wsamp_s)[0][k];
 464             FLOAT b1 = (*wsamp_s)[0][BLKSIZE_s-k]; /* k is never 0 */
 465             numre = (a1*b1);
 466             numim = (a1*a1-b1*b1)*(FLOAT)0.5;
 467             den = r1;
 468             r1 = sqrt(r1);
 469           } else {
 470             numre = 1;
 471             numim = 0;
 472             den = 1;
 473           }
 474         }
 475
 476
 477         { /* multiply by (x2,-y2) */
 478           r2 = gfc->energy_s[2][k];
 479           if( r2 != 0 ) {
 480             FLOAT a2 = (*wsamp_s)[2][k];
 481             FLOAT b2 = (*wsamp_s)[2][BLKSIZE_s-k];
 482
 483
 484             FLOAT tmp2 = (numim+numre)*(a2+b2)*(FLOAT)0.5;
 485             FLOAT tmp1 = -a2*numre+tmp2;
 486             numre =       -b2*numim+tmp2;
 487             numim = tmp1;
 488
 489             r2 = sqrt(r2);
 490             den *= r2;
 491           } else {
 492             /* do nothing */
 493           }
 494         }
 495
 496         { /* r-prime factor */
 497           FLOAT tmp = (2*r1-r2)/den;
 498           numre *= tmp;
 499           numim *= tmp;
 500         }
 501
 502         rn = sqrt(gfc->energy_s[1][k]);
 503         den=rn+fabs(2*r1-r2);
 504         if( den != 0 ) {
 505           FLOAT an = (*wsamp_s)[1][k];
 506           FLOAT bn = (*wsamp_s)[1][BLKSIZE_s-k];
 507           numre = (an+bn)*(FLOAT)0.5-numre;
 508           numim = (an-bn)*(FLOAT)0.5-numim;
 509           den = sqrt(numre*numre+numim*numim)/den;
 510         }
 511
 512         gfc->cw[j+1] = gfc->cw[j+2] = gfc->cw[j+3] = gfc->cw[j] = den;
 513       }
 514
 515 #if 0
 516     for ( j = 14; j < HBLKSIZE-4; j += 4 )
 517       {/* calculate energy from short ffts */
 518         FLOAT8 tot,ave;
 519         k = (j+2) / 4;
 520         for (tot=0, sblock=0; sblock < 3; sblock++)
 521           tot+=gfc->energy_s[sblock][k];
 522         ave = gfc->energy[j+1]+ gfc->energy[j+2]+ gfc->energy[j+3]+ gfc->energy[j];
 523         ave /= 4.;
 524         gfc->energy[j+1] = gfc->energy[j+2] = gfc->energy[j+3] =  gfc->energy[j]=tot;
 525       }
 526 #endif
 527
 528     /**********************************************************************
 529      *    Calculate the energy and the unpredictability in the threshold   *
 530      *    calculation partitions                                           *
 531      **********************************************************************/
 532
 533     b = 0;
 534     for (j = 0; j < gfc->cw_upper_index && gfc->numlines_l[b] && b<gfc->npart_l_orig; )
 535       {
 536         FLOAT8 ebb, cbb;
 537
 538         ebb = gfc->energy[j];
 539         cbb = gfc->energy[j] * gfc->cw[j];
 540         j++;
 541
 542         for (i = gfc->numlines_l[b] - 1; i > 0; i--)
 543           {
 544             ebb += gfc->energy[j];
 545             cbb += gfc->energy[j] * gfc->cw[j];
 546             j++;
 547           }
 548         eb[b] = ebb;
 549         cb[b] = cbb;
 550         b++;
 551       }
 552
 553     for (; b < gfc->npart_l_orig; b++ )
 554       {
 555         FLOAT8 ebb = gfc->energy[j++];
 556         assert(gfc->numlines_l[b]);
 557         for (i = gfc->numlines_l[b] - 1; i > 0; i--)
 558           {
 559             ebb += gfc->energy[j++];
 560           }
 561         eb[b] = ebb;
 562         cb[b] = ebb * 0.4;
 563       }
 564
 565     /**********************************************************************
 566      *      convolve the partitioned energy and unpredictability          *
 567      *      with the spreading function, s3_l[b][k](packed into s3_ll)    *
 568      ******************************************************************** */
 569     gfc->pe[chn] = 0;           /*  calculate percetual entropy */
 570     {
 571       int kk = 0;
 572       for ( b = 0;b < gfc->npart_l; b++ )
 573         {
 574           FLOAT8 tbb,ecb,ctb;
 575
 576           ecb = 0;
 577           ctb = 0;
 578
 579           for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ )
 580             {
 581               ecb += gfc->s3_ll[kk] * eb[k];    /* sprdngf for Layer III */
 582               ctb += gfc->s3_ll[kk] * cb[k];
 583               kk++;
 584             }
 585
 586           /* calculate the tonality of each threshold calculation partition
 587            * calculate the SNR in each threshold calculation partition
 588            * tonality = -0.299 - .43*log(ctb/ecb);
 589            * tonality = 0:           use NMT   (lots of masking)
 590            * tonality = 1:           use TMN   (little masking)
 591            */
 592
 593 /* ISO values */
 594 #define CONV1 (-.299)
 595 #define CONV2 (-.43)
 596
 597           tbb = ecb;
 598           if (tbb != 0)
 599             {
 600               tbb = ctb / tbb;
 601               if (tbb <= exp((1-CONV1)/CONV2))
 602                 { /* tonality near 1 */
 603                   tbb = exp(-LN_TO_LOG10 * TMN);
 604                 }
 605               else if (tbb >= exp((0-CONV1)/CONV2))
 606                 {/* tonality near 0 */
 607                   tbb = exp(-LN_TO_LOG10 * NMT);
 608                 }
 609               else
 610                 {
 611                   /* convert to tonality index */
 612                   /* tonality small:   tbb=1 */
 613                   /* tonality large:   tbb=-.299 */
 614                   tbb = CONV1 + CONV2*log(tbb);
 615                   tbb = exp(-LN_TO_LOG10 * ( TMN*tbb + (1-tbb)*NMT) );
 616                 }
 617             }
 618
 619           /* at this point, tbb represents the amount the spreading function
 620            * will be reduced.  The smaller the value, the less masking.
 621            * minval[] = 1 (0db)     says just use tbb.
 622            * minval[]= .01 (-20db)  says reduce spreading function by
 623            *                        at least 20db.
 624            */
 625
 626           tbb = Min(gfc->minval[b], tbb);
 627
 628             /* stabilize tonality estimation */
 629             if ( gfc->PSY->tonalityPatch ) {
 630                 if ( b > 5 )
 631                 {
 632                     FLOAT8 const x = 1.8699422;
 633                     FLOAT8 w = gfc->PSY->prvTonRed[b/2] * x;
 634                     if (tbb > w)
 635                         tbb = w;
 636                 }
 637                 gfc->PSY->prvTonRed[b] = tbb;
 638             }
 639
 640           ecb *= tbb;
 641
 642           /* long block pre-echo control.   */
 643           /* rpelev=2.0, rpelev2=16.0 */
 644           /* note: all surges in PE are because of this pre-echo formula
 645            * for thr[b].  If it this is not used, PE is always around 600
 646            */
 647           /* dont use long block pre-echo control if previous granule was
 648            * a short block.  This is to avoid the situation:
 649            * frame0:  quiet (very low masking)
 650            * frame1:  surge  (triggers short blocks)
 651            * frame2:  regular frame.  looks like pre-echo when compared to
 652            *          frame0, but all pre-echo was in frame1.
 653            */
 654           /* chn=0,1   L and R channels
 655              chn=2,3   S and M channels.
 656           */
 657
 658             if (gfc->blocktype_old[chn>1 ? chn-2 : chn] == SHORT_TYPE )
 659                 thr[b] = ecb; /* Min(ecb, rpelev*gfc->nb_1[chn][b]); */
 660             else
 661                 thr[b] = Min(ecb, Min(rpelev*gfc->nb_1[chn][b],rpelev2*gfc->nb_2[chn][b]) );
 662
 663           {
 664             FLOAT8 thrpe;
 665             thrpe = Max(thr[b],gfc->ATH->cb[b]);
 666             /*
 667               printf("%i thr=%e   ATH=%e  \n",b,thr[b],gfc->ATH->cb[b]);
 668             */
 669             if (thrpe < eb[b])
 670               gfc->pe[chn] -= gfc->numlines_l[b] * log(thrpe / eb[b]);
 671           }
 672
 673             if ( gfp->VBR == vbr_rh || gfp->VBR == vbr_mtrh ) {
 674                 thr[b] = Min(ecb, rpelev*gfc->nb_1[chn][b]);
 675                 if (gfc->blocktype_old[chn>1 ? chn-2 : chn] != SHORT_TYPE )
 676                     thr[b] = Min(thr[b], rpelev2*gfc->nb_2[chn][b]);
 677                 thr[b] = Max( thr[b], 1e-37 );
 678             }
 679
 680           gfc->nb_2[chn][b] = gfc->nb_1[chn][b];
 681           gfc->nb_1[chn][b] = ecb;
 682         }
 683     }
 684
 685     /***************************************************************
 686      * determine the block type (window type) based on L & R channels
 687      *
 688      ***************************************************************/
 689     {  /* compute PE for all 4 channels */
 690       FLOAT mn,mx,ma=0,mb=0,mc=0;
 691
 692       for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++)
 693         {
 694           ma += gfc->energy_s[0][j];
 695           mb += gfc->energy_s[1][j];
 696           mc += gfc->energy_s[2][j];
 697         }
 698       mn = Min(ma,mb);
 699       mn = Min(mn,mc);
 700       mx = Max(ma,mb);
 701       mx = Max(mx,mc);
 702
 703       /* bit allocation is based on pe.  */
 704       if (mx>mn) {
 705         FLOAT8 tmp = 400*log(mx/(1e-12+mn));
 706         if (tmp>gfc->pe[chn]) gfc->pe[chn]=tmp;
 707       }
 708
 709       /* block type is based just on L or R channel */
 710       if (chn<2) {
 711         uselongblock[chn] = 1;
 712
 713         /* tuned for t1.wav.  doesnt effect most other samples */
 714         if (gfc->pe[chn] > 3000)
 715             uselongblock[chn]=0;
 716
 717         if ( mx > 30*mn )
 718           {/* big surge of energy - always use short blocks */
 719             uselongblock[chn] = 0;
 720           }
 721         else if ((mx > 10*mn) && (gfc->pe[chn] > 1000))
 722           {/* medium surge, medium pe - use short blocks */
 723             uselongblock[chn] = 0;
 724           }
 725
 726         /* disable short blocks */
 727         if (gfp->short_blocks == short_block_dispensed)
 728           uselongblock[chn]=1;
 729         if (gfp->short_blocks == short_block_forced)
 730           uselongblock[chn]=0;
 731       }
 732     }
 733
 734 #if defined(HAVE_GTK)
 735     if (gfp->analysis) {
 736       FLOAT mn,mx,ma=0,mb=0,mc=0;
 737
 738       for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++)
 739       {
 740         ma += gfc->energy_s[0][j];
 741         mb += gfc->energy_s[1][j];
 742         mc += gfc->energy_s[2][j];
 743       }
 744       mn = Min(ma,mb);
 745       mn = Min(mn,mc);
 746       mx = Max(ma,mb);
 747       mx = Max(mx,mc);
 748
 749       gfc->pinfo->ers[gr_out][chn]=gfc->ers_save[chn];
 750       gfc->ers_save[chn]=(mx/(1e-12+mn));
 751       gfc->pinfo->pe[gr_out][chn]=gfc->pe_save[chn];
 752       gfc->pe_save[chn]=gfc->pe[chn];
 753     }
 754 #endif
 755
 756
 757     /***************************************************************
 758      * compute masking thresholds for both short and long blocks
 759      ***************************************************************/
 760     /* longblock threshold calculation (part 2) */
 761     for ( sb = 0; sb < NBPSY_l; sb++ )
 762       {
 763         FLOAT8 enn = gfc->w1_l[sb] * eb[gfc->bu_l[sb]] + gfc->w2_l[sb] * eb[gfc->bo_l[sb]];
 764         FLOAT8 thmm = gfc->w1_l[sb] *thr[gfc->bu_l[sb]] + gfc->w2_l[sb] * thr[gfc->bo_l[sb]];
 765
 766         for ( b = gfc->bu_l[sb]+1; b < gfc->bo_l[sb]; b++ )
 767           {
 768             enn  += eb[b];
 769             thmm += thr[b];
 770           }
 771
 772         gfc->en [chn].l[sb] = enn;
 773         gfc->thm[chn].l[sb] = thmm;
 774       }
 775
 776
 777     /* threshold calculation for short blocks */
 778     for ( sblock = 0; sblock < 3; sblock++ )
 779       {
 780         j = 0;
 781         for ( b = 0; b < gfc->npart_s_orig; b++ )
 782           {
 783             FLOAT ecb = gfc->energy_s[sblock][j++];
 784             for (i = 1 ; i<gfc->numlines_s[b]; ++i)
 785               {
 786                 ecb += gfc->energy_s[sblock][j++];
 787               }
 788             eb[b] = ecb;
 789           }
 790         {
 791           int kk = 0;
 792           for ( b = 0; b < gfc->npart_s; b++ )
 793             {
 794               FLOAT8 ecb = 0;
 795               for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ ) {
 796                 ecb += gfc->s3_ss[kk++] * eb[k];
 797               }
 798               ecb *= gfc->SNR_s[b];
 799
 800 /* 2001-07-13 */
 801         if ( gfp->VBR == vbr_off || gfp->VBR == vbr_abr ) {
 802             /* this looks like a BUG to me. robert */
 803             thr[b] = Max (1e-6, ecb);
 804         }
 805         else {
 806             thr[b] = Min( ecb, rpelev_s  * gfc->nb_s1[chn][b] );
 807             if (gfc->blocktype_old[chn>1 ? chn-2 : chn] == SHORT_TYPE ) {
 808                 thr[b] = Min( thr[b], rpelev2_s * gfc->nb_s2[chn][b] );
 809             }
 810             thr[b] = Max( thr[b], 1e-37 );
 811             gfc->nb_s2[chn][b] = gfc->nb_s1[chn][b];
 812             gfc->nb_s1[chn][b] = ecb;
 813         }
 814
 815             }
 816
 817           for ( sb = 0; sb < NBPSY_s; sb++ )
 818             {
 819               FLOAT8 enn  = gfc->w1_s[sb] * eb[gfc->bu_s[sb]] + gfc->w2_s[sb] * eb[gfc->bo_s[sb]];
 820               FLOAT8 thmm = gfc->w1_s[sb] *thr[gfc->bu_s[sb]] + gfc->w2_s[sb] * thr[gfc->bo_s[sb]];
 821
 822               for ( b = gfc->bu_s[sb]+1; b < gfc->bo_s[sb]; b++ ) {
 823                 enn  += eb[b];
 824                 thmm += thr[b];
 825               }
 826
 827               gfc->en [chn].s[sb][sblock] = enn;
 828               gfc->thm[chn].s[sb][sblock] = thmm;
 829             }
 830         }
 831       }
 832   } /* end loop over chn */
 833
 834
 835
 836   /***************************************************************
 837    * compute M/S thresholds
 838    ***************************************************************/
 839   /* compute M/S thresholds from Johnston & Ferreira 1992 ICASSP paper */
 840   if (gfp->mode == JOINT_STEREO) {
 841     FLOAT8 rside,rmid,mld;
 842     int chmid=2,chside=3;
 843
 844     for ( sb = 0; sb < NBPSY_l; sb++ ) {
 845       /* use this fix if L & R masking differs by 2db or less */
 846       /* if db = 10*log10(x2/x1) < 2 */
 847       /* if (x2 < 1.58*x1) { */
 848       if (gfc->thm[0].l[sb] <= 1.58*gfc->thm[1].l[sb]
 849           && gfc->thm[1].l[sb] <= 1.58*gfc->thm[0].l[sb]) {
 850
 851         mld = gfc->mld_l[sb]*gfc->en[chside].l[sb];
 852         rmid = Max(gfc->thm[chmid].l[sb], Min(gfc->thm[chside].l[sb],mld));
 853
 854         mld = gfc->mld_l[sb]*gfc->en[chmid].l[sb];
 855         rside = Max(gfc->thm[chside].l[sb],Min(gfc->thm[chmid].l[sb],mld));
 856
 857         gfc->thm[chmid].l[sb]=rmid;
 858         gfc->thm[chside].l[sb]=rside;
 859       }
 860     }
 861     for ( sb = 0; sb < NBPSY_s; sb++ ) {
 862       for ( sblock = 0; sblock < 3; sblock++ ) {
 863         if (gfc->thm[0].s[sb][sblock] <= 1.58*gfc->thm[1].s[sb][sblock]
 864             && gfc->thm[1].s[sb][sblock] <= 1.58*gfc->thm[0].s[sb][sblock]) {
 865
 866           mld = gfc->mld_s[sb]*gfc->en[chside].s[sb][sblock];
 867           rmid = Max(gfc->thm[chmid].s[sb][sblock],Min(gfc->thm[chside].s[sb][sblock],mld));
 868
 869           mld = gfc->mld_s[sb]*gfc->en[chmid].s[sb][sblock];
 870           rside = Max(gfc->thm[chside].s[sb][sblock],Min(gfc->thm[chmid].s[sb][sblock],mld));
 871
 872           gfc->thm[chmid].s[sb][sblock]=rmid;
 873           gfc->thm[chside].s[sb][sblock]=rside;
 874         }
 875       }
 876     }
 877   }
 878
 879
 880
 881   /***************************************************************
 882    * Adjust M/S maskings if user set "msfix"
 883    ***************************************************************/
 884   /* Naoki Shibata 2000 */
 885   if (numchn == 4 && gfp->msfix!=0) {
 886       FLOAT msfix = gfp->msfix;
 887
 888     for ( sb = 0; sb < NBPSY_l; sb++ )
 889       {
 890         FLOAT8 thmL,thmR,thmM,thmS,ath;
 891         ath  = (gfc->ATH->cb[(gfc->bu_l[sb] + gfc->bo_l[sb])/2])*pow(10,-gfp->ATHlower/10.0);
 892         thmL = Max(gfc->thm[0].l[sb],ath);
 893         thmR = Max(gfc->thm[1].l[sb],ath);
 894         thmM = Max(gfc->thm[2].l[sb],ath);
 895         thmS = Max(gfc->thm[3].l[sb],ath);
 896
 897         if (thmL*msfix < (thmM+thmS)/2) {
 898           FLOAT8 f = thmL*msfix / ((thmM+thmS)/2);
 899           thmM *= f;
 900           thmS *= f;
 901         }
 902         if (thmR*msfix < (thmM+thmS)/2) {
 903           FLOAT8 f = thmR*msfix / ((thmM+thmS)/2);
 904           thmM *= f;
 905           thmS *= f;
 906         }
 907
 908         gfc->thm[2].l[sb] = Min(thmM,gfc->thm[2].l[sb]);
 909         gfc->thm[3].l[sb] = Min(thmS,gfc->thm[3].l[sb]);
 910       }
 911
 912     for ( sb = 0; sb < NBPSY_s; sb++ ) {
 913       for ( sblock = 0; sblock < 3; sblock++ ) {
 914         FLOAT8 thmL,thmR,thmM,thmS,ath;
 915         ath  = (gfc->ATH->cb[(gfc->bu_s[sb] + gfc->bo_s[sb])/2])*pow(10,-gfp->ATHlower/10.0);
 916         thmL = Max(gfc->thm[0].s[sb][sblock],ath);
 917         thmR = Max(gfc->thm[1].s[sb][sblock],ath);
 918         thmM = Max(gfc->thm[2].s[sb][sblock],ath);
 919         thmS = Max(gfc->thm[3].s[sb][sblock],ath);
 920
 921         if (thmL*msfix < (thmM+thmS)/2) {
 922           FLOAT8 f = thmL*msfix / ((thmM+thmS)/2);
 923           thmM *= f;
 924           thmS *= f;
 925         }
 926         if (thmR*msfix < (thmM+thmS)/2) {
 927           FLOAT8 f = thmR*msfix / ((thmM+thmS)/2);
 928           thmM *= f;
 929           thmS *= f;
 930         }
 931
 932         gfc->thm[2].s[sb][sblock] = Min(gfc->thm[2].s[sb][sblock],thmM);
 933         gfc->thm[3].s[sb][sblock] = Min(gfc->thm[3].s[sb][sblock],thmS);
 934       }
 935     }
 936   }
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948   if (gfp->mode == JOINT_STEREO)  {
 949     /* determin ms_ratio from masking thresholds*/
 950     /* use ms_stereo (ms_ratio < .35) if average thresh. diff < 5 db */
 951     FLOAT8 db,x1,x2,sidetot=0,tot=0;
 952     for (sb= NBPSY_l/4 ; sb< NBPSY_l; sb ++ ) {
 953       x1 = Min(gfc->thm[0].l[sb],gfc->thm[1].l[sb]);
 954       x2 = Max(gfc->thm[0].l[sb],gfc->thm[1].l[sb]);
 955       /* thresholds difference in db */
 956       if (x2 >= 1000*x1)  db=3;
 957       else db = log10(x2/x1);
 958       /*  DEBUGF(gfc,"db = %f %e %e  \n",db,gfc->thm[0].l[sb],gfc->thm[1].l[sb]);*/
 959       sidetot += db;
 960       tot++;
 961     }
 962     ms_ratio_l= (sidetot/tot)*0.7; /* was .35*(sidetot/tot)/5.0*10 */
 963     ms_ratio_l = Min(ms_ratio_l,0.5);
 964
 965     sidetot=0; tot=0;
 966     for ( sblock = 0; sblock < 3; sblock++ )
 967       for ( sb = NBPSY_s/4; sb < NBPSY_s; sb++ ) {
 968         x1 = Min(gfc->thm[0].s[sb][sblock],gfc->thm[1].s[sb][sblock]);
 969         x2 = Max(gfc->thm[0].s[sb][sblock],gfc->thm[1].s[sb][sblock]);
 970         /* thresholds difference in db */
 971         if (x2 >= 1000*x1)  db=3;
 972         else db = log10(x2/x1);
 973         sidetot += db;
 974         tot++;
 975       }
 976     ms_ratio_s = (sidetot/tot)*0.7; /* was .35*(sidetot/tot)/5.0*10 */
 977     ms_ratio_s = Min(ms_ratio_s,.5);
 978   }
 979
 980   /***************************************************************
 981    * determine final block type
 982    ***************************************************************/
 983
 984   for (chn=0; chn<gfc->channels_out; chn++) {
 985     blocktype[chn] = NORM_TYPE;
 986   }
 987
 988
 989   if (gfp->short_blocks == short_block_coupled) {
 990       /* force both channels to use the same block type */
 991       /* this is necessary if the frame is to be encoded in ms_stereo.  */
 992       /* But even without ms_stereo, FhG  does this */
 993       int bothlong= (uselongblock[0] && uselongblock[1]);
 994       if (!bothlong) {
 995         uselongblock[0]=0;
 996         uselongblock[1]=0;
 997       }
 998   }
 999
1000
1001
1002   /* update the blocktype of the previous granule, since it depends on what
1003    * happend in this granule */
1004   for (chn=0; chn<gfc->channels_out; chn++) {
1005     if ( uselongblock[chn])
1006       {                         /* no attack : use long blocks */
1007         assert( gfc->blocktype_old[chn] != START_TYPE );
1008         switch( gfc->blocktype_old[chn] )
1009           {
1010           case NORM_TYPE:
1011           case STOP_TYPE:
1012             blocktype[chn] = NORM_TYPE;
1013             break;
1014           case SHORT_TYPE:
1015             blocktype[chn] = STOP_TYPE;
1016             break;
1017           }
1018       } else   {
1019         /* attack : use short blocks */
1020         blocktype[chn] = SHORT_TYPE;
1021         if ( gfc->blocktype_old[chn] == NORM_TYPE ) {
1022           gfc->blocktype_old[chn] = START_TYPE;
1023         }
1024         if ( gfc->blocktype_old[chn] == STOP_TYPE ) {
1025           gfc->blocktype_old[chn] = SHORT_TYPE ;
1026         }
1027       }
1028
1029     blocktype_d[chn] = gfc->blocktype_old[chn];  /* value returned to calling program */
1030     gfc->blocktype_old[chn] = blocktype[chn];    /* save for next call to l3psy_anal */
1031   }
1032
1033   if (blocktype_d[0]==2 && blocktype_d[1]==2)
1034     *ms_ratio = gfc->ms_ratio_s_old;
1035   else
1036     *ms_ratio = gfc->ms_ratio_l_old;
1037
1038   gfc->ms_ratio_s_old = ms_ratio_s;
1039   gfc->ms_ratio_l_old = ms_ratio_l;
1040
1041   /* we dont know the block type of this frame yet - assume long */
1042   *ms_ratio_next = ms_ratio_l;
1043
1044   return 0;
1045 }
1046
1047 /* addition of simultaneous masking   Naoki Shibata 2000/7 */
1048 inline static FLOAT8 mask_add(FLOAT8 m1,FLOAT8 m2,int k,int b, lame_internal_flags * const gfc)
1049 {
1050   static const FLOAT8 table1[] = {
1051     3.3246 *3.3246 ,3.23837*3.23837,3.15437*3.15437,3.00412*3.00412,2.86103*2.86103,2.65407*2.65407,2.46209*2.46209,2.284  *2.284  ,
1052     2.11879*2.11879,1.96552*1.96552,1.82335*1.82335,1.69146*1.69146,1.56911*1.56911,1.46658*1.46658,1.37074*1.37074,1.31036*1.31036,
1053     1.25264*1.25264,1.20648*1.20648,1.16203*1.16203,1.12765*1.12765,1.09428*1.09428,1.0659 *1.0659 ,1.03826*1.03826,1.01895*1.01895,
1054     1
1055   };
1056
1057   static const FLOAT8 table2[] = {
1058     1.33352*1.33352,1.35879*1.35879,1.38454*1.38454,1.39497*1.39497,1.40548*1.40548,1.3537 *1.3537 ,1.30382*1.30382,1.22321*1.22321,
1059     1.14758*1.14758
1060   };
1061
1062   static const FLOAT8 table3[] = {
1063     2.35364*2.35364,2.29259*2.29259,2.23313*2.23313,2.12675*2.12675,2.02545*2.02545,1.87894*1.87894,1.74303*1.74303,1.61695*1.61695,
1064     1.49999*1.49999,1.39148*1.39148,1.29083*1.29083,1.19746*1.19746,1.11084*1.11084,1.03826*1.03826
1065   };
1066
1067
1068   int i;
1069   FLOAT8 m;
1070
1071   if (m1 == 0) return m2;
1072
1073   if (b < 0) b = -b;
1074
1075   i = 10*log10(m2 / m1)/10*16;
1076   m = 10*log10((m1+m2)/gfc->ATH->cb[k]);
1077
1078   if (i < 0) i = -i;
1079
1080   if (b <= 3) {  /* approximately, 1 bark = 3 partitions */
1081     if (i > 8) return m1+m2;
1082     return (m1+m2)*table2[i];
1083   }
1084
1085   if (m<15) {
1086     if (m > 0) {
1087       FLOAT8 f=1.0,r;
1088       if (i > 24) return m1+m2;
1089       if (i > 13) f = 1; else f = table3[i];
1090       r = (m-0)/15;
1091       return (m1+m2)*(table1[i]*r+f*(1-r));
1092     }
1093     if (i > 13) return m1+m2;
1094     return (m1+m2)*table3[i];
1095   }
1096
1097   if (i > 24) return m1+m2;
1098   return (m1+m2)*table1[i];
1099 }
1100
1101 int L3psycho_anal_ns( lame_global_flags * gfp,
1102                     const sample_t *buffer[2], int gr_out,
1103                     FLOAT8 *ms_ratio,
1104                     FLOAT8 *ms_ratio_next,
1105                     III_psy_ratio masking_ratio[2][2],
1106                     III_psy_ratio masking_MS_ratio[2][2],
1107                     FLOAT8 percep_entropy[2],FLOAT8 percep_MS_entropy[2],
1108                     FLOAT8 energy[4],
1109                     int blocktype_d[2])
1110 {
1111 /* to get a good cache performance, one has to think about
1112  * the sequence, in which the variables are used.
1113  * (Note: these static variables have been moved to the gfc-> struct,
1114  * and their order in memory is layed out in util.h)
1115  */
1116   lame_internal_flags *gfc=gfp->internal_flags;
1117
1118   /* fft and energy calculation   */
1119   FLOAT (*wsamp_l)[BLKSIZE];
1120   FLOAT (*wsamp_s)[3][BLKSIZE_s];
1121
1122   /* convolution   */
1123   FLOAT8 eb[CBANDS],eb2[CBANDS];
1124   FLOAT8 thr[CBANDS];
1125
1126   FLOAT8 max[CBANDS],avg[CBANDS],tonality2[CBANDS];
1127
1128   /* ratios    */
1129   FLOAT8 ms_ratio_l=0,ms_ratio_s=0;
1130
1131   /* block type  */
1132   int blocktype[2],uselongblock[2];
1133
1134   /* usual variables like loop indices, etc..    */
1135   int numchn, chn;
1136   int   b, i, j, k;
1137   int   sb,sblock;
1138
1139   /* variables used for --nspsytune */
1140   int ns_attacks[4];
1141   FLOAT ns_hpfsmpl[4][576+576/3+NSFIRLEN];
1142   FLOAT pe_l[4],pe_s[4];
1143   FLOAT pcfact;
1144
1145
1146   if(gfc->psymodel_init==0) {
1147     psymodel_init(gfp);
1148     init_fft(gfc);
1149     gfc->psymodel_init=1;
1150   }
1151
1152
1153   numchn = gfc->channels_out;
1154   /* chn=2 and 3 = Mid and Side channels */
1155   if (gfp->mode == JOINT_STEREO) numchn=4;
1156
1157   if (gfp->VBR==vbr_off) pcfact = gfc->ResvMax == 0 ? 0 : ((FLOAT)gfc->ResvSize)/gfc->ResvMax*0.5;
1158   else if (gfp->VBR == vbr_rh  ||  gfp->VBR == vbr_mtrh  ||  gfp->VBR == vbr_mt) {
1159     static const FLOAT8 pcQns[10]={1.0,1.0,1.0,0.8,0.6,0.5,0.4,0.3,0.2,0.1};
1160     pcfact = pcQns[gfp->VBR_q];
1161   } else pcfact = 1;
1162
1163   /**********************************************************************
1164    *  Apply HPF of fs/4 to the input signal.
1165    *  This is used for attack detection / handling.
1166    **********************************************************************/
1167
1168   {
1169     static const FLOAT fircoef[] = {
1170       -8.65163e-18,-0.00851586,-6.74764e-18, 0.0209036,
1171       -3.36639e-17,-0.0438162 ,-1.54175e-17, 0.0931738,
1172       -5.52212e-17,-0.313819  , 0.5        ,-0.313819,
1173       -5.52212e-17, 0.0931738 ,-1.54175e-17,-0.0438162,
1174       -3.36639e-17, 0.0209036 ,-6.74764e-18,-0.00851586,
1175       -8.65163e-18,
1176     };
1177
1178     for(chn=0;chn<gfc->channels_out;chn++)
1179       {
1180         FLOAT firbuf[576+576/3+NSFIRLEN];
1181
1182         /* apply high pass filter of fs/4 */
1183
1184         for(i=-NSFIRLEN;i<576+576/3;i++)
1185           firbuf[i+NSFIRLEN] = buffer[chn][576-350+(i)];
1186
1187         for(i=0;i<576+576/3-NSFIRLEN;i++)
1188           {
1189             FLOAT sum = 0;
1190             for(j=0;j<NSFIRLEN;j++)
1191               sum += fircoef[j] * firbuf[i+j];
1192             ns_hpfsmpl[chn][i] = sum;
1193           }
1194         for(;i<576+576/3;i++)
1195           ns_hpfsmpl[chn][i] = 0;
1196       }
1197
1198     if (gfp->mode == JOINT_STEREO) {
1199       for(i=0;i<576+576/3;i++)
1200         {
1201           ns_hpfsmpl[2][i] = ns_hpfsmpl[0][i]+ns_hpfsmpl[1][i];
1202           ns_hpfsmpl[3][i] = ns_hpfsmpl[0][i]-ns_hpfsmpl[1][i];
1203         }
1204     }
1205   }
1206
1207
1208
1209   /* there is a one granule delay.  Copy maskings computed last call
1210    * into masking_ratio to return to calling program.
1211    */
1212   for (chn=0; chn<numchn; chn++) {
1213       for (i=0; i<numchn; ++i) {
1214           energy[chn]=gfc->tot_ener[chn];
1215       }
1216   }
1217
1218
1219   for (chn=0; chn<numchn; chn++) {
1220     /* there is a one granule delay.  Copy maskings computed last call
1221      * into masking_ratio to return to calling program.
1222      */
1223     pe_l[chn] = gfc->nsPsy.pe_l[chn];
1224     pe_s[chn] = gfc->nsPsy.pe_s[chn];
1225
1226     if (chn < 2) {
1227       /* LR maskings  */
1228       //percep_entropy            [chn]       = gfc -> pe  [chn];
1229       masking_ratio    [gr_out] [chn]  .en  = gfc -> en  [chn];
1230       masking_ratio    [gr_out] [chn]  .thm = gfc -> thm [chn];
1231     } else {
1232       /* MS maskings  */
1233       //percep_MS_entropy         [chn-2]     = gfc -> pe  [chn];
1234       masking_MS_ratio [gr_out] [chn-2].en  = gfc -> en  [chn];
1235       masking_MS_ratio [gr_out] [chn-2].thm = gfc -> thm [chn];
1236     }
1237   }
1238
1239   for (chn=0; chn<numchn; chn++) {
1240
1241     /**********************************************************************
1242      *  compute FFTs
1243      **********************************************************************/
1244
1245     wsamp_s = gfc->wsamp_S+(chn & 1);
1246     wsamp_l = gfc->wsamp_L+(chn & 1);
1247
1248     if (chn<2) {
1249       fft_long ( gfc, *wsamp_l, chn, buffer);
1250       fft_short( gfc, *wsamp_s, chn, buffer);
1251     }
1252
1253     /* FFT data for mid and side channel is derived from L & R */
1254
1255     if (chn == 2)
1256       {
1257         for (j = BLKSIZE-1; j >=0 ; --j)
1258           {
1259             FLOAT l = gfc->wsamp_L[0][j];
1260             FLOAT r = gfc->wsamp_L[1][j];
1261             gfc->wsamp_L[0][j] = (l+r)*(FLOAT)(SQRT2*0.5);
1262             gfc->wsamp_L[1][j] = (l-r)*(FLOAT)(SQRT2*0.5);
1263           }
1264         for (b = 2; b >= 0; --b)
1265           {
1266             for (j = BLKSIZE_s-1; j >= 0 ; --j)
1267               {
1268                 FLOAT l = gfc->wsamp_S[0][b][j];
1269                 FLOAT r = gfc->wsamp_S[1][b][j];
1270                 gfc->wsamp_S[0][b][j] = (l+r)*(FLOAT)(SQRT2*0.5);
1271                 gfc->wsamp_S[1][b][j] = (l-r)*(FLOAT)(SQRT2*0.5);
1272               }
1273           }
1274       }
1275
1276
1277     /**********************************************************************
1278      *  compute energies for each spectral line
1279      **********************************************************************/
1280
1281     /* long block */
1282
1283     gfc->energy[0]  = (*wsamp_l)[0];
1284     gfc->energy[0] *= gfc->energy[0];
1285
1286     gfc->tot_ener[chn] = gfc->energy[0]; /* sum total energy at nearly no extra cost */
1287
1288     for (j=BLKSIZE/2-1; j >= 0; --j)
1289     {
1290       FLOAT re = (*wsamp_l)[BLKSIZE/2-j];
1291       FLOAT im = (*wsamp_l)[BLKSIZE/2+j];
1292       gfc->energy[BLKSIZE/2-j] = (re * re + im * im) * (FLOAT)0.5;
1293
1294       if (BLKSIZE/2-j > 10)
1295         gfc->tot_ener[chn] += gfc->energy[BLKSIZE/2-j];
1296     }
1297
1298
1299     /* short block */
1300
1301     for (b = 2; b >= 0; --b)
1302     {
1303       gfc->energy_s[b][0]  = (*wsamp_s)[b][0];
1304       gfc->energy_s[b][0] *=  gfc->energy_s [b][0];
1305       for (j=BLKSIZE_s/2-1; j >= 0; --j)
1306       {
1307         FLOAT re = (*wsamp_s)[b][BLKSIZE_s/2-j];
1308         FLOAT im = (*wsamp_s)[b][BLKSIZE_s/2+j];
1309         gfc->energy_s[b][BLKSIZE_s/2-j] = (re * re + im * im) * (FLOAT)0.5;
1310       }
1311     }
1312
1313
1314     /* output data for analysis */
1315
1316 #if defined(HAVE_GTK)
1317     if (gfp->analysis) {
1318       for (j=0; j<HBLKSIZE ; j++) {
1319         gfc->pinfo->energy[gr_out][chn][j]=gfc->energy_save[chn][j];
1320         gfc->energy_save[chn][j]=gfc->energy[j];
1321       }
1322     }
1323 #endif
1324
1325
1326     /**********************************************************************
1327      * compute loudness approximation (used for ATH auto-level adjustment)
1328      **********************************************************************/
1329     if( gfp->athaa_loudapprox == 2 ) {
1330       if( chn < 2 ) {           /* no loudness for mid and side channels */
1331         gfc->loudness_sq[gr_out][chn] = gfc->loudness_sq_save[chn];
1332         gfc->loudness_sq_save[chn]
1333           = psycho_loudness_approx( gfc->energy, gfp);
1334       }
1335     }
1336
1337
1338
1339     /**********************************************************************
1340      *    Calculate the energy and the tonality of each partition.
1341      **********************************************************************/
1342
1343     for (b=0, j=0; b<gfc->npart_l_orig; b++)
1344       {
1345         FLOAT8 ebb,m,a;
1346
1347         ebb = gfc->energy[j];
1348         m = a = gfc->energy[j];
1349         j++;
1350
1351         for (i = gfc->numlines_l[b] - 1; i > 0; i--)
1352           {
1353             FLOAT8 el = gfc->energy[j];
1354             ebb += gfc->energy[j];
1355             a += el;
1356             m = m < el ? el : m;
1357             j++;
1358           }
1359         eb[b] = ebb;
1360         max[b] = m;
1361         avg[b] = a / gfc->numlines_l[b];
1362       }
1363
1364     j = 0;
1365     for (b=0; b < gfc->npart_l_orig; b++ )
1366       {
1367         int c1,c2;
1368         FLOAT8 m,a;
1369         tonality2[b] = 0;
1370         c1 = c2 = 0;
1371         m = a = 0;
1372         for(k=b-1;k<=b+1;k++)
1373           {
1374             if (k >= 0 && k < gfc->npart_l_orig) {
1375               c1++;
1376               c2 += gfc->numlines_l[k];
1377               a += avg[k];
1378               m = m < max[k] ? max[k] : m;
1379             }
1380           }
1381
1382         a /= c1;
1383         tonality2[b] = a == 0 ? 0 : (m / a - 1)/(c2-1);
1384       }
1385
1386     for (b=0; b < gfc->npart_l_orig; b++ )
1387       {
1388 #if 0
1389         static FLOAT8 tab[20] =
1390         {  0,  1,  2,  2,  2,  2,  2,  6,9.3,9.3,9.3,9.3,9.3,9.3,9.3,9.3,9.3,9.3,9.3,9.3};
1391
1392         static int init = 1;
1393         if (init) {
1394           int j;
1395           for(j=0;j<20;j++) {
1396             tab[j] = pow(10.0,-tab[j]/10.0);
1397           }
1398           init = 0;
1399         }
1400 #else
1401         static FLOAT8 tab[20] = {
1402           1,0.79433,0.63096,0.63096,0.63096,0.63096,0.63096,0.25119,0.11749,0.11749,
1403           0.11749,0.11749,0.11749,0.11749,0.11749,0.11749,0.11749,0.11749,0.11749,0.11749
1404         };
1405 #endif
1406
1407         int t = 20*tonality2[b];
1408         if (t > 19) t = 19;
1409         eb2[b] = eb[b] * tab[t];
1410       }
1411
1412
1413     /**********************************************************************
1414      *      convolve the partitioned energy and unpredictability
1415      *      with the spreading function, s3_l[b][k]
1416      ******************************************************************** */
1417     {
1418     int kk = 0;
1419     for ( b = 0;b < gfc->npart_l; b++ )
1420       {
1421         FLOAT8 ecb;
1422
1423         /****   convolve the partitioned energy with the spreading function   ****/
1424
1425         ecb = 0;
1426
1427 #if 1
1428         for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ )
1429           {
1430             ecb = mask_add(ecb,gfc->s3_ll[kk++] * eb2[k],k,k-b,gfc);
1431           }
1432
1433         ecb *= 0.158489319246111; // pow(10,-0.8)
1434 #endif
1435
1436 #if 0
1437         for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ )
1438           {
1439             ecb += gfc->s3_ll[kk++] * eb2[k];
1440           }
1441
1442         ecb *= 0.223872113856834; // pow(10,-0.65);
1443 #endif
1444
1445         /****   long block pre-echo control   ****/
1446
1447         /* dont use long block pre-echo control if previous granule was
1448          * a short block.  This is to avoid the situation:
1449          * frame0:  quiet (very low masking)
1450          * frame1:  surge  (triggers short blocks)
1451          * frame2:  regular frame.  looks like pre-echo when compared to
1452          *          frame0, but all pre-echo was in frame1.
1453          */
1454
1455         /* chn=0,1   L and R channels
1456            chn=2,3   S and M channels.
1457         */
1458
1459 #define NS_INTERP(x,y,r) (pow((x),(r))*pow((y),1-(r)))
1460
1461         if (gfc->blocktype_old[chn>1 ? chn-2 : chn] == SHORT_TYPE )
1462           thr[b] = ecb; /* Min(ecb, rpelev*gfc->nb_1[chn][b]); */
1463         else
1464           thr[b] = NS_INTERP(Min(ecb, Min(rpelev*gfc->nb_1[chn][b],rpelev2*gfc->nb_2[chn][b])),ecb,pcfact);
1465
1466         gfc->nb_2[chn][b] = gfc->nb_1[chn][b];
1467         gfc->nb_1[chn][b] = ecb;
1468       }
1469     }
1470
1471     /***************************************************************
1472      * determine the block type (window type)
1473      ***************************************************************/
1474
1475     {
1476       static int count=0;
1477       FLOAT en_subshort[12];
1478       FLOAT attack_intensity[12];
1479       int ns_uselongblock = 1;
1480
1481       /* calculate energies of each sub-shortblocks */
1482
1483       k = 0;
1484       for(i=0;i<12;i++)
1485         {
1486           en_subshort[i] = 0;
1487           for(j=0;j<576/9;j++)
1488             {
1489               en_subshort[i] += ns_hpfsmpl[chn][k] * ns_hpfsmpl[chn][k];
1490               k++;
1491             }
1492
1493           if (en_subshort[i] < 100) en_subshort[i] = 100;
1494         }
1495
1496       /* compare energies between sub-shortblocks */
1497
1498 #define NSATTACKTHRE 150
1499 #define NSATTACKTHRE_S 300
1500
1501       for(i=0;i<2;i++) {
1502         attack_intensity[i] = en_subshort[i] / gfc->nsPsy.last_en_subshort[chn][7+i];
1503       }
1504
1505       for(;i<12;i++) {
1506         attack_intensity[i] = en_subshort[i] / en_subshort[i-2];
1507       }
1508
1509       ns_attacks[0] = ns_attacks[1] = ns_attacks[2] = ns_attacks[3] = 0;
1510
1511       for(i=0;i<12;i++)
1512         {
1513           if (!ns_attacks[i/3] && attack_intensity[i] > (chn == 3 ? (gfc->presetTune.use ? gfc->presetTune.attackthre_s : NSATTACKTHRE_S)
1514                                                               : (gfc->presetTune.use ? gfc->presetTune.attackthre   : NSATTACKTHRE))) ns_attacks[i/3] = (i % 3)+1;
1515         }
1516
1517       if (ns_attacks[0] && gfc->nsPsy.last_attacks[chn][2]) ns_attacks[0] = 0;
1518       if (ns_attacks[1] && ns_attacks[0]) ns_attacks[1] = 0;
1519       if (ns_attacks[2] && ns_attacks[1]) ns_attacks[2] = 0;
1520       if (ns_attacks[3] && ns_attacks[2]) ns_attacks[3] = 0;
1521
1522       if (gfc->nsPsy.last_attacks[chn][2] == 3 ||
1523           ns_attacks[0] || ns_attacks[1] || ns_attacks[2] || ns_attacks[3]) ns_uselongblock = 0;
1524
1525       if (chn < 4) count++;
1526
1527       for(i=0;i<9;i++)
1528         {
1529           gfc->nsPsy.last_en_subshort[chn][i] = en_subshort[i];
1530           gfc->nsPsy.last_attack_intensity[chn][i] = attack_intensity[i];
1531         }
1532
1533       if (gfp->short_blocks == short_block_dispensed) {
1534         uselongblock[chn] = 1;
1535       }
1536       else if (gfp->short_blocks == short_block_forced) {
1537         uselongblock[chn] = 0;
1538       }
1539       else {
1540         if (chn < 2) {
1541           uselongblock[chn] = ns_uselongblock;
1542         } else {
1543           if (ns_uselongblock == 0) uselongblock[0] = uselongblock[1] = 0;
1544         }
1545       }
1546     }
1547
1548 #if defined(HAVE_GTK)
1549     if (gfp->analysis) {
1550       FLOAT mn,mx,ma=0,mb=0,mc=0;
1551
1552       for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++)
1553       {
1554         ma += gfc->energy_s[0][j];
1555         mb += gfc->energy_s[1][j];
1556         mc += gfc->energy_s[2][j];
1557       }
1558       mn = Min(ma,mb);
1559       mn = Min(mn,mc);
1560       mx = Max(ma,mb);
1561       mx = Max(mx,mc);
1562
1563       gfc->pinfo->ers[gr_out][chn]=gfc->ers_save[chn];
1564       gfc->ers_save[chn]=(mx/(1e-12+mn));
1565     }
1566 #endif
1567
1568
1569     /***************************************************************
1570      * compute masking thresholds for long blocks
1571      ***************************************************************/
1572
1573     for ( sb = 0; sb < NBPSY_l; sb++ )
1574       {
1575         FLOAT8 enn = gfc->w1_l[sb] * eb[gfc->bu_l[sb]] + gfc->w2_l[sb] * eb[gfc->bo_l[sb]];
1576         FLOAT8 thmm = gfc->w1_l[sb] *thr[gfc->bu_l[sb]] + gfc->w2_l[sb] * thr[gfc->bo_l[sb]];
1577
1578         for ( b = gfc->bu_l[sb]+1; b < gfc->bo_l[sb]; b++ )
1579           {
1580             enn  += eb[b];
1581             thmm += thr[b];
1582           }
1583
1584         gfc->en [chn].l[sb] = enn;
1585         gfc->thm[chn].l[sb] = thmm;
1586       }
1587
1588
1589     /***************************************************************
1590      * compute masking thresholds for short blocks
1591      ***************************************************************/
1592
1593     for ( sblock = 0; sblock < 3; sblock++ )
1594       {
1595         j = 0;
1596         for ( b = 0; b < gfc->npart_s_orig; b++ )
1597           {
1598             FLOAT ecb = gfc->energy_s[sblock][j++];
1599             for (i = 1 ; i<gfc->numlines_s[b]; ++i)
1600               {
1601                 ecb += gfc->energy_s[sblock][j++];
1602               }
1603             eb[b] = ecb;
1604           }
1605
1606         {
1607           int kk = 0;
1608           for ( b = 0; b < gfc->npart_s; b++ )
1609             {
1610               FLOAT8 ecb = 0;
1611               for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ )
1612                 {
1613                   ecb += gfc->s3_ss[kk++] * eb[k];
1614                 }
1615
1616 /* 2001-07-13 */
1617             /* this looks like a BUG  */
1618             thr[b] = Max (1e-6, ecb);
1619
1620             if (gfp->VBR == vbr_mtrh) {
1621                 thr[b] = Min( ecb, rpelev_s  * gfc->nb_s1[chn][b] );
1622                 if (gfc->blocktype_old[chn>1 ? chn-2 : chn] == SHORT_TYPE ) {
1623                    thr[b] = Min( thr[b], rpelev2_s * gfc->nb_s2[chn][b] );
1624                 }
1625                 thr[b] = Max( thr[b], 1e-37 );
1626                 gfc->nb_s2[chn][b] = gfc->nb_s1[chn][b];
1627                 gfc->nb_s1[chn][b] = ecb;
1628                 }
1629             }
1630         }
1631         for ( sb = 0; sb < NBPSY_s; sb++ )
1632           {
1633             FLOAT8 enn  = gfc->w1_s[sb] * eb[gfc->bu_s[sb]] + gfc->w2_s[sb] * eb[gfc->bo_s[sb]];
1634             FLOAT8 thmm = gfc->w1_s[sb] *thr[gfc->bu_s[sb]] + gfc->w2_s[sb] * thr[gfc->bo_s[sb]];
1635
1636             for ( b = gfc->bu_s[sb]+1; b < gfc->bo_s[sb]; b++ )
1637               {
1638                 enn  += eb[b];
1639                 thmm += thr[b];
1640               }
1641
1642             /****   short block pre-echo control   ****/
1643
1644 #define NS_PREECHO_ATT0 0.8
1645 #define NS_PREECHO_ATT1 0.6
1646 #define NS_PREECHO_ATT2 0.3
1647
1648             thmm *= NS_PREECHO_ATT0;
1649
1650             if (ns_attacks[sblock] >= 2) {
1651               if (sblock != 0) {
1652                 double p = NS_INTERP(gfc->thm[chn].s[sb][sblock-1],thmm,NS_PREECHO_ATT1*pcfact);
1653                 thmm = Min(thmm,p);
1654               } else {
1655                 double p = NS_INTERP(gfc->nsPsy.last_thm[chn][sb][2],thmm,NS_PREECHO_ATT1*pcfact);
1656                 thmm = Min(thmm,p);
1657               }
1658             } else if (ns_attacks[sblock+1] == 1) {
1659               if (sblock != 0) {
1660                 double p = NS_INTERP(gfc->thm[chn].s[sb][sblock-1],thmm,NS_PREECHO_ATT1*pcfact);
1661                 thmm = Min(thmm,p);
1662               } else {
1663                 double p = NS_INTERP(gfc->nsPsy.last_thm[chn][sb][2],thmm,NS_PREECHO_ATT1*pcfact);
1664                 thmm = Min(thmm,p);
1665               }
1666             }
1667
1668             if (ns_attacks[sblock] == 1) {
1669               double p = sblock == 0 ? gfc->nsPsy.last_thm[chn][sb][2] : gfc->thm[chn].s[sb][sblock-1];
1670               p = NS_INTERP(p,thmm,NS_PREECHO_ATT2*pcfact);
1671               thmm = Min(thmm,p);
1672             } else if ((sblock != 0 && ns_attacks[sblock-1] == 3) ||
1673                        (sblock == 0 && gfc->nsPsy.last_attacks[chn][2] == 3)) {
1674               double p = sblock <= 1 ? gfc->nsPsy.last_thm[chn][sb][sblock+1] : gfc->thm[chn].s[sb][0];
1675               p = NS_INTERP(p,thmm,NS_PREECHO_ATT2*pcfact);
1676               thmm = Min(thmm,p);
1677             }
1678
1679             gfc->en [chn].s[sb][sblock] = enn;
1680             gfc->thm[chn].s[sb][sblock] = thmm;
1681           }
1682       }
1683
1684
1685     /***************************************************************
1686      * save some values for analysis of the next granule
1687      ***************************************************************/
1688
1689     for ( sblock = 0; sblock < 3; sblock++ )
1690       {
1691         for ( sb = 0; sb < NBPSY_s; sb++ )
1692           {
1693             gfc->nsPsy.last_thm[chn][sb][sblock] = gfc->thm[chn].s[sb][sblock];
1694           }
1695       }
1696
1697     for(i=0;i<3;i++)
1698       gfc->nsPsy.last_attacks[chn][i] = ns_attacks[i];
1699
1700   } /* end loop over chn */
1701
1702
1703
1704   /***************************************************************
1705    * compute M/S thresholds
1706    ***************************************************************/
1707
1708   /* from Johnston & Ferreira 1992 ICASSP paper */
1709
1710   if ( numchn==4 /* mid/side and r/l */) {
1711     FLOAT8 rside,rmid,mld;
1712     int chmid=2,chside=3;
1713
1714     for ( sb = 0; sb < NBPSY_l; sb++ ) {
1715       /* use this fix if L & R masking differs by 2db or less */
1716       /* if db = 10*log10(x2/x1) < 2 */
1717       /* if (x2 < 1.58*x1) { */
1718       if (gfc->thm[0].l[sb] <= 1.58*gfc->thm[1].l[sb]
1719           && gfc->thm[1].l[sb] <= 1.58*gfc->thm[0].l[sb]) {
1720
1721         mld = gfc->mld_l[sb]*gfc->en[chside].l[sb];
1722         rmid = Max(gfc->thm[chmid].l[sb], Min(gfc->thm[chside].l[sb],mld));
1723
1724         mld = gfc->mld_l[sb]*gfc->en[chmid].l[sb];
1725         rside = Max(gfc->thm[chside].l[sb],Min(gfc->thm[chmid].l[sb],mld));
1726
1727         gfc->thm[chmid].l[sb]=rmid;
1728         gfc->thm[chside].l[sb]=rside;
1729       }
1730     }
1731     for ( sb = 0; sb < NBPSY_s; sb++ ) {
1732       for ( sblock = 0; sblock < 3; sblock++ ) {
1733         if (gfc->thm[0].s[sb][sblock] <= 1.58*gfc->thm[1].s[sb][sblock]
1734             && gfc->thm[1].s[sb][sblock] <= 1.58*gfc->thm[0].s[sb][sblock]) {
1735
1736           mld = gfc->mld_s[sb]*gfc->en[chside].s[sb][sblock];
1737           rmid = Max(gfc->thm[chmid].s[sb][sblock],Min(gfc->thm[chside].s[sb][sblock],mld));
1738
1739           mld = gfc->mld_s[sb]*gfc->en[chmid].s[sb][sblock];
1740           rside = Max(gfc->thm[chside].s[sb][sblock],Min(gfc->thm[chmid].s[sb][sblock],mld));
1741
1742           gfc->thm[chmid].s[sb][sblock]=rmid;
1743           gfc->thm[chside].s[sb][sblock]=rside;
1744         }
1745       }
1746     }
1747   }
1748
1749
1750   /* Naoki Shibata 2000 */
1751
1752 #define NS_MSFIX 3.5
1753
1754   if (numchn == 4) {
1755     FLOAT msfix = NS_MSFIX;
1756     if (gfc->nsPsy.safejoint) msfix = 1;
1757     if (gfp->msfix) msfix = gfp->msfix;
1758
1759     if (gfc->presetTune.use && gfc->ATH->adjust >=
1760                                        gfc->presetTune.athadjust_switch_level &&
1761                                                            gfc->presetTune.athadjust_msfix > 0)
1762                 msfix = gfc->presetTune.athadjust_msfix;
1763
1764     for ( sb = 0; sb < NBPSY_l; sb++ )
1765       {
1766         FLOAT8 thmL,thmR,thmM,thmS,ath;
1767         ath  = (gfc->ATH->cb[(gfc->bu_l[sb] + gfc->bo_l[sb])/2])*pow(10,-gfp->ATHlower/10.0);
1768         thmL = Max(gfc->thm[0].l[sb],ath);
1769         thmR = Max(gfc->thm[1].l[sb],ath);
1770         thmM = Max(gfc->thm[2].l[sb],ath);
1771         thmS = Max(gfc->thm[3].l[sb],ath);
1772
1773         if (thmL*msfix < (thmM+thmS)/2) {
1774           FLOAT8 f = thmL * (gfc->presetTune.use ? gfc->presetTune.ms_maskadjust : msfix) / ((thmM+thmS)/2);
1775           thmM *= f;
1776           thmS *= f;
1777         }
1778         if (thmR*msfix < (thmM+thmS)/2) {
1779           FLOAT8 f = thmR * (gfc->presetTune.use ? gfc->presetTune.ms_maskadjust : msfix) / ((thmM+thmS)/2);
1780           thmM *= f;
1781           thmS *= f;
1782         }
1783
1784         gfc->thm[2].l[sb] = Min(thmM,gfc->thm[2].l[sb]);
1785         gfc->thm[3].l[sb] = Min(thmS,gfc->thm[3].l[sb]);
1786       }
1787
1788     for ( sb = 0; sb < NBPSY_s; sb++ ) {
1789       for ( sblock = 0; sblock < 3; sblock++ ) {
1790         FLOAT8 thmL,thmR,thmM,thmS,ath;
1791         ath  = (gfc->ATH->cb[(gfc->bu_s[sb] + gfc->bo_s[sb])/2])*pow(10,-gfp->ATHlower/10.0);
1792         thmL = Max(gfc->thm[0].s[sb][sblock],ath);
1793         thmR = Max(gfc->thm[1].s[sb][sblock],ath);
1794         thmM = Max(gfc->thm[2].s[sb][sblock],ath);
1795         thmS = Max(gfc->thm[3].s[sb][sblock],ath);
1796
1797         if (thmL*msfix < (thmM+thmS)/2) {
1798           FLOAT8 f = thmL*msfix / ((thmM+thmS)/2);
1799           thmM *= f;
1800           thmS *= f;
1801         }
1802         if (thmR*msfix < (thmM+thmS)/2) {
1803           FLOAT8 f = thmR*msfix / ((thmM+thmS)/2);
1804           thmM *= f;
1805           thmS *= f;
1806         }
1807
1808         gfc->thm[2].s[sb][sblock] = Min(gfc->thm[2].s[sb][sblock],thmM);
1809         gfc->thm[3].s[sb][sblock] = Min(gfc->thm[3].s[sb][sblock],thmS);
1810       }
1811     }
1812   }
1813
1814   ms_ratio_l = 0;
1815   ms_ratio_s = 0;
1816
1817
1818   /***************************************************************
1819    * compute estimation of the amount of bit used in the granule
1820    ***************************************************************/
1821
1822   for(chn=0;chn<numchn;chn++)
1823     {
1824       {
1825         static FLOAT8 regcoef[] = {
1826           1124.23,10.0583,10.7484,7.29006,16.2714,6.2345,4.09743,3.05468,3.33196,2.54688,
1827           3.68168,5.83109,2.93817,-8.03277,-10.8458,8.48777,9.13182,2.05212,8.6674,50.3937,73.267,97.5664,0
1828         };
1829
1830         FLOAT8 msum = regcoef[0]/4;
1831         int sb;
1832
1833         for ( sb = 0; sb < NBPSY_l; sb++ )
1834           {
1835             FLOAT8 t;
1836
1837             if (gfc->thm[chn].l[sb]*gfc->masking_lower != 0 &&
1838                 gfc->en[chn].l[sb]/(gfc->thm[chn].l[sb]*gfc->masking_lower) > 1)
1839               t = log(gfc->en[chn].l[sb]/(gfc->thm[chn].l[sb]*gfc->masking_lower));
1840             else
1841               t = 0;
1842             msum += regcoef[sb+1] * t;
1843           }
1844
1845         gfc->nsPsy.pe_l[chn] = msum;
1846       }
1847
1848       {
1849         static FLOAT8 regcoef[] = {
1850           1236.28,0,0,0,0.434542,25.0738,0,0,0,19.5442,19.7486,60,100,0
1851         };
1852
1853         FLOAT8 msum = regcoef[0]/4;
1854         int sb,sblock;
1855
1856         for(sblock=0;sblock<3;sblock++)
1857           {
1858             for ( sb = 0; sb < NBPSY_s; sb++ )
1859               {
1860                 FLOAT8 t;
1861
1862                 if (gfc->thm[chn].s[sb][sblock] * gfc->masking_lower != 0 &&
1863                     gfc->en[chn].s[sb][sblock] / (gfc->thm[chn].s[sb][sblock] * gfc->masking_lower) > 1)
1864                   t = log(gfc->en[chn].s[sb][sblock] / (gfc->thm[chn].s[sb][sblock] * gfc->masking_lower));
1865                 else
1866                   t = 0;
1867                 msum += regcoef[sb+1] * t;
1868               }
1869           }
1870
1871         gfc->nsPsy.pe_s[chn] = msum;
1872       }
1873
1874       //gfc->pe[chn] -= 150;
1875     }
1876
1877   /***************************************************************
1878    * determine final block type
1879    ***************************************************************/
1880
1881   for (chn=0; chn<gfc->channels_out; chn++) {
1882     blocktype[chn] = NORM_TYPE;
1883   }
1884
1885   if (gfp->short_blocks == short_block_coupled) {
1886       /* force both channels to use the same block type */
1887       /* this is necessary if the frame is to be encoded in ms_stereo.  */
1888       /* But even without ms_stereo, FhG  does this */
1889       int bothlong= (uselongblock[0] && uselongblock[1]);
1890       if (!bothlong) {
1891         uselongblock[0]=0;
1892         uselongblock[1]=0;
1893       }
1894   }
1895
1896   /* update the blocktype of the previous granule, since it depends on what
1897    * happend in this granule */
1898   for (chn=0; chn<gfc->channels_out; chn++) {
1899     if ( uselongblock[chn])
1900       {                         /* no attack : use long blocks */
1901         assert( gfc->blocktype_old[chn] != START_TYPE );
1902         switch( gfc->blocktype_old[chn] )
1903           {
1904           case NORM_TYPE:
1905           case STOP_TYPE:
1906             blocktype[chn] = NORM_TYPE;
1907             break;
1908           case SHORT_TYPE:
1909             blocktype[chn] = STOP_TYPE;
1910             break;
1911           }
1912       } else   {
1913         /* attack : use short blocks */
1914         blocktype[chn] = SHORT_TYPE;
1915         if ( gfc->blocktype_old[chn] == NORM_TYPE ) {
1916           gfc->blocktype_old[chn] = START_TYPE;
1917         }
1918         if ( gfc->blocktype_old[chn] == STOP_TYPE ) {
1919           gfc->blocktype_old[chn] = SHORT_TYPE ;
1920         }
1921       }
1922
1923     blocktype_d[chn] = gfc->blocktype_old[chn];  /* value returned to calling program */
1924     gfc->blocktype_old[chn] = blocktype[chn];    /* save for next call to l3psy_anal */
1925
1926     if (gfc->presetTune.use) {
1927         if (blocktype_d[chn] != NORM_TYPE)
1928             gfc->presetTune.quantcomp_current = gfc->presetTune.quantcomp_type_s;
1929         else
1930             gfc->presetTune.quantcomp_current = gfp->experimentalX;
1931
1932         if (gfc->ATH->adjust >= gfc->presetTune.athadjust_switch_level &&
1933                                                      blocktype_d[chn] == NORM_TYPE &&
1934                                                                                  gfc->presetTune.quantcomp_alt_type > -1) {
1935                     gfc->presetTune.quantcomp_current = gfc->presetTune.quantcomp_alt_type;
1936                 }
1937     }
1938   }
1939
1940   /*********************************************************************
1941    * compute the value of PE to return (one granule delay)
1942    *********************************************************************/
1943   for(chn=0;chn<numchn;chn++)
1944     {
1945       if (chn < 2) {
1946         if (blocktype_d[chn] == SHORT_TYPE) {
1947           percep_entropy[chn] = pe_s[chn];
1948         } else {
1949           percep_entropy[chn] = pe_l[chn];
1950         }
1951         if (gfp->analysis) gfc->pinfo->pe[gr_out][chn] = percep_entropy[chn];
1952       } else {
1953         if (blocktype_d[0] == SHORT_TYPE) {
1954           percep_MS_entropy[chn-2] = pe_s[chn];
1955         } else {
1956           percep_MS_entropy[chn-2] = pe_l[chn];
1957         }
1958         if (gfp->analysis) gfc->pinfo->pe[gr_out][chn] = percep_MS_entropy[chn-2];
1959       }
1960     }
1961
1962   return 0;
1963 }
1964
1965
1966
1967
1968
1969 /*
1970  *   The spreading function.  Values returned in units of energy
1971  */
1972 FLOAT8 s3_func(FLOAT8 bark) {
1973
1974     FLOAT8 tempx,x,tempy,temp;
1975     tempx = bark;
1976     if (tempx>=0) tempx *= 3;
1977     else tempx *=1.5;
1978
1979     if(tempx>=0.5 && tempx<=2.5)
1980         {
1981             temp = tempx - 0.5;
1982             x = 8.0 * (temp*temp - 2.0 * temp);
1983         }
1984     else x = 0.0;
1985     tempx += 0.474;
1986     tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx);
1987
1988     if (tempy <= -60.0) return  0.0;
1989
1990     tempx = exp( (x + tempy)*LN_TO_LOG10 );
1991
1992     /* Normalization.  The spreading function should be normalized so that:
1993          +inf
1994            /
1995            |  s3 [ bark ]  d(bark)   =  1
1996            /
1997          -inf
1998     */
1999     tempx /= .6609193;
2000     return tempx;
2001
2002 }
2003
2004
2005
2006
2007
2008
2009
2010
2011 int L3para_read(lame_global_flags * gfp, FLOAT8 sfreq, int *numlines_l,int *numlines_s,
2012 FLOAT8 *minval,
2013 FLOAT8 s3_l[CBANDS][CBANDS], FLOAT8 s3_s[CBANDS][CBANDS],
2014 FLOAT8 *SNR,
2015 int *bu_l, int *bo_l, FLOAT8 *w1_l, FLOAT8 *w2_l,
2016 int *bu_s, int *bo_s, FLOAT8 *w1_s, FLOAT8 *w2_s,
2017 int *npart_l_orig,int *npart_l,int *npart_s_orig,int *npart_s)
2018 {
2019   lame_internal_flags *gfc=gfp->internal_flags;
2020
2021
2022   FLOAT8 bval_l[CBANDS], bval_s[CBANDS];
2023   FLOAT8 bval_l_width[CBANDS], bval_s_width[CBANDS];
2024   int  i,j;
2025   int partition[HBLKSIZE];
2026
2027
2028
2029   /* compute numlines, the number of spectral lines in each partition band */
2030   /* each partition band should be about DELBARK wide. */
2031   j=0;
2032   for(i=0;i<CBANDS;i++)
2033     {
2034       FLOAT8 ji, bark1,bark2;
2035       int k,j2;
2036
2037       j2 = j;
2038       j2 = Min(j2,BLKSIZE/2);
2039
2040       do {
2041         /* find smallest j2 >= j so that  (bark - bark_l[i-1]) < DELBARK */
2042         ji = j;
2043         bark1 = freq2bark(sfreq*ji/BLKSIZE);
2044
2045         ++j2;
2046         ji = j2;
2047         bark2  = freq2bark(sfreq*ji/BLKSIZE);
2048       } while ((bark2 - bark1) < DELBARK  && j2<=BLKSIZE/2);
2049
2050       for (k=j; k<j2; ++k)
2051         partition[k]=i;
2052       numlines_l[i]=(j2-j);
2053       j = j2;
2054       if (j > BLKSIZE/2) break;
2055     }
2056   *npart_l_orig = i+1;
2057   assert(*npart_l_orig <= CBANDS);
2058
2059   /* compute which partition bands are in which scalefactor bands */
2060   { int i1,i2,sfb,start,end;
2061     FLOAT8 freq1,freq2;
2062     for ( sfb = 0; sfb < SBMAX_l; sfb++ ) {
2063       start = gfc->scalefac_band.l[ sfb ];
2064       end   = gfc->scalefac_band.l[ sfb+1 ];
2065       freq1 = sfreq*(start-.5)/(2*576);
2066       freq2 = sfreq*(end-1+.5)/(2*576);
2067
2068       i1 = floor(.5 + BLKSIZE*freq1/sfreq);
2069       if (i1<0) i1=0;
2070       i2 = floor(.5 + BLKSIZE*freq2/sfreq);
2071       if (i2>BLKSIZE/2) i2=BLKSIZE/2;
2072
2073       //      DEBUGF(gfc,"longblock:  old: (%i,%i)  new: (%i,%i) %i %i \n",bu_l[sfb],bo_l[sfb],partition[i1],partition[i2],i1,i2);
2074
2075       w1_l[sfb]=.5;
2076       w2_l[sfb]=.5;
2077       bu_l[sfb]=partition[i1];
2078       bo_l[sfb]=partition[i2];
2079
2080     }
2081   }
2082
2083
2084   /* compute bark value and ATH of each critical band */
2085   j = 0;
2086   for ( i = 0; i < *npart_l_orig; i++ ) {
2087       int     k;
2088       FLOAT8  bark1,bark2;
2089       /* FLOAT8 mval,freq; */
2090
2091       // Calculating the medium bark scaled frequency of the spectral lines
2092       // from j ... j + numlines[i]-1  (=spectral lines in parition band i)
2093
2094       k         = numlines_l[i] - 1;
2095       bark1 = freq2bark(sfreq*(j+0)/BLKSIZE);
2096       bark2 = freq2bark(sfreq*(j+k)/BLKSIZE);
2097       bval_l[i] = .5*(bark1+bark2);
2098
2099       bark1 = freq2bark(sfreq*(j+0-.5)/BLKSIZE);
2100       bark2 = freq2bark(sfreq*(j+k+.5)/BLKSIZE);
2101       bval_l_width[i] = bark2-bark1;
2102
2103       gfc->ATH->cb [i] = 1.e37; // preinit for minimum search
2104       for (k=0; k < numlines_l[i]; k++, j++) {
2105         FLOAT8  freq = sfreq*j/(1000.0*BLKSIZE);
2106         FLOAT8  level;
2107         assert( freq <= 24 );              // or only '<'
2108         //      freq = Min(.1,freq);       // ATH below 100 Hz constant, not further climbing
2109         level  = ATHformula (freq*1000, gfp) - 20;   // scale to FFT units; returned value is in dB
2110         level  = pow ( 10., 0.1*level );   // convert from dB -> energy
2111         level *= numlines_l [i];
2112         if ( level < gfc->ATH->cb [i] )
2113             gfc->ATH->cb [i] = level;
2114       }
2115
2116
2117     }
2118
2119   /* MINVAL.  For low freq, the strength of the masking is limited by minval
2120    * this is an ISO MPEG1 thing, dont know if it is really needed */
2121   for(i=0;i<*npart_l_orig;i++){
2122     double x = (-20+bval_l[i]*20.0/10.0);
2123     if (bval_l[i]>10) x = 0;
2124     minval[i]=pow(10.0,x/10);
2125     gfc->PSY->prvTonRed[i] = minval[i];
2126   }
2127
2128
2129
2130
2131
2132
2133
2134   /************************************************************************/
2135   /* SHORT BLOCKS */
2136   /************************************************************************/
2137
2138   /* compute numlines */
2139   j=0;
2140   for(i=0;i<CBANDS;i++)
2141     {
2142       FLOAT8 ji, bark1,bark2;
2143       int k,j2;
2144
2145       j2 = j;
2146       j2 = Min(j2,BLKSIZE_s/2);
2147
2148       do {
2149         /* find smallest j2 >= j so that  (bark - bark_s[i-1]) < DELBARK */
2150         ji = j;
2151         bark1  = freq2bark(sfreq*ji/BLKSIZE_s);
2152
2153         ++j2;
2154         ji = j2;
2155         bark2  = freq2bark(sfreq*ji/BLKSIZE_s);
2156
2157       } while ((bark2 - bark1) < DELBARK  && j2<=BLKSIZE_s/2);
2158
2159       for (k=j; k<j2; ++k)
2160         partition[k]=i;
2161       numlines_s[i]=(j2-j);
2162       j = j2;
2163       if (j > BLKSIZE_s/2) break;
2164     }
2165   *npart_s_orig = i+1;
2166   assert(*npart_s_orig <= CBANDS);
2167
2168   /* compute which partition bands are in which scalefactor bands */
2169   { int i1,i2,sfb,start,end;
2170     FLOAT8 freq1,freq2;
2171     for ( sfb = 0; sfb < SBMAX_s; sfb++ ) {
2172       start = gfc->scalefac_band.s[ sfb ];
2173       end   = gfc->scalefac_band.s[ sfb+1 ];
2174       freq1 = sfreq*(start-.5)/(2*192);
2175       freq2 = sfreq*(end-1+.5)/(2*192);
2176
2177       i1 = floor(.5 + BLKSIZE_s*freq1/sfreq);
2178       if (i1<0) i1=0;
2179       i2 = floor(.5 + BLKSIZE_s*freq2/sfreq);
2180       if (i2>BLKSIZE_s/2) i2=BLKSIZE_s/2;
2181
2182       //DEBUGF(gfc,"shortblock: old: (%i,%i)  new: (%i,%i) %i %i \n",bu_s[sfb],bo_s[sfb], partition[i1],partition[i2],i1,i2);
2183
2184       w1_s[sfb]=.5;
2185       w2_s[sfb]=.5;
2186       bu_s[sfb]=partition[i1];
2187       bo_s[sfb]=partition[i2];
2188
2189     }
2190   }
2191
2192
2193
2194
2195
2196   /* compute bark values of each critical band */
2197   j = 0;
2198   for(i=0;i<*npart_s_orig;i++)
2199     {
2200       int     k;
2201       FLOAT8  bark1,bark2,snr;
2202       k    = numlines_s[i] - 1;
2203
2204       bark1 = freq2bark (sfreq*(j+0)/BLKSIZE_s);
2205       bark2 = freq2bark (sfreq*(j+k)/BLKSIZE_s);
2206       bval_s[i] = .5*(bark1+bark2);
2207
2208       bark1 = freq2bark (sfreq*(j+0-.5)/BLKSIZE_s);
2209       bark2 = freq2bark (sfreq*(j+k+.5)/BLKSIZE_s);
2210       bval_s_width[i] = bark2-bark1;
2211       j        += k+1;
2212
2213       /* SNR formula */
2214       if (bval_s[i]<13)
2215           snr=-8.25;
2216       else
2217           snr  = -4.5 * (bval_s[i]-13)/(24.0-13.0)  +
2218               -8.25*(bval_s[i]-24)/(13.0-24.0);
2219
2220       SNR[i]=pow(10.0,snr/10.0);
2221     }
2222
2223
2224
2225
2226
2227
2228   /************************************************************************
2229    * Now compute the spreading function, s[j][i], the value of the spread-*
2230    * ing function, centered at band j, for band i, store for later use    *
2231    ************************************************************************/
2232   /* i.e.: sum over j to spread into signal barkval=i
2233      NOTE: i and j are used opposite as in the ISO docs */
2234   for(i=0;i<*npart_l_orig;i++)    {
2235       for(j=0;j<*npart_l_orig;j++)      {
2236           s3_l[i][j]=s3_func(bval_l[i]-bval_l[j])*bval_l_width[j];
2237       }
2238   }
2239   for(i=0;i<*npart_s_orig;i++)     {
2240       for(j=0;j<*npart_s_orig;j++)      {
2241           s3_s[i][j]=s3_func(bval_s[i]-bval_s[j])*bval_s_width[j];
2242       }
2243   }
2244
2245
2246
2247
2248   /* compute: */
2249   /* npart_l_orig   = number of partition bands before convolution */
2250   /* npart_l  = number of partition bands after convolution */
2251
2252   *npart_l=bo_l[NBPSY_l-1]+1;
2253   *npart_s=bo_s[NBPSY_s-1]+1;
2254
2255   assert(*npart_l <= *npart_l_orig);
2256   assert(*npart_s <= *npart_s_orig);
2257
2258
2259     /* setup stereo demasking thresholds */
2260     /* formula reverse enginerred from plot in paper */
2261     for ( i = 0; i < NBPSY_s; i++ ) {
2262       FLOAT8 arg,mld;
2263       arg = freq2bark(sfreq*gfc->scalefac_band.s[i]/(2*192));
2264       arg = (Min(arg, 15.5)/15.5);
2265
2266       mld = 1.25*(1-cos(PI*arg))-2.5;
2267       gfc->mld_s[i] = pow(10.0,mld);
2268     }
2269     for ( i = 0; i < NBPSY_l; i++ ) {
2270       FLOAT8 arg,mld;
2271       arg = freq2bark(sfreq*gfc->scalefac_band.l[i]/(2*576));
2272       arg = (Min(arg, 15.5)/15.5);
2273
2274       mld = 1.25*(1-cos(PI*arg))-2.5;
2275       gfc->mld_l[i] = pow(10.0,mld);
2276     }
2277
2278 #define temporalmask_sustain_sec 0.01
2279
2280     /* setup temporal masking */
2281     gfc->decay = exp(-1.0*LOG10/(temporalmask_sustain_sec*sfreq/192.0));
2282
2283     return 0;
2284 }
2285
2286
2287
2288
2289
2290
2291
2292
2293
2294
2295
2296 int psymodel_init(lame_global_flags *gfp)
2297 {
2298     lame_internal_flags *gfc=gfp->internal_flags;
2299     int i,j,b,sb,k,samplerate;
2300
2301     FLOAT8 s3_s[CBANDS][CBANDS];
2302     FLOAT8 s3_l[CBANDS][CBANDS];
2303     int numberOfNoneZero;
2304
2305     samplerate = gfp->out_samplerate;
2306     gfc->ms_ener_ratio_old=.25;
2307     gfc->blocktype_old[0]=STOP_TYPE;
2308     gfc->blocktype_old[1]=STOP_TYPE;
2309     gfc->blocktype_old[0]=SHORT_TYPE;
2310     gfc->blocktype_old[1]=SHORT_TYPE;
2311
2312     for (i=0; i<4; ++i) {
2313       for (j=0; j<CBANDS; ++j) {
2314         gfc->nb_1[i][j]=1e20;
2315         gfc->nb_2[i][j]=1e20;
2316       }
2317       for ( sb = 0; sb < NBPSY_l; sb++ ) {
2318         gfc->en[i].l[sb] = 1e20;
2319         gfc->thm[i].l[sb] = 1e20;
2320       }
2321       for (j=0; j<3; ++j) {
2322         for ( sb = 0; sb < NBPSY_s; sb++ ) {
2323           gfc->en[i].s[sb][j] = 1e20;
2324           gfc->thm[i].s[sb][j] = 1e20;
2325         }
2326       }
2327     }
2328     for (i=0; i<4; ++i) {
2329       for (j=0; j<3; ++j) {
2330         for ( sb = 0; sb < NBPSY_s; sb++ ) {
2331           gfc->nsPsy.last_thm[i][sb][j] = 1e20;
2332         }
2333       }
2334     }
2335     for(i=0;i<4;i++) {
2336       for(j=0;j<9;j++)
2337         gfc->nsPsy.last_en_subshort[i][j] = 100;
2338       for(j=0;j<3;j++)
2339         gfc->nsPsy.last_attacks[i][j] = 0;
2340       gfc->nsPsy.pe_l[i] = gfc->nsPsy.pe_s[i] = 0;
2341     }
2342
2343
2344
2345
2346     /*  gfp->cwlimit = sfreq*j/1024.0;  */
2347     gfc->cw_lower_index=6;
2348     gfc->cw_upper_index = gfc->PSY->cwlimit*1024.0/gfp->out_samplerate;
2349     gfc->cw_upper_index=Min(HBLKSIZE-4,gfc->cw_upper_index);      /* j+3 < HBLKSIZE-1 */
2350     gfc->cw_upper_index=Max(6,gfc->cw_upper_index);
2351
2352     for ( j = 0; j < HBLKSIZE; j++ )
2353       gfc->cw[j] = 0.4f;
2354
2355
2356
2357     i=L3para_read( gfp,(FLOAT8) samplerate,gfc->numlines_l,gfc->numlines_s,
2358           gfc->minval,s3_l,s3_s,gfc->SNR_s,gfc->bu_l,
2359           gfc->bo_l,gfc->w1_l,gfc->w2_l, gfc->bu_s,gfc->bo_s,
2360           gfc->w1_s,gfc->w2_s,&gfc->npart_l_orig,&gfc->npart_l,
2361           &gfc->npart_s_orig,&gfc->npart_s );
2362     if (i!=0) return -1;
2363
2364
2365
2366     /* npart_l_orig   = number of partition bands before convolution */
2367     /* npart_l  = number of partition bands after convolution */
2368
2369     numberOfNoneZero = 0;
2370     for (i=0; i<gfc->npart_l; i++) {
2371       for (j = 0; j < gfc->npart_l_orig; j++) {
2372         if (s3_l[i][j] != 0.0)
2373           break;
2374       }
2375       gfc->s3ind[i][0] = j;
2376
2377       for (j = gfc->npart_l_orig - 1; j > 0; j--) {
2378         if (s3_l[i][j] != 0.0)
2379           break;
2380       }
2381       gfc->s3ind[i][1] = j;
2382       numberOfNoneZero += (gfc->s3ind[i][1] - gfc->s3ind[i][0] + 1);
2383     }
2384     gfc->s3_ll = malloc(sizeof(FLOAT8)*numberOfNoneZero);
2385     if (!gfc->s3_ll)
2386       return -1;
2387
2388     k = 0;
2389     for (i=0; i<gfc->npart_l; i++) {
2390       for (j = gfc->s3ind[i][0]; j <= gfc->s3ind[i][1]; j++) {
2391         gfc->s3_ll[k++] = s3_l[i][j];
2392       }
2393     }
2394
2395
2396
2397     numberOfNoneZero = 0;
2398     for (i=0; i<gfc->npart_s; i++) {
2399       for (j = 0; j < gfc->npart_s_orig; j++) {
2400         if (s3_s[i][j] != 0.0)
2401           break;
2402       }
2403       gfc->s3ind_s[i][0] = j;
2404
2405       for (j = gfc->npart_s_orig - 1; j > 0; j--) {
2406         if (s3_s[i][j] != 0.0)
2407           break;
2408       }
2409       gfc->s3ind_s[i][1] = j;
2410       numberOfNoneZero += (gfc->s3ind_s[i][1] - gfc->s3ind_s[i][0] + 1);
2411     }
2412     gfc->s3_ss = malloc(sizeof(FLOAT8)*numberOfNoneZero);
2413     if (!gfc->s3_ss)
2414       return -1;
2415
2416
2417
2418
2419     /*
2420       #include "debugscalefac.c"
2421     */
2422
2423
2424
2425     if (gfc->nsPsy.use) {
2426         /* long block spreading function normalization */
2427         for ( b = 0;b < gfc->npart_l; b++ ) {
2428             for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ ) {
2429                 // spreading function has been properly normalized by
2430                 // multiplying by DELBARK/.6609193 = .515.
2431                 // It looks like Naoki was
2432                 // way ahead of me and added this factor here!
2433                 // it is no longer needed.
2434                 //gfc->s3_l[b][k] *= 0.5;
2435             }
2436         }
2437         /* short block spreading function normalization */
2438         // no longer needs to be normalized, but nspsytune wants
2439         // SNR_s applied here istead of later to save CPU cycles
2440         for ( b = 0;b < gfc->npart_s; b++ ) {
2441             FLOAT8 norm=0;
2442             for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ ) {
2443                 norm += s3_s[b][k];
2444             }
2445             for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ ) {
2446                 s3_s[b][k] *= gfc->SNR_s[b] /* / norm */;
2447             }
2448         }
2449     }
2450
2451
2452
2453     if (gfc->nsPsy.use) {
2454 #if 1
2455         /* spread only from npart_l bands.  Normally, we use the spreading
2456          * function to convolve from npart_l_orig down to npart_l bands
2457          */
2458         for(b=0;b<gfc->npart_l;b++)
2459             if (gfc->s3ind[b][1] > gfc->npart_l-1) gfc->s3ind[b][1] = gfc->npart_l-1;
2460 #endif
2461     }
2462     k = 0;
2463     for (i=0; i<gfc->npart_s; i++) {
2464       for (j = gfc->s3ind_s[i][0]; j <= gfc->s3ind_s[i][1]; j++) {
2465         gfc->s3_ss[k++] = s3_s[i][j];
2466       }
2467     }
2468
2469
2470
2471                                 /* init. for loudness approx. -jd 2001 mar 27*/
2472     gfc->loudness_sq_save[0] = 0.0;
2473     gfc->loudness_sq_save[1] = 0.0;
2474
2475
2476
2477     return 0;
2478 }
2479
2480
2481
2482
2483
2484 /* psycho_loudness_approx
2485    jd - 2001 mar 12
2486 in:  energy   - BLKSIZE/2 elements of frequency magnitudes ^ 2
2487      gfp      - uses out_samplerate, ATHtype (also needed for ATHformula)
2488 returns: loudness^2 approximation, a positive value roughly tuned for a value
2489          of 1.0 for signals near clipping.
2490 notes:   When calibrated, feeding this function binary white noise at sample
2491          values +32767 or -32768 should return values that approach 3.
2492          ATHformula is used to approximate an equal loudness curve.
2493 future:  Data indicates that the shape of the equal loudness curve varies
2494          with intensity.  This function might be improved by using an equal
2495          loudness curve shaped for typical playback levels (instead of the
2496          ATH, that is shaped for the threshold).  A flexible realization might
2497          simply bend the existing ATH curve to achieve the desired shape.
2498          However, the potential gain may not be enough to justify an effort.
2499 */
2500 FLOAT
2501 psycho_loudness_approx( FLOAT *energy, lame_global_flags *gfp )
2502 {
2503   int i;
2504   static int eql_type = -1;
2505   static FLOAT eql_w[BLKSIZE/2];/* equal loudness weights (based on ATH) */
2506   const FLOAT vo_scale= 1./( 14752 ); /* tuned for output level */
2507                                       /* (sensitive to energy scale) */
2508   FLOAT loudness_power;
2509
2510   if( eql_type != gfp->ATHtype ) {
2511                                 /* compute equal loudness weights (eql_w) */
2512     FLOAT freq;
2513     FLOAT freq_inc = gfp->out_samplerate / (BLKSIZE);
2514     FLOAT eql_balance = 0.0;
2515     eql_type = gfp->ATHtype;
2516     freq = 0.0;
2517     for( i = 0; i < BLKSIZE/2; ++i ) {
2518       freq += freq_inc;
2519                                 /* convert ATH dB to relative power (not dB) */
2520                                 /*  to determine eql_w */
2521       eql_w[i] = 1. / pow( 10, ATHformula( freq, gfp ) / 10 );
2522       eql_balance += eql_w[i];
2523     }
2524     eql_balance = 1 / eql_balance;
2525     for( i = BLKSIZE/2; --i >= 0; ) { /* scale weights */
2526       eql_w[i] *= eql_balance;
2527     }
2528   }
2529
2530   loudness_power = 0.0;
2531   for( i = 0; i < BLKSIZE/2; ++i ) { /* apply weights to power in freq. bands*/
2532     loudness_power += energy[i] * eql_w[i];
2533   }
2534   loudness_power /= (BLKSIZE/2);
2535   loudness_power *= vo_scale * vo_scale;
2536
2537   return( loudness_power );
2538 }