libfaad2/sbr_hfgen.c

   1 /*
   2 ** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding
   3 ** Copyright (C) 2003-2004 M. Bakker, Ahead Software AG, http://www.nero.com
   4 **
   5 ** This program is free software; you can redistribute it and/or modify
   6 ** it under the terms of the GNU General Public License as published by
   7 ** the Free Software Foundation; either version 2 of the License, or
   8 ** (at your option) any later version.
   9 **
  10 ** This program is distributed in the hope that it will be useful,
  11 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 ** GNU General Public License for more details.
  14 **
  15 ** You should have received a copy of the GNU General Public License
  16 ** along with this program; if not, write to the Free Software
  17 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  18 **
  19 ** Any non-GPL usage of this software or parts of this software is strictly
  20 ** forbidden.
  21 **
  22 ** Commercial non-GPL licensing of this software is possible.
  23 ** For more info contact Ahead Software through Mpeg4AAClicense@nero.com.
  24 **
  25 ** Initially modified for use with MPlayer by Arpad Gereöffy on 2003/08/30
  26 ** $Id: sbr_hfgen.c,v 1.3 2004/06/02 22:59:03 diego Exp $
  27 ** detailed CVS changelog at http://www.mplayerhq.hu/cgi-bin/cvsweb.cgi/main/
  28 **/
  29
  30 /* High Frequency generation */
  31
  32 #include "common.h"
  33 #include "structs.h"
  34
  35 #ifdef SBR_DEC
  36
  37 #include "sbr_syntax.h"
  38 #include "sbr_hfgen.h"
  39 #include "sbr_fbt.h"
  40
  41
  42 /* static function declarations */
  43 static void calc_prediction_coef(sbr_info *sbr, qmf_t Xlow[MAX_NTSRHFG][32],
  44                                  complex_t *alpha_0, complex_t *alpha_1
  45 #ifdef SBR_LOW_POWER
  46                                  , real_t *rxx
  47 #endif
  48                                  );
  49 #ifdef SBR_LOW_POWER
  50 static void calc_aliasing_degree(sbr_info *sbr, real_t *rxx, real_t *deg);
  51 #endif
  52 static void calc_chirp_factors(sbr_info *sbr, uint8_t ch);
  53 static void patch_construction(sbr_info *sbr);
  54
  55
  56 void hf_generation(sbr_info *sbr, qmf_t Xlow[MAX_NTSRHFG][32],
  57                    qmf_t Xhigh[MAX_NTSRHFG][64]
  58 #ifdef SBR_LOW_POWER
  59                    ,real_t *deg
  60 #endif
  61                    ,uint8_t ch)
  62 {
  63     uint8_t l, i, x;
  64     ALIGN complex_t alpha_0[64], alpha_1[64];
  65 #ifdef SBR_LOW_POWER
  66     ALIGN real_t rxx[64];
  67 #endif
  68
  69     uint8_t offset = sbr->tHFAdj;
  70     uint8_t first = sbr->t_E[ch][0];
  71     uint8_t last = sbr->t_E[ch][sbr->L_E[ch]];
  72
  73 //    printf("%d %d\n", first, last);
  74
  75     calc_chirp_factors(sbr, ch);
  76
  77     for (i = first; i < last; i++)
  78     {
  79         memset(Xhigh[i + offset], 0, 64 * sizeof(qmf_t));
  80     }
  81
  82     if ((ch == 0) && (sbr->Reset))
  83         patch_construction(sbr);
  84
  85     /* calculate the prediction coefficients */
  86     calc_prediction_coef(sbr, Xlow, alpha_0, alpha_1
  87 #ifdef SBR_LOW_POWER
  88         , rxx
  89 #endif
  90         );
  91
  92 #ifdef SBR_LOW_POWER
  93     calc_aliasing_degree(sbr, rxx, deg);
  94 #endif
  95
  96     /* actual HF generation */
  97     for (i = 0; i < sbr->noPatches; i++)
  98     {
  99         for (x = 0; x < sbr->patchNoSubbands[i]; x++)
 100         {
 101             complex_t a0, a1;
 102             real_t bw, bw2;
 103             uint8_t q, p, k, g;
 104
 105             /* find the low and high band for patching */
 106             k = sbr->kx + x;
 107             for (q = 0; q < i; q++)
 108             {
 109                 k += sbr->patchNoSubbands[q];
 110             }
 111             p = sbr->patchStartSubband[i] + x;
 112
 113 #ifdef SBR_LOW_POWER
 114             if (x != 0 /*x < sbr->patchNoSubbands[i]-1*/)
 115                 deg[k] = deg[p];
 116             else
 117                 deg[k] = 0;
 118 #endif
 119
 120             g = sbr->table_map_k_to_g[k];
 121
 122             bw = sbr->bwArray[ch][g];
 123             bw2 = MUL_C(bw, bw);
 124
 125             /* do the patching */
 126             /* with or without filtering */
 127             if (bw2 > 0)
 128             {
 129                 RE(a0) = MUL_C(RE(alpha_0[p]), bw);
 130                 RE(a1) = MUL_C(RE(alpha_1[p]), bw2);
 131 #ifndef SBR_LOW_POWER
 132                 IM(a0) = MUL_C(IM(alpha_0[p]), bw);
 133                 IM(a1) = MUL_C(IM(alpha_1[p]), bw2);
 134 #endif
 135
 136                                 for (l = first; l < last; l++)
 137                 {
 138                     QMF_RE(Xhigh[l + offset][k]) = QMF_RE(Xlow[l + offset][p]);
 139 #ifndef SBR_LOW_POWER
 140                     QMF_IM(Xhigh[l + offset][k]) = QMF_IM(Xlow[l + offset][p]);
 141 #endif
 142
 143 #ifdef SBR_LOW_POWER
 144                     QMF_RE(Xhigh[l + offset][k]) += (
 145                         MUL_R(RE(a0), QMF_RE(Xlow[l - 1 + offset][p])) +
 146                         MUL_R(RE(a1), QMF_RE(Xlow[l - 2 + offset][p])));
 147 #else
 148                     QMF_RE(Xhigh[l + offset][k]) += (
 149                         RE(a0) * QMF_RE(Xlow[l - 1 + offset][p]) -
 150                         IM(a0) * QMF_IM(Xlow[l - 1 + offset][p]) +
 151                         RE(a1) * QMF_RE(Xlow[l - 2 + offset][p]) -
 152                         IM(a1) * QMF_IM(Xlow[l - 2 + offset][p]));
 153                     QMF_IM(Xhigh[l + offset][k]) += (
 154                         IM(a0) * QMF_RE(Xlow[l - 1 + offset][p]) +
 155                         RE(a0) * QMF_IM(Xlow[l - 1 + offset][p]) +
 156                         IM(a1) * QMF_RE(Xlow[l - 2 + offset][p]) +
 157                         RE(a1) * QMF_IM(Xlow[l - 2 + offset][p]));
 158 #endif
 159                 }
 160             } else {
 161                 for (l = first; l < last; l++)
 162                 {
 163                     QMF_RE(Xhigh[l + offset][k]) = QMF_RE(Xlow[l + offset][p]);
 164 #ifndef SBR_LOW_POWER
 165                     QMF_IM(Xhigh[l + offset][k]) = QMF_IM(Xlow[l + offset][p]);
 166 #endif
 167                 }
 168             }
 169         }
 170     }
 171
 172     if (sbr->Reset)
 173     {
 174         limiter_frequency_table(sbr);
 175     }
 176 }
 177
 178 typedef struct
 179 {
 180     complex_t r01;
 181     complex_t r02;
 182     complex_t r11;
 183     complex_t r12;
 184     complex_t r22;
 185     real_t det;
 186 } acorr_coef;
 187
 188 #define SBR_ABS(A) ((A) < 0) ? -(A) : (A)
 189
 190 #ifdef SBR_LOW_POWER
 191 static void auto_correlation(sbr_info *sbr, acorr_coef *ac,
 192                              qmf_t buffer[MAX_NTSRHFG][32],
 193                              uint8_t bd, uint8_t len)
 194 {
 195     real_t r01 = 0, r02 = 0, r11 = 0;
 196     int8_t j;
 197     uint8_t offset = sbr->tHFAdj;
 198     const real_t rel = 1 / (1 + 1e-6f);
 199
 200
 201     for (j = offset; j < len + offset; j++)
 202     {
 203         r01 += QMF_RE(buffer[j][bd]) * QMF_RE(buffer[j-1][bd]);
 204         r02 += QMF_RE(buffer[j][bd]) * QMF_RE(buffer[j-2][bd]);
 205         r11 += QMF_RE(buffer[j-1][bd]) * QMF_RE(buffer[j-1][bd]);
 206     }
 207     RE(ac->r12) = r01 -
 208         QMF_RE(buffer[len+offset-1][bd]) * QMF_RE(buffer[len+offset-2][bd]) +
 209         QMF_RE(buffer[offset-1][bd]) * QMF_RE(buffer[offset-2][bd]);
 210     RE(ac->r22) = r11 -
 211         QMF_RE(buffer[len+offset-2][bd]) * QMF_RE(buffer[len+offset-2][bd]) +
 212         QMF_RE(buffer[offset-2][bd]) * QMF_RE(buffer[offset-2][bd]);
 213     RE(ac->r01) = r01;
 214     RE(ac->r02) = r02;
 215     RE(ac->r11) = r11;
 216
 217     ac->det = MUL_R(RE(ac->r11), RE(ac->r22)) - MUL_C(MUL_R(RE(ac->r12), RE(ac->r12)), rel);
 218 }
 219 #else
 220 static void auto_correlation(sbr_info *sbr, acorr_coef *ac, qmf_t buffer[MAX_NTSRHFG][32],
 221                              uint8_t bd, uint8_t len)
 222 {
 223     real_t r01r = 0, r01i = 0, r02r = 0, r02i = 0, r11r = 0;
 224     const real_t rel = 1 / (1 + 1e-6f);
 225     int8_t j;
 226     uint8_t offset = sbr->tHFAdj;
 227
 228
 229     for (j = offset; j < len + offset; j++)
 230     {
 231         r01r += QMF_RE(buffer[j][bd]) * QMF_RE(buffer[j-1][bd]) +
 232             QMF_IM(buffer[j][bd]) * QMF_IM(buffer[j-1][bd]);
 233         r01i += QMF_IM(buffer[j][bd]) * QMF_RE(buffer[j-1][bd]) -
 234             QMF_RE(buffer[j][bd]) * QMF_IM(buffer[j-1][bd]);
 235         r02r += QMF_RE(buffer[j][bd]) * QMF_RE(buffer[j-2][bd]) +
 236             QMF_IM(buffer[j][bd]) * QMF_IM(buffer[j-2][bd]);
 237         r02i += QMF_IM(buffer[j][bd]) * QMF_RE(buffer[j-2][bd]) -
 238             QMF_RE(buffer[j][bd]) * QMF_IM(buffer[j-2][bd]);
 239         r11r += QMF_RE(buffer[j-1][bd]) * QMF_RE(buffer[j-1][bd]) +
 240             QMF_IM(buffer[j-1][bd]) * QMF_IM(buffer[j-1][bd]);
 241     }
 242
 243     RE(ac->r01) = r01r;
 244     IM(ac->r01) = r01i;
 245     RE(ac->r02) = r02r;
 246     IM(ac->r02) = r02i;
 247     RE(ac->r11) = r11r;
 248
 249     RE(ac->r12) = r01r -
 250         (QMF_RE(buffer[len+offset-1][bd]) * QMF_RE(buffer[len+offset-2][bd]) + QMF_IM(buffer[len+offset-1][bd]) * QMF_IM(buffer[len+offset-2][bd])) +
 251         (QMF_RE(buffer[offset-1][bd]) * QMF_RE(buffer[offset-2][bd]) + QMF_IM(buffer[offset-1][bd]) * QMF_IM(buffer[offset-2][bd]));
 252     IM(ac->r12) = r01i -
 253         (QMF_IM(buffer[len+offset-1][bd]) * QMF_RE(buffer[len+offset-2][bd]) - QMF_RE(buffer[len+offset-1][bd]) * QMF_IM(buffer[len+offset-2][bd])) +
 254         (QMF_IM(buffer[offset-1][bd]) * QMF_RE(buffer[offset-2][bd]) - QMF_RE(buffer[offset-1][bd]) * QMF_IM(buffer[offset-2][bd]));
 255     RE(ac->r22) = r11r -
 256         (QMF_RE(buffer[len+offset-2][bd]) * QMF_RE(buffer[len+offset-2][bd]) + QMF_IM(buffer[len+offset-2][bd]) * QMF_IM(buffer[len+offset-2][bd])) +
 257         (QMF_RE(buffer[offset-2][bd]) * QMF_RE(buffer[offset-2][bd]) + QMF_IM(buffer[offset-2][bd]) * QMF_IM(buffer[offset-2][bd]));
 258
 259     ac->det = RE(ac->r11) * RE(ac->r22) - rel * (RE(ac->r12) * RE(ac->r12) + IM(ac->r12) * IM(ac->r12));
 260 }
 261 #endif
 262
 263 /* calculate linear prediction coefficients using the covariance method */
 264 static void calc_prediction_coef(sbr_info *sbr, qmf_t Xlow[MAX_NTSRHFG][32],
 265                                  complex_t *alpha_0, complex_t *alpha_1
 266 #ifdef SBR_LOW_POWER
 267                                  , real_t *rxx
 268 #endif
 269                                  )
 270 {
 271     uint8_t k;
 272     real_t tmp;
 273     acorr_coef ac;
 274
 275     for (k = 1; k < sbr->f_master[0]; k++)
 276     {
 277         auto_correlation(sbr, &ac, Xlow, k, sbr->numTimeSlotsRate + 6);
 278
 279 #ifdef SBR_LOW_POWER
 280         if (ac.det == 0)
 281         {
 282             RE(alpha_1[k]) = 0;
 283         } else {
 284             tmp = MUL_R(RE(ac.r01), RE(ac.r12)) - MUL_R(RE(ac.r02), RE(ac.r11));
 285             RE(alpha_1[k]) = SBR_DIV(tmp, ac.det);
 286         }
 287
 288         if (RE(ac.r11) == 0)
 289         {
 290             RE(alpha_0[k]) = 0;
 291         } else {
 292             tmp = RE(ac.r01) + MUL_R(RE(alpha_1[k]), RE(ac.r12));
 293             RE(alpha_0[k]) = -SBR_DIV(tmp, RE(ac.r11));
 294         }
 295
 296         if ((RE(alpha_0[k]) >= REAL_CONST(4)) || (RE(alpha_1[k]) >= REAL_CONST(4)))
 297         {
 298             RE(alpha_0[k]) = REAL_CONST(0);
 299             RE(alpha_1[k]) = REAL_CONST(0);
 300         }
 301
 302         /* reflection coefficient */
 303         if (RE(ac.r11) == 0)
 304         {
 305             rxx[k] = REAL_CONST(0.0);
 306         } else {
 307             rxx[k] = -SBR_DIV(RE(ac.r01), RE(ac.r11));
 308             if (rxx[k] > REAL_CONST(1.0)) rxx[k] = REAL_CONST(1.0);
 309             if (rxx[k] < REAL_CONST(-1.0)) rxx[k] = REAL_CONST(-1.0);
 310         }
 311 #else
 312         if (ac.det == 0)
 313         {
 314             RE(alpha_1[k]) = 0;
 315             IM(alpha_1[k]) = 0;
 316         } else {
 317             tmp = REAL_CONST(1.0) / ac.det;
 318             RE(alpha_1[k]) = (RE(ac.r01) * RE(ac.r12) - IM(ac.r01) * IM(ac.r12) - RE(ac.r02) * RE(ac.r11)) * tmp;
 319             IM(alpha_1[k]) = (IM(ac.r01) * RE(ac.r12) + RE(ac.r01) * IM(ac.r12) - IM(ac.r02) * RE(ac.r11)) * tmp;
 320         }
 321
 322         if (RE(ac.r11) == 0)
 323         {
 324             RE(alpha_0[k]) = 0;
 325             IM(alpha_0[k]) = 0;
 326         } else {
 327             tmp = 1.0f / RE(ac.r11);
 328             RE(alpha_0[k]) = -(RE(ac.r01) + RE(alpha_1[k]) * RE(ac.r12) + IM(alpha_1[k]) * IM(ac.r12)) * tmp;
 329             IM(alpha_0[k]) = -(IM(ac.r01) + IM(alpha_1[k]) * RE(ac.r12) - RE(alpha_1[k]) * IM(ac.r12)) * tmp;
 330         }
 331
 332         if ((RE(alpha_0[k])*RE(alpha_0[k]) + IM(alpha_0[k])*IM(alpha_0[k]) >= 16) ||
 333             (RE(alpha_1[k])*RE(alpha_1[k]) + IM(alpha_1[k])*IM(alpha_1[k]) >= 16))
 334         {
 335             RE(alpha_0[k]) = 0;
 336             IM(alpha_0[k]) = 0;
 337             RE(alpha_1[k]) = 0;
 338             IM(alpha_1[k]) = 0;
 339         }
 340 #endif
 341     }
 342 }
 343
 344 #ifdef SBR_LOW_POWER
 345 static void calc_aliasing_degree(sbr_info *sbr, real_t *rxx, real_t *deg)
 346 {
 347     uint8_t k;
 348
 349     rxx[0] = REAL_CONST(0.0);
 350     deg[1] = REAL_CONST(0.0);
 351
 352     for (k = 2; k < sbr->k0; k++)
 353     {
 354         deg[k] = 0.0;
 355
 356         if ((k % 2 == 0) && (rxx[k] < REAL_CONST(0.0)))
 357         {
 358             if (rxx[k-1] < 0.0)
 359             {
 360                 deg[k] = REAL_CONST(1.0);
 361
 362                 if (rxx[k-2] > REAL_CONST(0.0))
 363                 {
 364                     deg[k-1] = REAL_CONST(1.0) - MUL_R(rxx[k-1], rxx[k-1]);
 365                 }
 366             } else if (rxx[k-2] > REAL_CONST(0.0)) {
 367                 deg[k]   = REAL_CONST(1.0) - MUL_R(rxx[k-1], rxx[k-1]);
 368             }
 369         }
 370
 371         if ((k % 2 == 1) && (rxx[k] > REAL_CONST(0.0)))
 372         {
 373             if (rxx[k-1] > REAL_CONST(0.0))
 374             {
 375                 deg[k] = REAL_CONST(1.0);
 376
 377                 if (rxx[k-2] < REAL_CONST(0.0))
 378                 {
 379                     deg[k-1] = REAL_CONST(1.0) - MUL_R(rxx[k-1], rxx[k-1]);
 380                 }
 381             } else if (rxx[k-2] < REAL_CONST(0.0)) {
 382                 deg[k] = REAL_CONST(1.0) - MUL_R(rxx[k-1], rxx[k-1]);
 383             }
 384         }
 385     }
 386 }
 387 #endif
 388
 389 /* FIXED POINT: bwArray = COEF */
 390 static real_t mapNewBw(uint8_t invf_mode, uint8_t invf_mode_prev)
 391 {
 392     switch (invf_mode)
 393     {
 394     case 1: /* LOW */
 395         if (invf_mode_prev == 0) /* NONE */
 396             return COEF_CONST(0.6);
 397         else
 398             return COEF_CONST(0.75);
 399
 400     case 2: /* MID */
 401         return COEF_CONST(0.9);
 402
 403     case 3: /* HIGH */
 404         return COEF_CONST(0.98);
 405
 406     default: /* NONE */
 407         if (invf_mode_prev == 1) /* LOW */
 408             return COEF_CONST(0.6);
 409         else
 410             return COEF_CONST(0.0);
 411     }
 412 }
 413
 414 /* FIXED POINT: bwArray = COEF */
 415 static void calc_chirp_factors(sbr_info *sbr, uint8_t ch)
 416 {
 417     uint8_t i;
 418
 419     for (i = 0; i < sbr->N_Q; i++)
 420     {
 421         sbr->bwArray[ch][i] = mapNewBw(sbr->bs_invf_mode[ch][i], sbr->bs_invf_mode_prev[ch][i]);
 422
 423         if (sbr->bwArray[ch][i] < sbr->bwArray_prev[ch][i])
 424             sbr->bwArray[ch][i] = MUL_F(sbr->bwArray[ch][i], FRAC_CONST(0.75)) + MUL_F(sbr->bwArray_prev[ch][i], FRAC_CONST(0.25));
 425         else
 426             sbr->bwArray[ch][i] = MUL_F(sbr->bwArray[ch][i], FRAC_CONST(0.90625)) + MUL_F(sbr->bwArray_prev[ch][i], FRAC_CONST(0.09375));
 427
 428         if (sbr->bwArray[ch][i] < COEF_CONST(0.015625))
 429             sbr->bwArray[ch][i] = COEF_CONST(0.0);
 430
 431         if (sbr->bwArray[ch][i] >= COEF_CONST(0.99609375))
 432             sbr->bwArray[ch][i] = COEF_CONST(0.99609375);
 433
 434         sbr->bwArray_prev[ch][i] = sbr->bwArray[ch][i];
 435         sbr->bs_invf_mode_prev[ch][i] = sbr->bs_invf_mode[ch][i];
 436     }
 437 }
 438
 439 static void patch_construction(sbr_info *sbr)
 440 {
 441     uint8_t i, k;
 442     uint8_t odd, sb;
 443     uint8_t msb = sbr->k0;
 444     uint8_t usb = sbr->kx;
 445     uint8_t goalSbTab[] = { 21, 23, 43, 46, 64, 85, 93, 128, 0, 0, 0 };
 446     /* (uint8_t)(2.048e6/sbr->sample_rate + 0.5); */
 447     uint8_t goalSb = goalSbTab[get_sr_index(sbr->sample_rate)];
 448
 449     sbr->noPatches = 0;
 450
 451     if (goalSb < (sbr->kx + sbr->M))
 452     {
 453         for (i = 0, k = 0; sbr->f_master[i] < goalSb; i++)
 454             k = i+1;
 455     } else {
 456         k = sbr->N_master;
 457     }
 458
 459     do
 460     {
 461         uint8_t j = k + 1;
 462
 463         do
 464         {
 465             j--;
 466
 467             sb = sbr->f_master[j];
 468             odd = (sb - 2 + sbr->k0) % 2;
 469         } while (sb > (sbr->k0 - 1 + msb - odd));
 470
 471         sbr->patchNoSubbands[sbr->noPatches] = max(sb - usb, 0);
 472         sbr->patchStartSubband[sbr->noPatches] = sbr->k0 - odd -
 473             sbr->patchNoSubbands[sbr->noPatches];
 474
 475         if (sbr->patchNoSubbands[sbr->noPatches] > 0)
 476         {
 477             usb = sb;
 478             msb = sb;
 479             sbr->noPatches++;
 480         } else {
 481             msb = sbr->kx;
 482         }
 483
 484         if (sbr->f_master[k] - sb < 3)
 485             k = sbr->N_master;
 486     } while (sb != (sbr->kx + sbr->M));
 487
 488     if ((sbr->patchNoSubbands[sbr->noPatches-1] < 3) && (sbr->noPatches > 1))
 489     {
 490         sbr->noPatches--;
 491     }
 492
 493     sbr->noPatches = min(sbr->noPatches, 5);
 494 }
 495
 496 #endif