src/lib/doslib/ext/faad/mdct.c

   1 /*
   2 ** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding
   3 ** Copyright (C) 2003-2005 M. Bakker, Nero 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 ** The "appropriate copyright message" mentioned in section 2c of the GPLv2
  23 ** must read: "Code from FAAD2 is copyright (c) Nero AG, www.nero.com"
  24 **
  25 ** Commercial non-GPL licensing of this software is possible.
  26 ** For more info contact Nero AG through Mpeg4AAClicense@nero.com.
  27 **
  28 ** $Id: mdct.c,v 1.47 2007/11/01 12:33:31 menno Exp $
  29 **/
  30
  31 /*
  32  * Fast (I)MDCT Implementation using (I)FFT ((Inverse) Fast Fourier Transform)
  33  * and consists of three steps: pre-(I)FFT complex multiplication, complex
  34  * (I)FFT, post-(I)FFT complex multiplication,
  35  *
  36  * As described in:
  37  *  P. Duhamel, Y. Mahieux, and J.P. Petit, "A Fast Algorithm for the
  38  *  Implementation of Filter Banks Based on 'Time Domain Aliasing
  39  *  Cancellation\92," IEEE Proc. on ICASSP\9191, 1991, pp. 2209-2212.
  40  *
  41  *
  42  * As of April 6th 2002 completely rewritten.
  43  * This (I)MDCT can now be used for any data size n, where n is divisible by 8.
  44  *
  45  */
  46
  47 #include "common.h"
  48 #include "structs.h"
  49
  50 #include <stdlib.h>
  51 #ifdef _WIN32_WCE
  52 #define assert(x)
  53 #else
  54 #include <assert.h>
  55 #endif
  56
  57 #include "cfft.h"
  58 #include "mdct.h"
  59 #include "mdct_tab.h"
  60
  61
  62 mdct_info *faad_mdct_init(uint16_t N)
  63 {
  64     mdct_info *mdct = (mdct_info*)faad_malloc(sizeof(mdct_info));
  65
  66     assert(N % 8 == 0);
  67
  68     mdct->N = N;
  69
  70     /* NOTE: For "small framelengths" in FIXED_POINT the coefficients need to be
  71      * scaled by sqrt("(nearest power of 2) > N" / N) */
  72
  73     /* RE(mdct->sincos[k]) = scale*(real_t)(cos(2.0*M_PI*(k+1./8.) / (real_t)N));
  74      * IM(mdct->sincos[k]) = scale*(real_t)(sin(2.0*M_PI*(k+1./8.) / (real_t)N)); */
  75     /* scale is 1 for fixed point, sqrt(N) for floating point */
  76     switch (N)
  77     {
  78     case 2048: mdct->sincos = (complex_t*)mdct_tab_2048; break;
  79     case 256:  mdct->sincos = (complex_t*)mdct_tab_256;  break;
  80 #ifdef LD_DEC
  81     case 1024: mdct->sincos = (complex_t*)mdct_tab_1024; break;
  82 #endif
  83 #ifdef ALLOW_SMALL_FRAMELENGTH
  84     case 1920: mdct->sincos = (complex_t*)mdct_tab_1920; break;
  85     case 240:  mdct->sincos = (complex_t*)mdct_tab_240;  break;
  86 #ifdef LD_DEC
  87     case 960:  mdct->sincos = (complex_t*)mdct_tab_960;  break;
  88 #endif
  89 #endif
  90 #ifdef SSR_DEC
  91     case 512:  mdct->sincos = (complex_t*)mdct_tab_512;  break;
  92     case 64:   mdct->sincos = (complex_t*)mdct_tab_64;   break;
  93 #endif
  94     }
  95
  96     /* initialise fft */
  97     mdct->cfft = cffti(N/4);
  98
  99 #ifdef PROFILE
 100     mdct->cycles = 0;
 101     mdct->fft_cycles = 0;
 102 #endif
 103
 104     return mdct;
 105 }
 106
 107 void faad_mdct_end(mdct_info *mdct)
 108 {
 109     if (mdct != NULL)
 110     {
 111 #ifdef PROFILE
 112         printf("MDCT[%.4d]:         %I64d cycles\n", mdct->N, mdct->cycles);
 113         printf("CFFT[%.4d]:         %I64d cycles\n", mdct->N/4, mdct->fft_cycles);
 114 #endif
 115
 116         cfftu(mdct->cfft);
 117
 118         faad_free(mdct);
 119     }
 120 }
 121
 122 void faad_imdct(mdct_info *mdct, real_t *X_in, real_t *X_out)
 123 {
 124     uint16_t k;
 125
 126     complex_t x;
 127 #ifdef ALLOW_SMALL_FRAMELENGTH
 128 #ifdef FIXED_POINT
 129     real_t scale, b_scale = 0;
 130 #endif
 131 #endif
 132     ALIGN complex_t Z1[512];
 133     complex_t *sincos = mdct->sincos;
 134
 135     uint16_t N  = mdct->N;
 136     uint16_t N2 = N >> 1;
 137     uint16_t N4 = N >> 2;
 138     uint16_t N8 = N >> 3;
 139
 140 #ifdef PROFILE
 141     int64_t count1, count2 = faad_get_ts();
 142 #endif
 143
 144 #ifdef ALLOW_SMALL_FRAMELENGTH
 145 #ifdef FIXED_POINT
 146     /* detect non-power of 2 */
 147     if (N & (N-1))
 148     {
 149         /* adjust scale for non-power of 2 MDCT */
 150         /* 2048/1920 */
 151         b_scale = 1;
 152         scale = COEF_CONST(1.0666666666666667);
 153     }
 154 #endif
 155 #endif
 156
 157     /* pre-IFFT complex multiplication */
 158     for (k = 0; k < N4; k++)
 159     {
 160         ComplexMult(&IM(Z1[k]), &RE(Z1[k]),
 161             X_in[2*k], X_in[N2 - 1 - 2*k], RE(sincos[k]), IM(sincos[k]));
 162     }
 163
 164 #ifdef PROFILE
 165     count1 = faad_get_ts();
 166 #endif
 167
 168     /* complex IFFT, any non-scaling FFT can be used here */
 169     cfftb(mdct->cfft, Z1);
 170
 171 #ifdef PROFILE
 172     count1 = faad_get_ts() - count1;
 173 #endif
 174
 175     /* post-IFFT complex multiplication */
 176     for (k = 0; k < N4; k++)
 177     {
 178         RE(x) = RE(Z1[k]);
 179         IM(x) = IM(Z1[k]);
 180         ComplexMult(&IM(Z1[k]), &RE(Z1[k]),
 181             IM(x), RE(x), RE(sincos[k]), IM(sincos[k]));
 182
 183 #ifdef ALLOW_SMALL_FRAMELENGTH
 184 #ifdef FIXED_POINT
 185         /* non-power of 2 MDCT scaling */
 186         if (b_scale)
 187         {
 188             RE(Z1[k]) = MUL_C(RE(Z1[k]), scale);
 189             IM(Z1[k]) = MUL_C(IM(Z1[k]), scale);
 190         }
 191 #endif
 192 #endif
 193     }
 194
 195     /* reordering */
 196     for (k = 0; k < N8; k+=2)
 197     {
 198         X_out[              2*k] =  IM(Z1[N8 +     k]);
 199         X_out[          2 + 2*k] =  IM(Z1[N8 + 1 + k]);
 200
 201         X_out[          1 + 2*k] = -RE(Z1[N8 - 1 - k]);
 202         X_out[          3 + 2*k] = -RE(Z1[N8 - 2 - k]);
 203
 204         X_out[N4 +          2*k] =  RE(Z1[         k]);
 205         X_out[N4 +    + 2 + 2*k] =  RE(Z1[     1 + k]);
 206
 207         X_out[N4 +      1 + 2*k] = -IM(Z1[N4 - 1 - k]);
 208         X_out[N4 +      3 + 2*k] = -IM(Z1[N4 - 2 - k]);
 209
 210         X_out[N2 +          2*k] =  RE(Z1[N8 +     k]);
 211         X_out[N2 +    + 2 + 2*k] =  RE(Z1[N8 + 1 + k]);
 212
 213         X_out[N2 +      1 + 2*k] = -IM(Z1[N8 - 1 - k]);
 214         X_out[N2 +      3 + 2*k] = -IM(Z1[N8 - 2 - k]);
 215
 216         X_out[N2 + N4 +     2*k] = -IM(Z1[         k]);
 217         X_out[N2 + N4 + 2 + 2*k] = -IM(Z1[     1 + k]);
 218
 219         X_out[N2 + N4 + 1 + 2*k] =  RE(Z1[N4 - 1 - k]);
 220         X_out[N2 + N4 + 3 + 2*k] =  RE(Z1[N4 - 2 - k]);
 221     }
 222
 223 #ifdef PROFILE
 224     count2 = faad_get_ts() - count2;
 225     mdct->fft_cycles += count1;
 226     mdct->cycles += (count2 - count1);
 227 #endif
 228 }
 229
 230 #ifdef LTP_DEC
 231 void faad_mdct(mdct_info *mdct, real_t *X_in, real_t *X_out)
 232 {
 233     uint16_t k;
 234
 235     complex_t x;
 236     ALIGN complex_t Z1[512];
 237     complex_t *sincos = mdct->sincos;
 238
 239     uint16_t N  = mdct->N;
 240     uint16_t N2 = N >> 1;
 241     uint16_t N4 = N >> 2;
 242     uint16_t N8 = N >> 3;
 243
 244 #ifndef FIXED_POINT
 245         real_t scale = REAL_CONST(N);
 246 #else
 247         real_t scale = REAL_CONST(4.0/N);
 248 #endif
 249
 250 #ifdef ALLOW_SMALL_FRAMELENGTH
 251 #ifdef FIXED_POINT
 252     /* detect non-power of 2 */
 253     if (N & (N-1))
 254     {
 255         /* adjust scale for non-power of 2 MDCT */
 256         /* *= sqrt(2048/1920) */
 257         scale = MUL_C(scale, COEF_CONST(1.0327955589886444));
 258     }
 259 #endif
 260 #endif
 261
 262     /* pre-FFT complex multiplication */
 263     for (k = 0; k < N8; k++)
 264     {
 265         uint16_t n = k << 1;
 266         RE(x) = X_in[N - N4 - 1 - n] + X_in[N - N4 +     n];
 267         IM(x) = X_in[    N4 +     n] - X_in[    N4 - 1 - n];
 268
 269         ComplexMult(&RE(Z1[k]), &IM(Z1[k]),
 270             RE(x), IM(x), RE(sincos[k]), IM(sincos[k]));
 271
 272         RE(Z1[k]) = MUL_R(RE(Z1[k]), scale);
 273         IM(Z1[k]) = MUL_R(IM(Z1[k]), scale);
 274
 275         RE(x) =  X_in[N2 - 1 - n] - X_in[        n];
 276         IM(x) =  X_in[N2 +     n] + X_in[N - 1 - n];
 277
 278         ComplexMult(&RE(Z1[k + N8]), &IM(Z1[k + N8]),
 279             RE(x), IM(x), RE(sincos[k + N8]), IM(sincos[k + N8]));
 280
 281         RE(Z1[k + N8]) = MUL_R(RE(Z1[k + N8]), scale);
 282         IM(Z1[k + N8]) = MUL_R(IM(Z1[k + N8]), scale);
 283     }
 284
 285     /* complex FFT, any non-scaling FFT can be used here  */
 286     cfftf(mdct->cfft, Z1);
 287
 288     /* post-FFT complex multiplication */
 289     for (k = 0; k < N4; k++)
 290     {
 291         uint16_t n = k << 1;
 292         ComplexMult(&RE(x), &IM(x),
 293             RE(Z1[k]), IM(Z1[k]), RE(sincos[k]), IM(sincos[k]));
 294
 295         X_out[         n] = -RE(x);
 296         X_out[N2 - 1 - n] =  IM(x);
 297         X_out[N2 +     n] = -IM(x);
 298         X_out[N  - 1 - n] =  RE(x);
 299     }
 300 }
 301 #endif