1
2 #define TONAL_WIDTH 0.1 /* When more than one tonal component is within
3 this width in Bark, the weaker one(s) are
4 eliminated */
5
6
7 #include "common.h"
8 #include "encoder.h"
9 #include "psycho_1.h"
10
11 /**********************************************************************
12 *
13 * This module implements the psychoacoustic model I for the
14 * MPEG encoder layer II. It uses simplified tonal and noise masking
15 * threshold analysis to generate SMR for the encoder bit allocation
16 * routine.
17 *
18 **********************************************************************/
19 #ifdef NEWP1TABLES
20
psycho_1_read_cbound(int lay,int freq,int * crit_band)21 int *psycho_1_read_cbound (int lay, int freq, int *crit_band)
22 /* this function reads in critical band boundaries */
23 {
24
25 #include "critband.h"
26 //static const int FirstCriticalBand[7][27] = {...
27 int *cbound;
28
29 int i, k;
30
31 if ((lay < 1) || (lay > 2)) {
32 printf ("Internal error (read_cbound())\n");
33 exit(99);
34 }
35 if ((freq < 0) || (freq > 6) || (freq == 3)) {
36 printf ("Internal error (read_cbound())\n");
37 exit(99);
38 }
39
40 *crit_band = SecondCriticalBand[freq][0];
41 cbound = (int *) mem_alloc (sizeof (int) * *crit_band, "cbound");
42 for (i = 0; i < *crit_band; i++) {
43 k = SecondCriticalBand[freq][i + 1];
44 if (k != 0) {
45 cbound[i] = k;
46 } else {
47 printf ("Internal error (read_cbound())\n");
48 exit(99);
49 }
50 }
51 return(cbound); // return the pointer to the table of critband boundaries.
52 }
53
psycho_1_read_freq_band(int * sub_size,g_ptr * ltg,int lay,int freq)54 void psycho_1_read_freq_band (int *sub_size, g_ptr *ltg, int lay, int freq)
55 {
56
57 #include "freqtable.h"
58
59 int i, k;
60
61 if ((freq < 0) || (freq > 6) || (freq == 3)) {
62 printf ("Internal error (read_freq_band())\n");
63 return;
64 }
65
66 /* read input for freq. subbands */
67
68 *sub_size = SecondFreqEntries[freq] + 1;
69 *ltg = (g_ptr) mem_alloc (sizeof (g_thres) * *sub_size, "ltg");
70 (*ltg)[0].line = 0; /* initialize global masking threshold */
71 (*ltg)[0].bark = 0.0;
72 (*ltg)[0].hear = 0.0;
73 for (i = 1; i < *sub_size; i++) {
74 k = SecondFreqSubband[freq][i - 1].line;
75 if (k != 0) {
76 (*ltg)[i].line = k;
77 (*ltg)[i].bark = SecondFreqSubband[freq][i - 1].bark;
78 (*ltg)[i].hear = SecondFreqSubband[freq][i - 1].hear;
79 } else {
80 printf ("Internal error (read_freq_band())\n");
81 return;
82 }
83 }
84 }
85
86 #else
read_crit_band(int lay,int freq)87 int read_crit_band (int lay, int freq)
88 {
89 int crit_band;
90 FILE *fp;
91 char r[16], t[80];
92
93 strcpy (r, "2cb1");
94 r[0] = (char) lay + '0';
95 r[3] = (char) freq + '0';
96 if (!(fp = OpenTableFile (r))) { /* check boundary values */
97 printf ("Please check %s boundary table\n", r);
98 exit (0);
99 }
100 fgets (t, 80, fp);
101 sscanf (t, "%d\n", &crit_band);
102 fclose (fp);
103 return (crit_band);
104 }
105
read_cbound(int lay,int freq,int crit_band,int * cbound)106 void read_cbound (int lay, int freq, int crit_band, int *cbound)
107 /* this function reads in critical band boundaries */
108 {
109 int i, j, k;
110 FILE *fp;
111 char r[16], t[80];
112
113 strcpy (r, "2cb1");
114 r[0] = (char) lay + '0';
115 r[3] = (char) freq + '0';
116 if (!(fp = OpenTableFile (r))) { /* check boundary values */
117 printf ("Please check %s boundary table\n", r);
118 exit (0);
119 }
120 fgets (t, 80, fp); /* skip input for critical bands */
121 sscanf (t, "%d\n", &k);
122 for (i = 0; i < crit_band; i++) {
123 fgets (t, 80, fp);
124 sscanf (t, "%d %d\n", &j, &k);
125 if (i == j)
126 cbound[j] = k;
127 else { /* error */
128 printf ("Please check index %d in cbound table %s\n", i, r);
129 exit (0);
130 }
131 }
132 fclose (fp);
133 }
134
read_freq_band(int * sub_size,g_ptr * ltg,int lay,int freq)135 void read_freq_band (int *sub_size, g_ptr * ltg, int lay, int freq)
136 /* this function reads in frequency bands and bark values */
137 {
138 int i, j, k;
139 double b, c;
140 FILE *fp;
141 char r[16], t[80];
142
143 strcpy (r, "2th1");
144 r[0] = (char) lay + '0';
145 r[3] = (char) freq + '0';
146 if (!(fp = OpenTableFile (r))) { /* check freq. values */
147 printf ("Please check frequency and cband table %s\n", r);
148 exit (0);
149 }
150 fgets (t, 80, fp); /* read input for freq. subbands */
151 sscanf (t, "%d\n", sub_size);
152 *ltg = (g_ptr) mem_alloc (sizeof (g_thres) * (*sub_size), "ltg");
153 (*ltg)[0].line = 0; /* initialize global masking threshold */
154 (*ltg)[0].bark = 0;
155 (*ltg)[0].hear = 0;
156 for (i = 1; i < (*sub_size); i++) { /* continue to read freq. subband */
157 fgets (t, 80, fp); /* and assign */
158 sscanf (t, "%d %d %lf %lf\n", &j, &k, &b, &c);
159 if (i == j) {
160 (*ltg)[j].line = k;
161 (*ltg)[j].bark = b;
162 (*ltg)[j].hear = c;
163 } else { /* error */
164 printf ("Please check index %d in freq-cb table %s\n", i, r);
165 exit (0);
166 }
167 }
168 fclose (fp);
169 }
170 #endif
171
172
make_map(int sub_size,mask * power,g_thres * ltg)173 void make_map (int sub_size, mask * power, g_thres * ltg)
174 /* this function calculates the global masking threshold */
175 {
176 int i, j;
177
178 for (i = 1; i < sub_size; i++)
179 for (j = ltg[i - 1].line; j <= ltg[i].line; j++)
180 power[j].map = i;
181 }
182
183 #ifdef DELETEME
add_db(double a,double b)184 double add_db (double a, double b)
185 {
186 a = pow (10.0, a / 10.0);
187 b = pow (10.0, b / 10.0);
188 return 10 * log10 (a + b);
189 }
190 #endif
191
192 #define DBTAB 1000
193 double dbtable[DBTAB];
194
psycho_1_init_add_db(void)195 void psycho_1_init_add_db (void)
196 {
197 int i;
198 double x;
199 for (i = 0; i < DBTAB; i++) {
200 x = (double) i / 10.0;
201 dbtable[i] = 10 * log10 (1 + pow (10.0, x / 10.0)) - x;
202 }
203 }
204
add_db(double a,double b)205 INLINE double add_db (double a, double b)
206 {
207 /* MFC - if the difference between a and b is large (>99), then just return the
208 largest one. (about 10% of the time)
209 - For differences between 0 and 99, return the largest value, but add
210 in a pre-calculated difference value.
211 - the value 99 was chosen arbitarily.
212 - maximum (a-b) i've seen is 572 */
213 FLOAT fdiff;
214 int idiff;
215 fdiff = (10.0 * (a - b));
216
217 if (fdiff > 990.0) {
218 return a;
219 }
220 if (fdiff < -990.0) {
221 return (b);
222 }
223
224 idiff = (int) fdiff;
225 if (idiff >= 0) {
226 return (a + dbtable[idiff]);
227 }
228
229 return (b + dbtable[-idiff]);
230 }
231
232 /****************************************************************
233 *
234 * Fast Fourier transform of the input samples.
235 *
236 ****************************************************************/
237
II_f_f_t(double * sample,mask * power)238 void II_f_f_t (double *sample, mask * power)
239
240 { /* this function calculates an */
241 /* FFT analysis for the freq. */
242 /* domain */
243 int i, j, k, ll, l = 0;
244 int ip, le, le1;
245 double t_r, t_i, u_r, u_i;
246 static int M, MM1, init = 0, N, NV2, NM1;
247 double *x_r, *x_i, *energy;
248 static int *rev;
249 static double *w_r, *w_i;
250
251 x_r = (double *) mem_alloc (sizeof (DFFT), "x_r");
252 x_i = (double *) mem_alloc (sizeof (DFFT), "x_i");
253 energy = (double *) mem_alloc (sizeof (DFFT), "energy");
254
255 if (!init) {
256 rev = (int *) mem_alloc (sizeof (IFFT), "rev");
257 w_r = (double *) mem_alloc (sizeof (D10), "w_r");
258 w_i = (double *) mem_alloc (sizeof (D10), "w_i");
259 M = 10;
260 MM1 = 9;
261 N = FFT_SIZE;
262 NV2 = FFT_SIZE >> 1;
263 NM1 = FFT_SIZE - 1;
264 for (ll = 0; ll < M; ll++) {
265 le = 1 << (M - ll);
266 le1 = le >> 1;
267 w_r[ll] = cos (PI / le1);
268 w_i[ll] = -sin (PI / le1);
269 }
270 for (i = 0; i < FFT_SIZE; rev[i] = l, i++)
271 for (j = 0, l = 0; j < 10; j++) {
272 k = (i >> j) & 1;
273 l |= (k << (9 - j) ); //MFC
274 }
275 init = 1;
276 }
277
278 for (i = 0; i < FFT_SIZE; i++) {
279 x_r[i] = sample[i];
280 x_i[i] = energy[i] = 0;
281 }
282 /*
283 memcpy ((char *) x_r, (char *) sample, sizeof (double) * FFT_SIZE);
284 */
285 for (ll = 0; ll < MM1; ll++) {
286 le = 1 << (M - ll);
287 le1 = le >> 1;
288 u_r = 1;
289 u_i = 0;
290 for (j = 0; j < le1; j++) {
291 for (i = j; i < N; i += le) {
292 ip = i + le1;
293 t_r = x_r[i] + x_r[ip];
294 t_i = x_i[i] + x_i[ip];
295 x_r[ip] = x_r[i] - x_r[ip];
296 x_i[ip] = x_i[i] - x_i[ip];
297 x_r[i] = t_r;
298 x_i[i] = t_i;
299 t_r = x_r[ip];
300 x_r[ip] = x_r[ip] * u_r - x_i[ip] * u_i;
301 x_i[ip] = x_i[ip] * u_r + t_r * u_i;
302 }
303 t_r = u_r;
304 u_r = u_r * w_r[ll] - u_i * w_i[ll];
305 u_i = u_i * w_r[ll] + t_r * w_i[ll];
306 }
307 }
308 for (i = 0; i < N; i += 2) {
309 ip = i + 1;
310 t_r = x_r[i] + x_r[ip];
311 t_i = x_i[i] + x_i[ip];
312 x_r[ip] = x_r[i] - x_r[ip];
313 x_i[ip] = x_i[i] - x_i[ip];
314 x_r[i] = t_r;
315 x_i[i] = t_i;
316 energy[i] = x_r[i] * x_r[i] + x_i[i] * x_i[i];
317 }
318 for (i = 0; i < FFT_SIZE; i++)
319 if (i < rev[i]) {
320 t_r = energy[i];
321 energy[i] = energy[rev[i]];
322 energy[rev[i]] = t_r;
323 }
324 for (i = 0; i < HAN_SIZE; i++) {
325 /* calculate power density spectrum */
326 if (energy[i] < 1E-20)
327 energy[i] = 1E-20;
328 /* power calculation corrected with a factor 4, both positive
329 and negative frequencies exist, 1992-11-06 shn */
330 power[i].x = 10 * log10 (energy[i] * 4.0) + POWERNORM;
331 power[i].next = STOP;
332 power[i].type = FALSE;
333 }
334
335 mem_free ((void **) &x_r);
336 mem_free ((void **) &x_i);
337 mem_free ((void **) &energy);
338
339 }
340
341 /****************************************************************
342 *
343 * Window the incoming audio signal.
344 *
345 ****************************************************************/
346
II_hann_win(double * sample)347 void II_hann_win (double *sample)
348 { /* this function calculates a */
349 /* Hann window for PCM (input) *//* samples for a 1024-pt. FFT */
350 register int i;
351 register double sqrt_8_over_3;
352 static int init = 0;
353 static double *window;
354
355 if (!init) {
356 /* calculate window function for the Fourier transform */
357 window = (double *) mem_alloc (sizeof (DFFT), "window");
358 sqrt_8_over_3 = pow (8.0 / 3.0, 0.5);
359 for (i = 0; i < FFT_SIZE; i++) {
360 /* Hann window formula */
361 window[i] =
362 sqrt_8_over_3 * 0.5 * (1 -
363 cos (2.0 * PI * i / (FFT_SIZE - 1))) / FFT_SIZE;
364 }
365 init = 1;
366 }
367
368 for (i = 0; i < FFT_SIZE; i++)
369 sample[i] *= window[i];
370 }
371
372 /*******************************************************************
373 *
374 * This function finds the maximum spectral component in each
375 * subband and return them to the encoder for time-domain threshold
376 * determination.
377 *
378 *******************************************************************/
379
II_pick_max(mask * power,double * spike)380 void II_pick_max (mask * power, double *spike)
381 {
382 double max;
383 int i, j;
384
385 for (i = 0; i < HAN_SIZE; spike[i >> 4] = max, i += 16) /* calculate the */
386 for (j = 0, max = DBMIN; j < 16; j++) /* maximum spectral */
387 max = (max > power[i + j].x) ? max : power[i + j].x; /* component in each */
388 } /* subband from bound */
389 /* 4-16 */
390
391 /****************************************************************
392 *
393 * This function labels the tonal component in the power
394 * spectrum.
395 *
396 ****************************************************************/
397
II_tonal_label(mask * power,int * tone)398 void II_tonal_label (mask * power, int *tone)
399 { /* this function extracts (tonal) */
400 /* sinusoidals from the spectrum */
401 int i, j, last = LAST, first, run, last_but_one = LAST; /* dpwe */
402 double max;
403
404 *tone = LAST;
405 for (i = 2; i < HAN_SIZE - 12; i++) {
406 if (power[i].x > power[i - 1].x && power[i].x >= power[i + 1].x) {
407 power[i].type = TONE;
408 power[i].next = LAST;
409 if (last != LAST)
410 power[last].next = i;
411 else
412 first = *tone = i;
413 last = i;
414 }
415 }
416 last = LAST;
417 first = *tone;
418 *tone = LAST;
419 while (first != LAST) { /* the conditions for the tonal */
420 if (first < 2 || first > 499)
421 run = 0; /* otherwise k+/-j will be out of bounds */
422 else if (first < 62)
423 run = 2; /* components in layer II, which */
424 else if (first < 126)
425 run = 3; /* are the boundaries for calc. */
426 else if (first < 254)
427 run = 6; /* the tonal components */
428 else
429 run = 12;
430 max = power[first].x - 7; /* after calculation of tonal */
431 for (j = 2; j <= run; j++) /* components, set to local max */
432 if (max < power[first - j].x || max < power[first + j].x) {
433 power[first].type = FALSE;
434 break;
435 }
436 if (power[first].type == TONE) { /* extract tonal components */
437 int help = first;
438 if (*tone == LAST)
439 *tone = first;
440 while ((power[help].next != LAST) && (power[help].next - first) <= run)
441 help = power[help].next;
442 help = power[help].next;
443 power[first].next = help;
444 if ((first - last) <= run) {
445 if (last_but_one != LAST)
446 power[last_but_one].next = first;
447 }
448 if (first > 1 && first < 255) { /* calculate the sum of the */
449 double tmp; /* powers of the components */
450 tmp = add_db (power[first - 1].x, power[first + 1].x);
451 power[first].x = add_db (power[first].x, tmp);
452 }
453 for (j = 1; j <= run; j++) {
454 power[first - j].x = power[first + j].x = DBMIN;
455 power[first - j].next = power[first + j].next = STOP;
456 power[first - j].type = power[first + j].type = FALSE;
457 }
458 last_but_one = last;
459 last = first;
460 first = power[first].next;
461 } else {
462 int ll;
463 if (last == LAST); /* *tone = power[first].next; dpwe */
464 else
465 power[last].next = power[first].next;
466 ll = first;
467 first = power[first].next;
468 power[ll].next = STOP;
469 }
470 }
471 }
472
473 /****************************************************************
474 *
475 * This function groups all the remaining non-tonal
476 * spectral lines into critical band where they are replaced by
477 * one single line.
478 *
479 ****************************************************************/
480
noise_label(int crit_band,int * cbound,mask * power,int * noise,g_thres * ltg)481 void noise_label (int crit_band, int *cbound, mask * power, int *noise,
482 g_thres * ltg)
483 {
484 int i, j, centre, last = LAST;
485 double index, weight, sum;
486 /* calculate the remaining spectral */
487 for (i = 0; i < crit_band - 1; i++) { /* lines for non-tonal components */
488 for (j = cbound[i], weight = 0.0, sum = DBMIN; j < cbound[i + 1]; j++) {
489 if (power[j].type != TONE) {
490 if (power[j].x != DBMIN) {
491 sum = add_db (power[j].x, sum);
492 weight +=
493 pow (10.0, power[j].x / 10.0) * (ltg[power[j].map].bark - i);
494 power[j].x = DBMIN;
495 }
496 } /* check to see if the spectral line is low dB, and if */
497 } /* so replace the center of the critical band, which is */
498 /* the center freq. of the noise component */
499 if (sum <= DBMIN)
500 centre = (cbound[i + 1] + cbound[i]) / 2;
501 else {
502 index = weight / pow (10.0, sum / 10.0);
503 centre = cbound[i] + (int) (index * (double) (cbound[i + 1] - cbound[i]));
504 } /* locate next non-tonal component until finished; */
505 /* add to list of non-tonal components */
506 if (power[centre].type == TONE)
507 centre++;
508 if (last == LAST)
509 *noise = centre;
510 else {
511 power[centre].next = LAST;
512 power[last].next = centre;
513 }
514 power[centre].x = sum;
515 power[centre].type = NOISE;
516 last = centre;
517 }
518 }
519
520 /****************************************************************
521 *
522 * This function reduces the number of noise and tonal
523 * component for further threshold analysis.
524 *
525 ****************************************************************/
526
subsampling(mask * power,g_thres * ltg,int * tone,int * noise)527 void subsampling (mask * power, g_thres * ltg, int *tone, int *noise)
528 {
529 int i, old;
530
531 i = *tone;
532 old = STOP;
533 /* calculate tonal components for */
534 while (i != LAST) {
535 /* reduction of spectral lines */
536 if (power[i].x < ltg[power[i].map].hear) {
537 power[i].type = FALSE;
538 power[i].x = DBMIN;
539 if (old == STOP)
540 *tone = power[i].next;
541 else
542 power[old].next = power[i].next;
543 } else
544 old = i;
545 i = power[i].next;
546 }
547 i = *noise;
548 old = STOP;
549 /* calculate non-tonal components for */
550 while (i != LAST) {
551 /* reduction of spectral lines */
552 if (power[i].x < ltg[power[i].map].hear) {
553 power[i].type = FALSE;
554 power[i].x = DBMIN;
555 if (old == STOP)
556 *noise = power[i].next;
557 else
558 power[old].next = power[i].next;
559 } else
560 old = i;
561 i = power[i].next;
562 }
563 i = *tone;
564 old = STOP;
565 while (i != LAST) {
566 /* if more than one */
567 if (power[i].next == LAST)
568 break; /* tonal component */
569 if (ltg[power[power[i].next].map].bark - /* is less than .5 */
570 ltg[power[i].map].bark < TONAL_WIDTH) {
571 /* bark, take the */
572 if (power[power[i].next].x > power[i].x) {
573 /* maximum */
574 if (old == STOP)
575 *tone = power[i].next;
576 else
577 power[old].next = power[i].next;
578 power[i].type = FALSE;
579 power[i].x = DBMIN;
580 } else {
581 power[power[i].next].type = FALSE;
582 power[power[i].next].x = DBMIN;
583 power[i].next = power[power[i].next].next;
584 }
585 } else
586 old = i;
587 i = power[i].next;
588 }
589 }
590
591
592 /* ----------------------------------------------------------------------------
593 The masking function parameters are here set according to the parameters in
594 the IRT real time implementation. The constant definitions are for convenience.
595 1993-07-23 shn
596 ---------------------------------------------------------------------------- */
597
598 #define AV_TONAL_K -9.0 /* Masking index, tonal, constant part [dB] */
599 #define AV_NOISE_K -5.0 /* Masking index, noisy, constant part [dB] */
600 #define AV_TONAL_DZ -0.3 /* Masking index, tonal, CBR dependence [dB/Bark] */
601 #define AV_NOISE_DZ -0.3 /* Masking index, noisy, CBR dependence [dB/Bark] */
602
603 #define LOW_LIM_1 -1.0 /* 1st lower slope from 0 to LOW_LIM_1 [Bark] */
604 #define LOW_LIM_2 -3.0 /* 2nd lower slope from LOW_LIM_1 to LOW_LIM_2 [Bark] */
605
606 #define LOW_DZ_K_1 6.0 /* 1st lower slope, constant part [dB/Bark] */
607 #define LOW_DZ_SPL_1 0.4 /* 1st lower slope, SPL dependence [dB/(Bark*dB)] */
608 #define LOW_DZ_MIN_1 17.0 /* 1st lower slope, minimum value [dB/Bark] */
609 #define LOW_DZ_2 17.0 /* 2nd lower slope [dB/Bark] */
610
611 #define UP_LIM_1 1.0 /* 1st upper slope from 0 to UP_LIM_1 [Bark] */
612 #define UP_LIM_2 8.0 /* 2nd upper slope from UP_LIM_1 to UP_LIM_2 [Bark] */
613
614 #define UP_DZ_1 -18.0 /* 1st upper slope, constant part [dB/Bark] */
615 #define UP_SPL_1 0.0 /* 1st upper slope, SPL dependence [dB/(Bark*dB)] */
616 #define UP_DZ_2 -17.0 /* 2nd upper slope, constant part [dB/Bark] */
617 #define UP_SPL_2 -0.1 /* 2nd upper slope, SPL dependence [dB/(Bark*dB)] */
618
619 #define H_THR_OFFSET -12.0 /* Hearing threshold offset [dB] */
620 #define H_THR_OS_BR 0 /*96*/ /* Minimum datarate for offset, [kbit/s per channel] */
621
622 #define MASK_ADD 2.0 /* Addition of maskers [dB] */
623 #define QUIET_ADD 3.0 /* Addition of masker and threshold in quiet [dB] */
624
625 /****************************************************************
626 *
627 * This function calculates the individual threshold and
628 * sum with the quiet threshold to find the global threshold.
629 *
630 ****************************************************************/
631
threshold(int sub_size,mask * power,g_thres * ltg,int * tone,int * noise,int bit_rate)632 void threshold (int sub_size, mask * power, g_thres * ltg, int *tone,
633 int *noise, int bit_rate)
634 {
635 int k, t;
636 double z, dz, spl, vf=0, tmps;
637
638 for (k = 1; k < sub_size; k++) { /* Target frequencies */
639 ltg[k].x = DBMIN;
640 t = *tone; /* calculate individual masking threshold */
641 while (t != LAST) { /* for tonal components, t, to find LTG */
642 z = ltg[power[t].map].bark; /* critical band rate of masker */
643 dz = ltg[k].bark - z; /* distance of bark value */
644 spl = power[t].x; /* sound pressure level of masker */
645
646 if (dz >= LOW_LIM_2 && dz < UP_LIM_2) {
647 tmps = spl + AV_TONAL_K + AV_TONAL_DZ * z;
648
649 /* masking function for lower & upper slopes */
650 if (LOW_LIM_2 <= dz && dz < LOW_LIM_1)
651 if (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1 > LOW_DZ_MIN_1)
652 vf = LOW_DZ_2 * (dz - LOW_LIM_1) +
653 (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1) * LOW_LIM_1;
654 else
655 vf = LOW_DZ_2 * (dz - LOW_LIM_1) + LOW_DZ_MIN_1 * LOW_LIM_1;
656 else if (LOW_LIM_1 <= dz && dz < 0)
657 if (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1 > LOW_DZ_MIN_1)
658 vf = (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1) * dz;
659 else
660 vf = LOW_DZ_MIN_1 * dz;
661 else if (0 <= dz && dz < UP_LIM_1)
662 vf = (UP_DZ_1 * dz);
663 else if (UP_LIM_1 <= dz && dz < UP_LIM_2)
664 vf = (dz - UP_LIM_1) * (UP_DZ_2 - UP_SPL_2 * spl) +
665 UP_DZ_1 * UP_LIM_1;
666 tmps += vf;
667 ltg[k].x = non_lin_add (ltg[k].x, tmps, MASK_ADD);
668 }
669 t = power[t].next;
670 } /* while */
671
672 t = *noise; /* calculate individual masking threshold */
673 while (t != LAST) { /* for non-tonal components, t, to find LTG */
674 z = ltg[power[t].map].bark; /* critical band rate of masker */
675 dz = ltg[k].bark - z; /* distance of bark value */
676 spl = power[t].x; /* sound pressure level of masker */
677
678 if (dz >= LOW_LIM_2 && dz < UP_LIM_2) {
679 tmps = spl + AV_NOISE_K + AV_NOISE_DZ * z;
680
681 /* masking function for lower & upper slopes */
682 if (LOW_LIM_2 <= dz && dz < LOW_LIM_1)
683 if (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1 > LOW_DZ_MIN_1)
684 vf = LOW_DZ_2 * (dz - LOW_LIM_1) +
685 (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1) * LOW_LIM_1;
686 else
687 vf = LOW_DZ_2 * (dz - LOW_LIM_1) + LOW_DZ_MIN_1 * LOW_LIM_1;
688 else if (LOW_LIM_1 <= dz && dz < 0)
689 if (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1 > LOW_DZ_MIN_1)
690 vf = (LOW_DZ_SPL_1 * spl + LOW_DZ_K_1) * dz;
691 else
692 vf = LOW_DZ_MIN_1 * dz;
693 else if (0 <= dz && dz < UP_LIM_1)
694 vf = (UP_DZ_1 * dz);
695 else if (UP_LIM_1 <= dz && dz < UP_LIM_2)
696 vf = (dz - UP_LIM_1) * (UP_DZ_2 - UP_SPL_2 * spl) +
697 UP_DZ_1 * UP_LIM_1;
698 tmps += vf;
699 ltg[k].x = non_lin_add (ltg[k].x, tmps, MASK_ADD);
700 }
701 t = power[t].next;
702 }
703
704 if (bit_rate < H_THR_OS_BR)
705 ltg[k].x = non_lin_add (ltg[k].hear, ltg[k].x, QUIET_ADD);
706 else
707 ltg[k].x = non_lin_add (ltg[k].hear + H_THR_OFFSET, ltg[k].x, QUIET_ADD);
708
709 } /* for */
710 fflush (stderr);
711 }
712
713
714
715 /* --------------------------------------------------------------
716 non_lin_add2
717 A flexible addition function for levels.
718 Input: a,b: the levels to be added.
719 c: the number of dB increase when a and b are equal.
720 Common values for c are 3.01 (power addition)
721 and 6.02 (voltage addition).
722 10.0/(10*log10(2)) = 3.3219
723 Function added 1993-04-14 Soren H. Nielsen
724 -------------------------------------------------------------- */
non_lin_add(double a,double b,double c)725 double non_lin_add (double a, double b, double c)
726 {
727 //double theirs, mine;
728 c *= 3.3219;
729
730 //mine = c/10.0 * add_db(10.0*a/c,10.0*b/c);
731
732 a = pow (10.0, a / c);
733 b = pow (10.0, b / c);
734 //theirs = c * log10(a+b);
735 //fprintf(stdout,"[%f %f]\n",theirs, mine);
736
737 return (c * log10 (a + b));
738 }
739
740
741 /****************************************************************
742 *
743 * This function finds the minimum masking threshold and
744 * return the value to the encoder.
745 *
746 ****************************************************************/
747
II_minimum_mask(int sub_size,g_thres * ltg,double * ltmin,int sblimit)748 void II_minimum_mask (int sub_size, g_thres * ltg, double *ltmin, int sblimit)
749 {
750 double min;
751 int i, j;
752
753
754 j = 1;
755 for (i = 0; i < sblimit; i++)
756 if (j >= sub_size - 1) /* check subband limit, and */
757 ltmin[i] = ltg[sub_size - 1].hear; /* calculate the minimum masking */
758 else { /* level of LTMIN for each subband */
759 min = ltg[j].x;
760 while (ltg[j].line >> 4 == i && j < sub_size) {
761 if (min > ltg[j].x)
762 min = ltg[j].x;
763 j++;
764 }
765 ltmin[i] = min;
766 }
767 }
768
769 /*****************************************************************
770 *
771 * This procedure is called in musicin to pick out the
772 * smaller of the scalefactor or threshold.
773 *
774 *****************************************************************/
775
II_smr(double * ltmin,double * smr,double * spike,double * scale,int sblimit,int l,int m)776 void II_smr (double *ltmin, double *smr, double *spike, double *scale,
777 int sblimit, int l, int m)
778 {
779 int i;
780 double max;
781
782 for (i = l; i < m; i++) {
783 /* determine the signal */
784 max = 20 * log10 (scale[i] * 32768) - 10; /* level for each subband */
785 if (spike[i] > max)
786 max = spike[i]; /* for the maximum scale */
787 max -= ltmin[i]; /* factors */
788 smr[i] = max;
789 }
790 }
791
792 /****************************************************************
793 *
794 * This procedure calls all the necessary functions to
795 * complete the psychoacoustic analysis.
796 *
797 ****************************************************************/
798
psycho_1(double (* buffer)[1152],double (* scale)[32],double (* ltmin)[32],frame_params * fr_ps,double (* smr)[32],double (* spiki)[32],int aiff)799 void psycho_1 (double (*buffer)[1152],
800 double (*scale)[32],
801 double (*ltmin)[32],
802 frame_params * fr_ps,
803 double (*smr)[32], double (*spiki)[32], int aiff)
804 {
805 layer *info = fr_ps->header;
806 int stereo = fr_ps->stereo;
807 int stereomc = fr_ps->stereomc;
808 int stereoaug = fr_ps->stereoaug;
809 int sblimit = fr_ps->sblimit;
810 int k, i, tone = 0, noise = 0;
811 static char init = 0;
812 static int off[12] = { 256, 256, 256, 256, 256, 256, 256, 256, 256, 256, 256, 256 }; /* max 7 MC channels + 5 compatible channels */
813 double *sample;
814 DSBL *spike;
815 static D1408 *fft_buf;
816 static mask_ptr power;
817 static g_ptr ltg;
818 int j, l, z, q;
819 static int crit_band;
820 static int *cbound;
821 static int sub_size;
822
823 sample = (double *) mem_alloc (sizeof (DFFT), "sample");
824 spike = (DSBL *) mem_alloc (sizeof (D12SBL), "spike");
825
826 if (!init) {
827 psycho_1_init_add_db();
828 /* bands, bark values, and mapping */
829 /* changed 5 to 7 for matricing 8/10/93,SR */
830 /* changed 7 to 12 for aug matricing 8/11/96,FdB */
831 fft_buf = (D1408 *) mem_alloc ((long) sizeof (D1408) * 12, "fft_buf");
832 power = (mask_ptr) mem_alloc (sizeof (mask) * HAN_SIZE, "power");
833 /* call functions for critical boundaries, freq. */
834 #ifdef NEWP1TABLES
835 cbound = psycho_1_read_cbound (info->lay, info->sampling_frequency, &crit_band);
836 psycho_1_read_freq_band (&sub_size, <g, info->lay, info->sampling_frequency);
837 #else
838 crit_band = read_crit_band (info->lay, info->sampling_frequency);
839 cbound = (int *) mem_alloc (sizeof (int) * crit_band, "cbound");
840 read_cbound (info->lay, info->sampling_frequency, crit_band, cbound);
841 read_freq_band (&sub_size, <g, info->lay, info->sampling_frequency);
842 #endif
843 make_map (sub_size, power, ltg);
844 for (i = 0; i < 1408; i++)
845 fft_buf[0][i] = fft_buf[1][i] = fft_buf[2][i] =
846 fft_buf[3][i] = fft_buf[4][i] = fft_buf[5][i] =
847 fft_buf[6][i] = fft_buf[7][i] = fft_buf[8][i] =
848 fft_buf[9][i] = fft_buf[10][i] = fft_buf[11][i] = 0;
849 init = 1;
850 }
851
852 if (aiff != 1) {
853 j = 0;
854 l = 2;
855 } else {
856 j = 0;
857 if (stereoaug == 2)
858 l = 12; /*fr_ps->stereo + fr_ps->stereomc + fr_ps->stereoaug + 5 compatible */
859 else
860 l = 7; /*fr_ps->stereo + fr_ps->stereomc + 2 compatible */
861 }
862
863 for (k = j; k < l; k++) {
864 for (i = 0; i < 1152; i++)
865 fft_buf[k][(i + off[k]) % 1408] = (double) buffer[k][i] / SCALE;
866 for (i = 0; i < FFT_SIZE; i++)
867 sample[i] = fft_buf[k][(i + 1216 + off[k]) % 1408];
868
869 off[k] += 1152;
870 off[k] %= 1408;
871 /* call functions for windowing PCM samples, */
872 II_hann_win (sample); /* location of spectral components in each */
873 for (i = 0; i < HAN_SIZE; i++)
874 power[i].x = DBMIN; /* subband with labeling */
875 II_f_f_t (sample, power); /* locate remaining non- */
876
877 if (fr_ps->header->center == 3 && k == 2) {
878 /* set center to 0, 9/2/93,SR */
879 /* add to Left and Right ? WtK */
880 for (z = 184; z < HAN_SIZE; z++)
881 power[z].x = -103.670; /* DBMIN + 96.330; */
882 }
883 II_pick_max (power, spike[k]); /* tonal sinusoidals, */
884
885 #ifdef PRINTOUT
886 if (verbosity >= 3) {
887 printf ("\nChannel %d", k);
888 printf ("\nSignal value per subband, from the FFT:\n");
889 for (i = 0; i < sblimit; i++)
890 printf ("%5.1f dB ", spike[k][i]);
891 printf ("\nMax. signal peak per subband, SCF SPL:\n");
892 for (i = 0; i < sblimit; i++) /* from [II_smr] determine the SCF SPL */
893 printf ("%5.1f dB ", 20 * log10 (scale[k][i] * 32768));
894 fflush (stdout);
895 }
896 #endif
897
898 II_tonal_label (power, &tone); /* reduce noise & tonal components , find */
899 noise_label (crit_band, cbound, power, &noise, ltg);
900
901 #ifdef PRINTOUT
902 if (verbosity >= 3) {
903 printf ("\nMaskers before sorting, FFT based levels:\n");
904 for (i = 0; i < 511; i++) {
905 if ((power[i].type == NOISE) && (power[i].x > -200))
906 printf ("N:%3u %5.1f dB ", i, power[i].x);
907 if ((power[i].type == TONE) && (power[i].x > -200))
908 printf ("T:%3u %5.1f dB ", i, power[i].x);
909 }
910 printf ("tone = %d noise = %d \n", tone, noise);
911 fflush (stdout);
912 }
913 #endif
914
915 subsampling (power, ltg, &tone, &noise); /* global & minimal */
916
917 #ifdef PRINTOUT
918 if (verbosity >= 3) {
919 printf ("\nMaskers after sorting:\n");
920 for (i = 0; i < 511; i++) {
921 if ((power[i].type == NOISE) && (power[i].x > -200))
922 printf ("N:%3u %5.1f dB ", i, power[i].x);
923 if ((power[i].type == TONE) && (power[i].x > -200))
924 printf ("T:%3u %5.1f dB ", i, power[i].x);
925 }
926 fflush (stdout);
927 }
928 #endif
929
930 threshold (sub_size, power, ltg, &tone, &noise, bitrate[info->lay - 1][info->bitrate_index] / (stereo + stereomc + stereoaug)); /*to-mask ratio *//* 21/03/1995 JMZ BUG ???!!!??? */
931 II_minimum_mask (sub_size, ltg, <min[k][0], sblimit);
932
933 #ifdef PRINTOUT
934 if (verbosity >= 3) {
935 printf ("\nMinimum masking threshold:\n");
936 for (i = 0; i < sblimit; i++)
937 printf ("%5.1f dB ", ltmin[k][i]);
938 fflush (stdout);
939 }
940 #endif
941
942 for (i = 0; i < SBLIMIT; i++)
943 spiki[k][i] = spike[k][i];
944
945 i = 0;
946 q = sblimit;
947 II_smr (ltmin[k], smr[k], spike[k], scale[k], sblimit, i, q);
948 }
949
950 mem_free ((void **) &sample);
951 mem_free ((void **) &spike);
952 } /*II_Psycho_One */
953
II_Psycho_One_ml(double (* buffer)[1152],double (* scale)[32],double (* ltmin)[32],frame_params * fr_ps,double (* smr)[32],double (* spiki)[32])954 void II_Psycho_One_ml (double (*buffer)[1152],
955 double (*scale)[32],
956 double (*ltmin)[32],
957 frame_params * fr_ps,
958 double (*smr)[32], double (*spiki)[32]
959 )
960 {
961 layer *info = fr_ps->header;
962 int n_ml_ch = info->multiling_ch;
963 int sblimit_ml = fr_ps->sblimit_ml;
964 int k, i, tone = 0, noise = 0;
965 static char init = 0;
966 static int off[7] = { 256, 256, 256, 256, 256, 256, 256 }; /* max 7 ML channels */
967 double *sample;
968 DSBL *spike;
969 static D1408 *fft_buf;
970 static mask_ptr power;
971 static g_ptr ltg_ml;
972 static int crit_band_ml;
973 static int *cbound_ml;
974 static int sub_size_ml;
975 int q;
976
977
978 sample = (double *) mem_alloc (sizeof (DFFT), "sample");
979 spike = (DSBL *) mem_alloc (sizeof (D7SBL), "spike");
980
981 if (!init) {
982 psycho_1_init_add_db(); /* MFC: This is already done in the normal psycho_1 call */
983 /* bands, bark values, and mapping */
984 fft_buf = (D1408 *) mem_alloc ((long) sizeof (D1408) * 7, "fft_buf");
985 power = (mask_ptr) mem_alloc (sizeof (mask) * HAN_SIZE, "power");
986 /* call functions for critical boundaries, freq. */
987 if (info->multiling_fs == 0) {
988 /* ML channels are the same frequency as the main channels */
989 #ifdef NEWP1TABLES
990 cbound_ml = psycho_1_read_cbound(info->lay, info->sampling_frequency, &crit_band_ml);
991 psycho_1_read_freq_band(&sub_size_ml, <g_ml, info->lay,info->sampling_frequency);
992 #else
993 crit_band_ml = read_crit_band (info->lay, info->sampling_frequency);
994 cbound_ml = (int *) mem_alloc (sizeof (int) * crit_band_ml, "cbound_ml");
995 read_cbound (info->lay, info->sampling_frequency, crit_band_ml,
996 cbound_ml);
997 read_freq_band (&sub_size_ml, <g_ml, info->lay,
998 info->sampling_frequency);
999 #endif
1000 /* values are equal to those for the mc audio data */
1001 } else {
1002 /* ML channels are half the frequency of the main channels - LSF */
1003 #ifdef NEWP1TABLES
1004 cbound_ml = psycho_1_read_cbound(info->lay, info->sampling_frequency+4, &crit_band_ml);
1005 psycho_1_read_freq_band(&sub_size_ml, <g_ml, info->lay,info->sampling_frequency+4);
1006 #else
1007 crit_band_ml = read_crit_band (info->lay, info->sampling_frequency + 4);
1008 cbound_ml = (int *) mem_alloc (sizeof (int) * crit_band_ml, "cbound_ml");
1009 read_cbound (info->lay, info->sampling_frequency + 4, crit_band_ml,
1010 cbound_ml);
1011 read_freq_band (&sub_size_ml, <g_ml, info->lay,
1012 info->sampling_frequency + 4);
1013 #endif
1014 }
1015 make_map (sub_size_ml, power, ltg_ml);
1016 for (i = 0; i < 1408; i++)
1017 fft_buf[0][i] = fft_buf[1][i] = fft_buf[2][i] =
1018 fft_buf[3][i] = fft_buf[4][i] = fft_buf[5][i] = fft_buf[6][i] = 0;
1019 init = 1;
1020 }
1021
1022 if (n_ml_ch > 0) {
1023 for (k = 0; k < n_ml_ch; k++) {
1024 for (i = 0; i < 1152; i++)
1025 fft_buf[k][(i + off[k]) % 1408] = (double) buffer[7 + k][i] / SCALE;
1026 for (i = 0; i < FFT_SIZE; i++)
1027 sample[i] = fft_buf[k][(i + 1216 + off[k]) % 1408];
1028
1029 off[k] += 1152;
1030 off[k] %= 1408;
1031
1032 /* call functions for windowing PCM samples, */
1033 II_hann_win (sample); /* location of spectral components in each */
1034 for (i = 0; i < HAN_SIZE; i++)
1035 power[i].x = DBMIN; /* subband with labeling */
1036 II_f_f_t (sample, power); /* locate remaining non- */
1037
1038 II_pick_max (power, spike[k]);
1039 /* tonal sinusoidals, */
1040
1041 II_tonal_label (power, &tone); /* reduce noise & tonal components , find */
1042 noise_label (crit_band_ml, cbound_ml, power, &noise, ltg_ml);
1043
1044 subsampling (power, ltg_ml, &tone, &noise); /* global & minimal */
1045
1046 threshold (sub_size_ml, power, ltg_ml, &tone, &noise, bitrate[info->lay - 1][info->bitrate_index] / 2); /* to-mask ratio */
1047 /* threshold, and sgnl- */
1048 /* fprintf(stderr, "sblimit_ml : %d\n", sblimit_ml); fflush(stderr); */
1049 II_minimum_mask (sub_size_ml, ltg_ml, <min[7 + k][0], sblimit_ml);
1050
1051 for (i = 0; i < SBLIMIT; i++)
1052 spiki[7 + k][i] = spike[k][i];
1053
1054 i = 0;
1055 q = sblimit_ml;
1056 II_smr (ltmin[7 + k], smr[7 + k], spike[k], scale[7 + k], sblimit_ml, i,
1057 q);
1058
1059 } /* k-loop */
1060 }
1061 /*n_ml_ch>0 */
1062 mem_free ((void **) &sample);
1063 mem_free ((void **) &spike);
1064 } /* II_psycho_One_ml */
1065
1066