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