From cc2160de5cc05dc3aa77f3a34358e66f6193c8c9 Mon Sep 17 00:00:00 2001 From: MichaÅ‚ CichoÅ„ Date: Tue, 26 Jun 2012 20:35:30 +0200 Subject: Add support for AAC. --- faad2/src/libfaad/mdct.c | 301 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 301 insertions(+) create mode 100644 faad2/src/libfaad/mdct.c (limited to 'faad2/src/libfaad/mdct.c') diff --git a/faad2/src/libfaad/mdct.c b/faad2/src/libfaad/mdct.c new file mode 100644 index 0000000..6c4f584 --- /dev/null +++ b/faad2/src/libfaad/mdct.c @@ -0,0 +1,301 @@ +/* +** 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" + +#include +#ifdef _WIN32_WCE +#define assert(x) +#else +#include +#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, 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_DEC +void 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 -- cgit v1.2.3