libSBRenc/src/resampler.cpp

   1
   2 /* -----------------------------------------------------------------------------------------------------------
   3 Software License for The Fraunhofer FDK AAC Codec Library for Android
   4
   5 © Copyright  1995 - 2013 Fraunhofer-Gesellschaft zur Förderung der angewandten Forschung e.V.
   6   All rights reserved.
   7
   8  1.    INTRODUCTION
   9 The Fraunhofer FDK AAC Codec Library for Android ("FDK AAC Codec") is software that implements
  10 the MPEG Advanced Audio Coding ("AAC") encoding and decoding scheme for digital audio.
  11 This FDK AAC Codec software is intended to be used on a wide variety of Android devices.
  12
  13 AAC's HE-AAC and HE-AAC v2 versions are regarded as today's most efficient general perceptual
  14 audio codecs. AAC-ELD is considered the best-performing full-bandwidth communications codec by
  15 independent studies and is widely deployed. AAC has been standardized by ISO and IEC as part
  16 of the MPEG specifications.
  17
  18 Patent licenses for necessary patent claims for the FDK AAC Codec (including those of Fraunhofer)
  19 may be obtained through Via Licensing (www.vialicensing.com) or through the respective patent owners
  20 individually for the purpose of encoding or decoding bit streams in products that are compliant with
  21 the ISO/IEC MPEG audio standards. Please note that most manufacturers of Android devices already license
  22 these patent claims through Via Licensing or directly from the patent owners, and therefore FDK AAC Codec
  23 software may already be covered under those patent licenses when it is used for those licensed purposes only.
  24
  25 Commercially-licensed AAC software libraries, including floating-point versions with enhanced sound quality,
  26 are also available from Fraunhofer. Users are encouraged to check the Fraunhofer website for additional
  27 applications information and documentation.
  28
  29 2.    COPYRIGHT LICENSE
  30
  31 Redistribution and use in source and binary forms, with or without modification, are permitted without
  32 payment of copyright license fees provided that you satisfy the following conditions:
  33
  34 You must retain the complete text of this software license in redistributions of the FDK AAC Codec or
  35 your modifications thereto in source code form.
  36
  37 You must retain the complete text of this software license in the documentation and/or other materials
  38 provided with redistributions of the FDK AAC Codec or your modifications thereto in binary form.
  39 You must make available free of charge copies of the complete source code of the FDK AAC Codec and your
  40 modifications thereto to recipients of copies in binary form.
  41
  42 The name of Fraunhofer may not be used to endorse or promote products derived from this library without
  43 prior written permission.
  44
  45 You may not charge copyright license fees for anyone to use, copy or distribute the FDK AAC Codec
  46 software or your modifications thereto.
  47
  48 Your modified versions of the FDK AAC Codec must carry prominent notices stating that you changed the software
  49 and the date of any change. For modified versions of the FDK AAC Codec, the term
  50 "Fraunhofer FDK AAC Codec Library for Android" must be replaced by the term
  51 "Third-Party Modified Version of the Fraunhofer FDK AAC Codec Library for Android."
  52
  53 3.    NO PATENT LICENSE
  54
  55 NO EXPRESS OR IMPLIED LICENSES TO ANY PATENT CLAIMS, including without limitation the patents of Fraunhofer,
  56 ARE GRANTED BY THIS SOFTWARE LICENSE. Fraunhofer provides no warranty of patent non-infringement with
  57 respect to this software.
  58
  59 You may use this FDK AAC Codec software or modifications thereto only for purposes that are authorized
  60 by appropriate patent licenses.
  61
  62 4.    DISCLAIMER
  63
  64 This FDK AAC Codec software is provided by Fraunhofer on behalf of the copyright holders and contributors
  65 "AS IS" and WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES, including but not limited to the implied warranties
  66 of merchantability and fitness for a particular purpose. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
  67 CONTRIBUTORS BE LIABLE for any direct, indirect, incidental, special, exemplary, or consequential damages,
  68 including but not limited to procurement of substitute goods or services; loss of use, data, or profits,
  69 or business interruption, however caused and on any theory of liability, whether in contract, strict
  70 liability, or tort (including negligence), arising in any way out of the use of this software, even if
  71 advised of the possibility of such damage.
  72
  73 5.    CONTACT INFORMATION
  74
  75 Fraunhofer Institute for Integrated Circuits IIS
  76 Attention: Audio and Multimedia Departments - FDK AAC LL
  77 Am Wolfsmantel 33
  78 91058 Erlangen, Germany
  79
  80 www.iis.fraunhofer.de/amm
  81 amm-info@iis.fraunhofer.de
  82 ----------------------------------------------------------------------------------------------------------- */
  83
  84 /*!
  85   \file
  86   \brief  FDK resampler tool box:
  87   \author M. Werner
  88 */
  89
  90 #include "resampler.h"
  91
  92 #include "genericStds.h"
  93
  94
  95 /**************************************************************************/
  96 /*                   BIQUAD Filter Specifications                         */
  97 /**************************************************************************/
  98
  99 #define B1 0
 100 #define B2 1
 101 #define A1 2
 102 #define A2 3
 103
 104 #define BQC(x) FL2FXCONST_SGL(x/2)
 105
 106
 107 struct FILTER_PARAM {
 108   const FIXP_SGL *coeffa;    /*! SOS matrix One row/section. Scaled using BQC(). Order of coefficients: B1,B2,A1,A2. B0=A0=1.0 */
 109   FIXP_DBL g;                /*! overall gain */
 110   int Wc;                    /*! normalized passband bandwidth at input samplerate * 1000 */
 111   int noCoeffs;              /*! number of filter coeffs */
 112   int delay;                 /*! delay in samples at input samplerate */
 113 };
 114
 115 #define BIQUAD_COEFSTEP 4
 116
 117 /**
 118  *\brief Low Pass
 119  Wc = 0,5, order 30, Stop Band -96dB. Wc criteria is "almost 0dB passband", not the usual -3db gain point.
 120  [b,a]=cheby2(30,96,0.505)
 121  [sos,g]=tf2sos(b,a)
 122  bandwidth 0.48
 123  */
 124 static const FIXP_SGL sos48[] = {
 125  BQC(1.98941075681938),   BQC(0.999999996890811), BQC(0.863264527201963),     BQC( 0.189553799960663),
 126  BQC(1.90733804822445),   BQC(1.00000001736189),  BQC(0.836321575841691),     BQC( 0.203505809266564),
 127  BQC(1.75616665495325),   BQC(0.999999946079721), BQC(0.784699225121588),     BQC( 0.230471265506986),
 128  BQC(1.55727745512726),   BQC(1.00000011737815),  BQC(0.712515423588351),     BQC( 0.268752723900498),
 129  BQC(1.33407591943643),   BQC(0.999999795953228), BQC(0.625059117330989),     BQC( 0.316194685288965),
 130  BQC(1.10689898412458),   BQC(1.00000035057114),  BQC(0.52803514366398),      BQC( 0.370517843224669),
 131  BQC(0.89060371078454),   BQC(0.999999343962822), BQC(0.426920462165257),     BQC( 0.429608200207746),
 132  BQC(0.694438261209433),  BQC( 1.0000008629792),  BQC(0.326530699561716),     BQC( 0.491714450654174),
 133  BQC(0.523237800935322),  BQC(1.00000101349782),  BQC(0.230829556274851),     BQC( 0.555559034843281),
 134  BQC(0.378631165929563),  BQC(0.99998986482665),  BQC(0.142906422036095),     BQC( 0.620338874442411),
 135  BQC(0.260786911308437),  BQC(1.00003261460178),  BQC(0.0651008576256505),    BQC( 0.685759923926262),
 136  BQC(0.168409429188098),  BQC(0.999933049695828), BQC(-0.000790067789975562), BQC( 0.751905896602325),
 137  BQC(0.100724533818628),  BQC(1.00009472669872),  BQC(-0.0533772830257041),   BQC( 0.81930744384525),
 138  BQC(0.0561434357867363), BQC(0.999911636304276), BQC(-0.0913550299236405),   BQC( 0.88883625875915),
 139  BQC(0.0341680678662057), BQC(1.00003667508676),  BQC(-0.113405185536697),    BQC( 0.961756638268446)
 140 };
 141
 142 #ifdef RS_BIQUAD_SCATTERGAIN
 143 static const FIXP_DBL g48 = FL2FXCONST_DBL(0.67436532061161992682404480717671 - 0.001);
 144 #else
 145 static const FIXP_DBL g48 = FL2FXCONST_DBL(0.002712866530047) - (FIXP_DBL)0x8000;
 146 #endif
 147
 148 static const struct FILTER_PARAM param_set48 = {
 149   sos48,
 150   g48,
 151   480,
 152   15,
 153   4 /* LF 2 */
 154 };
 155
 156 /**
 157  *\brief Low Pass
 158  Wc = 0,5, order 24, Stop Band -96dB. Wc criteria is "almost 0dB passband", not the usual -3db gain point.
 159  [b,a]=cheby2(24,96,0.5)
 160  [sos,g]=tf2sos(b,a)
 161  bandwidth 0.45
 162  */
 163 static const FIXP_SGL sos45[] = {
 164   BQC(1.982962601444),     BQC(1.00000000007504),  BQC(0.646113303737836),   BQC( 0.10851149979981),
 165   BQC(1.85334094281111),   BQC(0.999999999677192), BQC(0.612073220102006),   BQC( 0.130022141698044),
 166   BQC(1.62541051415425),   BQC(1.00000000080398),  BQC(0.547879702855959),   BQC( 0.171165825133192),
 167   BQC(1.34554656923247),   BQC(0.9999999980169),   BQC(0.460373914508491),   BQC( 0.228677463376354),
 168   BQC(1.05656568503116),   BQC(1.00000000569363),  BQC(0.357891894038287),   BQC( 0.298676843912185),
 169   BQC(0.787967587877312),  BQC(0.999999984415017), BQC(0.248826893211877),   BQC( 0.377441803512978),
 170   BQC(0.555480971120497),  BQC(1.00000003583307),  BQC(0.140614263345315),   BQC( 0.461979302213679),
 171   BQC(0.364986207070964),  BQC(0.999999932084303), BQC(0.0392669446074516),  BQC( 0.55033451180825),
 172   BQC(0.216827267631558),  BQC(1.00000010534682),  BQC(-0.0506232228865103), BQC( 0.641691581560946),
 173   BQC(0.108951672277119),  BQC(0.999999871167516), BQC(-0.125584840183225),  BQC( 0.736367748771803),
 174   BQC(0.0387988607229035), BQC(1.00000011205574),  BQC(-0.182814849097974),  BQC( 0.835802108714964),
 175   BQC(0.0042866175809225), BQC(0.999999954830813), BQC(-0.21965740617151),   BQC( 0.942623047782363)
 176 };
 177
 178 #ifdef RS_BIQUAD_SCATTERGAIN
 179 static const FIXP_DBL g45 = FL2FXCONST_DBL(0.60547428891341319051142629706723 - 0.001);
 180 #else
 181 static const FIXP_DBL g45 = FL2FXCONST_DBL(0.00242743980909524) - (FIXP_DBL)0x8000;
 182 #endif
 183
 184 static const struct FILTER_PARAM param_set45 = {
 185   sos45,
 186   g45,
 187   450,
 188   12,
 189   4 /* LF 2 */
 190 };
 191
 192 /*
 193  Created by Octave 2.1.73, Mon Oct 13 17:31:32 2008 CEST
 194  Wc = 0,5, order 16, Stop Band -96dB damping.
 195  [b,a]=cheby2(16,96,0.5)
 196  [sos,g]=tf2sos(b,a)
 197  bandwidth = 0.41
 198  */
 199
 200 static const FIXP_SGL sos41[] =
 201 {
 202   BQC(1.96193625292),       BQC(0.999999999999964), BQC(0.169266178786789),   BQC(0.0128823300475907),
 203   BQC(1.68913437662092),    BQC(1.00000000000053),  BQC(0.124751503206552),   BQC(0.0537472273950989),
 204   BQC(1.27274692366017),    BQC(0.999999999995674), BQC(0.0433108625178357),  BQC(0.131015753236317),
 205   BQC(0.85214175088395),    BQC(1.00000000001813),  BQC(-0.0625658152550408), BQC(0.237763778993806),
 206   BQC(0.503841579939009),   BQC(0.999999999953223), BQC(-0.179176128722865),  BQC(0.367475236424474),
 207   BQC(0.249990711986162),   BQC(1.00000000007952),  BQC(-0.294425165824676),  BQC(0.516594857170212),
 208   BQC(0.087971668680286),   BQC(0.999999999915528), BQC(-0.398956566777928),  BQC(0.686417767801123),
 209   BQC(0.00965373325350294), BQC(1.00000000003744),  BQC(-0.48579173764817),   BQC(0.884931534239068)
 210 };
 211
 212 #ifdef RS_BIQUAD_SCATTERGAIN
 213 static const FIXP_DBL g41 = FL2FXCONST_DBL(0.44578514476476679750811222123569);
 214 #else
 215 static const FIXP_DBL g41 = FL2FXCONST_DBL(0.00155956951169248);
 216 #endif
 217
 218 static const struct FILTER_PARAM param_set41 = {
 219   sos41,
 220   g41,
 221   410,
 222   8,
 223   5 /* LF 3 */
 224 };
 225
 226 /*
 227  # Created by Octave 2.1.73, Mon Oct 13 17:55:33 2008 CEST
 228  Wc = 0,5, order 12, Stop Band -96dB damping.
 229  [b,a]=cheby2(12,96,0.5);
 230  [sos,g]=tf2sos(b,a)
 231 */
 232 static const FIXP_SGL sos35[] =
 233 {
 234   BQC(1.93299325235762),   BQC(0.999999999999985), BQC(-0.140733187246596), BQC(0.0124139497836062),
 235   BQC(1.4890416764109),    BQC(1.00000000000011),  BQC(-0.198215402588504), BQC(0.0746730616584138),
 236   BQC(0.918450161309795),  BQC(0.999999999999619), BQC(-0.30133912791941),  BQC(0.192276468839529),
 237   BQC(0.454877024246818),  BQC(1.00000000000086),  BQC(-0.432337328809815), BQC(0.356852933642815),
 238   BQC(0.158017147118507),  BQC(0.999999999998876), BQC(-0.574817494249777), BQC(0.566380436970833),
 239   BQC(0.0171834649478749), BQC(1.00000000000055),  BQC(-0.718581178041165), BQC(0.83367484487889)
 240 };
 241
 242 #ifdef RS_BIQUAD_SCATTERGAIN
 243 static const FIXP_DBL g35 = FL2FXCONST_DBL(0.34290853574973898694521267606792);
 244 #else
 245 static const FIXP_DBL g35 = FL2FXCONST_DBL(0.00162580994125131);
 246 #endif
 247
 248 static const struct FILTER_PARAM param_set35 = {
 249   sos35,
 250   g35,
 251   350,
 252   6,
 253   4
 254 };
 255
 256 /*
 257  # Created by Octave 2.1.73, Mon Oct 13 18:15:38 2008 CEST
 258  Wc = 0,5, order 8, Stop Band -96dB damping.
 259  [b,a]=cheby2(8,96,0.5);
 260  [sos,g]=tf2sos(b,a)
 261 */
 262 static const FIXP_SGL sos25[] =
 263 {
 264   BQC(1.85334094301225),   BQC(1.0),               BQC(-0.702127214212663), BQC(0.132452403998767),
 265   BQC(1.056565682167),     BQC(0.999999999999997), BQC(-0.789503667880785), BQC(0.236328693569128),
 266   BQC(0.364986307455489),  BQC(0.999999999999996), BQC(-0.955191189843375), BQC(0.442966457936379),
 267   BQC(0.0387985751642125), BQC(1.0),               BQC(-1.19817786088084),  BQC(0.770493895456328)
 268 };
 269
 270 #ifdef RS_BIQUAD_SCATTERGAIN
 271 static const FIXP_DBL g25 = FL2FXCONST_DBL(0.17533917408936346960080259950471);
 272 #else
 273 static const FIXP_DBL g25 = FL2FXCONST_DBL(0.000945182835294559);
 274 #endif
 275
 276 static const struct FILTER_PARAM param_set25 = {
 277   sos25,
 278   g25,
 279   250,
 280   4,
 281   5
 282 };
 283
 284 /* Must be sorted in descending order */
 285 static const struct FILTER_PARAM *const filter_paramSet[] = {
 286   &param_set48,
 287   &param_set45,
 288   &param_set41,
 289   &param_set35,
 290   &param_set25
 291 };
 292
 293
 294 /**************************************************************************/
 295 /*                         Resampler Functions                            */
 296 /**************************************************************************/
 297
 298
 299 /*!
 300   \brief   Reset downsampler instance and clear delay lines
 301
 302   \return  success of operation
 303 */
 304
 305 INT FDKaacEnc_InitDownsampler(DOWNSAMPLER *DownSampler, /*!< pointer to downsampler instance */
 306                               int Wc,                   /*!< normalized cutoff freq * 1000*  */
 307                               int ratio)                /*!< downsampler ratio (only 2 supported at the momment) */
 308
 309 {
 310   UINT i;
 311   const struct FILTER_PARAM *currentSet=NULL;
 312
 313   FDK_ASSERT(ratio == 2);
 314   FDKmemclear(DownSampler->downFilter.states, sizeof(DownSampler->downFilter.states));
 315   DownSampler->downFilter.ptr   =   0;
 316
 317   /*
 318     find applicable parameter set
 319   */
 320   currentSet = filter_paramSet[0];
 321   for(i=1;i<sizeof(filter_paramSet)/sizeof(struct FILTER_PARAM *);i++){
 322     if (filter_paramSet[i]->Wc <= Wc) {
 323       break;
 324     }
 325     currentSet = filter_paramSet[i];
 326   }
 327
 328   DownSampler->downFilter.coeffa = currentSet->coeffa;
 329
 330
 331   DownSampler->downFilter.gain = currentSet->g;
 332   FDK_ASSERT(currentSet->noCoeffs <= MAXNR_SECTIONS*2);
 333
 334   DownSampler->downFilter.noCoeffs = currentSet->noCoeffs;
 335   DownSampler->delay = currentSet->delay;
 336   DownSampler->downFilter.Wc = currentSet->Wc;
 337
 338   DownSampler->ratio =   ratio;
 339   DownSampler->pending = ratio-1;
 340   return(1);
 341 }
 342
 343
 344 /*!
 345   \brief   faster simple folding operation
 346            Filter:
 347            H(z) = A(z)/B(z)
 348            with
 349            A(z) = a[0]*z^0 + a[1]*z^1 + a[2]*z^2 ... a[n]*z^n
 350
 351   \return  filtered value
 352 */
 353
 354 static inline INT_PCM AdvanceFilter(LP_FILTER *downFilter,  /*!< pointer to iir filter instance */
 355                                      INT_PCM  *pInput,          /*!< input of filter                */
 356                                      int downRatio,
 357                                      int inStride)
 358 {
 359   INT_PCM output;
 360   int i, n;
 361
 362
 363 #ifdef RS_BIQUAD_SCATTERGAIN
 364 #define BIQUAD_SCALE 3
 365 #else
 366 #define BIQUAD_SCALE 12
 367 #endif
 368
 369   FIXP_DBL y = FL2FXCONST_DBL(0.0f);
 370   FIXP_DBL input;
 371
 372   for (n=0; n<downRatio; n++)
 373   {
 374     FIXP_BQS (*states)[2] = downFilter->states;
 375     const FIXP_SGL *coeff = downFilter->coeffa;
 376     int s1,s2;
 377
 378     s1 = downFilter->ptr;
 379     s2 = s1 ^ 1;
 380
 381 #if (SAMPLE_BITS == 16)
 382     input =  ((FIXP_DBL)pInput[n*inStride]) << (DFRACT_BITS-SAMPLE_BITS-BIQUAD_SCALE);
 383 #elif (SAMPLE_BITS == 32)
 384     input =  pInput[n*inStride] >> BIQUAD_SCALE;
 385 #else
 386 #error NOT IMPLEMENTED
 387 #endif
 388
 389 #ifndef RS_BIQUAD_SCATTERGAIN /* Merged Direct form I */
 390
 391     FIXP_BQS state1, state2, state1b, state2b;
 392
 393     state1 = states[0][s1];
 394     state2 = states[0][s2];
 395
 396     /* Loop over sections */
 397     for (i=0; i<downFilter->noCoeffs; i++)
 398     {
 399       FIXP_DBL state0;
 400
 401       /* Load merged states (from next section) */
 402       state1b = states[i+1][s1];
 403       state2b = states[i+1][s2];
 404
 405       state0 = input  + fMult(state1, coeff[B1]) + fMult(state2, coeff[B2]);
 406       y      = state0 - fMult(state1b, coeff[A1]) - fMult(state2b, coeff[A2]);
 407
 408       /* Store new feed forward merge state */
 409       states[i+1][s2] = y<<1;
 410       /* Store new feed backward state */
 411       states[i][s2] = input<<1;
 412
 413       /* Feedback output to next section. */
 414       input = y;
 415
 416       /* Transfer merged states */
 417       state1 = state1b;
 418       state2 = state2b;
 419
 420       /* Step to next coef set */
 421       coeff += BIQUAD_COEFSTEP;
 422     }
 423     downFilter->ptr ^= 1;
 424   }
 425   /* Apply global gain */
 426   y = fMult(y, downFilter->gain);
 427
 428 #else /* Direct form II */
 429
 430     /* Loop over sections */
 431     for (i=0; i<downFilter->noCoeffs; i++)
 432     {
 433       FIXP_BQS state1, state2;
 434       FIXP_DBL state0;
 435
 436       /* Load states */
 437       state1 = states[i][s1];
 438       state2 = states[i][s2];
 439
 440       state0 = input  - fMult(state1, coeff[A1]) - fMult(state2, coeff[A2]);
 441       y = state0      + fMult(state1, coeff[B1]) + fMult(state2, coeff[B2]);
 442       /* Apply scattered gain */
 443       y = fMult(y, downFilter->gain);
 444
 445       /* Store new state in normalized form */
 446 #ifdef RS_BIQUAD_STATES16
 447       /* Do not saturate any state value ! The result would be unacceptable. Rounding makes SNR around 10dB better. */
 448       states[i][s2] = (FIXP_BQS)(LONG)((state0 + (FIXP_DBL)(1<<(DFRACT_BITS-FRACT_BITS-2))) >> (DFRACT_BITS-FRACT_BITS-1));
 449 #else
 450       states[i][s2] = state0<<1;
 451 #endif
 452
 453       /* Feedback output to next section. */
 454       input=y;
 455
 456       /* Step to next coef set */
 457       coeff += BIQUAD_COEFSTEP;
 458     }
 459     downFilter->ptr ^= 1;
 460   }
 461
 462 #endif
 463
 464   /* Apply final gain/scaling to output */
 465 #if (SAMPLE_BITS == 16)
 466   output = (INT_PCM) SATURATE_RIGHT_SHIFT(y+(FIXP_DBL)(1<<(DFRACT_BITS-SAMPLE_BITS-BIQUAD_SCALE-1)), DFRACT_BITS-SAMPLE_BITS-BIQUAD_SCALE, SAMPLE_BITS);
 467   //output = (INT_PCM) SATURATE_RIGHT_SHIFT(y, DFRACT_BITS-SAMPLE_BITS-BIQUAD_SCALE, SAMPLE_BITS);
 468 #else
 469   output = SATURATE_LEFT_SHIFT(y, BIQUAD_SCALE, SAMPLE_BITS);
 470 #endif
 471
 472
 473   return output;
 474 }
 475
 476
 477
 478
 479 /*!
 480   \brief   FDKaacEnc_Downsample numInSamples of type INT_PCM
 481            Returns number of output samples in numOutSamples
 482
 483   \return  success of operation
 484 */
 485
 486 INT FDKaacEnc_Downsample(DOWNSAMPLER *DownSampler,  /*!< pointer to downsampler instance */
 487                          INT_PCM *inSamples,        /*!< pointer to input samples */
 488                          INT numInSamples,          /*!< number  of input samples  */
 489                          INT inStride,              /*!< increment of input samples */
 490                          INT_PCM *outSamples,       /*!< pointer to output samples */
 491                          INT *numOutSamples,        /*!< pointer tp number of output samples */
 492                          INT outStride              /*!< increment of output samples */
 493                          )
 494 {
 495     INT i;
 496     *numOutSamples=0;
 497
 498     for(i=0; i<numInSamples; i+=DownSampler->ratio)
 499     {
 500       *outSamples = AdvanceFilter(&(DownSampler->downFilter), &inSamples[i*inStride], DownSampler->ratio, inStride);
 501       outSamples += outStride;
 502     }
 503     *numOutSamples = numInSamples/DownSampler->ratio;
 504
 505     return 0;
 506 }
 507