1 /*
2 * Musepack audio compression
3 * Copyright (c) 2005-2009, The Musepack Development Team
4 * Copyright (C) 1999-2004 Buschmann/Klemm/Piecha/Wolf
5 *
6 * This library is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 2.1 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 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with this library; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 */
20
21 /*
22 * Prediction
23 * Short-Block-detection with smooth inset
24 * revise CalcMSThreshold
25 * /dev/audio for Windows too
26 * revise PNS/IS
27 * CVS with smoother inset
28 * several files per call
29 * revise ANS with changing SCFs
30
31 * No IS
32 * PNS estimation very rough, also IS should be used to reduce data rate in the side channel
33 * ANS problems at Frame boundaries when resolution changes
34 * ANS problems at Subframe boundaries when SCF changes
35 * CVS+ with smoother transition
36
37 ----------------------------------------
38
39 Optimize Tabelle[18] (use second table)
40 CVS+
41
42 - ANS is disregarded during the search for the best Res
43 - ANS messes up if res changes (each 36 samples) and/or SCF changes (each 12 samples)
44 - PNS not in difference signal
45
46 - implement IS in decoder
47 - Experimental Quantizer with complete energy preservation
48 - 1D, calculated
49 - 2D, calculated
50 - 2D, manually modified, coeffs set to 1.f
51
52 */
53
54 #include <string.h>
55 #include <stdio.h>
56
57 #include "libmpcpsy.h"
58 #include <mpc/datatypes.h>
59 #include <mpc/minimax.h>
60 #include <mpc/mpcmath.h>
61
62 // psy_tab.c
63 extern const float iw [PART_LONG]; // inverse partition-width for long
64 extern const float iw_short [PART_SHORT]; // inverse partition-width for short
65 extern const int wl [PART_LONG]; // w_low for long
66 extern const int wl_short [PART_SHORT]; // w_low for short
67 extern const int wh [PART_LONG]; // w_high for long
68 extern const int wh_short [PART_SHORT]; // w_high for short
69 extern float MinVal [PART_LONG]; // minimum quality that's adapted to the model, minval for long
70 extern float Loudness [PART_LONG]; // weighting factors for loudness calculation
71 extern float SPRD [PART_LONG] [PART_LONG]; // tabulated spreading function
72 extern float O_MAX;
73 extern float O_MIN;
74 extern float FAC1;
75 extern float FAC2; // constants for offset calculation
76 extern float partLtq [PART_LONG]; // threshold in quiet (partitions)
77 extern float invLtq [PART_LONG]; // inverse threshold in quiet (partitions, long)
78 extern float fftLtq [512]; // threshold in quiet (FFT)
79
80 // ans.c
81 extern float ANSspec_L [MAX_ANS_LINES];
82 extern float ANSspec_R [MAX_ANS_LINES]; // L/R-masking threshold for ANS
83 extern float ANSspec_M [MAX_ANS_LINES];
84 extern float ANSspec_S [MAX_ANS_LINES]; // M/S-masking threshold for ANS
85
86 void Init_Psychoakustiktabellen ( PsyModel* );
87 int CVD2048 ( PsyModel*, const float*, int* );
88
89 // Antialiasing for calculation of the subband power
90 const float Butfly [7] = { 0.5f, 0.2776f, 0.1176f, 0.0361f, 0.0075f, 0.000948f, 0.0000598f };
91
92 // Antialiasing for calculation of the masking thresholds
93 const float InvButfly [7] = { 2.f, 3.6023f, 8.5034f, 27.701f, 133.33f, 1054.852f, 16722.408f };
94
95 /* V A R I A B L E S */
96
97 static float a [PART_LONG];
98 static float b [PART_LONG];
99 static float c [PART_LONG];
100 static float d [PART_LONG]; // Integrations for tmpMask
101 static float Xsave_L [3 * 512];
102 static float Xsave_R [3 * 512]; // FFT-Amplitudes L/R
103 static float Ysave_L [3 * 512];
104 static float Ysave_R [3 * 512]; // FFT-Phases L/R
105 static float T_L [PART_LONG];
106 static float T_R [PART_LONG]; // time-constants for tmpMask
107 static float pre_erg_L[2][PART_SHORT];
108 static float pre_erg_R[2][PART_SHORT]; // Preecho-control short
109 static float PreThr_L [PART_LONG];
110 static float PreThr_R [PART_LONG]; // for Pre-Echo-control L/R
111 static float tmp_Mask_L [PART_LONG];
112 static float tmp_Mask_R [PART_LONG]; // for Post-Masking L/R
113 static int Vocal_L [MAX_CVD_LINE + 4];
114 static int Vocal_R [MAX_CVD_LINE + 4]; // FFT-Line belongs to harmonic?
115
116 /* F U N C T I O N S */
117
118 // fft_routines.c
119 void Init_FFT ( PsyModel* );
120 void PowSpec256 ( const float*, float* );
121 void PowSpec1024 ( const float*, float* );
122 void PowSpec2048 ( const float*, float* );
123 void PolarSpec1024 ( const float*, float*, float* );
124 void Cepstrum2048 ( float* cep, const int );
125
126 void Init_ANS ( void );
127
128 // Resets Arrays
Init_Psychoakustik(PsyModel * m)129 void Init_Psychoakustik ( PsyModel* m)
130 {
131 int i;
132
133 // initializing arrays with zero
134 memset ( Xsave_L, 0, sizeof Xsave_L );
135 memset ( Xsave_R, 0, sizeof Xsave_R );
136 memset ( Ysave_L, 0, sizeof Ysave_L );
137 memset ( Ysave_R, 0, sizeof Ysave_R );
138 memset ( a, 0, sizeof a );
139 memset ( b, 0, sizeof b );
140 memset ( c, 0, sizeof c );
141 memset ( d, 0, sizeof d );
142 memset ( T_L, 0, sizeof T_L );
143 memset ( T_R, 0, sizeof T_R );
144 memset ( Vocal_L, 0, sizeof Vocal_L );
145 memset ( Vocal_R, 0, sizeof Vocal_R );
146
147 m->SampleFreq = 0.;
148 m->BandWidth = 0.;
149 m->KBD1 = 2.;
150 m->KBD2 = -1.;
151 m->Ltq_offset = 0;
152 m->Ltq_max = 0;
153 m->EarModelFlag = 0;
154 m->PNS = 0.;
155 m->CombPenalities = -1;
156
157 // generate FFT lookup-tables with largest FFT-size of 1024
158 Init_FFT (m);
159
160 Init_ANS ();
161
162 Init_Psychoakustiktabellen (m);
163
164 // setting pre-echo variables to Ltq
165 for ( i = 0; i < PART_LONG; i++ ) {
166 pre_erg_L [0][i/3] = pre_erg_R [0][i/3] =
167 pre_erg_L [1][i/3] = pre_erg_R [1][i/3] =
168 tmp_Mask_L [i] = tmp_Mask_R [i] =
169 PreThr_L [i] = PreThr_R [i] = partLtq [i];
170 }
171
172 return;
173 }
174
175
176 // VBRmode 1: Adjustment of all SMRs via a factor (offset of SMRoffset dB)
177 // VBRmode 2: SMRs have a minimum of minSMR dB
178 static void
RaiseSMR_Signal(const int MaxBand,float * signal,float tmp)179 RaiseSMR_Signal ( const int MaxBand, float* signal, float tmp )
180 {
181 int Band;
182 float z = 0.;
183
184 for ( Band = MaxBand; Band >= 0; Band-- ) {
185 if ( z < signal [Band] ) z = signal [Band];
186 if ( z > tmp ) z = tmp;
187 if ( signal [Band] < z ) signal [Band] = z;
188 }
189 }
190
191
RaiseSMR(PsyModel * m,const int MaxBand,SMRTyp * smr)192 void RaiseSMR (PsyModel* m, const int MaxBand, SMRTyp* smr )
193 {
194 float tmp = POW10 ( 0.1 * m->minSMR );
195
196 RaiseSMR_Signal ( MaxBand, smr->L, tmp );
197 RaiseSMR_Signal ( MaxBand, smr->R, tmp );
198 RaiseSMR_Signal ( MaxBand, smr->M, tmp );
199 RaiseSMR_Signal ( MaxBand, smr->S, 0.5 * tmp );
200
201 return;
202 }
203
204 // input : *smr
205 // output: *smr, *ms, *x (only the entries for L/R contain relevant data)
206 // Check if either M/S- or L/R-coding has a lower perceptual entropy
207 // Choose the better mode, copy the appropriate data into the
208 // arrays that belong to L and R and set the ms-Flag accordingly.
209 void
MS_LR_Entscheidung(const int MaxBand,unsigned char * ms,SMRTyp * smr,SubbandFloatTyp * x)210 MS_LR_Entscheidung ( const int MaxBand, unsigned char* ms, SMRTyp* smr, SubbandFloatTyp* x )
211 {
212 int Band;
213 int n;
214 float PE_MS;
215 float PE_LR;
216 float tmpM;
217 float tmpS;
218 float* l;
219 float* r;
220
221
222 for ( Band = 0; Band <= MaxBand; Band++ ) { // calculate perceptual entropy
223 PE_LR = PE_MS = 1.f;
224 if (smr->L[Band] > 1.) PE_LR *= smr->L[Band];
225 if (smr->R[Band] > 1.) PE_LR *= smr->R[Band];
226 if (smr->M[Band] > 1.) PE_MS *= smr->M[Band];
227 if (smr->S[Band] > 1.) PE_MS *= smr->S[Band];
228
229 if ( PE_MS < PE_LR ) {
230 ms[Band] = 1;
231
232 // calculate M/S-signal and copies it to L/R-array
233 l = x[Band].L;
234 r = x[Band].R;
235 for ( n = 0; n < 36; n++, l++, r++ ) {
236 tmpM = (*l + *r) * 0.5f;
237 tmpS = (*l - *r) * 0.5f;
238 *l = tmpM;
239 *r = tmpS;
240 }
241
242 // copy M/S - SMR to L/R-fields
243 smr->L[Band] = smr->M[Band];
244 smr->R[Band] = smr->S[Band];
245 }
246 else {
247 ms[Band] = 0;
248 }
249 }
250
251 return;
252 }
253
254 // input : FFT-spectrums *spec0 und *spec1
255 // output: energy in the individual subbands *erg0 and *erg1
256 // With Butfly[], you can calculate the results of aliasing during calculation
257 // of subband energy from the FFT-spectrums.
258 static void
SubbandEnergy(const int MaxBand,float * erg0,float * erg1,const float * spec0,const float * spec1)259 SubbandEnergy ( const int MaxBand,
260 float* erg0,
261 float* erg1,
262 const float* spec0,
263 const float* spec1 )
264 {
265 int n;
266 int k;
267 int alias;
268 float tmp0;
269 float tmp1;
270
271
272 // Is this here correct for FFT-based data or is this calculation rule only for MDCTs???
273
274 for ( k = 0; k <= MaxBand; k++ ) { // subband index
275 tmp0 = tmp1 = 0.f;
276 for ( n = 0; n < 16; n++, spec0++, spec1++ ) { // spectral index
277 tmp0 += *spec0;
278 tmp1 += *spec1;
279
280 // Consideration of Aliasing between the subbands
281 if ( n < +sizeof(Butfly)/sizeof(*Butfly) && k != 0 ) {
282 alias = -1 - (n<<1);
283 tmp0 += Butfly [n] * (spec0[alias] - *spec0);
284 tmp1 += Butfly [n] * (spec1[alias] - *spec1);
285 }
286 else if ( n > 15-sizeof(Butfly)/sizeof(*Butfly) && k != 31 ) {
287 alias = 31 - (n<<1);
288 tmp0 += Butfly [15-n] * (spec0[alias] - *spec0);
289 tmp1 += Butfly [15-n] * (spec1[alias] - *spec1);
290 }
291 }
292 *erg0++ = tmp0;
293 *erg1++ = tmp1;
294 }
295
296 return;
297 }
298
299 // input : FFT-Spectrums *spec0 and *spec1
300 // output: energy in the individual partitions *erg0 and *erg1
301 static void
PartitionEnergy(float * erg0,float * erg1,const float * spec0,const float * spec1)302 PartitionEnergy ( float* erg0,
303 float* erg1,
304 const float* spec0,
305 const float* spec1 )
306 {
307 unsigned int n;
308 unsigned int k;
309 float e0;
310 float e1;
311
312 n = 0;
313
314 for ( ; n < 23; n++ ) { // 11 or 23
315 k = wh[n] - wl[n];
316 e0 = *spec0++;
317 e1 = *spec1++;
318 while ( k-- ) {
319 e0 += *spec0++;
320 e1 += *spec1++;
321 }
322 *erg0++ = e0;
323 *erg1++ = e1;
324 }
325
326 for ( ; n < 48; n++ ) { // 37 ... 46, 48, 57
327 k = wh[n] - wl[n];
328 e0 = sqrt (*spec0++);
329 e1 = sqrt (*spec1++);
330 while ( k-- ) {
331 e0 += sqrt (*spec0++);
332 e1 += sqrt (*spec1++);
333 }
334 *erg0++ = e0*e0 * iw[n];
335 *erg1++ = e1*e1 * iw[n];
336 }
337
338 for ( ; n < PART_LONG; n++ ) {
339 k = wh[n] - wl[n];
340 e0 = *spec0++;
341 e1 = *spec1++;
342 while ( k-- ) {
343 e0 += *spec0++;
344 e1 += *spec1++;
345 }
346 *erg0++ = e0;
347 *erg1++ = e1;
348 }
349 }
350
351
352 // input : FFT-Spectrums *spec0, *spec1 and unpredictability *cw0 and *cw1
353 // output: weighted energy in the individual partitions *erg0, *erg1
354 static void
WeightedPartitionEnergy(float * erg0,float * erg1,const float * spec0,const float * spec1,const float * cw0,const float * cw1)355 WeightedPartitionEnergy ( float* erg0,
356 float* erg1,
357 const float* spec0,
358 const float* spec1,
359 const float* cw0,
360 const float* cw1 )
361 {
362 unsigned int n;
363 unsigned int k;
364 float e0;
365 float e1;
366
367 n = 0;
368
369 for ( ; n < 23; n++ ) {
370 e0 = *spec0++ * *cw0++;
371 e1 = *spec1++ * *cw1++;
372 k = wh[n] - wl[n];
373 while ( k-- ) {
374 e0 += *spec0++ * *cw0++;
375 e1 += *spec1++ * *cw1++;
376 }
377 *erg0++ = e0;
378 *erg1++ = e1;
379 }
380
381 for ( ; n < 48; n++ ) {
382 e0 = sqrt (*spec0++ * *cw0++);
383 e1 = sqrt (*spec1++ * *cw1++);
384 k = wh[n] - wl[n];
385 while ( k-- ) {
386 e0 += sqrt (*spec0++ * *cw0++);
387 e1 += sqrt (*spec1++ * *cw1++);
388 }
389 *erg0++ = e0*e0 * iw[n];
390 *erg1++ = e1*e1 * iw[n];
391 }
392
393 for ( ; n < PART_LONG; n++ ) {
394 e0 = *spec0++ * *cw0++;
395 e1 = *spec1++ * *cw1++;
396 k = wh[n] - wl[n];
397 while ( k-- ) {
398 e0 += *spec0++ * *cw0++;
399 e1 += *spec1++ * *cw1++;
400 }
401 *erg0++ = e0;
402 *erg1++ = e1;
403 }
404 }
405
406 // input : masking thresholds, first half of the arrays *shaped0 and *shaped1
407 // output: masking thresholds, second half of the arrays *shaped0 and *shaped1
408 // Considering the result of aliasing via InvButfly[]
409 // The input *thr0, *thr1 is gathered via address calculation from *shaped0, *shaped1
410
411 static void
AdaptThresholds(const int MaxLine,float * shaped0,float * shaped1)412 AdaptThresholds ( const int MaxLine, float* shaped0, float* shaped1 )
413 {
414 int n;
415 int mod;
416 int alias;
417 float tmp;
418 const float* invb = InvButfly;
419 const float* thr0 = shaped0 - 512;
420 const float* thr1 = shaped1 - 512;
421 float tmp0;
422 float tmp1;
423
424
425 // should be able to optimize it with coasting. [ 9 ] + n * [ 7 + 7 + 2 ] + [ 7 ]
426 // Schleife Schl Schl Ausr Schleife
427 for ( n = 0; n < MaxLine; n++, thr0++, thr1++ ) {
428 mod = n & 15; // n%16
429 tmp0 = *thr0;
430 tmp1 = *thr1;
431
432 if ( mod < +sizeof(InvButfly)/sizeof(*InvButfly) && n > 12 ) {
433 alias = -1 - (mod<<1);
434 tmp = thr0[alias] * invb[mod];
435 if ( tmp < tmp0 ) tmp0 = tmp;
436 tmp = thr1[alias] * invb[mod];
437 if ( tmp < tmp1 ) tmp1 = tmp;
438 }
439 else if ( mod > 15-sizeof(InvButfly)/sizeof(*InvButfly) && n < 499 ) {
440 alias = 31 - (mod<<1);
441 tmp = thr0[alias] * invb[15-mod];
442 if ( tmp < tmp0 ) tmp0 = tmp;
443 tmp = thr1[alias] * invb[15-mod];
444 if ( tmp < tmp1 ) tmp1 = tmp;
445 }
446 *shaped0++ = tmp0;
447 *shaped1++ = tmp1;
448 }
449
450 return;
451 }
452
453
454 // input : current spectrum in the form of power *spec and phase *phase,
455 // the last two earlier spectrums are at position
456 // 512 and 1024 of the corresponding Input-Arrays.
457 // Array *vocal, which can mark an FFT_Linie as harmonic
458 // output: current amplitude *amp and unpredictability *cw
459 static void
CalcUnpred(PsyModel * m,const int MaxLine,const float * spec,const float * phase,const int * vocal,float * amp0,float * phs0,float * cw)460 CalcUnpred (PsyModel* m,
461 const int MaxLine,
462 const float* spec,
463 const float* phase,
464 const int* vocal,
465 float* amp0,
466 float* phs0,
467 float* cw )
468 {
469 int n;
470 float amp;
471 float tmp;
472 #define amp1 ((amp0) + 512) // amp[ 512...1023] contains data of frame-1
473 #define amp2 ((amp0) + 1024) // amp[1024...1535] contains data of frame-2
474 #define phs1 ((phs0) + 512) // phs[ 512...1023] contains data of frame-1
475 #define phs2 ((phs0) + 1024) // phs[1024...1535] contains data of frame-2
476
477
478 for ( n = 0; n < MaxLine; n++ ) {
479 tmp = COSF ((phs0[n] = phase[n]) - 2*phs1[n] + phs2[n]); // copy phase to output-array, predict phase and calculate predictive error
480 amp0[n] = SQRTF (spec[n]); // calculate and set amplitude
481 amp = 2*amp1[n] - amp2[n]; // predict amplitude
482
483 // calculate unpredictability
484 cw[n] = SQRTF (spec[n] + amp * (amp - 2*amp0[n] * tmp)) / (amp0[n] + FABS(amp));
485 }
486
487 // postprocessing of harmonic FFT-lines (*cw is set to CVD_UNPRED)
488 if ( m->CVD_used && vocal != NULL ) {
489 for ( n = 0; n < MAX_CVD_LINE; n++, cw++, vocal++ )
490 if ( *vocal != 0 && *cw > CVD_UNPRED * 0.01 * *vocal )
491 *cw = CVD_UNPRED * 0.01 * *vocal;
492 }
493
494 return;
495 }
496 #undef amp1
497 #undef amp2
498 #undef phs1
499 #undef phs2
500
501
502 // input : Energy *erg, calibrated energy *werg
503 // output: spread energy *res, spread weighted energy *wres
504 // SPRD describes the spreading function as calculated in psy_tab.c
505 static void
SpreadingSignal(const float * erg,const float * werg,float * res,float * wres)506 SpreadingSignal ( const float* erg, const float* werg, float* res,
507 float* wres )
508 {
509 int n;
510 int k;
511 int start;
512 int stop;
513 const float* sprd;
514 float e;
515 float ew;
516
517
518 for (k=0; k<PART_LONG; ++k, ++erg, ++werg) { // Source (masking partition)
519 start = maxi(k-5, 0); // minimum affected partition
520 stop = mini(k+7, PART_LONG-1); // maximum affected partition
521 sprd = SPRD[k] + start; // load vector
522 e = *erg;
523 ew = *werg;
524
525 for (n=start; n<=stop; ++n, ++sprd) {
526 res [n] += *sprd * e; // spreading signal
527 wres[n] += *sprd * ew; // spreading weighted signal
528 }
529 }
530
531 return;
532 }
533
534 // input : spread weighted energy *werg, spread energy *erg
535 // output: masking threshold *erg after applying the tonality-offset
536 static void
ApplyTonalityOffset(float * erg0,float * erg1,const float * werg0,const float * werg1)537 ApplyTonalityOffset ( float* erg0, float* erg1, const float* werg0, const float* werg1 )
538 {
539 int n;
540 float Offset;
541 float quot;
542
543
544 // calculation of the masked threshold in the partition range
545 for ( n = 0; n < PART_LONG; n++ ) {
546 quot = *werg0++ / *erg0;
547 if (quot <= 0.05737540597f) Offset = O_MAX;
548 else if (quot < 0.5871011603f ) Offset = FAC1 * POW (quot, FAC2);
549 else Offset = O_MIN;
550 *erg0++ *= iw[n] * minf(MinVal[n], Offset);
551
552 quot = *werg1++ / *erg1;
553 if (quot <= 0.05737540597f) Offset = O_MAX;
554 else if (quot < 0.5871011603f ) Offset = FAC1 * POW (quot, FAC2);
555 else Offset = O_MIN;
556 *erg1++ *= iw[n] * minf(MinVal[n], Offset);
557 }
558
559 return;
560 }
561
562 // input: previous loudness *loud, energies *erg, threshold in quiet *adapted_ltq
563 // output: tracked loudness *loud, adapted threshold in quiet <Return value>
564 static float
AdaptLtq(PsyModel * m,const float * erg0,const float * erg1)565 AdaptLtq ( PsyModel* m, const float* erg0, const float* erg1 )
566 {
567 static float loud = 0.f;
568 float* weight = Loudness;
569 float sum = 0.f;
570 int n;
571
572 // calculate loudness
573 for ( n = 0; n < PART_LONG; n++ )
574 sum += (*erg0++ + *erg1++) * *weight++;
575
576 // Utilization of the time constants (fast drop of Ltq T=5, slow rise of Ltq T=20)
577 //loud = (sum < loud) ? (4 * sum + loud)*0.2f : (19 * loud + sum)*0.05f;
578 loud = 0.98 * loud + 0.02 * (0.5 * sum);
579
580 // calculate dynamic offset for threshold in quiet, 0...+20 dB, at 96 dB loudness, an offset of 20 dB is assumed
581 return 1.f + m->varLtq * loud * 5.023772e-08f;
582 }
583
584 // input : simultaneous masking threshold *frqthr,
585 // previous masking threshold *tmpthr,
586 // Integrations *a (short-time) and *b (long-time)
587 // output: tracked Integrations *a and *b, time constant *tau
588 static void
CalcTemporalThreshold(float * a,float * b,float * tau,float * frqthr,float * tmpthr)589 CalcTemporalThreshold ( float* a, float* b, float* tau, float* frqthr, float* tmpthr )
590 {
591 int n;
592 float tmp;
593
594
595 for ( n = 0; n < PART_LONG; n++ ) {
596 // following calculations relative to threshold in quiet
597 frqthr[n] *= invLtq[n];
598 tmpthr[n] *= invLtq[n];
599
600 // new post-masking 'tmp' via time constant tau, if old post-masking > Ltq (=1)
601 tmp = tmpthr[n] > 1.f ? POW ( tmpthr[n], tau[n] ) : 1.f;
602
603 // calculate time constant for post-masking in next frame,
604 // if new time constant has to be calculated (new tmpMask < frqMask)
605 a[n] += 0.5f * (frqthr[n] - a[n]); // short time integrator
606 b[n] += 0.15f * (frqthr[n] - b[n]); // long time integrator
607 if (tmp < frqthr[n])
608 tau[n] = a[n] <= b[n] ? 0.8f : 0.2f + b[n] / a[n] * 0.6f;
609
610 // use post-masking of (Re-Normalization)
611 tmpthr[n] = maxf (frqthr[n], tmp) * partLtq[n];
612 }
613
614 return;
615 }
616
617 // input : L/R-Masking thresholds in Partitions *thrL, *thrR
618 // L/R-Subband energies *ergL, *ergR
619 // M/S-Subband energies *ergM, *ergS
620 // output: M/S-Masking thresholds in Partitions *thrM, *thrS
621 static void
CalcMSThreshold(PsyModel * m,const float * const ergL,const float * const ergR,const float * const ergM,const float * const ergS,float * const thrL,float * const thrR,float * const thrM,float * const thrS)622 CalcMSThreshold ( PsyModel* m,
623 const float* const ergL,
624 const float* const ergR,
625 const float* const ergM,
626 const float* const ergS,
627 float* const thrL,
628 float* const thrR,
629 float* const thrM,
630 float* const thrS )
631 {
632 int n;
633 float norm;
634 float tmp;
635
636 // All hardcoded numbers here should be pulled from somewhere,
637 // the "4.", the -2 dB, the 0.0625 and the 0.9375, as well as all bands where this is done
638
639 for ( n = 0; n < PART_LONG; n++ ) {
640 // estimate M/S thresholds out of L/R thresholds and M/S and L/R energies
641 thrS[n] = thrM[n] = maxf (ergM[n], ergS[n]) / maxf (ergL[n], ergR[n]) * minf (thrL[n], thrR[n]);
642
643 switch ( m->MS_Channelmode ) { // preserve 'near-mid' signal components
644 case 3:
645 if ( n > 0 ) {
646 double ratioMS = ergM[n] > ergS[n] ? ergS[n] / ergM[n] : ergM[n] / ergS[n];
647 double ratioLR = ergL[n] > ergR[n] ? ergR[n] / ergL[n] : ergL[n] / ergR[n];
648 if ( ratioMS < ratioLR ) { // MS
649 if ( ergM[n] > ergS[n] )
650 thrS[n] = thrL[n] = thrR[n] = 1.e18f;
651 else
652 thrM[n] = thrL[n] = thrR[n] = 1.e18f;
653 }
654 else { // LR
655 if ( ergL[n] > ergR[n] )
656 thrR[n] = thrM[n] = thrS[n] = 1.e18f;
657 else
658 thrL[n] = thrM[n] = thrS[n] = 1.e18f;
659 }
660 }
661 break;
662 case 4:
663 if ( n > 0 ) {
664 double ratioMS = ergM[n] > ergS[n] ? ergS[n] / ergM[n] : ergM[n] / ergS[n];
665 double ratioLR = ergL[n] > ergR[n] ? ergR[n] / ergL[n] : ergL[n] / ergR[n];
666 if ( ratioMS < ratioLR ) { // MS
667 if ( ergM[n] > ergS[n] )
668 thrS[n] = 1.e18f;
669 else
670 thrM[n] = 1.e18f;
671 }
672 else { // LR
673 if ( ergL[n] > ergR[n] )
674 thrR[n] = 1.e18f;
675 else
676 thrL[n] = 1.e18f;
677 }
678 }
679 break;
680 case 5:
681 thrS[n] *= 2.; // +3 dB
682 break;
683 case 6:
684 break;
685 default:
686 fprintf ( stderr, "Unknown stereo mode\n");
687 case 10:
688 if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB
689 norm = 0.70794578f * iw[n]; // -1.5 dB * iwidth
690 if ( ergM[n] > ergS[n] ) {
691 tmp = ergS[n] * norm;
692 if ( thrS[n] > tmp )
693 thrS[n] = MS2SPAT1 * thrS[n] + (1.f-MS2SPAT1) * tmp; // raises masking threshold by up to 3 dB
694 } else if ( ergS[n] > ergM[n] ) {
695 tmp = ergM[n] * norm;
696 if ( thrM[n] > tmp )
697 thrM[n] = MS2SPAT1 * thrM[n] + (1.f-MS2SPAT1) * tmp;
698 }
699 }
700 break;
701 case 11:
702 if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB
703 norm = 0.63095734f * iw[n]; // -2.0 dB * iwidth
704 if ( ergM[n] > ergS[n] ) {
705 tmp = ergS[n] * norm;
706 if ( thrS[n] > tmp )
707 thrS[n] = MS2SPAT2 * thrS[n] + (1.f-MS2SPAT2) * tmp; // raises masking threshold by up to 6 dB
708 } else if ( ergS[n] > ergM[n] ) {
709 tmp = ergM[n] * norm;
710 if ( thrM[n] > tmp )
711 thrM[n] = MS2SPAT2 * thrM[n] + (1.f-MS2SPAT2) * tmp;
712 }
713 }
714 break;
715 case 12:
716 if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB
717 norm = 0.56234133f * iw[n]; // -2.5 dB * iwidth
718 if ( ergM[n] > ergS[n] ) {
719 tmp = ergS[n] * norm;
720 if ( thrS[n] > tmp )
721 thrS[n] = MS2SPAT3 * thrS[n] + (1.f-MS2SPAT3) * tmp; // raises masking threshold by up to 9 dB
722 } else if ( ergS[n] > ergM[n] ) {
723 tmp = ergM[n] * norm;
724 if ( thrM[n] > tmp )
725 thrM[n] = MS2SPAT3 * thrM[n] + (1.f-MS2SPAT3) * tmp;
726 }
727 }
728 break;
729 case 13:
730 if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB
731 norm = 0.50118723f * iw[n]; // -3.0 dB * iwidth
732 if ( ergM[n] > ergS[n] ) {
733 tmp = ergS[n] * norm;
734 if ( thrS[n] > tmp )
735 thrS[n] = MS2SPAT4 * thrS[n] + (1.f-MS2SPAT4) * tmp; // raises masking threshold by up to 12 dB
736 } else if ( ergS[n] > ergM[n] ) {
737 tmp = ergM[n] * norm;
738 if ( thrM[n] > tmp )
739 thrM[n] = MS2SPAT4 * thrM[n] + (1.f-MS2SPAT4) * tmp;
740 }
741 }
742 break;
743 case 15:
744 if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB
745 norm = 0.50118723f * iw[n]; // -3.0 dB * iwidth
746 if ( ergM[n] > ergS[n] ) {
747 tmp = ergS[n] * norm;
748 if ( thrS[n] > tmp )
749 thrS[n] = tmp; // raises masking threshold by up to +oo dB an
750 } else if ( ergS[n] > ergM[n] ) {
751 tmp = ergM[n] * norm;
752 if ( thrM[n] > tmp )
753 thrM[n] = tmp;
754 }
755 }
756 break;
757 case 22:
758 if ( 4. * ergL[n] > ergR[n] && ergL[n] < 4. * ergR[n] ) {// Energy between both channels differs by less than 6 dB
759 norm = 0.56234133f * iw[n]; // -2.5 dB * iwidth
760 if ( ergM[n] > ergS[n] ) {
761 tmp = ergS[n] * norm;
762 if ( thrS[n] > tmp )
763 thrS[n] = maxf (tmp, ergM[n]*iw[n]*0.025); // +/- 1.414
764 } else if ( ergS[n] > ergM[n] ) {
765 tmp = ergM[n] * norm;
766 if ( thrM[n] > tmp )
767 thrM[n] = maxf (tmp, ergS[n]*iw[n]*0.025); // +/- 1.414
768 }
769 }
770 break;
771 }
772 }
773
774 return;
775 }
776
777 // input : Masking thresholds in Partitions *partThr0, *partThr1
778 // level of threshold in quiet *ltq in FFT-resolution
779 // output: Masking thresholds in FFT-resolution *thr0, *thr1
780 // inline, because it's called 4x
781 static void
ApplyLtq(float * thr0,float * thr1,const float * partThr0,const float * partThr1,const float AdaptedLTQ,int MSflag)782 ApplyLtq ( float* thr0,
783 float* thr1,
784 const float* partThr0,
785 const float* partThr1,
786 const float AdaptedLTQ,
787 int MSflag )
788 {
789 int n, k;
790 float ms, ltq, tmp, tmpThr0, tmpThr1;
791
792 ms = AdaptedLTQ * (MSflag ? 0.125f : 0.25f);
793 for( n = 0; n < PART_LONG; n++ )
794 {
795 tmpThr0 = sqrt(partThr0[n]);
796 tmpThr1 = sqrt(partThr1[n]);
797 for ( k = wl[n]; k <= wh[n]; k++, thr0++, thr1++ )
798 {
799 // threshold in quiet (Partition)
800 // Applies a much more gentle ATH rolloff + 6 dB more dynamic
801 ltq = sqrt (ms * fftLtq [k]);
802 tmp = tmpThr0 + ltq;
803 *thr0 = tmp * tmp;
804 tmp = tmpThr1 + ltq;
805 *thr1 = tmp * tmp;
806 }
807 }
808 }
809
810 // input : Subband energies *erg0, *erg1
811 // Masking thresholds in FFT-resolution *thr0, *thr1
812 // output: SMR per Subband *smr0, *smr1
813 static void
CalculateSMR(const int MaxBand,const float * erg0,const float * erg1,const float * thr0,const float * thr1,float * smr0,float * smr1)814 CalculateSMR ( const int MaxBand,
815 const float* erg0,
816 const float* erg1,
817 const float* thr0,
818 const float* thr1,
819 float* smr0,
820 float* smr1 )
821 {
822 int n;
823 int k;
824 float tmp0;
825 float tmp1;
826
827 // calculation of the masked thresholds in the subbands
828 for (n = 0; n <= MaxBand; n++ ) {
829 tmp0 = *thr0++;
830 tmp1 = *thr1++;
831 for (k=1; k<16; ++k, ++thr0, ++thr1) {
832 if (*thr0 < tmp0) tmp0 = *thr0;
833 if (*thr1 < tmp1) tmp1 = *thr1;
834 }
835 *smr0++ = 0.0625f * *erg0++ / tmp0;
836 *smr1++ = 0.0625f * *erg1++ / tmp1;
837 }
838
839 return;
840 }
841
842 // input : energy spectrums erg[4][128] (4 delayed FFTs)
843 // Energy of the last short block *preerg in short partitions
844 // PreechoFac declares allowed traved of the masking threshold
845 // output: masking threshold *thr in short partitions
846 // Energy of the last short block *preerg in short partitions
847 static void
CalcShortThreshold(PsyModel * m,const float erg[4][128],const float ShortThr,float * thr,float old_erg[2][PART_SHORT],int * transient)848 CalcShortThreshold ( PsyModel* m,
849 const float erg [4] [128],
850 const float ShortThr,
851 float* thr,
852 float old_erg [2][PART_SHORT],
853 int* transient )
854 {
855 const int* index_lo = wl_short; // lower FFT-index
856 const int* index_hi = wh_short; // upper FFT-index
857 const float* iwidth = iw_short; // inverse partition-width
858 int k;
859 int n;
860 int l;
861 float new_erg;
862 float th, TransDetect = m->TransDetect;
863 const float* ep;
864
865 for ( k = 0; k < PART_SHORT; k++ ) {
866 transient [k] = 0;
867 th = old_erg [0][k];
868 for ( n = 0; n < 4; n++ ) {
869 ep = erg[n] + index_lo [k];
870 l = index_hi [k] - index_lo [k];
871
872 new_erg = *ep++;
873 while (l--)
874 new_erg += *ep++; // e = Short_Partition-energy in piece n
875
876 if ( new_erg > old_erg [0][k] ) { // bigger than the old?
877
878 if ( new_erg > old_erg [0][k] * TransDetect ||
879 new_erg > old_erg [1][k] * TransDetect*2 ) // is signal transient?
880 transient [k] = 1;
881 }
882 else {
883 th = minf ( th, new_erg ); // assume short threshold = engr*PreechoFac
884 }
885
886 old_erg [1][k] = old_erg [0][k];
887 old_erg [0][k] = new_erg; // save the current one
888 }
889 thr [k] = th * ShortThr * *iwidth++; // pull out and multiply only when transient[k]=1
890 }
891
892 return;
893 }
894
895 // input : previous simultaneous masking threshold *preThr,
896 // current simultaneous masking threshold *simThr
897 // output: update of *preThr for next call,
898 // current masking threshold *partThr
899 static void
PreechoControl(float * partThr0,float * preThr0,const float * simThr0,float * partThr1,float * preThr1,const float * simThr1)900 PreechoControl ( float* partThr0,
901 float* preThr0,
902 const float* simThr0,
903 float* partThr1,
904 float* preThr1,
905 const float* simThr1 )
906 {
907 int n;
908
909 for ( n = 0; n < PART_LONG; n++ ) {
910 *partThr0++ = minf ( *simThr0, *preThr0 * PREFAC_LONG);
911 *partThr1++ = minf ( *simThr1, *preThr1 * PREFAC_LONG);
912 *preThr0++ = *simThr0++;
913 *preThr1++ = *simThr1++;
914 }
915 return;
916 }
917
918
919 void
TransientenCalc(int * T,const int * TL,const int * TR)920 TransientenCalc ( int* T,
921 const int* TL,
922 const int* TR )
923 {
924 int i;
925 int x1;
926 int x2;
927
928 memset ( T, 0, 32*sizeof(*T) );
929
930 for ( i = 0; i < PART_SHORT; i++ )
931 if ( TL[i] || TR[i] ) {
932 x1 = wl_short[i] >> 2;
933 x2 = wh_short[i] >> 2;
934 while ( x1 <= x2 )
935 T [x1++] = 1;
936 }
937 }
938
939
940 // input : PCM-Data *data
941 // output: SMRs for the input data
942 SMRTyp
Psychoakustisches_Modell(PsyModel * m,const int MaxBand,const PCMDataTyp * data,int * TransientL,int * TransientR)943 Psychoakustisches_Modell ( PsyModel* m,
944 const int MaxBand,
945 const PCMDataTyp* data,
946 int* TransientL,
947 int* TransientR )
948 {
949 float Xi_L[32], Xi_R[32]; // acoustic pressure per Subband L/R
950 float Xi_M[32], Xi_S[32]; // acoustic pressure per Subband M/S
951 float cw_L[512], cw_R[512]; // unpredictability (only L/R)
952 float erg0[512], erg1[512]; // holds energy spectrum of long FFT
953 float phs0[512], phs1[512]; // holds phase spectrum of long FFT
954 float Thr_L[2*512], Thr_R[2*512]; // masking thresholds L/R, second half for triangle swap
955 float Thr_M[2*512], Thr_S[2*512]; // masking thresholds M/S, second half for triangle swap
956 float F_256[4][128]; // holds energies of short FFTs (L/R only)
957 float Xerg[1024]; // holds energy spectrum of very long FFT
958 float Ls_L[PART_LONG], Ls_R[PART_LONG]; // acoustic pressure in Partition L/R
959 float Ls_M[PART_LONG], Ls_S[PART_LONG]; // acoustic pressure per each partition M/S
960 float PartThr_L[PART_LONG], PartThr_R[PART_LONG]; // masking thresholds L/R (Partition)
961 float PartThr_M[PART_LONG], PartThr_S[PART_LONG]; // masking thresholds M/S (Partition)
962 float sim_Mask_L[PART_LONG], sim_Mask_R[PART_LONG]; // simultaneous masking (only L/R)
963 float clow_L[PART_LONG], clow_R[PART_LONG]; // spread, weighted energy (only L/R)
964 float cLs_L[PART_LONG], cLs_R[PART_LONG]; // weighted partition energy (only L/R)
965 float shortThr_L[PART_SHORT],shortThr_R[PART_SHORT]; // threshold for short FFT (only L/R)
966 int n;
967 int MaxLine = (MaxBand+1)*16; // set FFT-resolution according to MaxBand
968 SMRTyp SMR0;
969 SMRTyp SMR1; // holds SMR's for first and second Analysis
970 int isvoc_L = 0;
971 int isvoc_R = 0;
972 float factorLTQ = 1.f; // Offset after variable LTQ
973
974 // 'ClearVocalDetection'-Process
975 if ( m->CVD_used ) {
976 memset ( Vocal_L, 0, sizeof Vocal_L );
977 memset ( Vocal_R, 0, sizeof Vocal_R );
978
979 // left channel
980 PowSpec2048 ( &data->L[0], Xerg );
981 isvoc_L = CVD2048 ( m, Xerg, Vocal_L );
982 // right channel
983 PowSpec2048 ( &data->R[0], Xerg );
984 isvoc_R = CVD2048 ( m, Xerg, Vocal_R );
985 }
986
987 // calculation of the spectral energy via FFT
988 PolarSpec1024 ( &data->L[0], erg0, phs0 ); // left
989 PolarSpec1024 ( &data->R[0], erg1, phs1 ); // right
990
991 // calculation of the acoustic pressures per each subband for L/R-signals
992 SubbandEnergy ( MaxBand, Xi_L, Xi_R, erg0, erg1 );
993
994 // calculation of the acoustic pressures per each partition
995 PartitionEnergy ( Ls_L, Ls_R, erg0, erg1 );
996
997 // calculate the predictability of the signal
998 // left
999 memmove ( Xsave_L+512, Xsave_L, 1024*sizeof(float) );
1000 memmove ( Ysave_L+512, Ysave_L, 1024*sizeof(float) );
1001 CalcUnpred ( m, MaxLine, erg0, phs0, isvoc_L ? Vocal_L : NULL, Xsave_L, Ysave_L, cw_L );
1002 // right
1003 memmove ( Xsave_R+512, Xsave_R, 1024*sizeof(float) );
1004 memmove ( Ysave_R+512, Ysave_R, 1024*sizeof(float) );
1005 CalcUnpred ( m, MaxLine, erg1, phs1, isvoc_R ? Vocal_R : NULL, Xsave_R, Ysave_R, cw_R );
1006
1007 // calculation of the weighted acoustic pressures per each partition
1008 WeightedPartitionEnergy ( cLs_L, cLs_R, erg0, erg1, cw_L, cw_R );
1009
1010 // Spreading Signal & weighted unpredictability-signal
1011 // left
1012 memset ( clow_L , 0, sizeof clow_L );
1013 memset ( sim_Mask_L, 0, sizeof sim_Mask_L );
1014 SpreadingSignal ( Ls_L, cLs_L, sim_Mask_L, clow_L );
1015 // right
1016 memset ( clow_R , 0, sizeof clow_R );
1017 memset ( sim_Mask_R, 0, sizeof sim_Mask_R );
1018 SpreadingSignal ( Ls_R, cLs_R, sim_Mask_R, clow_R );
1019
1020 // Offset depending on tonality
1021 ApplyTonalityOffset ( sim_Mask_L, sim_Mask_R, clow_L, clow_R );
1022
1023 // handling of transient signals
1024 // calculate four short FFTs (left)
1025 PowSpec256 ( &data->L[ 0+SHORTFFT_OFFSET], F_256[0] );
1026 PowSpec256 ( &data->L[144+SHORTFFT_OFFSET], F_256[1] );
1027 PowSpec256 ( &data->L[288+SHORTFFT_OFFSET], F_256[2] );
1028 PowSpec256 ( &data->L[432+SHORTFFT_OFFSET], F_256[3] );
1029 // calculate short Threshold
1030 CalcShortThreshold ( m, F_256, m->ShortThr, shortThr_L, pre_erg_L, TransientL );
1031
1032 // calculate four short FFTs (right)
1033 PowSpec256 ( &data->R[ 0+SHORTFFT_OFFSET], F_256[0] );
1034 PowSpec256 ( &data->R[144+SHORTFFT_OFFSET], F_256[1] );
1035 PowSpec256 ( &data->R[288+SHORTFFT_OFFSET], F_256[2] );
1036 PowSpec256 ( &data->R[432+SHORTFFT_OFFSET], F_256[3] );
1037 // calculate short Threshold
1038 CalcShortThreshold ( m, F_256, m->ShortThr, shortThr_R, pre_erg_R, TransientR );
1039
1040 // dynamic adjustment of the threshold in quiet to the loudness of the current sequence
1041 if ( m->varLtq > 0. )
1042 factorLTQ = AdaptLtq (m, Ls_L, Ls_R );
1043
1044 // utilization of the temporal post-masking
1045 if ( m->tmpMask_used ) {
1046 CalcTemporalThreshold ( a, b, T_L, sim_Mask_L, tmp_Mask_L );
1047 CalcTemporalThreshold ( c, d, T_R, sim_Mask_R, tmp_Mask_R );
1048 memcpy ( sim_Mask_L, tmp_Mask_L, sizeof sim_Mask_L );
1049 memcpy ( sim_Mask_R, tmp_Mask_R, sizeof sim_Mask_R );
1050 }
1051
1052 // transient signal?
1053 for ( n = 0; n < PART_SHORT; n++ ) {
1054 if ( TransientL [n] ) {
1055 sim_Mask_L [3*n ] = minf ( sim_Mask_L [3*n ], shortThr_L [n] );
1056 sim_Mask_L [3*n+1] = minf ( sim_Mask_L [3*n+1], shortThr_L [n] );
1057 sim_Mask_L [3*n+2] = minf ( sim_Mask_L [3*n+2], shortThr_L [n] );
1058 }
1059 if ( TransientR[n] ) {
1060 sim_Mask_R [3*n ] = minf ( sim_Mask_R [3*n ], shortThr_R [n] );
1061 sim_Mask_R [3*n+1] = minf ( sim_Mask_R [3*n+1], shortThr_R [n] );
1062 sim_Mask_R [3*n+2] = minf ( sim_Mask_R [3*n+2], shortThr_R [n] );
1063 }
1064 }
1065
1066 // Pre-Echo control
1067 PreechoControl ( PartThr_L,PreThr_L, sim_Mask_L, PartThr_R, PreThr_R, sim_Mask_R );
1068
1069 // utilization of the threshold in quiet
1070 ApplyLtq ( Thr_L, Thr_R, PartThr_L, PartThr_R, factorLTQ, 0 );
1071
1072 // Consideration of aliasing between the subbands (noise is smeared)
1073 // In: Thr[0..511], Out: Thr[512...1023]
1074 AdaptThresholds ( MaxLine, Thr_L+512, Thr_R+512 );
1075 memmove ( Thr_L, Thr_L+512, 512*sizeof(float) );
1076 memmove ( Thr_R, Thr_R+512, 512*sizeof(float) );
1077
1078 // calculation of the Signal-to-Mask-Ratio
1079 CalculateSMR ( MaxBand, Xi_L, Xi_R, Thr_L, Thr_R, SMR0.L, SMR0.R );
1080
1081 /***************************************************************************************/
1082 /***************************************************************************************/
1083 if ( m->MS_Channelmode > 0 ) {
1084 // calculation of the spectral energy via FFT
1085 PowSpec1024 ( &data->M[0], erg0 ); // mid
1086 PowSpec1024 ( &data->S[0], erg1 ); // side
1087
1088 // calculation of the acoustic pressures per each subband for M/S-signals
1089 SubbandEnergy ( MaxBand, Xi_M, Xi_S, erg0, erg1 );
1090
1091 // calculation of the acoustic pressures per each partition
1092 PartitionEnergy ( Ls_M, Ls_S, erg0, erg1 );
1093
1094 // calculate masking thresholds for M/S
1095 CalcMSThreshold ( m, Ls_L, Ls_R, Ls_M, Ls_S, PartThr_L, PartThr_R, PartThr_M, PartThr_S );
1096 ApplyLtq ( Thr_M, Thr_S, PartThr_M, PartThr_S, factorLTQ, 1 );
1097
1098 // Consideration of aliasing between the subbands (noise is smeared)
1099 // In: Thr[0..511], Out: Thr[512...1023]
1100 AdaptThresholds ( MaxLine, Thr_M+512, Thr_S+512 );
1101 memmove ( Thr_M, Thr_M+512, 512*sizeof(float) );
1102 memmove ( Thr_S, Thr_S+512, 512*sizeof(float) );
1103
1104 // calculation of the Signal-to-Mask-Ratio
1105 CalculateSMR ( MaxBand, Xi_M, Xi_S, Thr_M, Thr_S, SMR0.M, SMR0.S );
1106 }
1107
1108 if ( m->NS_Order > 0 ) { // providing the Noise Shaping thresholds
1109 memcpy ( ANSspec_L, Thr_L, sizeof ANSspec_L );
1110 memcpy ( ANSspec_R, Thr_R, sizeof ANSspec_R );
1111 memcpy ( ANSspec_M, Thr_M, sizeof ANSspec_M );
1112 memcpy ( ANSspec_S, Thr_S, sizeof ANSspec_S );
1113 }
1114 /***************************************************************************************/
1115 /***************************************************************************************/
1116
1117 //
1118 //-------- second model calculation via shifted FFT ------------------------
1119 //
1120 // calculation of the spectral power via FFT
1121 PolarSpec1024 ( &data->L[576], erg0, phs0 ); // left
1122 PolarSpec1024 ( &data->R[576], erg1, phs1 ); // right
1123
1124 // calculation of the acoustic pressures per each subband for L/R-signals
1125 SubbandEnergy ( MaxBand, Xi_L, Xi_R, erg0, erg1 );
1126
1127 // calculation of the acoustic pressures per each partition
1128 PartitionEnergy ( Ls_L, Ls_R, erg0, erg1 );
1129
1130 // calculate the predictability of the signal
1131 // left
1132 memmove ( Xsave_L+512, Xsave_L, 1024*sizeof(float) );
1133 memmove ( Ysave_L+512, Ysave_L, 1024*sizeof(float) );
1134 CalcUnpred ( m, MaxLine, erg0, phs0, isvoc_L ? Vocal_L : NULL, Xsave_L, Ysave_L, cw_L );
1135 // right
1136 memmove ( Xsave_R+512, Xsave_R, 1024*sizeof(float) );
1137 memmove ( Ysave_R+512, Ysave_R, 1024*sizeof(float) );
1138 CalcUnpred ( m, MaxLine, erg1, phs1, isvoc_R ? Vocal_R : NULL, Xsave_R, Ysave_R, cw_R );
1139
1140 // calculation of the weighted acoustic pressure per each partition
1141 WeightedPartitionEnergy ( cLs_L, cLs_R, erg0, erg1, cw_L, cw_R );
1142
1143 // Spreading Signal & weighted unpredictability-signal
1144 // left
1145 memset ( clow_L , 0, sizeof clow_L );
1146 memset ( sim_Mask_L, 0, sizeof sim_Mask_L );
1147 SpreadingSignal ( Ls_L, cLs_L, sim_Mask_L, clow_L );
1148 // right
1149 memset ( clow_R , 0, sizeof clow_R );
1150 memset ( sim_Mask_R, 0, sizeof sim_Mask_R );
1151 SpreadingSignal ( Ls_R, cLs_R, sim_Mask_R, clow_R );
1152
1153 // Offset depending on tonality
1154 ApplyTonalityOffset ( sim_Mask_L, sim_Mask_R, clow_L, clow_R );
1155
1156 // Handling of transient signals
1157 // calculate four short FFTs (left)
1158 PowSpec256 ( &data->L[ 576+SHORTFFT_OFFSET], F_256[0] );
1159 PowSpec256 ( &data->L[ 720+SHORTFFT_OFFSET], F_256[1] );
1160 PowSpec256 ( &data->L[ 864+SHORTFFT_OFFSET], F_256[2] );
1161 PowSpec256 ( &data->L[1008+SHORTFFT_OFFSET], F_256[3] );
1162 // calculate short Threshold
1163 CalcShortThreshold ( m, F_256, m->ShortThr, shortThr_L, pre_erg_L, TransientL );
1164
1165 // calculate four short FFTs (right)
1166 PowSpec256 ( &data->R[ 576+SHORTFFT_OFFSET], F_256[0] );
1167 PowSpec256 ( &data->R[ 720+SHORTFFT_OFFSET], F_256[1] );
1168 PowSpec256 ( &data->R[ 864+SHORTFFT_OFFSET], F_256[2] );
1169 PowSpec256 ( &data->R[1008+SHORTFFT_OFFSET], F_256[3] );
1170 // calculate short Threshold
1171 CalcShortThreshold ( m, F_256, m->ShortThr, shortThr_R, pre_erg_R, TransientR );
1172
1173 // dynamic adjustment of threshold in quiet to loudness of the current sequence
1174 if ( m->varLtq > 0. )
1175 factorLTQ = AdaptLtq ( m, Ls_L, Ls_R );
1176
1177 // utilization of temporal post-masking
1178 if (m->tmpMask_used) {
1179 CalcTemporalThreshold ( a, b, T_L, sim_Mask_L, tmp_Mask_L );
1180 CalcTemporalThreshold ( c, d, T_R, sim_Mask_R, tmp_Mask_R );
1181 memcpy ( sim_Mask_L, tmp_Mask_L, sizeof sim_Mask_L );
1182 memcpy ( sim_Mask_R, tmp_Mask_R, sizeof sim_Mask_R );
1183 }
1184
1185 // transient signal?
1186 for ( n = 0; n < PART_SHORT; n++ ) {
1187 if ( TransientL[n] ) {
1188 sim_Mask_L [3*n ] = minf ( sim_Mask_L [3*n ], shortThr_L [n] );
1189 sim_Mask_L [3*n+1] = minf ( sim_Mask_L [3*n+1], shortThr_L [n] );
1190 sim_Mask_L [3*n+2] = minf ( sim_Mask_L [3*n+2], shortThr_L [n] );
1191 }
1192 if ( TransientR[n] ) {
1193 sim_Mask_R [3*n ] = minf ( sim_Mask_R [3*n ], shortThr_R [n] );
1194 sim_Mask_R [3*n+1] = minf ( sim_Mask_R [3*n+1], shortThr_R [n] );
1195 sim_Mask_R [3*n+2] = minf ( sim_Mask_R [3*n+2], shortThr_R [n] );
1196 }
1197 }
1198
1199 // Pre-Echo control
1200 PreechoControl ( PartThr_L, PreThr_L, sim_Mask_L, PartThr_R, PreThr_R, sim_Mask_R );
1201
1202 // utilization of threshold in quiet
1203 ApplyLtq ( Thr_L, Thr_R, PartThr_L, PartThr_R, factorLTQ, 0 );
1204
1205 // Consideration of aliasing between the subbands (noise is smeared)
1206 // In: Thr[0..511], Out: Thr[512...1023]
1207 AdaptThresholds ( MaxLine, Thr_L+512, Thr_R+512 );
1208 memmove ( Thr_L, Thr_L+512, 512*sizeof(float) );
1209 memmove ( Thr_R, Thr_R+512, 512*sizeof(float) );
1210
1211 // calculation of the Signal-to-Mask-Ratio
1212 CalculateSMR ( MaxBand, Xi_L, Xi_R, Thr_L, Thr_R, SMR1.L, SMR1.R );
1213
1214 /***************************************************************************************/
1215 /***************************************************************************************/
1216 if ( m->MS_Channelmode > 0 ) {
1217 // calculation of the spectral energy via FFT
1218 PowSpec1024 ( &data->M[576], erg0 ); // mid
1219 PowSpec1024 ( &data->S[576], erg1 ); // side
1220
1221 // calculation of the acoustic pressure per each subband for M/S-signals
1222 SubbandEnergy ( MaxBand, Xi_M, Xi_S, erg0, erg1 );
1223
1224 // calculation of the acoustic pressure per each partition
1225 PartitionEnergy ( Ls_M, Ls_S, erg0, erg1 );
1226
1227 // calculate masking thresholds for M/S
1228 CalcMSThreshold ( m, Ls_L, Ls_R, Ls_M, Ls_S, PartThr_L, PartThr_R, PartThr_M, PartThr_S );
1229 ApplyLtq ( Thr_M, Thr_S, PartThr_M, PartThr_S, factorLTQ, 1 );
1230
1231 // Consideration of aliasing between the subbands (noise is smeared)
1232 // In: Thr[0..511], Out: Thr[512...1023]
1233 AdaptThresholds ( MaxLine, Thr_M+512, Thr_S+512 );
1234 memmove ( Thr_M, Thr_M+512, 512*sizeof(float) );
1235 memmove ( Thr_S, Thr_S+512, 512*sizeof(float) );
1236
1237 // calculation of the Signal-to-Mask-Ratio
1238 CalculateSMR ( MaxBand, Xi_M, Xi_S, Thr_M, Thr_S, SMR1.M, SMR1.S );
1239 }
1240 /***************************************************************************************/
1241 /***************************************************************************************/
1242
1243 if ( m->NS_Order > 0 ) {
1244 for ( n = 0; n < MAX_ANS_LINES; n++ ) { // providing Noise Shaping thresholds
1245 ANSspec_L [n] = minf ( ANSspec_L [n], Thr_L [n] );
1246 ANSspec_R [n] = minf ( ANSspec_R [n], Thr_R [n] );
1247 ANSspec_M [n] = minf ( ANSspec_M [n], Thr_M [n] );
1248 ANSspec_S [n] = minf ( ANSspec_S [n], Thr_S [n] );
1249 }
1250 }
1251
1252 for ( n = 0; n <= MaxBand; n++ ) { // choose 'worst case'-SMR from shifted analysis windows
1253 SMR0.L[n] = maxf ( SMR0.L[n], SMR1.L[n] );
1254 SMR0.R[n] = maxf ( SMR0.R[n], SMR1.R[n] );
1255 SMR0.M[n] = maxf ( SMR0.M[n], SMR1.M[n] );
1256 SMR0.S[n] = maxf ( SMR0.S[n], SMR1.S[n] );
1257 }
1258 return SMR0;
1259 }
1260