1 /* -----------------------------------------------------------------------------
2 Software License for The Fraunhofer FDK AAC Codec Library for Android
3 
4 © Copyright  1995 - 2019 Fraunhofer-Gesellschaft zur Förderung der angewandten
5 Forschung e.V. All rights reserved.
6 
7  1.    INTRODUCTION
8 The Fraunhofer FDK AAC Codec Library for Android ("FDK AAC Codec") is software
9 that implements the MPEG Advanced Audio Coding ("AAC") encoding and decoding
10 scheme for digital audio. This FDK AAC Codec software is intended to be used on
11 a wide variety of Android devices.
12 
13 AAC's HE-AAC and HE-AAC v2 versions are regarded as today's most efficient
14 general perceptual audio codecs. AAC-ELD is considered the best-performing
15 full-bandwidth communications codec by independent studies and is widely
16 deployed. AAC has been standardized by ISO and IEC as part of the MPEG
17 specifications.
18 
19 Patent licenses for necessary patent claims for the FDK AAC Codec (including
20 those of Fraunhofer) may be obtained through Via Licensing
21 (www.vialicensing.com) or through the respective patent owners individually for
22 the purpose of encoding or decoding bit streams in products that are compliant
23 with the ISO/IEC MPEG audio standards. Please note that most manufacturers of
24 Android devices already license these patent claims through Via Licensing or
25 directly from the patent owners, and therefore FDK AAC Codec software may
26 already be covered under those patent licenses when it is used for those
27 licensed purposes only.
28 
29 Commercially-licensed AAC software libraries, including floating-point versions
30 with enhanced sound quality, are also available from Fraunhofer. Users are
31 encouraged to check the Fraunhofer website for additional applications
32 information and documentation.
33 
34 2.    COPYRIGHT LICENSE
35 
36 Redistribution and use in source and binary forms, with or without modification,
37 are permitted without payment of copyright license fees provided that you
38 satisfy the following conditions:
39 
40 You must retain the complete text of this software license in redistributions of
41 the FDK AAC Codec or your modifications thereto in source code form.
42 
43 You must retain the complete text of this software license in the documentation
44 and/or other materials provided with redistributions of the FDK AAC Codec or
45 your modifications thereto in binary form. You must make available free of
46 charge copies of the complete source code of the FDK AAC Codec and your
47 modifications thereto to recipients of copies in binary form.
48 
49 The name of Fraunhofer may not be used to endorse or promote products derived
50 from this library without prior written permission.
51 
52 You may not charge copyright license fees for anyone to use, copy or distribute
53 the FDK AAC Codec software or your modifications thereto.
54 
55 Your modified versions of the FDK AAC Codec must carry prominent notices stating
56 that you changed the software and the date of any change. For modified versions
57 of the FDK AAC Codec, the term "Fraunhofer FDK AAC Codec Library for Android"
58 must be replaced by the term "Third-Party Modified Version of the Fraunhofer FDK
59 AAC Codec Library for Android."
60 
61 3.    NO PATENT LICENSE
62 
63 NO EXPRESS OR IMPLIED LICENSES TO ANY PATENT CLAIMS, including without
64 limitation the patents of Fraunhofer, ARE GRANTED BY THIS SOFTWARE LICENSE.
65 Fraunhofer provides no warranty of patent non-infringement with respect to this
66 software.
67 
68 You may use this FDK AAC Codec software or modifications thereto only for
69 purposes that are authorized by appropriate patent licenses.
70 
71 4.    DISCLAIMER
72 
73 This FDK AAC Codec software is provided by Fraunhofer on behalf of the copyright
74 holders and contributors "AS IS" and WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES,
75 including but not limited to the implied warranties of merchantability and
76 fitness for a particular purpose. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
77 CONTRIBUTORS BE LIABLE for any direct, indirect, incidental, special, exemplary,
78 or consequential damages, including but not limited to procurement of substitute
79 goods or services; loss of use, data, or profits, or business interruption,
80 however caused and on any theory of liability, whether in contract, strict
81 liability, or tort (including negligence), arising in any way out of the use of
82 this software, even if advised of the possibility of such damage.
83 
84 5.    CONTACT INFORMATION
85 
86 Fraunhofer Institute for Integrated Circuits IIS
87 Attention: Audio and Multimedia Departments - FDK AAC LL
88 Am Wolfsmantel 33
89 91058 Erlangen, Germany
90 
91 www.iis.fraunhofer.de/amm
92 amm-info@iis.fraunhofer.de
93 ----------------------------------------------------------------------------- */
94 
95 /**************************** AAC encoder library ******************************
96 
97    Author(s):   A. Horndasch (code originally from lwr) / Josef Hoepfl (FDK)
98 
99    Description: intensity stereo processing
100 
101 *******************************************************************************/
102 
103 #include "intensity.h"
104 
105 #include "interface.h"
106 #include "psy_configuration.h"
107 #include "psy_const.h"
108 #include "qc_main.h"
109 #include "bit_cnt.h"
110 
111 /* only set an IS seed it left/right channel correlation is above IS_CORR_THRESH
112  */
113 #define IS_CORR_THRESH FL2FXCONST_DBL(0.95f)
114 
115 /* when expanding the IS region to more SFBs only accept an error that is
116  * not more than IS_TOTAL_ERROR_THRESH overall and
117  * not more than IS_LOCAL_ERROR_THRESH for the current SFB */
118 #define IS_TOTAL_ERROR_THRESH FL2FXCONST_DBL(0.04f)
119 #define IS_LOCAL_ERROR_THRESH FL2FXCONST_DBL(0.01f)
120 
121 /* the maximum allowed change of the intensity direction (unit: IS scale) -
122  * scaled with factor 0.25 - */
123 #define IS_DIRECTION_DEVIATION_THRESH_SF 2
124 #define IS_DIRECTION_DEVIATION_THRESH \
125   FL2FXCONST_DBL(2.0f / (1 << IS_DIRECTION_DEVIATION_THRESH_SF))
126 
127 /* IS regions need to have a minimal percentage of the overall loudness, e.g.
128  * 0.06 == 6% */
129 #define IS_REGION_MIN_LOUDNESS FL2FXCONST_DBL(0.1f)
130 
131 /* only perform IS if IS_MIN_SFBS neighboring SFBs can be processed */
132 #define IS_MIN_SFBS 6
133 
134 /* only do IS if
135  * if IS_LEFT_RIGHT_RATIO_THRESH < sfbEnergyLeft[sfb]/sfbEnergyRight[sfb] < 1 /
136  * IS_LEFT_RIGHT_RATIO_THRESH
137  * -> no IS if the panning angle is not far from the middle, MS will do */
138 /* this is equivalent to a scale of +/-1.02914634566 */
139 #define IS_LEFT_RIGHT_RATIO_THRESH FL2FXCONST_DBL(0.7f)
140 
141 /* scalefactor of realScale */
142 #define REAL_SCALE_SF 1
143 
144 /* scalefactor overallLoudness */
145 #define OVERALL_LOUDNESS_SF 6
146 
147 /* scalefactor for sum over max samples per goup */
148 #define MAX_SFB_PER_GROUP_SF 6
149 
150 /* scalefactor for sum of mdct spectrum */
151 #define MDCT_SPEC_SF 6
152 
153 typedef struct {
154   FIXP_DBL corr_thresh; /*!< Only set an IS seed it left/right channel
155                            correlation is above corr_thresh */
156 
157   FIXP_DBL total_error_thresh; /*!< When expanding the IS region to more SFBs
158                                   only accept an error that is not more than
159                                   'total_error_thresh' overall. */
160 
161   FIXP_DBL local_error_thresh; /*!< When expanding the IS region to more SFBs
162                                   only accept an error that is not more than
163                                   'local_error_thresh' for the current SFB. */
164 
165   FIXP_DBL direction_deviation_thresh; /*!< The maximum allowed change of the
166                                           intensity direction (unit: IS scale)
167                                         */
168 
169   FIXP_DBL is_region_min_loudness; /*!< IS regions need to have a minimal
170                                       percentage of the overall loudness, e.g.
171                                       0.06 == 6% */
172 
173   INT min_is_sfbs; /*!< Only perform IS if 'min_is_sfbs' neighboring SFBs can be
174                       processed */
175 
176   FIXP_DBL left_right_ratio_threshold; /*!< No IS if the panning angle is not
177                                           far from the middle, MS will do */
178 
179 } INTENSITY_PARAMETERS;
180 
181 /*****************************************************************************
182 
183     functionname: calcSfbMaxScale
184 
185     description:  Calc max value in scalefactor band
186 
187     input:        *mdctSpectrum
188                    l1
189                    l2
190 
191     output:       none
192 
193     returns:      scalefactor
194 
195 *****************************************************************************/
calcSfbMaxScale(const FIXP_DBL * mdctSpectrum,const INT l1,const INT l2)196 static INT calcSfbMaxScale(const FIXP_DBL *mdctSpectrum, const INT l1,
197                            const INT l2) {
198   INT i;
199   INT sfbMaxScale;
200   FIXP_DBL maxSpc;
201 
202   maxSpc = FL2FXCONST_DBL(0.0);
203   for (i = l1; i < l2; i++) {
204     FIXP_DBL tmp = fixp_abs((FIXP_DBL)mdctSpectrum[i]);
205     maxSpc = fixMax(maxSpc, tmp);
206   }
207   sfbMaxScale = (maxSpc == FL2FXCONST_DBL(0.0)) ? (DFRACT_BITS - 2)
208                                                 : CntLeadingZeros(maxSpc) - 1;
209 
210   return sfbMaxScale;
211 }
212 
213 /*****************************************************************************
214 
215     functionname: FDKaacEnc_initIsParams
216 
217     description:  Initialization of intensity parameters
218 
219     input:        isParams
220 
221     output:       isParams
222 
223     returns:      none
224 
225 *****************************************************************************/
FDKaacEnc_initIsParams(INTENSITY_PARAMETERS * isParams)226 static void FDKaacEnc_initIsParams(INTENSITY_PARAMETERS *isParams) {
227   isParams->corr_thresh = IS_CORR_THRESH;
228   isParams->total_error_thresh = IS_TOTAL_ERROR_THRESH;
229   isParams->local_error_thresh = IS_LOCAL_ERROR_THRESH;
230   isParams->direction_deviation_thresh = IS_DIRECTION_DEVIATION_THRESH;
231   isParams->is_region_min_loudness = IS_REGION_MIN_LOUDNESS;
232   isParams->min_is_sfbs = IS_MIN_SFBS;
233   isParams->left_right_ratio_threshold = IS_LEFT_RIGHT_RATIO_THRESH;
234 }
235 
236 /*****************************************************************************
237 
238     functionname: FDKaacEnc_prepareIntensityDecision
239 
240     description:  Prepares intensity decision
241 
242     input:        sfbEnergyLeft
243                   sfbEnergyRight
244                   sfbEnergyLdDataLeft
245                   sfbEnergyLdDataRight
246                   mdctSpectrumLeft
247                   sfbEnergyLdDataRight
248                   isParams
249 
250     output:       hrrErr            scale: none
251                   isMask            scale: none
252                   realScale         scale: LD_DATA_SHIFT + REAL_SCALE_SF
253                   normSfbLoudness   scale: none
254 
255     returns:      none
256 
257 *****************************************************************************/
FDKaacEnc_prepareIntensityDecision(const FIXP_DBL * sfbEnergyLeft,const FIXP_DBL * sfbEnergyRight,const FIXP_DBL * sfbEnergyLdDataLeft,const FIXP_DBL * sfbEnergyLdDataRight,const FIXP_DBL * mdctSpectrumLeft,const FIXP_DBL * mdctSpectrumRight,const INTENSITY_PARAMETERS * isParams,FIXP_DBL * hrrErr,INT * isMask,FIXP_DBL * realScale,FIXP_DBL * normSfbLoudness,const INT sfbCnt,const INT sfbPerGroup,const INT maxSfbPerGroup,const INT * sfbOffset)258 static void FDKaacEnc_prepareIntensityDecision(
259     const FIXP_DBL *sfbEnergyLeft, const FIXP_DBL *sfbEnergyRight,
260     const FIXP_DBL *sfbEnergyLdDataLeft, const FIXP_DBL *sfbEnergyLdDataRight,
261     const FIXP_DBL *mdctSpectrumLeft, const FIXP_DBL *mdctSpectrumRight,
262     const INTENSITY_PARAMETERS *isParams, FIXP_DBL *hrrErr, INT *isMask,
263     FIXP_DBL *realScale, FIXP_DBL *normSfbLoudness, const INT sfbCnt,
264     const INT sfbPerGroup, const INT maxSfbPerGroup, const INT *sfbOffset) {
265   INT j, sfb, sfboffs;
266   INT grpCounter;
267 
268   /* temporary variables to compute loudness */
269   FIXP_DBL overallLoudness[MAX_NO_OF_GROUPS];
270 
271   /* temporary variables to compute correlation */
272   FIXP_DBL channelCorr[MAX_GROUPED_SFB];
273   FIXP_DBL ml, mr;
274   FIXP_DBL prod_lr;
275   FIXP_DBL square_l, square_r;
276   FIXP_DBL tmp_l, tmp_r;
277   FIXP_DBL inv_n;
278 
279   FDKmemclear(channelCorr, MAX_GROUPED_SFB * sizeof(FIXP_DBL));
280   FDKmemclear(normSfbLoudness, MAX_GROUPED_SFB * sizeof(FIXP_DBL));
281   FDKmemclear(overallLoudness, MAX_NO_OF_GROUPS * sizeof(FIXP_DBL));
282   FDKmemclear(realScale, MAX_GROUPED_SFB * sizeof(FIXP_DBL));
283 
284   for (grpCounter = 0, sfboffs = 0; sfboffs < sfbCnt;
285        sfboffs += sfbPerGroup, grpCounter++) {
286     overallLoudness[grpCounter] = FL2FXCONST_DBL(0.0f);
287     for (sfb = 0; sfb < maxSfbPerGroup; sfb++) {
288       INT sL, sR, s;
289       FIXP_DBL isValue = sfbEnergyLdDataLeft[sfb + sfboffs] -
290                          sfbEnergyLdDataRight[sfb + sfboffs];
291 
292       /* delimitate intensity scale value to representable range */
293       realScale[sfb + sfboffs] = fixMin(
294           FL2FXCONST_DBL(60.f / (1 << (REAL_SCALE_SF + LD_DATA_SHIFT))),
295           fixMax(FL2FXCONST_DBL(-60.f / (1 << (REAL_SCALE_SF + LD_DATA_SHIFT))),
296                  isValue));
297 
298       sL = fixMax(0, (CntLeadingZeros(sfbEnergyLeft[sfb + sfboffs]) - 1));
299       sR = fixMax(0, (CntLeadingZeros(sfbEnergyRight[sfb + sfboffs]) - 1));
300       s = (fixMin(sL, sR) >> 2) << 2;
301       normSfbLoudness[sfb + sfboffs] =
302           sqrtFixp(sqrtFixp(((sfbEnergyLeft[sfb + sfboffs] << s) >> 1) +
303                             ((sfbEnergyRight[sfb + sfboffs] << s) >> 1))) >>
304           (s >> 2);
305 
306       overallLoudness[grpCounter] +=
307           normSfbLoudness[sfb + sfboffs] >> OVERALL_LOUDNESS_SF;
308       /* don't do intensity if
309        * - panning angle is too close to the middle or
310        * - one channel is non-existent or
311        * - if it is dual mono */
312       if ((sfbEnergyLeft[sfb + sfboffs] >=
313            fMult(isParams->left_right_ratio_threshold,
314                  sfbEnergyRight[sfb + sfboffs])) &&
315           (fMult(isParams->left_right_ratio_threshold,
316                  sfbEnergyLeft[sfb + sfboffs]) <=
317            sfbEnergyRight[sfb + sfboffs])) {
318         /* this will prevent post processing from considering this SFB for
319          * merging */
320         hrrErr[sfb + sfboffs] = FL2FXCONST_DBL(1.0 / 8.0);
321       }
322     }
323   }
324 
325   for (grpCounter = 0, sfboffs = 0; sfboffs < sfbCnt;
326        sfboffs += sfbPerGroup, grpCounter++) {
327     INT invOverallLoudnessSF;
328     FIXP_DBL invOverallLoudness;
329 
330     if (overallLoudness[grpCounter] == FL2FXCONST_DBL(0.0)) {
331       invOverallLoudness = FL2FXCONST_DBL(0.0);
332       invOverallLoudnessSF = 0;
333     } else {
334       invOverallLoudness =
335           fDivNorm((FIXP_DBL)MAXVAL_DBL, overallLoudness[grpCounter],
336                    &invOverallLoudnessSF);
337       invOverallLoudnessSF =
338           invOverallLoudnessSF - OVERALL_LOUDNESS_SF +
339           1; /* +1: compensate fMultDiv2() in subsequent loop */
340     }
341     invOverallLoudnessSF = fixMin(
342         fixMax(invOverallLoudnessSF, -(DFRACT_BITS - 1)), DFRACT_BITS - 1);
343 
344     for (sfb = 0; sfb < maxSfbPerGroup; sfb++) {
345       FIXP_DBL tmp;
346 
347       tmp = fMultDiv2((normSfbLoudness[sfb + sfboffs] >> OVERALL_LOUDNESS_SF)
348                           << OVERALL_LOUDNESS_SF,
349                       invOverallLoudness);
350 
351       normSfbLoudness[sfb + sfboffs] = scaleValue(tmp, invOverallLoudnessSF);
352 
353       channelCorr[sfb + sfboffs] = FL2FXCONST_DBL(0.0f);
354 
355       /* max width of scalefactorband is 96; width's are always even */
356       /* inv_n is scaled with factor 2 to compensate fMultDiv2() in subsequent
357        * loops */
358       inv_n = GetInvInt(
359           (sfbOffset[sfb + sfboffs + 1] - sfbOffset[sfb + sfboffs]) >> 1);
360 
361       if (inv_n > FL2FXCONST_DBL(0.0f)) {
362         INT s, sL, sR;
363 
364         /* correlation := Pearson's product-moment coefficient */
365         /* compute correlation between channels and check if it is over
366          * threshold */
367         ml = FL2FXCONST_DBL(0.0f);
368         mr = FL2FXCONST_DBL(0.0f);
369         prod_lr = FL2FXCONST_DBL(0.0f);
370         square_l = FL2FXCONST_DBL(0.0f);
371         square_r = FL2FXCONST_DBL(0.0f);
372 
373         sL = calcSfbMaxScale(mdctSpectrumLeft, sfbOffset[sfb + sfboffs],
374                              sfbOffset[sfb + sfboffs + 1]);
375         sR = calcSfbMaxScale(mdctSpectrumRight, sfbOffset[sfb + sfboffs],
376                              sfbOffset[sfb + sfboffs + 1]);
377         s = fixMin(sL, sR);
378 
379         for (j = sfbOffset[sfb + sfboffs]; j < sfbOffset[sfb + sfboffs + 1];
380              j++) {
381           ml += fMultDiv2((mdctSpectrumLeft[j] << s),
382                           inv_n);  // scaled with mdctScale - s + inv_n
383           mr += fMultDiv2((mdctSpectrumRight[j] << s),
384                           inv_n);  // scaled with mdctScale - s + inv_n
385         }
386         ml = fMultDiv2(ml, inv_n);  // scaled with mdctScale - s + inv_n
387         mr = fMultDiv2(mr, inv_n);  // scaled with mdctScale - s + inv_n
388 
389         for (j = sfbOffset[sfb + sfboffs]; j < sfbOffset[sfb + sfboffs + 1];
390              j++) {
391           tmp_l = fMultDiv2((mdctSpectrumLeft[j] << s), inv_n) -
392                   ml;  // scaled with mdctScale - s + inv_n
393           tmp_r = fMultDiv2((mdctSpectrumRight[j] << s), inv_n) -
394                   mr;  // scaled with mdctScale - s + inv_n
395 
396           prod_lr += fMultDiv2(
397               tmp_l, tmp_r);  // scaled with 2*(mdctScale - s + inv_n) + 1
398           square_l +=
399               fPow2Div2(tmp_l);  // scaled with 2*(mdctScale - s + inv_n) + 1
400           square_r +=
401               fPow2Div2(tmp_r);  // scaled with 2*(mdctScale - s + inv_n) + 1
402         }
403         prod_lr = prod_lr << 1;    // scaled with 2*(mdctScale - s + inv_n)
404         square_l = square_l << 1;  // scaled with 2*(mdctScale - s + inv_n)
405         square_r = square_r << 1;  // scaled with 2*(mdctScale - s + inv_n)
406 
407         if (square_l > FL2FXCONST_DBL(0.0f) &&
408             square_r > FL2FXCONST_DBL(0.0f)) {
409           INT channelCorrSF = 0;
410 
411           /* local scaling of square_l and square_r is compensated after sqrt
412            * calculation */
413           sL = fixMax(0, (CntLeadingZeros(square_l) - 1));
414           sR = fixMax(0, (CntLeadingZeros(square_r) - 1));
415           s = ((sL + sR) >> 1) << 1;
416           sL = fixMin(sL, s);
417           sR = s - sL;
418           tmp = fMult(square_l << sL, square_r << sR);
419           tmp = sqrtFixp(tmp);
420 
421           FDK_ASSERT(tmp > FL2FXCONST_DBL(0.0f));
422 
423           /* numerator and denominator have the same scaling */
424           if (prod_lr < FL2FXCONST_DBL(0.0f)) {
425             channelCorr[sfb + sfboffs] =
426                 -(fDivNorm(-prod_lr, tmp, &channelCorrSF));
427 
428           } else {
429             channelCorr[sfb + sfboffs] =
430                 (fDivNorm(prod_lr, tmp, &channelCorrSF));
431           }
432           channelCorrSF = fixMin(
433               fixMax((channelCorrSF + ((sL + sR) >> 1)), -(DFRACT_BITS - 1)),
434               DFRACT_BITS - 1);
435 
436           if (channelCorrSF < 0) {
437             channelCorr[sfb + sfboffs] =
438                 channelCorr[sfb + sfboffs] >> (-channelCorrSF);
439           } else {
440             /* avoid overflows due to limited computational accuracy */
441             if (fAbs(channelCorr[sfb + sfboffs]) >
442                 (((FIXP_DBL)MAXVAL_DBL) >> channelCorrSF)) {
443               if (channelCorr[sfb + sfboffs] < FL2FXCONST_DBL(0.0f))
444                 channelCorr[sfb + sfboffs] = -(FIXP_DBL)MAXVAL_DBL;
445               else
446                 channelCorr[sfb + sfboffs] = (FIXP_DBL)MAXVAL_DBL;
447             } else {
448               channelCorr[sfb + sfboffs] = channelCorr[sfb + sfboffs]
449                                            << channelCorrSF;
450             }
451           }
452         }
453       }
454 
455       /* for post processing: hrrErr is the error in terms of (too little)
456        * correlation weighted with the loudness of the SFB; SFBs with small
457        * hrrErr can be merged */
458       if (hrrErr[sfb + sfboffs] == FL2FXCONST_DBL(1.0 / 8.0)) {
459         continue;
460       }
461 
462       hrrErr[sfb + sfboffs] =
463           fMultDiv2((FL2FXCONST_DBL(0.25f) - (channelCorr[sfb + sfboffs] >> 2)),
464                     normSfbLoudness[sfb + sfboffs]);
465 
466       /* set IS mask/vector to 1, if correlation is high enough */
467       if (fAbs(channelCorr[sfb + sfboffs]) >= isParams->corr_thresh) {
468         isMask[sfb + sfboffs] = 1;
469       }
470     }
471   }
472 }
473 
474 /*****************************************************************************
475 
476     functionname: FDKaacEnc_finalizeIntensityDecision
477 
478     description:  Finalizes intensity decision
479 
480     input:        isParams          scale: none
481                   hrrErr            scale: none
482                   realIsScale       scale: LD_DATA_SHIFT + REAL_SCALE_SF
483                   normSfbLoudness   scale: none
484 
485     output:       isMask            scale: none
486 
487     returns:      none
488 
489 *****************************************************************************/
FDKaacEnc_finalizeIntensityDecision(const FIXP_DBL * hrrErr,INT * isMask,const FIXP_DBL * realIsScale,const FIXP_DBL * normSfbLoudness,const INTENSITY_PARAMETERS * isParams,const INT sfbCnt,const INT sfbPerGroup,const INT maxSfbPerGroup)490 static void FDKaacEnc_finalizeIntensityDecision(
491     const FIXP_DBL *hrrErr, INT *isMask, const FIXP_DBL *realIsScale,
492     const FIXP_DBL *normSfbLoudness, const INTENSITY_PARAMETERS *isParams,
493     const INT sfbCnt, const INT sfbPerGroup, const INT maxSfbPerGroup) {
494   INT sfb, sfboffs, j;
495   FIXP_DBL isScaleLast = FL2FXCONST_DBL(0.0f);
496   INT isStartValueFound = 0;
497 
498   for (sfboffs = 0; sfboffs < sfbCnt; sfboffs += sfbPerGroup) {
499     INT startIsSfb = 0;
500     INT inIsBlock = 0;
501     INT currentIsSfbCount = 0;
502     FIXP_DBL overallHrrError = FL2FXCONST_DBL(0.0f);
503     FIXP_DBL isRegionLoudness = FL2FXCONST_DBL(0.0f);
504 
505     for (sfb = 0; sfb < maxSfbPerGroup; sfb++) {
506       if (isMask[sfboffs + sfb] == 1) {
507         if (currentIsSfbCount == 0) {
508           startIsSfb = sfboffs + sfb;
509         }
510         if (isStartValueFound == 0) {
511           isScaleLast = realIsScale[sfboffs + sfb];
512           isStartValueFound = 1;
513         }
514         inIsBlock = 1;
515         currentIsSfbCount++;
516         overallHrrError += hrrErr[sfboffs + sfb] >> (MAX_SFB_PER_GROUP_SF - 3);
517         isRegionLoudness +=
518             normSfbLoudness[sfboffs + sfb] >> MAX_SFB_PER_GROUP_SF;
519       } else {
520         /* based on correlation, IS should not be used
521          * -> use it anyway, if overall error is below threshold
522          *    and if local error does not exceed threshold
523          * otherwise: check if there are enough IS SFBs
524          */
525         if (inIsBlock) {
526           overallHrrError +=
527               hrrErr[sfboffs + sfb] >> (MAX_SFB_PER_GROUP_SF - 3);
528           isRegionLoudness +=
529               normSfbLoudness[sfboffs + sfb] >> MAX_SFB_PER_GROUP_SF;
530 
531           if ((hrrErr[sfboffs + sfb] < (isParams->local_error_thresh >> 3)) &&
532               (overallHrrError <
533                (isParams->total_error_thresh >> MAX_SFB_PER_GROUP_SF))) {
534             currentIsSfbCount++;
535             /* overwrite correlation based decision */
536             isMask[sfboffs + sfb] = 1;
537           } else {
538             inIsBlock = 0;
539           }
540         }
541       }
542       /* check for large direction deviation */
543       if (inIsBlock) {
544         if (fAbs(isScaleLast - realIsScale[sfboffs + sfb]) <
545             (isParams->direction_deviation_thresh >>
546              (REAL_SCALE_SF + LD_DATA_SHIFT -
547               IS_DIRECTION_DEVIATION_THRESH_SF))) {
548           isScaleLast = realIsScale[sfboffs + sfb];
549         } else {
550           isMask[sfboffs + sfb] = 0;
551           inIsBlock = 0;
552           currentIsSfbCount--;
553         }
554       }
555 
556       if (currentIsSfbCount > 0 && (!inIsBlock || sfb == maxSfbPerGroup - 1)) {
557         /* not enough SFBs -> do not use IS */
558         if (currentIsSfbCount < isParams->min_is_sfbs ||
559             (isRegionLoudness<isParams->is_region_min_loudness>>
560              MAX_SFB_PER_GROUP_SF)) {
561           for (j = startIsSfb; j <= sfboffs + sfb; j++) {
562             isMask[j] = 0;
563           }
564           isScaleLast = FL2FXCONST_DBL(0.0f);
565           isStartValueFound = 0;
566           for (j = 0; j < startIsSfb; j++) {
567             if (isMask[j] != 0) {
568               isScaleLast = realIsScale[j];
569               isStartValueFound = 1;
570             }
571           }
572         }
573         currentIsSfbCount = 0;
574         overallHrrError = FL2FXCONST_DBL(0.0f);
575         isRegionLoudness = FL2FXCONST_DBL(0.0f);
576       }
577     }
578   }
579 }
580 
581 /*****************************************************************************
582 
583     functionname: FDKaacEnc_IntensityStereoProcessing
584 
585     description:  Intensity stereo processing tool
586 
587     input:        sfbEnergyLeft
588                   sfbEnergyRight
589                   mdctSpectrumLeft
590                   mdctSpectrumRight
591                   sfbThresholdLeft
592                   sfbThresholdRight
593                   sfbSpreadEnLeft
594                   sfbSpreadEnRight
595                   sfbEnergyLdDataLeft
596                   sfbEnergyLdDataRight
597 
598     output:       isBook
599                   isScale
600                   pnsData->pnsFlag
601                   msDigest                 zeroed from start to sfbCnt
602                   msMask                   zeroed from start to sfbCnt
603                   mdctSpectrumRight        zeroed where isBook!=0
604                   sfbEnergyRight           zeroed where isBook!=0
605                   sfbSpreadEnRight       zeroed where isBook!=0
606                   sfbThresholdRight        zeroed where isBook!=0
607                   sfbEnergyLdDataRight     FL2FXCONST_DBL(-1.0) where isBook!=0
608                   sfbThresholdLdDataRight  FL2FXCONST_DBL(-0.515625f) where
609 isBook!=0
610 
611     returns:      none
612 
613 *****************************************************************************/
FDKaacEnc_IntensityStereoProcessing(FIXP_DBL * sfbEnergyLeft,FIXP_DBL * sfbEnergyRight,FIXP_DBL * mdctSpectrumLeft,FIXP_DBL * mdctSpectrumRight,FIXP_DBL * sfbThresholdLeft,FIXP_DBL * sfbThresholdRight,FIXP_DBL * sfbThresholdLdDataRight,FIXP_DBL * sfbSpreadEnLeft,FIXP_DBL * sfbSpreadEnRight,FIXP_DBL * sfbEnergyLdDataLeft,FIXP_DBL * sfbEnergyLdDataRight,INT * msDigest,INT * msMask,const INT sfbCnt,const INT sfbPerGroup,const INT maxSfbPerGroup,const INT * sfbOffset,const INT allowIS,INT * isBook,INT * isScale,PNS_DATA * RESTRICT pnsData[2])614 void FDKaacEnc_IntensityStereoProcessing(
615     FIXP_DBL *sfbEnergyLeft, FIXP_DBL *sfbEnergyRight,
616     FIXP_DBL *mdctSpectrumLeft, FIXP_DBL *mdctSpectrumRight,
617     FIXP_DBL *sfbThresholdLeft, FIXP_DBL *sfbThresholdRight,
618     FIXP_DBL *sfbThresholdLdDataRight, FIXP_DBL *sfbSpreadEnLeft,
619     FIXP_DBL *sfbSpreadEnRight, FIXP_DBL *sfbEnergyLdDataLeft,
620     FIXP_DBL *sfbEnergyLdDataRight, INT *msDigest, INT *msMask,
621     const INT sfbCnt, const INT sfbPerGroup, const INT maxSfbPerGroup,
622     const INT *sfbOffset, const INT allowIS, INT *isBook, INT *isScale,
623     PNS_DATA *RESTRICT pnsData[2]) {
624   INT sfb, sfboffs, j;
625   FIXP_DBL scale;
626   FIXP_DBL lr;
627   FIXP_DBL hrrErr[MAX_GROUPED_SFB];
628   FIXP_DBL normSfbLoudness[MAX_GROUPED_SFB];
629   FIXP_DBL realIsScale[MAX_GROUPED_SFB];
630   INTENSITY_PARAMETERS isParams;
631   INT isMask[MAX_GROUPED_SFB];
632 
633   FDKmemclear((void *)isBook, sfbCnt * sizeof(INT));
634   FDKmemclear((void *)isMask, sfbCnt * sizeof(INT));
635   FDKmemclear((void *)realIsScale, sfbCnt * sizeof(FIXP_DBL));
636   FDKmemclear((void *)isScale, sfbCnt * sizeof(INT));
637   FDKmemclear((void *)hrrErr, sfbCnt * sizeof(FIXP_DBL));
638 
639   if (!allowIS) return;
640 
641   FDKaacEnc_initIsParams(&isParams);
642 
643   /* compute / set the following values per SFB:
644    * - left/right ratio between channels
645    * - normalized loudness
646    *   + loudness == average of energy in channels to 0.25
647    *   + normalization: division by sum of all SFB loudnesses
648    * - isMask (is set to 0 if channels are the same or one is 0)
649    */
650   FDKaacEnc_prepareIntensityDecision(
651       sfbEnergyLeft, sfbEnergyRight, sfbEnergyLdDataLeft, sfbEnergyLdDataRight,
652       mdctSpectrumLeft, mdctSpectrumRight, &isParams, hrrErr, isMask,
653       realIsScale, normSfbLoudness, sfbCnt, sfbPerGroup, maxSfbPerGroup,
654       sfbOffset);
655 
656   FDKaacEnc_finalizeIntensityDecision(hrrErr, isMask, realIsScale,
657                                       normSfbLoudness, &isParams, sfbCnt,
658                                       sfbPerGroup, maxSfbPerGroup);
659 
660   for (sfb = 0; sfb < sfbCnt; sfb += sfbPerGroup) {
661     for (sfboffs = 0; sfboffs < maxSfbPerGroup; sfboffs++) {
662       INT sL, sR;
663       FIXP_DBL inv_n;
664       INT mdct_spec_sf = MDCT_SPEC_SF;
665 
666       msMask[sfb + sfboffs] = 0;
667       if (isMask[sfb + sfboffs] == 0) {
668         continue;
669       }
670 
671       if ((sfbEnergyLeft[sfb + sfboffs] < sfbThresholdLeft[sfb + sfboffs]) &&
672           (fMult(FL2FXCONST_DBL(1.0f / 1.5f), sfbEnergyRight[sfb + sfboffs]) >
673            sfbThresholdRight[sfb + sfboffs])) {
674         continue;
675       }
676       /* NEW: if there is a big-enough IS region, switch off PNS */
677       if (pnsData[0]) {
678         if (pnsData[0]->pnsFlag[sfb + sfboffs]) {
679           pnsData[0]->pnsFlag[sfb + sfboffs] = 0;
680         }
681         if (pnsData[1]->pnsFlag[sfb + sfboffs]) {
682           pnsData[1]->pnsFlag[sfb + sfboffs] = 0;
683         }
684       }
685 
686       if (sfbOffset[sfb + sfboffs + 1] - sfbOffset[sfb + sfboffs] >
687           1 << mdct_spec_sf) {
688         mdct_spec_sf++; /* This is for rare cases where the number of bins in a
689                            scale factor band is > 64 */
690       }
691 
692       inv_n = GetInvInt(
693           (sfbOffset[sfb + sfboffs + 1] - sfbOffset[sfb + sfboffs]) >>
694           1);  // scaled with 2 to compensate fMultDiv2() in subsequent loop
695       sL = calcSfbMaxScale(mdctSpectrumLeft, sfbOffset[sfb + sfboffs],
696                            sfbOffset[sfb + sfboffs + 1]);
697       sR = calcSfbMaxScale(mdctSpectrumRight, sfbOffset[sfb + sfboffs],
698                            sfbOffset[sfb + sfboffs + 1]);
699 
700       lr = FL2FXCONST_DBL(0.0f);
701       for (j = sfbOffset[sfb + sfboffs]; j < sfbOffset[sfb + sfboffs + 1]; j++)
702         lr += fMultDiv2(
703             fMultDiv2(mdctSpectrumLeft[j] << sL, mdctSpectrumRight[j] << sR),
704             inv_n);
705       lr = lr << 1;
706 
707       if (lr < FL2FXCONST_DBL(0.0f)) {
708         /* This means OUT OF phase intensity stereo, cf. standard */
709         INT s0, s1, s2;
710         FIXP_DBL tmp, d, ed = FL2FXCONST_DBL(0.0f);
711 
712         s0 = fixMin(sL, sR);
713         for (j = sfbOffset[sfb + sfboffs]; j < sfbOffset[sfb + sfboffs + 1];
714              j++) {
715           d = ((mdctSpectrumLeft[j] << s0) >> 1) -
716               ((mdctSpectrumRight[j] << s0) >> 1);
717           ed += fMultDiv2(d, d) >> (mdct_spec_sf - 1);
718         }
719         msMask[sfb + sfboffs] = 1;
720         tmp = fDivNorm(sfbEnergyLeft[sfb + sfboffs], ed, &s1);
721         s2 = (s1) + (2 * s0) - 2 - mdct_spec_sf;
722         if (s2 & 1) {
723           tmp = tmp >> 1;
724           s2 = s2 + 1;
725         }
726         s2 = (s2 >> 1) + 1;  // +1 compensate fMultDiv2() in subsequent loop
727         s2 = fixMin(fixMax(s2, -(DFRACT_BITS - 1)), (DFRACT_BITS - 1));
728         scale = sqrtFixp(tmp);
729         if (s2 < 0) {
730           s2 = -s2;
731           for (j = sfbOffset[sfb + sfboffs]; j < sfbOffset[sfb + sfboffs + 1];
732                j++) {
733             mdctSpectrumLeft[j] = (fMultDiv2(mdctSpectrumLeft[j], scale) -
734                                    fMultDiv2(mdctSpectrumRight[j], scale)) >>
735                                   s2;
736             mdctSpectrumRight[j] = FL2FXCONST_DBL(0.0f);
737           }
738         } else {
739           for (j = sfbOffset[sfb + sfboffs]; j < sfbOffset[sfb + sfboffs + 1];
740                j++) {
741             mdctSpectrumLeft[j] = (fMultDiv2(mdctSpectrumLeft[j], scale) -
742                                    fMultDiv2(mdctSpectrumRight[j], scale))
743                                   << s2;
744             mdctSpectrumRight[j] = FL2FXCONST_DBL(0.0f);
745           }
746         }
747       } else {
748         /* This means IN phase intensity stereo, cf. standard */
749         INT s0, s1, s2;
750         FIXP_DBL tmp, s, es = FL2FXCONST_DBL(0.0f);
751 
752         s0 = fixMin(sL, sR);
753         for (j = sfbOffset[sfb + sfboffs]; j < sfbOffset[sfb + sfboffs + 1];
754              j++) {
755           s = ((mdctSpectrumLeft[j] << s0) >> 1) +
756               ((mdctSpectrumRight[j] << s0) >> 1);
757           es += fMultDiv2(s, s) >>
758                 (mdct_spec_sf -
759                  1);  // scaled 2*(mdctScale - s0 + 1) + mdct_spec_sf
760         }
761         msMask[sfb + sfboffs] = 0;
762         tmp = fDivNorm(sfbEnergyLeft[sfb + sfboffs], es, &s1);
763         s2 = (s1) + (2 * s0) - 2 - mdct_spec_sf;
764         if (s2 & 1) {
765           tmp = tmp >> 1;
766           s2 = s2 + 1;
767         }
768         s2 = (s2 >> 1) + 1;  // +1 compensate fMultDiv2() in subsequent loop
769         s2 = fixMin(fixMax(s2, -(DFRACT_BITS - 1)), (DFRACT_BITS - 1));
770         scale = sqrtFixp(tmp);
771         if (s2 < 0) {
772           s2 = -s2;
773           for (j = sfbOffset[sfb + sfboffs]; j < sfbOffset[sfb + sfboffs + 1];
774                j++) {
775             mdctSpectrumLeft[j] = (fMultDiv2(mdctSpectrumLeft[j], scale) +
776                                    fMultDiv2(mdctSpectrumRight[j], scale)) >>
777                                   s2;
778             mdctSpectrumRight[j] = FL2FXCONST_DBL(0.0f);
779           }
780         } else {
781           for (j = sfbOffset[sfb + sfboffs]; j < sfbOffset[sfb + sfboffs + 1];
782                j++) {
783             mdctSpectrumLeft[j] = (fMultDiv2(mdctSpectrumLeft[j], scale) +
784                                    fMultDiv2(mdctSpectrumRight[j], scale))
785                                   << s2;
786             mdctSpectrumRight[j] = FL2FXCONST_DBL(0.0f);
787           }
788         }
789       }
790 
791       isBook[sfb + sfboffs] = CODE_BOOK_IS_IN_PHASE_NO;
792 
793       if (realIsScale[sfb + sfboffs] < FL2FXCONST_DBL(0.0f)) {
794         isScale[sfb + sfboffs] =
795             (INT)(((realIsScale[sfb + sfboffs] >> 1) -
796                    FL2FXCONST_DBL(
797                        0.5f / (1 << (REAL_SCALE_SF + LD_DATA_SHIFT + 1)))) >>
798                   (DFRACT_BITS - 1 - REAL_SCALE_SF - LD_DATA_SHIFT - 1)) +
799             1;
800       } else {
801         isScale[sfb + sfboffs] =
802             (INT)(((realIsScale[sfb + sfboffs] >> 1) +
803                    FL2FXCONST_DBL(
804                        0.5f / (1 << (REAL_SCALE_SF + LD_DATA_SHIFT + 1)))) >>
805                   (DFRACT_BITS - 1 - REAL_SCALE_SF - LD_DATA_SHIFT - 1));
806       }
807 
808       sfbEnergyRight[sfb + sfboffs] = FL2FXCONST_DBL(0.0f);
809       sfbEnergyLdDataRight[sfb + sfboffs] = FL2FXCONST_DBL(-1.0f);
810       sfbThresholdRight[sfb + sfboffs] = FL2FXCONST_DBL(0.0f);
811       sfbThresholdLdDataRight[sfb + sfboffs] = FL2FXCONST_DBL(-0.515625f);
812       sfbSpreadEnRight[sfb + sfboffs] = FL2FXCONST_DBL(0.0f);
813 
814       *msDigest = MS_SOME;
815     }
816   }
817 }
818