diff src/Input/aac/libfaad2/mdct.c @ 0:13389e613d67 trunk

[svn] - initial import of audacious-plugins tree (lots to do)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Input/aac/libfaad2/mdct.c	Mon Sep 18 01:11:49 2006 -0700
@@ -0,0 +1,298 @@
+/*
+** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding
+** Copyright (C) 2003-2004 M. Bakker, Ahead Software 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+**
+** Any non-GPL usage of this software or parts of this software is strictly
+** forbidden.
+**
+** Commercial non-GPL licensing of this software is possible.
+** For more info contact Ahead Software through [e-mail].
+**
+** $Id: mdct.c,v 1.43 2004/09/04 14:56:28 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 <stdlib.h>
+#ifdef _WIN32_WCE
+#define assert(x)
+#else
+#include <assert.h>
+#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