/*** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding** Copyright (C) 2003-2005 M. Bakker, Nero AG, http://www.nero.com**** This program is free software; you can redistribute it and/or modify** it under the terms of the GNU General Public License as published by** the Free Software Foundation; either version 2 of the License, or** (at your option) any later version.**** This program is distributed in the hope that it will be useful,** but WITHOUT ANY WARRANTY; without even the implied warranty of** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the** GNU General Public License for more details.**** You should have received a copy of the GNU General Public License** along with this program; if not, write to the Free Software** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.**** Any non-GPL usage of this software or parts of this software is strictly** forbidden.**** The "appropriate copyright message" mentioned in section 2c of the GPLv2** must read: "Code from FAAD2 is copyright (c) Nero AG, www.nero.com"**** Commercial non-GPL licensing of this software is possible.** For more info contact Nero AG through Mpeg4AAClicense@nero.com.**** $Id: mdct.c,v 1.47 2007/11/01 12:33:31 menno Exp $**/ /* * Fast (I)MDCT Implementation using (I)FFT ((Inverse) Fast Fourier Transform) * and consists of three steps: pre-(I)FFT complex multiplication, complex * (I)FFT, post-(I)FFT complex multiplication, * * As described in: *  P. Duhamel, Y. Mahieux, and J.P. Petit, "A Fast Algorithm for the *  Implementation of Filter Banks Based on 'Time Domain Aliasing *  Cancellation'," IEEE Proc. on ICASSP'91, 1991, pp. 2209-2212. * * * As of April 6th 2002 completely rewritten. * This (I)MDCT can now be used for any data size n, where n is divisible by 8. * */ #include "common.h"#include "structs.h" #ifndef __plan9__#include <stdlib.h>#ifdef _WIN32_WCE#define assert(x)#else#include <assert.h>#endif#endif #include "cfft.h"#include "mdct.h"#include "mdct_tab.h"  mdct_info *faad_mdct_init(uint16_t N){    mdct_info *mdct = (mdct_info*)faad_malloc(sizeof(mdct_info));     assert(N % 8 == 0);     mdct->N = N;     /* NOTE: For "small framelengths" in FIXED_POINT the coefficients need to be     * scaled by sqrt("(nearest power of 2) > N" / N) */     /* RE(mdct->sincos[k]) = scale*(real_t)(cos(2.0*M_PI*(k+1./8.) / (real_t)N));     * IM(mdct->sincos[k]) = scale*(real_t)(sin(2.0*M_PI*(k+1./8.) / (real_t)N)); */    /* scale is 1 for fixed point, sqrt(N) for floating point */    switch (N)    {    case 2048: mdct->sincos = (complex_t*)mdct_tab_2048; break;    case 256:  mdct->sincos = (complex_t*)mdct_tab_256;  break;#ifdef LD_DEC    case 1024: mdct->sincos = (complex_t*)mdct_tab_1024; break;#endif#ifdef ALLOW_SMALL_FRAMELENGTH    case 1920: mdct->sincos = (complex_t*)mdct_tab_1920; break;    case 240:  mdct->sincos = (complex_t*)mdct_tab_240;  break;#ifdef LD_DEC    case 960:  mdct->sincos = (complex_t*)mdct_tab_960;  break;#endif#endif#ifdef SSR_DEC    case 512:  mdct->sincos = (complex_t*)mdct_tab_512;  break;    case 64:   mdct->sincos = (complex_t*)mdct_tab_64;   break;#endif    }     /* initialise fft */    mdct->cfft = cffti(N/4); #ifdef PROFILE    mdct->cycles = 0;    mdct->fft_cycles = 0;#endif     return mdct;} void faad_mdct_end(mdct_info *mdct){    if (mdct != NULL)    {#ifdef PROFILE        printf("MDCT[%.4d]:         %I64d cycles\n", mdct->N, mdct->cycles);        printf("CFFT[%.4d]:         %I64d cycles\n", mdct->N/4, mdct->fft_cycles);#endif         cfftu(mdct->cfft);         faad_free(mdct);    }} void faad_imdct(mdct_info *mdct, real_t *X_in, real_t *X_out){    uint16_t k;     complex_t x;#ifdef ALLOW_SMALL_FRAMELENGTH#ifdef FIXED_POINT    real_t scale = 0, b_scale = 0;#endif#endif    ALIGN complex_t Z1[512];    complex_t *sincos = mdct->sincos;     uint16_t N  = mdct->N;    uint16_t N2 = N >> 1;    uint16_t N4 = N >> 2;    uint16_t N8 = N >> 3; #ifdef PROFILE    int64_t count1, count2 = faad_get_ts();#endif #ifdef ALLOW_SMALL_FRAMELENGTH#ifdef FIXED_POINT    /* detect non-power of 2 */    if (N & (N-1))    {        /* adjust scale for non-power of 2 MDCT */        /* 2048/1920 */        b_scale = 1;        scale = COEF_CONST(1.0666666666666667);    }#endif#endif     /* pre-IFFT complex multiplication */    for (k = 0; k < N4; k++)    {        ComplexMult(&IM(Z1[k]), &RE(Z1[k]),            X_in[2*k], X_in[N2 - 1 - 2*k], RE(sincos[k]), IM(sincos[k]));    } #ifdef PROFILE    count1 = faad_get_ts();#endif     /* complex IFFT, any non-scaling FFT can be used here */    cfftb(mdct->cfft, Z1); #ifdef PROFILE    count1 = faad_get_ts() - count1;#endif     /* post-IFFT complex multiplication */    for (k = 0; k < N4; k++)    {        RE(x) = RE(Z1[k]);        IM(x) = IM(Z1[k]);        ComplexMult(&IM(Z1[k]), &RE(Z1[k]),            IM(x), RE(x), RE(sincos[k]), IM(sincos[k])); #ifdef ALLOW_SMALL_FRAMELENGTH#ifdef FIXED_POINT        /* non-power of 2 MDCT scaling */        if (b_scale)        {            RE(Z1[k]) = MUL_C(RE(Z1[k]), scale);            IM(Z1[k]) = MUL_C(IM(Z1[k]), scale);        }#endif#endif    }     /* reordering */    for (k = 0; k < N8; k+=2)    {        X_out[              2*k] =  IM(Z1[N8 +     k]);        X_out[          2 + 2*k] =  IM(Z1[N8 + 1 + k]);         X_out[          1 + 2*k] = -RE(Z1[N8 - 1 - k]);        X_out[          3 + 2*k] = -RE(Z1[N8 - 2 - k]);         X_out[N4 +          2*k] =  RE(Z1[         k]);        X_out[N4 +    + 2 + 2*k] =  RE(Z1[     1 + k]);         X_out[N4 +      1 + 2*k] = -IM(Z1[N4 - 1 - k]);        X_out[N4 +      3 + 2*k] = -IM(Z1[N4 - 2 - k]);         X_out[N2 +          2*k] =  RE(Z1[N8 +     k]);        X_out[N2 +    + 2 + 2*k] =  RE(Z1[N8 + 1 + k]);         X_out[N2 +      1 + 2*k] = -IM(Z1[N8 - 1 - k]);        X_out[N2 +      3 + 2*k] = -IM(Z1[N8 - 2 - k]);         X_out[N2 + N4 +     2*k] = -IM(Z1[         k]);        X_out[N2 + N4 + 2 + 2*k] = -IM(Z1[     1 + k]);         X_out[N2 + N4 + 1 + 2*k] =  RE(Z1[N4 - 1 - k]);        X_out[N2 + N4 + 3 + 2*k] =  RE(Z1[N4 - 2 - k]);    } #ifdef PROFILE    count2 = faad_get_ts() - count2;    mdct->fft_cycles += count1;    mdct->cycles += (count2 - count1);#endif} #ifdef LTP_DECvoid faad_mdct(mdct_info *mdct, real_t *X_in, real_t *X_out){    uint16_t k;     complex_t x;    ALIGN complex_t Z1[512];    complex_t *sincos = mdct->sincos;     uint16_t N  = mdct->N;    uint16_t N2 = N >> 1;    uint16_t N4 = N >> 2;    uint16_t N8 = N >> 3; #ifndef FIXED_POINT	real_t scale = REAL_CONST(N);#else	real_t scale = REAL_CONST(4.0/N);#endif #ifdef ALLOW_SMALL_FRAMELENGTH#ifdef FIXED_POINT    /* detect non-power of 2 */    if (N & (N-1))    {        /* adjust scale for non-power of 2 MDCT */        /* *= sqrt(2048/1920) */        scale = MUL_C(scale, COEF_CONST(1.0327955589886444));    }#endif#endif     /* pre-FFT complex multiplication */    for (k = 0; k < N8; k++)    {        uint16_t n = k << 1;        RE(x) = X_in[N - N4 - 1 - n] + X_in[N - N4 +     n];        IM(x) = X_in[    N4 +     n] - X_in[    N4 - 1 - n];         ComplexMult(&RE(Z1[k]), &IM(Z1[k]),            RE(x), IM(x), RE(sincos[k]), IM(sincos[k]));         RE(Z1[k]) = MUL_R(RE(Z1[k]), scale);        IM(Z1[k]) = MUL_R(IM(Z1[k]), scale);         RE(x) =  X_in[N2 - 1 - n] - X_in[        n];        IM(x) =  X_in[N2 +     n] + X_in[N - 1 - n];         ComplexMult(&RE(Z1[k + N8]), &IM(Z1[k + N8]),            RE(x), IM(x), RE(sincos[k + N8]), IM(sincos[k + N8]));         RE(Z1[k + N8]) = MUL_R(RE(Z1[k + N8]), scale);        IM(Z1[k + N8]) = MUL_R(IM(Z1[k + N8]), scale);    }     /* complex FFT, any non-scaling FFT can be used here  */    cfftf(mdct->cfft, Z1);     /* post-FFT complex multiplication */    for (k = 0; k < N4; k++)    {        uint16_t n = k << 1;        ComplexMult(&RE(x), &IM(x),            RE(Z1[k]), IM(Z1[k]), RE(sincos[k]), IM(sincos[k]));         X_out[         n] = -RE(x);        X_out[N2 - 1 - n] =  IM(x);        X_out[N2 +     n] = -IM(x);        X_out[N  - 1 - n] =  RE(x);    }}#endif