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