4 * Copyright (c) 1999 Mark Taylor
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.
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.
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.
22 /* $Id: psymodel.c,v 1.1 2002/04/28 17:30:23 kramm Exp $ */
29 This routine computes the psycho acoustics, delayed by
32 Input: buffer of PCM data (1024 samples).
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
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.
51 barks: a non-linear frequency scale. Mapping from frequency to
52 barks is given by freq2bark()
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
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.
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).
76 The general outline is as follows:
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.
85 Each partition band is considiered a "masker". The strength
86 of the i'th masker in band j is given by:
88 s3(bark(i)-bark(j))*strength(i)
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.
96 s3() is the "spreading function". It is given by a formula
97 determined via listening tests.
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."
103 masking(j) = sum_over_i s3(i-j)*strength(i) = s3 o strength
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.
108 Note: instead of a simple convolution, LAME also has the
109 option of using "additive masking"
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.
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.
125 Other data computed by the psy-model:
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
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
142 #include "config_static.h"
146 #include "psymodel.h"
162 /* size of each partition band, in barks: */
167 /* AAC values, results in more masking over MP3 values */
176 #define NBPSY_l (SBMAX_l)
177 #define NBPSY_s (SBMAX_s)
181 #define LN_TO_LOG10 (M_LN10/10)
183 #define LN_TO_LOG10 0.2302585093
187 psycho_loudness_approx( FLOAT *energy, lame_global_flags *gfp );
194 L3psycho_anal. Compute psycho acoustics.
196 Data returned to the calling program must be delayed by one
199 This is done in two places.
200 If we do not need to know the blocktype, the copying
201 can be done here at the top of the program: we copy the data for
202 the last granule (computed during the last call) before it is
203 overwritten with the new data. It looks like this:
205 0. static psymodel_data
206 1. calling_program_data = psymodel_data
207 2. compute psymodel_data
209 For data which needs to know the blocktype, the copying must be
210 done at the end of this loop, and the old values must be saved:
212 0. static psymodel_data_old
213 1. compute psymodel_data
214 2. compute possible block type of this granule
215 3. compute final block type of previous granule based on #2.
216 4. calling_program_data = psymodel_data_old
217 5. psymodel_data_old = psymodel_data
219 int L3psycho_anal( lame_global_flags * gfp,
220 const sample_t *buffer[2], int gr_out,
222 FLOAT8 *ms_ratio_next,
223 III_psy_ratio masking_ratio[2][2],
224 III_psy_ratio masking_MS_ratio[2][2],
225 FLOAT8 percep_entropy[2],FLOAT8 percep_MS_entropy[2],
229 /* to get a good cache performance, one has to think about
230 * the sequence, in which the variables are used.
231 * (Note: these static variables have been moved to the gfc-> struct,
232 * and their order in memory is layed out in util.h)
234 lame_internal_flags *gfc=gfp->internal_flags;
237 /* fft and energy calculation */
238 FLOAT (*wsamp_l)[BLKSIZE];
239 FLOAT (*wsamp_s)[3][BLKSIZE_s];
247 FLOAT8 ms_ratio_l=0,ms_ratio_s=0;
250 int blocktype[2],uselongblock[2];
252 /* usual variables like loop indices, etc.. */
258 if(gfc->psymodel_init==0) {
261 gfc->psymodel_init=1;
263 for (chn = 0; chn < 4; ++chn )
264 for (b = 0; b < CBANDS; ++b )
265 gfc->nb_s1[chn][b] = gfc->nb_s2[chn][b] = 1.0;
273 numchn = gfc->channels_out;
274 /* chn=2 and 3 = Mid and Side channels */
275 if (gfp->mode == JOINT_STEREO) numchn=4;
277 for (chn=0; chn<numchn; chn++) {
278 for (i=0; i<numchn; ++i) {
279 energy[chn]=gfc->tot_ener[chn];
283 for (chn=0; chn<numchn; chn++) {
285 /* there is a one granule delay. Copy maskings computed last call
286 * into masking_ratio to return to calling program.
290 percep_entropy [chn] = gfc -> pe [chn];
291 masking_ratio [gr_out] [chn] .en = gfc -> en [chn];
292 masking_ratio [gr_out] [chn] .thm = gfc -> thm [chn];
295 percep_MS_entropy [chn-2] = gfc -> pe [chn];
296 masking_MS_ratio [gr_out] [chn-2].en = gfc -> en [chn];
297 masking_MS_ratio [gr_out] [chn-2].thm = gfc -> thm [chn];
302 /**********************************************************************
304 **********************************************************************/
305 wsamp_s = gfc->wsamp_S+(chn & 1);
306 wsamp_l = gfc->wsamp_L+(chn & 1);
308 fft_long ( gfc, *wsamp_l, chn, buffer);
309 fft_short( gfc, *wsamp_s, chn, buffer);
311 /* FFT data for mid and side channel is derived from L & R */
314 for (j = BLKSIZE-1; j >=0 ; --j)
316 FLOAT l = gfc->wsamp_L[0][j];
317 FLOAT r = gfc->wsamp_L[1][j];
318 gfc->wsamp_L[0][j] = (l+r)*(FLOAT)(SQRT2*0.5);
319 gfc->wsamp_L[1][j] = (l-r)*(FLOAT)(SQRT2*0.5);
321 for (b = 2; b >= 0; --b)
323 for (j = BLKSIZE_s-1; j >= 0 ; --j)
325 FLOAT l = gfc->wsamp_S[0][b][j];
326 FLOAT r = gfc->wsamp_S[1][b][j];
327 gfc->wsamp_S[0][b][j] = (l+r)*(FLOAT)(SQRT2*0.5);
328 gfc->wsamp_S[1][b][j] = (l-r)*(FLOAT)(SQRT2*0.5);
334 /**********************************************************************
336 **********************************************************************/
340 gfc->energy[0] = (*wsamp_l)[0];
341 gfc->energy[0] *= gfc->energy[0];
343 gfc->tot_ener[chn] = gfc->energy[0]; /* sum total energy at nearly no extra cost */
345 for (j=BLKSIZE/2-1; j >= 0; --j)
347 FLOAT re = (*wsamp_l)[BLKSIZE/2-j];
348 FLOAT im = (*wsamp_l)[BLKSIZE/2+j];
349 gfc->energy[BLKSIZE/2-j] = (re * re + im * im) * (FLOAT)0.5;
351 if (BLKSIZE/2-j > 10)
352 gfc->tot_ener[chn] += gfc->energy[BLKSIZE/2-j];
354 for (b = 2; b >= 0; --b)
356 gfc->energy_s[b][0] = (*wsamp_s)[b][0];
357 gfc->energy_s[b][0] *= gfc->energy_s [b][0];
358 for (j=BLKSIZE_s/2-1; j >= 0; --j)
360 FLOAT re = (*wsamp_s)[b][BLKSIZE_s/2-j];
361 FLOAT im = (*wsamp_s)[b][BLKSIZE_s/2+j];
362 gfc->energy_s[b][BLKSIZE_s/2-j] = (re * re + im * im) * (FLOAT)0.5;
366 #if defined(HAVE_GTK)
368 for (j=0; j<HBLKSIZE ; j++) {
369 gfc->pinfo->energy[gr_out][chn][j]=gfc->energy_save[chn][j];
370 gfc->energy_save[chn][j]=gfc->energy[j];
377 /**********************************************************************
378 * compute loudness approximation (used for ATH auto-level adjustment)
379 **********************************************************************/
380 if( gfp->athaa_loudapprox == 2 ) {
381 if( chn < 2 ) { /* no loudness for mid and side channels */
382 gfc->loudness_sq[gr_out][chn] = gfc->loudness_sq_save[chn];
383 gfc->loudness_sq_save[chn]
384 = psycho_loudness_approx( gfc->energy, gfp);
390 /**********************************************************************
391 * compute unpredicatability of first six spectral lines *
392 **********************************************************************/
393 for ( j = 0; j < gfc->cw_lower_index; j++ )
394 { /* calculate unpredictability measure cw */
398 FLOAT numre, numim, den;
400 a2 = gfc-> ax_sav[chn][1][j];
401 b2 = gfc-> bx_sav[chn][1][j];
402 r2 = gfc-> rx_sav[chn][1][j];
403 a1 = gfc-> ax_sav[chn][1][j] = gfc-> ax_sav[chn][0][j];
404 b1 = gfc-> bx_sav[chn][1][j] = gfc-> bx_sav[chn][0][j];
405 r1 = gfc-> rx_sav[chn][1][j] = gfc-> rx_sav[chn][0][j];
406 an = gfc-> ax_sav[chn][0][j] = (*wsamp_l)[j];
407 bn = gfc-> bx_sav[chn][0][j] = j==0 ? (*wsamp_l)[0] : (*wsamp_l)[BLKSIZE-j];
408 rn = gfc-> rx_sav[chn][0][j] = sqrt(gfc->energy[j]);
410 { /* square (x1,y1) */
413 numim = (a1*a1-b1*b1)*(FLOAT)0.5;
422 { /* multiply by (x2,-y2) */
424 FLOAT tmp2 = (numim+numre)*(a2+b2)*(FLOAT)0.5;
425 FLOAT tmp1 = -a2*numre+tmp2;
426 numre = -b2*numim+tmp2;
434 { /* r-prime factor */
435 FLOAT tmp = (2*r1-r2)/den;
439 den=rn+fabs(2*r1-r2);
441 numre = (an+bn)*(FLOAT)0.5-numre;
442 numim = (an-bn)*(FLOAT)0.5-numim;
443 den = sqrt(numre*numre+numim*numim)/den;
450 /**********************************************************************
451 * compute unpredicatibility of next 200 spectral lines *
452 **********************************************************************/
453 for ( j = gfc->cw_lower_index; j < gfc->cw_upper_index; j += 4 )
454 {/* calculate unpredictability measure cw */
456 FLOAT numre, numim, den;
460 { /* square (x1,y1) */
461 r1 = gfc->energy_s[0][k];
463 FLOAT a1 = (*wsamp_s)[0][k];
464 FLOAT b1 = (*wsamp_s)[0][BLKSIZE_s-k]; /* k is never 0 */
466 numim = (a1*a1-b1*b1)*(FLOAT)0.5;
477 { /* multiply by (x2,-y2) */
478 r2 = gfc->energy_s[2][k];
480 FLOAT a2 = (*wsamp_s)[2][k];
481 FLOAT b2 = (*wsamp_s)[2][BLKSIZE_s-k];
484 FLOAT tmp2 = (numim+numre)*(a2+b2)*(FLOAT)0.5;
485 FLOAT tmp1 = -a2*numre+tmp2;
486 numre = -b2*numim+tmp2;
496 { /* r-prime factor */
497 FLOAT tmp = (2*r1-r2)/den;
502 rn = sqrt(gfc->energy_s[1][k]);
503 den=rn+fabs(2*r1-r2);
505 FLOAT an = (*wsamp_s)[1][k];
506 FLOAT bn = (*wsamp_s)[1][BLKSIZE_s-k];
507 numre = (an+bn)*(FLOAT)0.5-numre;
508 numim = (an-bn)*(FLOAT)0.5-numim;
509 den = sqrt(numre*numre+numim*numim)/den;
512 gfc->cw[j+1] = gfc->cw[j+2] = gfc->cw[j+3] = gfc->cw[j] = den;
516 for ( j = 14; j < HBLKSIZE-4; j += 4 )
517 {/* calculate energy from short ffts */
520 for (tot=0, sblock=0; sblock < 3; sblock++)
521 tot+=gfc->energy_s[sblock][k];
522 ave = gfc->energy[j+1]+ gfc->energy[j+2]+ gfc->energy[j+3]+ gfc->energy[j];
524 gfc->energy[j+1] = gfc->energy[j+2] = gfc->energy[j+3] = gfc->energy[j]=tot;
528 /**********************************************************************
529 * Calculate the energy and the unpredictability in the threshold *
530 * calculation partitions *
531 **********************************************************************/
534 for (j = 0; j < gfc->cw_upper_index && gfc->numlines_l[b] && b<gfc->npart_l_orig; )
538 ebb = gfc->energy[j];
539 cbb = gfc->energy[j] * gfc->cw[j];
542 for (i = gfc->numlines_l[b] - 1; i > 0; i--)
544 ebb += gfc->energy[j];
545 cbb += gfc->energy[j] * gfc->cw[j];
553 for (; b < gfc->npart_l_orig; b++ )
555 FLOAT8 ebb = gfc->energy[j++];
556 assert(gfc->numlines_l[b]);
557 for (i = gfc->numlines_l[b] - 1; i > 0; i--)
559 ebb += gfc->energy[j++];
565 /**********************************************************************
566 * convolve the partitioned energy and unpredictability *
567 * with the spreading function, s3_l[b][k](packed into s3_ll) *
568 ******************************************************************** */
569 gfc->pe[chn] = 0; /* calculate percetual entropy */
572 for ( b = 0;b < gfc->npart_l; b++ )
579 for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ )
581 ecb += gfc->s3_ll[kk] * eb[k]; /* sprdngf for Layer III */
582 ctb += gfc->s3_ll[kk] * cb[k];
586 /* calculate the tonality of each threshold calculation partition
587 * calculate the SNR in each threshold calculation partition
588 * tonality = -0.299 - .43*log(ctb/ecb);
589 * tonality = 0: use NMT (lots of masking)
590 * tonality = 1: use TMN (little masking)
594 #define CONV1 (-.299)
601 if (tbb <= exp((1-CONV1)/CONV2))
602 { /* tonality near 1 */
603 tbb = exp(-LN_TO_LOG10 * TMN);
605 else if (tbb >= exp((0-CONV1)/CONV2))
606 {/* tonality near 0 */
607 tbb = exp(-LN_TO_LOG10 * NMT);
611 /* convert to tonality index */
612 /* tonality small: tbb=1 */
613 /* tonality large: tbb=-.299 */
614 tbb = CONV1 + CONV2*log(tbb);
615 tbb = exp(-LN_TO_LOG10 * ( TMN*tbb + (1-tbb)*NMT) );
619 /* at this point, tbb represents the amount the spreading function
620 * will be reduced. The smaller the value, the less masking.
621 * minval[] = 1 (0db) says just use tbb.
622 * minval[]= .01 (-20db) says reduce spreading function by
626 tbb = Min(gfc->minval[b], tbb);
628 /* stabilize tonality estimation */
629 if ( gfc->PSY->tonalityPatch ) {
632 FLOAT8 const x = 1.8699422;
633 FLOAT8 w = gfc->PSY->prvTonRed[b/2] * x;
637 gfc->PSY->prvTonRed[b] = tbb;
642 /* long block pre-echo control. */
643 /* rpelev=2.0, rpelev2=16.0 */
644 /* note: all surges in PE are because of this pre-echo formula
645 * for thr[b]. If it this is not used, PE is always around 600
647 /* dont use long block pre-echo control if previous granule was
648 * a short block. This is to avoid the situation:
649 * frame0: quiet (very low masking)
650 * frame1: surge (triggers short blocks)
651 * frame2: regular frame. looks like pre-echo when compared to
652 * frame0, but all pre-echo was in frame1.
654 /* chn=0,1 L and R channels
655 chn=2,3 S and M channels.
658 if (gfc->blocktype_old[chn>1 ? chn-2 : chn] == SHORT_TYPE )
659 thr[b] = ecb; /* Min(ecb, rpelev*gfc->nb_1[chn][b]); */
661 thr[b] = Min(ecb, Min(rpelev*gfc->nb_1[chn][b],rpelev2*gfc->nb_2[chn][b]) );
665 thrpe = Max(thr[b],gfc->ATH->cb[b]);
667 printf("%i thr=%e ATH=%e \n",b,thr[b],gfc->ATH->cb[b]);
670 gfc->pe[chn] -= gfc->numlines_l[b] * log(thrpe / eb[b]);
673 if ( gfp->VBR == vbr_rh || gfp->VBR == vbr_mtrh ) {
674 thr[b] = Min(ecb, rpelev*gfc->nb_1[chn][b]);
675 if (gfc->blocktype_old[chn>1 ? chn-2 : chn] != SHORT_TYPE )
676 thr[b] = Min(thr[b], rpelev2*gfc->nb_2[chn][b]);
677 thr[b] = Max( thr[b], 1e-37 );
680 gfc->nb_2[chn][b] = gfc->nb_1[chn][b];
681 gfc->nb_1[chn][b] = ecb;
685 /***************************************************************
686 * determine the block type (window type) based on L & R channels
688 ***************************************************************/
689 { /* compute PE for all 4 channels */
690 FLOAT mn,mx,ma=0,mb=0,mc=0;
692 for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++)
694 ma += gfc->energy_s[0][j];
695 mb += gfc->energy_s[1][j];
696 mc += gfc->energy_s[2][j];
703 /* bit allocation is based on pe. */
705 FLOAT8 tmp = 400*log(mx/(1e-12+mn));
706 if (tmp>gfc->pe[chn]) gfc->pe[chn]=tmp;
709 /* block type is based just on L or R channel */
711 uselongblock[chn] = 1;
713 /* tuned for t1.wav. doesnt effect most other samples */
714 if (gfc->pe[chn] > 3000)
718 {/* big surge of energy - always use short blocks */
719 uselongblock[chn] = 0;
721 else if ((mx > 10*mn) && (gfc->pe[chn] > 1000))
722 {/* medium surge, medium pe - use short blocks */
723 uselongblock[chn] = 0;
726 /* disable short blocks */
727 if (gfp->short_blocks == short_block_dispensed)
729 if (gfp->short_blocks == short_block_forced)
734 #if defined(HAVE_GTK)
736 FLOAT mn,mx,ma=0,mb=0,mc=0;
738 for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++)
740 ma += gfc->energy_s[0][j];
741 mb += gfc->energy_s[1][j];
742 mc += gfc->energy_s[2][j];
749 gfc->pinfo->ers[gr_out][chn]=gfc->ers_save[chn];
750 gfc->ers_save[chn]=(mx/(1e-12+mn));
751 gfc->pinfo->pe[gr_out][chn]=gfc->pe_save[chn];
752 gfc->pe_save[chn]=gfc->pe[chn];
757 /***************************************************************
758 * compute masking thresholds for both short and long blocks
759 ***************************************************************/
760 /* longblock threshold calculation (part 2) */
761 for ( sb = 0; sb < NBPSY_l; sb++ )
763 FLOAT8 enn = gfc->w1_l[sb] * eb[gfc->bu_l[sb]] + gfc->w2_l[sb] * eb[gfc->bo_l[sb]];
764 FLOAT8 thmm = gfc->w1_l[sb] *thr[gfc->bu_l[sb]] + gfc->w2_l[sb] * thr[gfc->bo_l[sb]];
766 for ( b = gfc->bu_l[sb]+1; b < gfc->bo_l[sb]; b++ )
772 gfc->en [chn].l[sb] = enn;
773 gfc->thm[chn].l[sb] = thmm;
777 /* threshold calculation for short blocks */
778 for ( sblock = 0; sblock < 3; sblock++ )
781 for ( b = 0; b < gfc->npart_s_orig; b++ )
783 FLOAT ecb = gfc->energy_s[sblock][j++];
784 for (i = 1 ; i<gfc->numlines_s[b]; ++i)
786 ecb += gfc->energy_s[sblock][j++];
792 for ( b = 0; b < gfc->npart_s; b++ )
795 for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ ) {
796 ecb += gfc->s3_ss[kk++] * eb[k];
798 ecb *= gfc->SNR_s[b];
801 if ( gfp->VBR == vbr_off || gfp->VBR == vbr_abr ) {
802 /* this looks like a BUG to me. robert */
803 thr[b] = Max (1e-6, ecb);
806 thr[b] = Min( ecb, rpelev_s * gfc->nb_s1[chn][b] );
807 if (gfc->blocktype_old[chn>1 ? chn-2 : chn] == SHORT_TYPE ) {
808 thr[b] = Min( thr[b], rpelev2_s * gfc->nb_s2[chn][b] );
810 thr[b] = Max( thr[b], 1e-37 );
811 gfc->nb_s2[chn][b] = gfc->nb_s1[chn][b];
812 gfc->nb_s1[chn][b] = ecb;
817 for ( sb = 0; sb < NBPSY_s; sb++ )
819 FLOAT8 enn = gfc->w1_s[sb] * eb[gfc->bu_s[sb]] + gfc->w2_s[sb] * eb[gfc->bo_s[sb]];
820 FLOAT8 thmm = gfc->w1_s[sb] *thr[gfc->bu_s[sb]] + gfc->w2_s[sb] * thr[gfc->bo_s[sb]];
822 for ( b = gfc->bu_s[sb]+1; b < gfc->bo_s[sb]; b++ ) {
827 gfc->en [chn].s[sb][sblock] = enn;
828 gfc->thm[chn].s[sb][sblock] = thmm;
832 } /* end loop over chn */
836 /***************************************************************
837 * compute M/S thresholds
838 ***************************************************************/
839 /* compute M/S thresholds from Johnston & Ferreira 1992 ICASSP paper */
840 if (gfp->mode == JOINT_STEREO) {
841 FLOAT8 rside,rmid,mld;
842 int chmid=2,chside=3;
844 for ( sb = 0; sb < NBPSY_l; sb++ ) {
845 /* use this fix if L & R masking differs by 2db or less */
846 /* if db = 10*log10(x2/x1) < 2 */
847 /* if (x2 < 1.58*x1) { */
848 if (gfc->thm[0].l[sb] <= 1.58*gfc->thm[1].l[sb]
849 && gfc->thm[1].l[sb] <= 1.58*gfc->thm[0].l[sb]) {
851 mld = gfc->mld_l[sb]*gfc->en[chside].l[sb];
852 rmid = Max(gfc->thm[chmid].l[sb], Min(gfc->thm[chside].l[sb],mld));
854 mld = gfc->mld_l[sb]*gfc->en[chmid].l[sb];
855 rside = Max(gfc->thm[chside].l[sb],Min(gfc->thm[chmid].l[sb],mld));
857 gfc->thm[chmid].l[sb]=rmid;
858 gfc->thm[chside].l[sb]=rside;
861 for ( sb = 0; sb < NBPSY_s; sb++ ) {
862 for ( sblock = 0; sblock < 3; sblock++ ) {
863 if (gfc->thm[0].s[sb][sblock] <= 1.58*gfc->thm[1].s[sb][sblock]
864 && gfc->thm[1].s[sb][sblock] <= 1.58*gfc->thm[0].s[sb][sblock]) {
866 mld = gfc->mld_s[sb]*gfc->en[chside].s[sb][sblock];
867 rmid = Max(gfc->thm[chmid].s[sb][sblock],Min(gfc->thm[chside].s[sb][sblock],mld));
869 mld = gfc->mld_s[sb]*gfc->en[chmid].s[sb][sblock];
870 rside = Max(gfc->thm[chside].s[sb][sblock],Min(gfc->thm[chmid].s[sb][sblock],mld));
872 gfc->thm[chmid].s[sb][sblock]=rmid;
873 gfc->thm[chside].s[sb][sblock]=rside;
881 /***************************************************************
882 * Adjust M/S maskings if user set "msfix"
883 ***************************************************************/
884 /* Naoki Shibata 2000 */
885 if (numchn == 4 && gfp->msfix!=0) {
886 FLOAT msfix = gfp->msfix;
888 for ( sb = 0; sb < NBPSY_l; sb++ )
890 FLOAT8 thmL,thmR,thmM,thmS,ath;
891 ath = (gfc->ATH->cb[(gfc->bu_l[sb] + gfc->bo_l[sb])/2])*pow(10,-gfp->ATHlower/10.0);
892 thmL = Max(gfc->thm[0].l[sb],ath);
893 thmR = Max(gfc->thm[1].l[sb],ath);
894 thmM = Max(gfc->thm[2].l[sb],ath);
895 thmS = Max(gfc->thm[3].l[sb],ath);
897 if (thmL*msfix < (thmM+thmS)/2) {
898 FLOAT8 f = thmL*msfix / ((thmM+thmS)/2);
902 if (thmR*msfix < (thmM+thmS)/2) {
903 FLOAT8 f = thmR*msfix / ((thmM+thmS)/2);
908 gfc->thm[2].l[sb] = Min(thmM,gfc->thm[2].l[sb]);
909 gfc->thm[3].l[sb] = Min(thmS,gfc->thm[3].l[sb]);
912 for ( sb = 0; sb < NBPSY_s; sb++ ) {
913 for ( sblock = 0; sblock < 3; sblock++ ) {
914 FLOAT8 thmL,thmR,thmM,thmS,ath;
915 ath = (gfc->ATH->cb[(gfc->bu_s[sb] + gfc->bo_s[sb])/2])*pow(10,-gfp->ATHlower/10.0);
916 thmL = Max(gfc->thm[0].s[sb][sblock],ath);
917 thmR = Max(gfc->thm[1].s[sb][sblock],ath);
918 thmM = Max(gfc->thm[2].s[sb][sblock],ath);
919 thmS = Max(gfc->thm[3].s[sb][sblock],ath);
921 if (thmL*msfix < (thmM+thmS)/2) {
922 FLOAT8 f = thmL*msfix / ((thmM+thmS)/2);
926 if (thmR*msfix < (thmM+thmS)/2) {
927 FLOAT8 f = thmR*msfix / ((thmM+thmS)/2);
932 gfc->thm[2].s[sb][sblock] = Min(gfc->thm[2].s[sb][sblock],thmM);
933 gfc->thm[3].s[sb][sblock] = Min(gfc->thm[3].s[sb][sblock],thmS);
948 if (gfp->mode == JOINT_STEREO) {
949 /* determin ms_ratio from masking thresholds*/
950 /* use ms_stereo (ms_ratio < .35) if average thresh. diff < 5 db */
951 FLOAT8 db,x1,x2,sidetot=0,tot=0;
952 for (sb= NBPSY_l/4 ; sb< NBPSY_l; sb ++ ) {
953 x1 = Min(gfc->thm[0].l[sb],gfc->thm[1].l[sb]);
954 x2 = Max(gfc->thm[0].l[sb],gfc->thm[1].l[sb]);
955 /* thresholds difference in db */
956 if (x2 >= 1000*x1) db=3;
957 else db = log10(x2/x1);
958 /* DEBUGF(gfc,"db = %f %e %e \n",db,gfc->thm[0].l[sb],gfc->thm[1].l[sb]);*/
962 ms_ratio_l= (sidetot/tot)*0.7; /* was .35*(sidetot/tot)/5.0*10 */
963 ms_ratio_l = Min(ms_ratio_l,0.5);
966 for ( sblock = 0; sblock < 3; sblock++ )
967 for ( sb = NBPSY_s/4; sb < NBPSY_s; sb++ ) {
968 x1 = Min(gfc->thm[0].s[sb][sblock],gfc->thm[1].s[sb][sblock]);
969 x2 = Max(gfc->thm[0].s[sb][sblock],gfc->thm[1].s[sb][sblock]);
970 /* thresholds difference in db */
971 if (x2 >= 1000*x1) db=3;
972 else db = log10(x2/x1);
976 ms_ratio_s = (sidetot/tot)*0.7; /* was .35*(sidetot/tot)/5.0*10 */
977 ms_ratio_s = Min(ms_ratio_s,.5);
980 /***************************************************************
981 * determine final block type
982 ***************************************************************/
984 for (chn=0; chn<gfc->channels_out; chn++) {
985 blocktype[chn] = NORM_TYPE;
989 if (gfp->short_blocks == short_block_coupled) {
990 /* force both channels to use the same block type */
991 /* this is necessary if the frame is to be encoded in ms_stereo. */
992 /* But even without ms_stereo, FhG does this */
993 int bothlong= (uselongblock[0] && uselongblock[1]);
1002 /* update the blocktype of the previous granule, since it depends on what
1003 * happend in this granule */
1004 for (chn=0; chn<gfc->channels_out; chn++) {
1005 if ( uselongblock[chn])
1006 { /* no attack : use long blocks */
1007 assert( gfc->blocktype_old[chn] != START_TYPE );
1008 switch( gfc->blocktype_old[chn] )
1012 blocktype[chn] = NORM_TYPE;
1015 blocktype[chn] = STOP_TYPE;
1019 /* attack : use short blocks */
1020 blocktype[chn] = SHORT_TYPE;
1021 if ( gfc->blocktype_old[chn] == NORM_TYPE ) {
1022 gfc->blocktype_old[chn] = START_TYPE;
1024 if ( gfc->blocktype_old[chn] == STOP_TYPE ) {
1025 gfc->blocktype_old[chn] = SHORT_TYPE ;
1029 blocktype_d[chn] = gfc->blocktype_old[chn]; /* value returned to calling program */
1030 gfc->blocktype_old[chn] = blocktype[chn]; /* save for next call to l3psy_anal */
1033 if (blocktype_d[0]==2 && blocktype_d[1]==2)
1034 *ms_ratio = gfc->ms_ratio_s_old;
1036 *ms_ratio = gfc->ms_ratio_l_old;
1038 gfc->ms_ratio_s_old = ms_ratio_s;
1039 gfc->ms_ratio_l_old = ms_ratio_l;
1041 /* we dont know the block type of this frame yet - assume long */
1042 *ms_ratio_next = ms_ratio_l;
1047 /* addition of simultaneous masking Naoki Shibata 2000/7 */
1048 inline static FLOAT8 mask_add(FLOAT8 m1,FLOAT8 m2,int k,int b, lame_internal_flags * const gfc)
1050 static const FLOAT8 table1[] = {
1051 3.3246 *3.3246 ,3.23837*3.23837,3.15437*3.15437,3.00412*3.00412,2.86103*2.86103,2.65407*2.65407,2.46209*2.46209,2.284 *2.284 ,
1052 2.11879*2.11879,1.96552*1.96552,1.82335*1.82335,1.69146*1.69146,1.56911*1.56911,1.46658*1.46658,1.37074*1.37074,1.31036*1.31036,
1053 1.25264*1.25264,1.20648*1.20648,1.16203*1.16203,1.12765*1.12765,1.09428*1.09428,1.0659 *1.0659 ,1.03826*1.03826,1.01895*1.01895,
1057 static const FLOAT8 table2[] = {
1058 1.33352*1.33352,1.35879*1.35879,1.38454*1.38454,1.39497*1.39497,1.40548*1.40548,1.3537 *1.3537 ,1.30382*1.30382,1.22321*1.22321,
1062 static const FLOAT8 table3[] = {
1063 2.35364*2.35364,2.29259*2.29259,2.23313*2.23313,2.12675*2.12675,2.02545*2.02545,1.87894*1.87894,1.74303*1.74303,1.61695*1.61695,
1064 1.49999*1.49999,1.39148*1.39148,1.29083*1.29083,1.19746*1.19746,1.11084*1.11084,1.03826*1.03826
1071 if (m1 == 0) return m2;
1075 i = 10*log10(m2 / m1)/10*16;
1076 m = 10*log10((m1+m2)/gfc->ATH->cb[k]);
1080 if (b <= 3) { /* approximately, 1 bark = 3 partitions */
1081 if (i > 8) return m1+m2;
1082 return (m1+m2)*table2[i];
1088 if (i > 24) return m1+m2;
1089 if (i > 13) f = 1; else f = table3[i];
1091 return (m1+m2)*(table1[i]*r+f*(1-r));
1093 if (i > 13) return m1+m2;
1094 return (m1+m2)*table3[i];
1097 if (i > 24) return m1+m2;
1098 return (m1+m2)*table1[i];
1101 int L3psycho_anal_ns( lame_global_flags * gfp,
1102 const sample_t *buffer[2], int gr_out,
1104 FLOAT8 *ms_ratio_next,
1105 III_psy_ratio masking_ratio[2][2],
1106 III_psy_ratio masking_MS_ratio[2][2],
1107 FLOAT8 percep_entropy[2],FLOAT8 percep_MS_entropy[2],
1111 /* to get a good cache performance, one has to think about
1112 * the sequence, in which the variables are used.
1113 * (Note: these static variables have been moved to the gfc-> struct,
1114 * and their order in memory is layed out in util.h)
1116 lame_internal_flags *gfc=gfp->internal_flags;
1118 /* fft and energy calculation */
1119 FLOAT (*wsamp_l)[BLKSIZE];
1120 FLOAT (*wsamp_s)[3][BLKSIZE_s];
1123 FLOAT8 eb[CBANDS],eb2[CBANDS];
1126 FLOAT8 max[CBANDS],avg[CBANDS],tonality2[CBANDS];
1129 FLOAT8 ms_ratio_l=0,ms_ratio_s=0;
1132 int blocktype[2],uselongblock[2];
1134 /* usual variables like loop indices, etc.. */
1139 /* variables used for --nspsytune */
1141 FLOAT ns_hpfsmpl[4][576+576/3+NSFIRLEN];
1142 FLOAT pe_l[4],pe_s[4];
1146 if(gfc->psymodel_init==0) {
1149 gfc->psymodel_init=1;
1153 numchn = gfc->channels_out;
1154 /* chn=2 and 3 = Mid and Side channels */
1155 if (gfp->mode == JOINT_STEREO) numchn=4;
1157 if (gfp->VBR==vbr_off) pcfact = gfc->ResvMax == 0 ? 0 : ((FLOAT)gfc->ResvSize)/gfc->ResvMax*0.5;
1158 else if (gfp->VBR == vbr_rh || gfp->VBR == vbr_mtrh || gfp->VBR == vbr_mt) {
1159 static const FLOAT8 pcQns[10]={1.0,1.0,1.0,0.8,0.6,0.5,0.4,0.3,0.2,0.1};
1160 pcfact = pcQns[gfp->VBR_q];
1163 /**********************************************************************
1164 * Apply HPF of fs/4 to the input signal.
1165 * This is used for attack detection / handling.
1166 **********************************************************************/
1169 static const FLOAT fircoef[] = {
1170 -8.65163e-18,-0.00851586,-6.74764e-18, 0.0209036,
1171 -3.36639e-17,-0.0438162 ,-1.54175e-17, 0.0931738,
1172 -5.52212e-17,-0.313819 , 0.5 ,-0.313819,
1173 -5.52212e-17, 0.0931738 ,-1.54175e-17,-0.0438162,
1174 -3.36639e-17, 0.0209036 ,-6.74764e-18,-0.00851586,
1178 for(chn=0;chn<gfc->channels_out;chn++)
1180 FLOAT firbuf[576+576/3+NSFIRLEN];
1182 /* apply high pass filter of fs/4 */
1184 for(i=-NSFIRLEN;i<576+576/3;i++)
1185 firbuf[i+NSFIRLEN] = buffer[chn][576-350+(i)];
1187 for(i=0;i<576+576/3-NSFIRLEN;i++)
1190 for(j=0;j<NSFIRLEN;j++)
1191 sum += fircoef[j] * firbuf[i+j];
1192 ns_hpfsmpl[chn][i] = sum;
1194 for(;i<576+576/3;i++)
1195 ns_hpfsmpl[chn][i] = 0;
1198 if (gfp->mode == JOINT_STEREO) {
1199 for(i=0;i<576+576/3;i++)
1201 ns_hpfsmpl[2][i] = ns_hpfsmpl[0][i]+ns_hpfsmpl[1][i];
1202 ns_hpfsmpl[3][i] = ns_hpfsmpl[0][i]-ns_hpfsmpl[1][i];
1209 /* there is a one granule delay. Copy maskings computed last call
1210 * into masking_ratio to return to calling program.
1212 for (chn=0; chn<numchn; chn++) {
1213 for (i=0; i<numchn; ++i) {
1214 energy[chn]=gfc->tot_ener[chn];
1219 for (chn=0; chn<numchn; chn++) {
1220 /* there is a one granule delay. Copy maskings computed last call
1221 * into masking_ratio to return to calling program.
1223 pe_l[chn] = gfc->nsPsy.pe_l[chn];
1224 pe_s[chn] = gfc->nsPsy.pe_s[chn];
1228 //percep_entropy [chn] = gfc -> pe [chn];
1229 masking_ratio [gr_out] [chn] .en = gfc -> en [chn];
1230 masking_ratio [gr_out] [chn] .thm = gfc -> thm [chn];
1233 //percep_MS_entropy [chn-2] = gfc -> pe [chn];
1234 masking_MS_ratio [gr_out] [chn-2].en = gfc -> en [chn];
1235 masking_MS_ratio [gr_out] [chn-2].thm = gfc -> thm [chn];
1239 for (chn=0; chn<numchn; chn++) {
1241 /**********************************************************************
1243 **********************************************************************/
1245 wsamp_s = gfc->wsamp_S+(chn & 1);
1246 wsamp_l = gfc->wsamp_L+(chn & 1);
1249 fft_long ( gfc, *wsamp_l, chn, buffer);
1250 fft_short( gfc, *wsamp_s, chn, buffer);
1253 /* FFT data for mid and side channel is derived from L & R */
1257 for (j = BLKSIZE-1; j >=0 ; --j)
1259 FLOAT l = gfc->wsamp_L[0][j];
1260 FLOAT r = gfc->wsamp_L[1][j];
1261 gfc->wsamp_L[0][j] = (l+r)*(FLOAT)(SQRT2*0.5);
1262 gfc->wsamp_L[1][j] = (l-r)*(FLOAT)(SQRT2*0.5);
1264 for (b = 2; b >= 0; --b)
1266 for (j = BLKSIZE_s-1; j >= 0 ; --j)
1268 FLOAT l = gfc->wsamp_S[0][b][j];
1269 FLOAT r = gfc->wsamp_S[1][b][j];
1270 gfc->wsamp_S[0][b][j] = (l+r)*(FLOAT)(SQRT2*0.5);
1271 gfc->wsamp_S[1][b][j] = (l-r)*(FLOAT)(SQRT2*0.5);
1277 /**********************************************************************
1278 * compute energies for each spectral line
1279 **********************************************************************/
1283 gfc->energy[0] = (*wsamp_l)[0];
1284 gfc->energy[0] *= gfc->energy[0];
1286 gfc->tot_ener[chn] = gfc->energy[0]; /* sum total energy at nearly no extra cost */
1288 for (j=BLKSIZE/2-1; j >= 0; --j)
1290 FLOAT re = (*wsamp_l)[BLKSIZE/2-j];
1291 FLOAT im = (*wsamp_l)[BLKSIZE/2+j];
1292 gfc->energy[BLKSIZE/2-j] = (re * re + im * im) * (FLOAT)0.5;
1294 if (BLKSIZE/2-j > 10)
1295 gfc->tot_ener[chn] += gfc->energy[BLKSIZE/2-j];
1301 for (b = 2; b >= 0; --b)
1303 gfc->energy_s[b][0] = (*wsamp_s)[b][0];
1304 gfc->energy_s[b][0] *= gfc->energy_s [b][0];
1305 for (j=BLKSIZE_s/2-1; j >= 0; --j)
1307 FLOAT re = (*wsamp_s)[b][BLKSIZE_s/2-j];
1308 FLOAT im = (*wsamp_s)[b][BLKSIZE_s/2+j];
1309 gfc->energy_s[b][BLKSIZE_s/2-j] = (re * re + im * im) * (FLOAT)0.5;
1314 /* output data for analysis */
1316 #if defined(HAVE_GTK)
1317 if (gfp->analysis) {
1318 for (j=0; j<HBLKSIZE ; j++) {
1319 gfc->pinfo->energy[gr_out][chn][j]=gfc->energy_save[chn][j];
1320 gfc->energy_save[chn][j]=gfc->energy[j];
1326 /**********************************************************************
1327 * compute loudness approximation (used for ATH auto-level adjustment)
1328 **********************************************************************/
1329 if( gfp->athaa_loudapprox == 2 ) {
1330 if( chn < 2 ) { /* no loudness for mid and side channels */
1331 gfc->loudness_sq[gr_out][chn] = gfc->loudness_sq_save[chn];
1332 gfc->loudness_sq_save[chn]
1333 = psycho_loudness_approx( gfc->energy, gfp);
1339 /**********************************************************************
1340 * Calculate the energy and the tonality of each partition.
1341 **********************************************************************/
1343 for (b=0, j=0; b<gfc->npart_l_orig; b++)
1347 ebb = gfc->energy[j];
1348 m = a = gfc->energy[j];
1351 for (i = gfc->numlines_l[b] - 1; i > 0; i--)
1353 FLOAT8 el = gfc->energy[j];
1354 ebb += gfc->energy[j];
1356 m = m < el ? el : m;
1361 avg[b] = a / gfc->numlines_l[b];
1365 for (b=0; b < gfc->npart_l_orig; b++ )
1372 for(k=b-1;k<=b+1;k++)
1374 if (k >= 0 && k < gfc->npart_l_orig) {
1376 c2 += gfc->numlines_l[k];
1378 m = m < max[k] ? max[k] : m;
1383 tonality2[b] = a == 0 ? 0 : (m / a - 1)/(c2-1);
1386 for (b=0; b < gfc->npart_l_orig; b++ )
1389 static FLOAT8 tab[20] =
1390 { 0, 1, 2, 2, 2, 2, 2, 6,9.3,9.3,9.3,9.3,9.3,9.3,9.3,9.3,9.3,9.3,9.3,9.3};
1392 static int init = 1;
1396 tab[j] = pow(10.0,-tab[j]/10.0);
1401 static FLOAT8 tab[20] = {
1402 1,0.79433,0.63096,0.63096,0.63096,0.63096,0.63096,0.25119,0.11749,0.11749,
1403 0.11749,0.11749,0.11749,0.11749,0.11749,0.11749,0.11749,0.11749,0.11749,0.11749
1407 int t = 20*tonality2[b];
1409 eb2[b] = eb[b] * tab[t];
1413 /**********************************************************************
1414 * convolve the partitioned energy and unpredictability
1415 * with the spreading function, s3_l[b][k]
1416 ******************************************************************** */
1419 for ( b = 0;b < gfc->npart_l; b++ )
1423 /**** convolve the partitioned energy with the spreading function ****/
1428 for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ )
1430 ecb = mask_add(ecb,gfc->s3_ll[kk++] * eb2[k],k,k-b,gfc);
1433 ecb *= 0.158489319246111; // pow(10,-0.8)
1437 for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ )
1439 ecb += gfc->s3_ll[kk++] * eb2[k];
1442 ecb *= 0.223872113856834; // pow(10,-0.65);
1445 /**** long block pre-echo control ****/
1447 /* dont use long block pre-echo control if previous granule was
1448 * a short block. This is to avoid the situation:
1449 * frame0: quiet (very low masking)
1450 * frame1: surge (triggers short blocks)
1451 * frame2: regular frame. looks like pre-echo when compared to
1452 * frame0, but all pre-echo was in frame1.
1455 /* chn=0,1 L and R channels
1456 chn=2,3 S and M channels.
1459 #define NS_INTERP(x,y,r) (pow((x),(r))*pow((y),1-(r)))
1461 if (gfc->blocktype_old[chn>1 ? chn-2 : chn] == SHORT_TYPE )
1462 thr[b] = ecb; /* Min(ecb, rpelev*gfc->nb_1[chn][b]); */
1464 thr[b] = NS_INTERP(Min(ecb, Min(rpelev*gfc->nb_1[chn][b],rpelev2*gfc->nb_2[chn][b])),ecb,pcfact);
1466 gfc->nb_2[chn][b] = gfc->nb_1[chn][b];
1467 gfc->nb_1[chn][b] = ecb;
1471 /***************************************************************
1472 * determine the block type (window type)
1473 ***************************************************************/
1477 FLOAT en_subshort[12];
1478 FLOAT attack_intensity[12];
1479 int ns_uselongblock = 1;
1481 /* calculate energies of each sub-shortblocks */
1487 for(j=0;j<576/9;j++)
1489 en_subshort[i] += ns_hpfsmpl[chn][k] * ns_hpfsmpl[chn][k];
1493 if (en_subshort[i] < 100) en_subshort[i] = 100;
1496 /* compare energies between sub-shortblocks */
1498 #define NSATTACKTHRE 150
1499 #define NSATTACKTHRE_S 300
1502 attack_intensity[i] = en_subshort[i] / gfc->nsPsy.last_en_subshort[chn][7+i];
1506 attack_intensity[i] = en_subshort[i] / en_subshort[i-2];
1509 ns_attacks[0] = ns_attacks[1] = ns_attacks[2] = ns_attacks[3] = 0;
1513 if (!ns_attacks[i/3] && attack_intensity[i] > (chn == 3 ? (gfc->presetTune.use ? gfc->presetTune.attackthre_s : NSATTACKTHRE_S)
1514 : (gfc->presetTune.use ? gfc->presetTune.attackthre : NSATTACKTHRE))) ns_attacks[i/3] = (i % 3)+1;
1517 if (ns_attacks[0] && gfc->nsPsy.last_attacks[chn][2]) ns_attacks[0] = 0;
1518 if (ns_attacks[1] && ns_attacks[0]) ns_attacks[1] = 0;
1519 if (ns_attacks[2] && ns_attacks[1]) ns_attacks[2] = 0;
1520 if (ns_attacks[3] && ns_attacks[2]) ns_attacks[3] = 0;
1522 if (gfc->nsPsy.last_attacks[chn][2] == 3 ||
1523 ns_attacks[0] || ns_attacks[1] || ns_attacks[2] || ns_attacks[3]) ns_uselongblock = 0;
1525 if (chn < 4) count++;
1529 gfc->nsPsy.last_en_subshort[chn][i] = en_subshort[i];
1530 gfc->nsPsy.last_attack_intensity[chn][i] = attack_intensity[i];
1533 if (gfp->short_blocks == short_block_dispensed) {
1534 uselongblock[chn] = 1;
1536 else if (gfp->short_blocks == short_block_forced) {
1537 uselongblock[chn] = 0;
1541 uselongblock[chn] = ns_uselongblock;
1543 if (ns_uselongblock == 0) uselongblock[0] = uselongblock[1] = 0;
1548 #if defined(HAVE_GTK)
1549 if (gfp->analysis) {
1550 FLOAT mn,mx,ma=0,mb=0,mc=0;
1552 for ( j = HBLKSIZE_s/2; j < HBLKSIZE_s; j ++)
1554 ma += gfc->energy_s[0][j];
1555 mb += gfc->energy_s[1][j];
1556 mc += gfc->energy_s[2][j];
1563 gfc->pinfo->ers[gr_out][chn]=gfc->ers_save[chn];
1564 gfc->ers_save[chn]=(mx/(1e-12+mn));
1569 /***************************************************************
1570 * compute masking thresholds for long blocks
1571 ***************************************************************/
1573 for ( sb = 0; sb < NBPSY_l; sb++ )
1575 FLOAT8 enn = gfc->w1_l[sb] * eb[gfc->bu_l[sb]] + gfc->w2_l[sb] * eb[gfc->bo_l[sb]];
1576 FLOAT8 thmm = gfc->w1_l[sb] *thr[gfc->bu_l[sb]] + gfc->w2_l[sb] * thr[gfc->bo_l[sb]];
1578 for ( b = gfc->bu_l[sb]+1; b < gfc->bo_l[sb]; b++ )
1584 gfc->en [chn].l[sb] = enn;
1585 gfc->thm[chn].l[sb] = thmm;
1589 /***************************************************************
1590 * compute masking thresholds for short blocks
1591 ***************************************************************/
1593 for ( sblock = 0; sblock < 3; sblock++ )
1596 for ( b = 0; b < gfc->npart_s_orig; b++ )
1598 FLOAT ecb = gfc->energy_s[sblock][j++];
1599 for (i = 1 ; i<gfc->numlines_s[b]; ++i)
1601 ecb += gfc->energy_s[sblock][j++];
1608 for ( b = 0; b < gfc->npart_s; b++ )
1611 for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ )
1613 ecb += gfc->s3_ss[kk++] * eb[k];
1617 /* this looks like a BUG */
1618 thr[b] = Max (1e-6, ecb);
1620 if (gfp->VBR == vbr_mtrh) {
1621 thr[b] = Min( ecb, rpelev_s * gfc->nb_s1[chn][b] );
1622 if (gfc->blocktype_old[chn>1 ? chn-2 : chn] == SHORT_TYPE ) {
1623 thr[b] = Min( thr[b], rpelev2_s * gfc->nb_s2[chn][b] );
1625 thr[b] = Max( thr[b], 1e-37 );
1626 gfc->nb_s2[chn][b] = gfc->nb_s1[chn][b];
1627 gfc->nb_s1[chn][b] = ecb;
1631 for ( sb = 0; sb < NBPSY_s; sb++ )
1633 FLOAT8 enn = gfc->w1_s[sb] * eb[gfc->bu_s[sb]] + gfc->w2_s[sb] * eb[gfc->bo_s[sb]];
1634 FLOAT8 thmm = gfc->w1_s[sb] *thr[gfc->bu_s[sb]] + gfc->w2_s[sb] * thr[gfc->bo_s[sb]];
1636 for ( b = gfc->bu_s[sb]+1; b < gfc->bo_s[sb]; b++ )
1642 /**** short block pre-echo control ****/
1644 #define NS_PREECHO_ATT0 0.8
1645 #define NS_PREECHO_ATT1 0.6
1646 #define NS_PREECHO_ATT2 0.3
1648 thmm *= NS_PREECHO_ATT0;
1650 if (ns_attacks[sblock] >= 2) {
1652 double p = NS_INTERP(gfc->thm[chn].s[sb][sblock-1],thmm,NS_PREECHO_ATT1*pcfact);
1655 double p = NS_INTERP(gfc->nsPsy.last_thm[chn][sb][2],thmm,NS_PREECHO_ATT1*pcfact);
1658 } else if (ns_attacks[sblock+1] == 1) {
1660 double p = NS_INTERP(gfc->thm[chn].s[sb][sblock-1],thmm,NS_PREECHO_ATT1*pcfact);
1663 double p = NS_INTERP(gfc->nsPsy.last_thm[chn][sb][2],thmm,NS_PREECHO_ATT1*pcfact);
1668 if (ns_attacks[sblock] == 1) {
1669 double p = sblock == 0 ? gfc->nsPsy.last_thm[chn][sb][2] : gfc->thm[chn].s[sb][sblock-1];
1670 p = NS_INTERP(p,thmm,NS_PREECHO_ATT2*pcfact);
1672 } else if ((sblock != 0 && ns_attacks[sblock-1] == 3) ||
1673 (sblock == 0 && gfc->nsPsy.last_attacks[chn][2] == 3)) {
1674 double p = sblock <= 1 ? gfc->nsPsy.last_thm[chn][sb][sblock+1] : gfc->thm[chn].s[sb][0];
1675 p = NS_INTERP(p,thmm,NS_PREECHO_ATT2*pcfact);
1679 gfc->en [chn].s[sb][sblock] = enn;
1680 gfc->thm[chn].s[sb][sblock] = thmm;
1685 /***************************************************************
1686 * save some values for analysis of the next granule
1687 ***************************************************************/
1689 for ( sblock = 0; sblock < 3; sblock++ )
1691 for ( sb = 0; sb < NBPSY_s; sb++ )
1693 gfc->nsPsy.last_thm[chn][sb][sblock] = gfc->thm[chn].s[sb][sblock];
1698 gfc->nsPsy.last_attacks[chn][i] = ns_attacks[i];
1700 } /* end loop over chn */
1704 /***************************************************************
1705 * compute M/S thresholds
1706 ***************************************************************/
1708 /* from Johnston & Ferreira 1992 ICASSP paper */
1710 if ( numchn==4 /* mid/side and r/l */) {
1711 FLOAT8 rside,rmid,mld;
1712 int chmid=2,chside=3;
1714 for ( sb = 0; sb < NBPSY_l; sb++ ) {
1715 /* use this fix if L & R masking differs by 2db or less */
1716 /* if db = 10*log10(x2/x1) < 2 */
1717 /* if (x2 < 1.58*x1) { */
1718 if (gfc->thm[0].l[sb] <= 1.58*gfc->thm[1].l[sb]
1719 && gfc->thm[1].l[sb] <= 1.58*gfc->thm[0].l[sb]) {
1721 mld = gfc->mld_l[sb]*gfc->en[chside].l[sb];
1722 rmid = Max(gfc->thm[chmid].l[sb], Min(gfc->thm[chside].l[sb],mld));
1724 mld = gfc->mld_l[sb]*gfc->en[chmid].l[sb];
1725 rside = Max(gfc->thm[chside].l[sb],Min(gfc->thm[chmid].l[sb],mld));
1727 gfc->thm[chmid].l[sb]=rmid;
1728 gfc->thm[chside].l[sb]=rside;
1731 for ( sb = 0; sb < NBPSY_s; sb++ ) {
1732 for ( sblock = 0; sblock < 3; sblock++ ) {
1733 if (gfc->thm[0].s[sb][sblock] <= 1.58*gfc->thm[1].s[sb][sblock]
1734 && gfc->thm[1].s[sb][sblock] <= 1.58*gfc->thm[0].s[sb][sblock]) {
1736 mld = gfc->mld_s[sb]*gfc->en[chside].s[sb][sblock];
1737 rmid = Max(gfc->thm[chmid].s[sb][sblock],Min(gfc->thm[chside].s[sb][sblock],mld));
1739 mld = gfc->mld_s[sb]*gfc->en[chmid].s[sb][sblock];
1740 rside = Max(gfc->thm[chside].s[sb][sblock],Min(gfc->thm[chmid].s[sb][sblock],mld));
1742 gfc->thm[chmid].s[sb][sblock]=rmid;
1743 gfc->thm[chside].s[sb][sblock]=rside;
1750 /* Naoki Shibata 2000 */
1752 #define NS_MSFIX 3.5
1755 FLOAT msfix = NS_MSFIX;
1756 if (gfc->nsPsy.safejoint) msfix = 1;
1757 if (gfp->msfix) msfix = gfp->msfix;
1759 if (gfc->presetTune.use && gfc->ATH->adjust >=
1760 gfc->presetTune.athadjust_switch_level &&
1761 gfc->presetTune.athadjust_msfix > 0)
1762 msfix = gfc->presetTune.athadjust_msfix;
1764 for ( sb = 0; sb < NBPSY_l; sb++ )
1766 FLOAT8 thmL,thmR,thmM,thmS,ath;
1767 ath = (gfc->ATH->cb[(gfc->bu_l[sb] + gfc->bo_l[sb])/2])*pow(10,-gfp->ATHlower/10.0);
1768 thmL = Max(gfc->thm[0].l[sb],ath);
1769 thmR = Max(gfc->thm[1].l[sb],ath);
1770 thmM = Max(gfc->thm[2].l[sb],ath);
1771 thmS = Max(gfc->thm[3].l[sb],ath);
1773 if (thmL*msfix < (thmM+thmS)/2) {
1774 FLOAT8 f = thmL * (gfc->presetTune.use ? gfc->presetTune.ms_maskadjust : msfix) / ((thmM+thmS)/2);
1778 if (thmR*msfix < (thmM+thmS)/2) {
1779 FLOAT8 f = thmR * (gfc->presetTune.use ? gfc->presetTune.ms_maskadjust : msfix) / ((thmM+thmS)/2);
1784 gfc->thm[2].l[sb] = Min(thmM,gfc->thm[2].l[sb]);
1785 gfc->thm[3].l[sb] = Min(thmS,gfc->thm[3].l[sb]);
1788 for ( sb = 0; sb < NBPSY_s; sb++ ) {
1789 for ( sblock = 0; sblock < 3; sblock++ ) {
1790 FLOAT8 thmL,thmR,thmM,thmS,ath;
1791 ath = (gfc->ATH->cb[(gfc->bu_s[sb] + gfc->bo_s[sb])/2])*pow(10,-gfp->ATHlower/10.0);
1792 thmL = Max(gfc->thm[0].s[sb][sblock],ath);
1793 thmR = Max(gfc->thm[1].s[sb][sblock],ath);
1794 thmM = Max(gfc->thm[2].s[sb][sblock],ath);
1795 thmS = Max(gfc->thm[3].s[sb][sblock],ath);
1797 if (thmL*msfix < (thmM+thmS)/2) {
1798 FLOAT8 f = thmL*msfix / ((thmM+thmS)/2);
1802 if (thmR*msfix < (thmM+thmS)/2) {
1803 FLOAT8 f = thmR*msfix / ((thmM+thmS)/2);
1808 gfc->thm[2].s[sb][sblock] = Min(gfc->thm[2].s[sb][sblock],thmM);
1809 gfc->thm[3].s[sb][sblock] = Min(gfc->thm[3].s[sb][sblock],thmS);
1818 /***************************************************************
1819 * compute estimation of the amount of bit used in the granule
1820 ***************************************************************/
1822 for(chn=0;chn<numchn;chn++)
1825 static FLOAT8 regcoef[] = {
1826 1124.23,10.0583,10.7484,7.29006,16.2714,6.2345,4.09743,3.05468,3.33196,2.54688,
1827 3.68168,5.83109,2.93817,-8.03277,-10.8458,8.48777,9.13182,2.05212,8.6674,50.3937,73.267,97.5664,0
1830 FLOAT8 msum = regcoef[0]/4;
1833 for ( sb = 0; sb < NBPSY_l; sb++ )
1837 if (gfc->thm[chn].l[sb]*gfc->masking_lower != 0 &&
1838 gfc->en[chn].l[sb]/(gfc->thm[chn].l[sb]*gfc->masking_lower) > 1)
1839 t = log(gfc->en[chn].l[sb]/(gfc->thm[chn].l[sb]*gfc->masking_lower));
1842 msum += regcoef[sb+1] * t;
1845 gfc->nsPsy.pe_l[chn] = msum;
1849 static FLOAT8 regcoef[] = {
1850 1236.28,0,0,0,0.434542,25.0738,0,0,0,19.5442,19.7486,60,100,0
1853 FLOAT8 msum = regcoef[0]/4;
1856 for(sblock=0;sblock<3;sblock++)
1858 for ( sb = 0; sb < NBPSY_s; sb++ )
1862 if (gfc->thm[chn].s[sb][sblock] * gfc->masking_lower != 0 &&
1863 gfc->en[chn].s[sb][sblock] / (gfc->thm[chn].s[sb][sblock] * gfc->masking_lower) > 1)
1864 t = log(gfc->en[chn].s[sb][sblock] / (gfc->thm[chn].s[sb][sblock] * gfc->masking_lower));
1867 msum += regcoef[sb+1] * t;
1871 gfc->nsPsy.pe_s[chn] = msum;
1874 //gfc->pe[chn] -= 150;
1877 /***************************************************************
1878 * determine final block type
1879 ***************************************************************/
1881 for (chn=0; chn<gfc->channels_out; chn++) {
1882 blocktype[chn] = NORM_TYPE;
1885 if (gfp->short_blocks == short_block_coupled) {
1886 /* force both channels to use the same block type */
1887 /* this is necessary if the frame is to be encoded in ms_stereo. */
1888 /* But even without ms_stereo, FhG does this */
1889 int bothlong= (uselongblock[0] && uselongblock[1]);
1896 /* update the blocktype of the previous granule, since it depends on what
1897 * happend in this granule */
1898 for (chn=0; chn<gfc->channels_out; chn++) {
1899 if ( uselongblock[chn])
1900 { /* no attack : use long blocks */
1901 assert( gfc->blocktype_old[chn] != START_TYPE );
1902 switch( gfc->blocktype_old[chn] )
1906 blocktype[chn] = NORM_TYPE;
1909 blocktype[chn] = STOP_TYPE;
1913 /* attack : use short blocks */
1914 blocktype[chn] = SHORT_TYPE;
1915 if ( gfc->blocktype_old[chn] == NORM_TYPE ) {
1916 gfc->blocktype_old[chn] = START_TYPE;
1918 if ( gfc->blocktype_old[chn] == STOP_TYPE ) {
1919 gfc->blocktype_old[chn] = SHORT_TYPE ;
1923 blocktype_d[chn] = gfc->blocktype_old[chn]; /* value returned to calling program */
1924 gfc->blocktype_old[chn] = blocktype[chn]; /* save for next call to l3psy_anal */
1926 if (gfc->presetTune.use) {
1927 if (blocktype_d[chn] != NORM_TYPE)
1928 gfc->presetTune.quantcomp_current = gfc->presetTune.quantcomp_type_s;
1930 gfc->presetTune.quantcomp_current = gfp->experimentalX;
1932 if (gfc->ATH->adjust >= gfc->presetTune.athadjust_switch_level &&
1933 blocktype_d[chn] == NORM_TYPE &&
1934 gfc->presetTune.quantcomp_alt_type > -1) {
1935 gfc->presetTune.quantcomp_current = gfc->presetTune.quantcomp_alt_type;
1940 /*********************************************************************
1941 * compute the value of PE to return (one granule delay)
1942 *********************************************************************/
1943 for(chn=0;chn<numchn;chn++)
1946 if (blocktype_d[chn] == SHORT_TYPE) {
1947 percep_entropy[chn] = pe_s[chn];
1949 percep_entropy[chn] = pe_l[chn];
1951 if (gfp->analysis) gfc->pinfo->pe[gr_out][chn] = percep_entropy[chn];
1953 if (blocktype_d[0] == SHORT_TYPE) {
1954 percep_MS_entropy[chn-2] = pe_s[chn];
1956 percep_MS_entropy[chn-2] = pe_l[chn];
1958 if (gfp->analysis) gfc->pinfo->pe[gr_out][chn] = percep_MS_entropy[chn-2];
1970 * The spreading function. Values returned in units of energy
1972 FLOAT8 s3_func(FLOAT8 bark) {
1974 FLOAT8 tempx,x,tempy,temp;
1976 if (tempx>=0) tempx *= 3;
1979 if(tempx>=0.5 && tempx<=2.5)
1982 x = 8.0 * (temp*temp - 2.0 * temp);
1986 tempy = 15.811389 + 7.5*tempx - 17.5*sqrt(1.0+tempx*tempx);
1988 if (tempy <= -60.0) return 0.0;
1990 tempx = exp( (x + tempy)*LN_TO_LOG10 );
1992 /* Normalization. The spreading function should be normalized so that:
1995 | s3 [ bark ] d(bark) = 1
2011 int L3para_read(lame_global_flags * gfp, FLOAT8 sfreq, int *numlines_l,int *numlines_s,
2013 FLOAT8 s3_l[CBANDS][CBANDS], FLOAT8 s3_s[CBANDS][CBANDS],
2015 int *bu_l, int *bo_l, FLOAT8 *w1_l, FLOAT8 *w2_l,
2016 int *bu_s, int *bo_s, FLOAT8 *w1_s, FLOAT8 *w2_s,
2017 int *npart_l_orig,int *npart_l,int *npart_s_orig,int *npart_s)
2019 lame_internal_flags *gfc=gfp->internal_flags;
2022 FLOAT8 bval_l[CBANDS], bval_s[CBANDS];
2023 FLOAT8 bval_l_width[CBANDS], bval_s_width[CBANDS];
2025 int partition[HBLKSIZE];
2029 /* compute numlines, the number of spectral lines in each partition band */
2030 /* each partition band should be about DELBARK wide. */
2032 for(i=0;i<CBANDS;i++)
2034 FLOAT8 ji, bark1,bark2;
2038 j2 = Min(j2,BLKSIZE/2);
2041 /* find smallest j2 >= j so that (bark - bark_l[i-1]) < DELBARK */
2043 bark1 = freq2bark(sfreq*ji/BLKSIZE);
2047 bark2 = freq2bark(sfreq*ji/BLKSIZE);
2048 } while ((bark2 - bark1) < DELBARK && j2<=BLKSIZE/2);
2050 for (k=j; k<j2; ++k)
2052 numlines_l[i]=(j2-j);
2054 if (j > BLKSIZE/2) break;
2056 *npart_l_orig = i+1;
2057 assert(*npart_l_orig <= CBANDS);
2059 /* compute which partition bands are in which scalefactor bands */
2060 { int i1,i2,sfb,start,end;
2062 for ( sfb = 0; sfb < SBMAX_l; sfb++ ) {
2063 start = gfc->scalefac_band.l[ sfb ];
2064 end = gfc->scalefac_band.l[ sfb+1 ];
2065 freq1 = sfreq*(start-.5)/(2*576);
2066 freq2 = sfreq*(end-1+.5)/(2*576);
2068 i1 = floor(.5 + BLKSIZE*freq1/sfreq);
2070 i2 = floor(.5 + BLKSIZE*freq2/sfreq);
2071 if (i2>BLKSIZE/2) i2=BLKSIZE/2;
2073 // DEBUGF(gfc,"longblock: old: (%i,%i) new: (%i,%i) %i %i \n",bu_l[sfb],bo_l[sfb],partition[i1],partition[i2],i1,i2);
2077 bu_l[sfb]=partition[i1];
2078 bo_l[sfb]=partition[i2];
2084 /* compute bark value and ATH of each critical band */
2086 for ( i = 0; i < *npart_l_orig; i++ ) {
2089 /* FLOAT8 mval,freq; */
2091 // Calculating the medium bark scaled frequency of the spectral lines
2092 // from j ... j + numlines[i]-1 (=spectral lines in parition band i)
2094 k = numlines_l[i] - 1;
2095 bark1 = freq2bark(sfreq*(j+0)/BLKSIZE);
2096 bark2 = freq2bark(sfreq*(j+k)/BLKSIZE);
2097 bval_l[i] = .5*(bark1+bark2);
2099 bark1 = freq2bark(sfreq*(j+0-.5)/BLKSIZE);
2100 bark2 = freq2bark(sfreq*(j+k+.5)/BLKSIZE);
2101 bval_l_width[i] = bark2-bark1;
2103 gfc->ATH->cb [i] = 1.e37; // preinit for minimum search
2104 for (k=0; k < numlines_l[i]; k++, j++) {
2105 FLOAT8 freq = sfreq*j/(1000.0*BLKSIZE);
2107 assert( freq <= 24 ); // or only '<'
2108 // freq = Min(.1,freq); // ATH below 100 Hz constant, not further climbing
2109 level = ATHformula (freq*1000, gfp) - 20; // scale to FFT units; returned value is in dB
2110 level = pow ( 10., 0.1*level ); // convert from dB -> energy
2111 level *= numlines_l [i];
2112 if ( level < gfc->ATH->cb [i] )
2113 gfc->ATH->cb [i] = level;
2119 /* MINVAL. For low freq, the strength of the masking is limited by minval
2120 * this is an ISO MPEG1 thing, dont know if it is really needed */
2121 for(i=0;i<*npart_l_orig;i++){
2122 double x = (-20+bval_l[i]*20.0/10.0);
2123 if (bval_l[i]>10) x = 0;
2124 minval[i]=pow(10.0,x/10);
2125 gfc->PSY->prvTonRed[i] = minval[i];
2134 /************************************************************************/
2136 /************************************************************************/
2138 /* compute numlines */
2140 for(i=0;i<CBANDS;i++)
2142 FLOAT8 ji, bark1,bark2;
2146 j2 = Min(j2,BLKSIZE_s/2);
2149 /* find smallest j2 >= j so that (bark - bark_s[i-1]) < DELBARK */
2151 bark1 = freq2bark(sfreq*ji/BLKSIZE_s);
2155 bark2 = freq2bark(sfreq*ji/BLKSIZE_s);
2157 } while ((bark2 - bark1) < DELBARK && j2<=BLKSIZE_s/2);
2159 for (k=j; k<j2; ++k)
2161 numlines_s[i]=(j2-j);
2163 if (j > BLKSIZE_s/2) break;
2165 *npart_s_orig = i+1;
2166 assert(*npart_s_orig <= CBANDS);
2168 /* compute which partition bands are in which scalefactor bands */
2169 { int i1,i2,sfb,start,end;
2171 for ( sfb = 0; sfb < SBMAX_s; sfb++ ) {
2172 start = gfc->scalefac_band.s[ sfb ];
2173 end = gfc->scalefac_band.s[ sfb+1 ];
2174 freq1 = sfreq*(start-.5)/(2*192);
2175 freq2 = sfreq*(end-1+.5)/(2*192);
2177 i1 = floor(.5 + BLKSIZE_s*freq1/sfreq);
2179 i2 = floor(.5 + BLKSIZE_s*freq2/sfreq);
2180 if (i2>BLKSIZE_s/2) i2=BLKSIZE_s/2;
2182 //DEBUGF(gfc,"shortblock: old: (%i,%i) new: (%i,%i) %i %i \n",bu_s[sfb],bo_s[sfb], partition[i1],partition[i2],i1,i2);
2186 bu_s[sfb]=partition[i1];
2187 bo_s[sfb]=partition[i2];
2196 /* compute bark values of each critical band */
2198 for(i=0;i<*npart_s_orig;i++)
2201 FLOAT8 bark1,bark2,snr;
2202 k = numlines_s[i] - 1;
2204 bark1 = freq2bark (sfreq*(j+0)/BLKSIZE_s);
2205 bark2 = freq2bark (sfreq*(j+k)/BLKSIZE_s);
2206 bval_s[i] = .5*(bark1+bark2);
2208 bark1 = freq2bark (sfreq*(j+0-.5)/BLKSIZE_s);
2209 bark2 = freq2bark (sfreq*(j+k+.5)/BLKSIZE_s);
2210 bval_s_width[i] = bark2-bark1;
2217 snr = -4.5 * (bval_s[i]-13)/(24.0-13.0) +
2218 -8.25*(bval_s[i]-24)/(13.0-24.0);
2220 SNR[i]=pow(10.0,snr/10.0);
2228 /************************************************************************
2229 * Now compute the spreading function, s[j][i], the value of the spread-*
2230 * ing function, centered at band j, for band i, store for later use *
2231 ************************************************************************/
2232 /* i.e.: sum over j to spread into signal barkval=i
2233 NOTE: i and j are used opposite as in the ISO docs */
2234 for(i=0;i<*npart_l_orig;i++) {
2235 for(j=0;j<*npart_l_orig;j++) {
2236 s3_l[i][j]=s3_func(bval_l[i]-bval_l[j])*bval_l_width[j];
2239 for(i=0;i<*npart_s_orig;i++) {
2240 for(j=0;j<*npart_s_orig;j++) {
2241 s3_s[i][j]=s3_func(bval_s[i]-bval_s[j])*bval_s_width[j];
2249 /* npart_l_orig = number of partition bands before convolution */
2250 /* npart_l = number of partition bands after convolution */
2252 *npart_l=bo_l[NBPSY_l-1]+1;
2253 *npart_s=bo_s[NBPSY_s-1]+1;
2255 assert(*npart_l <= *npart_l_orig);
2256 assert(*npart_s <= *npart_s_orig);
2259 /* setup stereo demasking thresholds */
2260 /* formula reverse enginerred from plot in paper */
2261 for ( i = 0; i < NBPSY_s; i++ ) {
2263 arg = freq2bark(sfreq*gfc->scalefac_band.s[i]/(2*192));
2264 arg = (Min(arg, 15.5)/15.5);
2266 mld = 1.25*(1-cos(PI*arg))-2.5;
2267 gfc->mld_s[i] = pow(10.0,mld);
2269 for ( i = 0; i < NBPSY_l; i++ ) {
2271 arg = freq2bark(sfreq*gfc->scalefac_band.l[i]/(2*576));
2272 arg = (Min(arg, 15.5)/15.5);
2274 mld = 1.25*(1-cos(PI*arg))-2.5;
2275 gfc->mld_l[i] = pow(10.0,mld);
2278 #define temporalmask_sustain_sec 0.01
2280 /* setup temporal masking */
2281 gfc->decay = exp(-1.0*LOG10/(temporalmask_sustain_sec*sfreq/192.0));
2296 int psymodel_init(lame_global_flags *gfp)
2298 lame_internal_flags *gfc=gfp->internal_flags;
2299 int i,j,b,sb,k,samplerate;
2301 FLOAT8 s3_s[CBANDS][CBANDS];
2302 FLOAT8 s3_l[CBANDS][CBANDS];
2303 int numberOfNoneZero;
2305 samplerate = gfp->out_samplerate;
2306 gfc->ms_ener_ratio_old=.25;
2307 gfc->blocktype_old[0]=STOP_TYPE;
2308 gfc->blocktype_old[1]=STOP_TYPE;
2309 gfc->blocktype_old[0]=SHORT_TYPE;
2310 gfc->blocktype_old[1]=SHORT_TYPE;
2312 for (i=0; i<4; ++i) {
2313 for (j=0; j<CBANDS; ++j) {
2314 gfc->nb_1[i][j]=1e20;
2315 gfc->nb_2[i][j]=1e20;
2317 for ( sb = 0; sb < NBPSY_l; sb++ ) {
2318 gfc->en[i].l[sb] = 1e20;
2319 gfc->thm[i].l[sb] = 1e20;
2321 for (j=0; j<3; ++j) {
2322 for ( sb = 0; sb < NBPSY_s; sb++ ) {
2323 gfc->en[i].s[sb][j] = 1e20;
2324 gfc->thm[i].s[sb][j] = 1e20;
2328 for (i=0; i<4; ++i) {
2329 for (j=0; j<3; ++j) {
2330 for ( sb = 0; sb < NBPSY_s; sb++ ) {
2331 gfc->nsPsy.last_thm[i][sb][j] = 1e20;
2337 gfc->nsPsy.last_en_subshort[i][j] = 100;
2339 gfc->nsPsy.last_attacks[i][j] = 0;
2340 gfc->nsPsy.pe_l[i] = gfc->nsPsy.pe_s[i] = 0;
2346 /* gfp->cwlimit = sfreq*j/1024.0; */
2347 gfc->cw_lower_index=6;
2348 gfc->cw_upper_index = gfc->PSY->cwlimit*1024.0/gfp->out_samplerate;
2349 gfc->cw_upper_index=Min(HBLKSIZE-4,gfc->cw_upper_index); /* j+3 < HBLKSIZE-1 */
2350 gfc->cw_upper_index=Max(6,gfc->cw_upper_index);
2352 for ( j = 0; j < HBLKSIZE; j++ )
2357 i=L3para_read( gfp,(FLOAT8) samplerate,gfc->numlines_l,gfc->numlines_s,
2358 gfc->minval,s3_l,s3_s,gfc->SNR_s,gfc->bu_l,
2359 gfc->bo_l,gfc->w1_l,gfc->w2_l, gfc->bu_s,gfc->bo_s,
2360 gfc->w1_s,gfc->w2_s,&gfc->npart_l_orig,&gfc->npart_l,
2361 &gfc->npart_s_orig,&gfc->npart_s );
2362 if (i!=0) return -1;
2366 /* npart_l_orig = number of partition bands before convolution */
2367 /* npart_l = number of partition bands after convolution */
2369 numberOfNoneZero = 0;
2370 for (i=0; i<gfc->npart_l; i++) {
2371 for (j = 0; j < gfc->npart_l_orig; j++) {
2372 if (s3_l[i][j] != 0.0)
2375 gfc->s3ind[i][0] = j;
2377 for (j = gfc->npart_l_orig - 1; j > 0; j--) {
2378 if (s3_l[i][j] != 0.0)
2381 gfc->s3ind[i][1] = j;
2382 numberOfNoneZero += (gfc->s3ind[i][1] - gfc->s3ind[i][0] + 1);
2384 gfc->s3_ll = malloc(sizeof(FLOAT8)*numberOfNoneZero);
2389 for (i=0; i<gfc->npart_l; i++) {
2390 for (j = gfc->s3ind[i][0]; j <= gfc->s3ind[i][1]; j++) {
2391 gfc->s3_ll[k++] = s3_l[i][j];
2397 numberOfNoneZero = 0;
2398 for (i=0; i<gfc->npart_s; i++) {
2399 for (j = 0; j < gfc->npart_s_orig; j++) {
2400 if (s3_s[i][j] != 0.0)
2403 gfc->s3ind_s[i][0] = j;
2405 for (j = gfc->npart_s_orig - 1; j > 0; j--) {
2406 if (s3_s[i][j] != 0.0)
2409 gfc->s3ind_s[i][1] = j;
2410 numberOfNoneZero += (gfc->s3ind_s[i][1] - gfc->s3ind_s[i][0] + 1);
2412 gfc->s3_ss = malloc(sizeof(FLOAT8)*numberOfNoneZero);
2420 #include "debugscalefac.c"
2425 if (gfc->nsPsy.use) {
2426 /* long block spreading function normalization */
2427 for ( b = 0;b < gfc->npart_l; b++ ) {
2428 for ( k = gfc->s3ind[b][0]; k <= gfc->s3ind[b][1]; k++ ) {
2429 // spreading function has been properly normalized by
2430 // multiplying by DELBARK/.6609193 = .515.
2431 // It looks like Naoki was
2432 // way ahead of me and added this factor here!
2433 // it is no longer needed.
2434 //gfc->s3_l[b][k] *= 0.5;
2437 /* short block spreading function normalization */
2438 // no longer needs to be normalized, but nspsytune wants
2439 // SNR_s applied here istead of later to save CPU cycles
2440 for ( b = 0;b < gfc->npart_s; b++ ) {
2442 for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ ) {
2445 for ( k = gfc->s3ind_s[b][0]; k <= gfc->s3ind_s[b][1]; k++ ) {
2446 s3_s[b][k] *= gfc->SNR_s[b] /* / norm */;
2453 if (gfc->nsPsy.use) {
2455 /* spread only from npart_l bands. Normally, we use the spreading
2456 * function to convolve from npart_l_orig down to npart_l bands
2458 for(b=0;b<gfc->npart_l;b++)
2459 if (gfc->s3ind[b][1] > gfc->npart_l-1) gfc->s3ind[b][1] = gfc->npart_l-1;
2463 for (i=0; i<gfc->npart_s; i++) {
2464 for (j = gfc->s3ind_s[i][0]; j <= gfc->s3ind_s[i][1]; j++) {
2465 gfc->s3_ss[k++] = s3_s[i][j];
2471 /* init. for loudness approx. -jd 2001 mar 27*/
2472 gfc->loudness_sq_save[0] = 0.0;
2473 gfc->loudness_sq_save[1] = 0.0;
2484 /* psycho_loudness_approx
2486 in: energy - BLKSIZE/2 elements of frequency magnitudes ^ 2
2487 gfp - uses out_samplerate, ATHtype (also needed for ATHformula)
2488 returns: loudness^2 approximation, a positive value roughly tuned for a value
2489 of 1.0 for signals near clipping.
2490 notes: When calibrated, feeding this function binary white noise at sample
2491 values +32767 or -32768 should return values that approach 3.
2492 ATHformula is used to approximate an equal loudness curve.
2493 future: Data indicates that the shape of the equal loudness curve varies
2494 with intensity. This function might be improved by using an equal
2495 loudness curve shaped for typical playback levels (instead of the
2496 ATH, that is shaped for the threshold). A flexible realization might
2497 simply bend the existing ATH curve to achieve the desired shape.
2498 However, the potential gain may not be enough to justify an effort.
2501 psycho_loudness_approx( FLOAT *energy, lame_global_flags *gfp )
2504 static int eql_type = -1;
2505 static FLOAT eql_w[BLKSIZE/2];/* equal loudness weights (based on ATH) */
2506 const FLOAT vo_scale= 1./( 14752 ); /* tuned for output level */
2507 /* (sensitive to energy scale) */
2508 FLOAT loudness_power;
2510 if( eql_type != gfp->ATHtype ) {
2511 /* compute equal loudness weights (eql_w) */
2513 FLOAT freq_inc = gfp->out_samplerate / (BLKSIZE);
2514 FLOAT eql_balance = 0.0;
2515 eql_type = gfp->ATHtype;
2517 for( i = 0; i < BLKSIZE/2; ++i ) {
2519 /* convert ATH dB to relative power (not dB) */
2520 /* to determine eql_w */
2521 eql_w[i] = 1. / pow( 10, ATHformula( freq, gfp ) / 10 );
2522 eql_balance += eql_w[i];
2524 eql_balance = 1 / eql_balance;
2525 for( i = BLKSIZE/2; --i >= 0; ) { /* scale weights */
2526 eql_w[i] *= eql_balance;
2530 loudness_power = 0.0;
2531 for( i = 0; i < BLKSIZE/2; ++i ) { /* apply weights to power in freq. bands*/
2532 loudness_power += energy[i] * eql_w[i];
2534 loudness_power /= (BLKSIZE/2);
2535 loudness_power *= vo_scale * vo_scale;
2537 return( loudness_power );