initial commit
[swftools.git] / 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:54:55 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 }