- Main Page
- Namespaces
- Classes
- Files
- Directories

- File List
- File Members

# mdct.c

    00001 /*
    00002 ** FAAD2 - Freeware Advanced Audio (AAC) Decoder including SBR decoding
    00003 ** Copyright (C) 2003-2005 M. Bakker, Ahead Software AG, http://www.nero.com
    00004 **
    00005 ** This program is free software; you can redistribute it and/or modify
    00006 ** it under the terms of the GNU General Public License as published by
    00007 ** the Free Software Foundation; either version 2 of the License, or
    00008 ** (at your option) any later version.
    00009 **
    00010 ** This program is distributed in the hope that it will be useful,
    00011 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
    00012 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    00013 ** GNU General Public License for more details.
    00014 **
    00015 ** You should have received a copy of the GNU General Public License
    00016 ** along with this program; if not, write to the Free Software
    00017 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
    00018 **
    00019 ** Any non-GPL usage of this software or parts of this software is strictly
    00020 ** forbidden.
    00021 **
    00022 ** Software using this code must display the following message visibly in the
    00023 ** software:
    00024 ** "FAAD2 AAC/HE-AAC/HE-AACv2/DRM decoder (c) Ahead Software, www.nero.com"
    00025 ** in, for example, the about-box or help/startup screen.
    00026 **
    00027 ** Commercial non-GPL licensing of this software is possible.
    00028 ** For more info contact Ahead Software through Mpeg4AAClicense@nero.com.
    00029 **
    00030 ** $Id: mdct.c,v 1.2 2005/11/01 21:41:43 gabest Exp $
    00031 **/
    00032
    00033 /*
    00034  * Fast (I)MDCT Implementation using (I)FFT ((Inverse) Fast Fourier Transform)
    00035  * and consists of three steps: pre-(I)FFT complex multiplication, complex
    00036  * (I)FFT, post-(I)FFT complex multiplication,
    00037  *
    00038  * As described in:
    00039  *  P. Duhamel, Y. Mahieux, and J.P. Petit, "A Fast Algorithm for the
    00040  *  Implementation of Filter Banks Based on 'Time Domain Aliasing
    00041  *  Cancellation’," IEEE Proc. on ICASSP‘91, 1991, pp. 2209-2212.
    00042  *
    00043  *
    00044  * As of April 6th 2002 completely rewritten.
    00045  * This (I)MDCT can now be used for any data size n, where n is divisible by 8.
    00046  *
    00047  */
    00048
    00049 #include "common.h"
    00050 #include "structs.h"
    00051
    00052 #include <stdlib.h>
    00053 #ifdef _WIN32_WCE
    00054 #define assert(x)
    00055 #else
    00056 #include <assert.h>
    00057 #endif
    00058
    00059 #include "cfft.h"
    00060 #include "mdct.h"
    00061 #include "mdct_tab.h"
    00062
    00063
    00064 mdct_info *faad_mdct_init(uint16_t N)
    00065 {
    00066     mdct_info *mdct = (mdct_info*)faad_malloc(sizeof(mdct_info));
    00067
    00068     assert(N % 8 == 0);
    00069
    00070     mdct->N = N;
    00071
    00072     /* NOTE: For "small framelengths" in FIXED_POINT the coefficients need to be
    00073      * scaled by sqrt("(nearest power of 2) > N" / N) */
    00074
    00075     /* RE(mdct->sincos[k]) = scale*(real_t)(cos(2.0*M_PI*(k+1./8.) / (real_t)N));
    00076      * IM(mdct->sincos[k]) = scale*(real_t)(sin(2.0*M_PI*(k+1./8.) / (real_t)N)); */
    00077     /* scale is 1 for fixed point, sqrt(N) for floating point */
    00078     switch (N)
    00079     {
    00080     case 2048: mdct->sincos = (complex_t*)mdct_tab_2048; break;
    00081     case 256:  mdct->sincos = (complex_t*)mdct_tab_256;  break;
    00082 #ifdef LD_DEC
    00083     case 1024: mdct->sincos = (complex_t*)mdct_tab_1024; break;
    00084 #endif
    00085 #ifdef ALLOW_SMALL_FRAMELENGTH
    00086     case 1920: mdct->sincos = (complex_t*)mdct_tab_1920; break;
    00087     case 240:  mdct->sincos = (complex_t*)mdct_tab_240;  break;
    00088 #ifdef LD_DEC
    00089     case 960:  mdct->sincos = (complex_t*)mdct_tab_960;  break;
    00090 #endif
    00091 #endif
    00092 #ifdef SSR_DEC
    00093     case 512:  mdct->sincos = (complex_t*)mdct_tab_512;  break;
    00094     case 64:   mdct->sincos = (complex_t*)mdct_tab_64;   break;
    00095 #endif
    00096     }
    00097
    00098     /* initialise fft */
    00099     mdct->cfft = cffti(N/4);
    00100
    00101 #ifdef PROFILE
    00102     mdct->cycles = 0;
    00103     mdct->fft_cycles = 0;
    00104 #endif
    00105
    00106     return mdct;
    00107 }
    00108
    00109 void faad_mdct_end(mdct_info *mdct)
    00110 {
    00111     if (mdct != NULL)
    00112     {
    00113 #ifdef PROFILE
    00114         printf("MDCT[%.4d]:         %I64d cycles\n", mdct->N, mdct->cycles);
    00115         printf("CFFT[%.4d]:         %I64d cycles\n", mdct->N/4, mdct->fft_cycles);
    00116 #endif
    00117
    00118         cfftu(mdct->cfft);
    00119
    00120         faad_free(mdct);
    00121     }
    00122 }
    00123
    00124 void faad_imdct(mdct_info *mdct, real_t *X_in, real_t *X_out)
    00125 {
    00126     uint16_t k;
    00127
    00128     complex_t x;
    00129 #ifdef ALLOW_SMALL_FRAMELENGTH
    00130 #ifdef FIXED_POINT
    00131     real_t scale, b_scale = 0;
    00132 #endif
    00133 #endif
    00134     ALIGN complex_t Z1[512];
    00135     complex_t *sincos = mdct->sincos;
    00136
    00137     uint16_t N  = mdct->N;
    00138     uint16_t N2 = N >> 1;
    00139     uint16_t N4 = N >> 2;
    00140     uint16_t N8 = N >> 3;
    00141
    00142 #ifdef PROFILE
    00143     int64_t count1, count2 = faad_get_ts();
    00144 #endif
    00145
    00146 #ifdef ALLOW_SMALL_FRAMELENGTH
    00147 #ifdef FIXED_POINT
    00148     /* detect non-power of 2 */
    00149     if (N & (N-1))
    00150     {
    00151         /* adjust scale for non-power of 2 MDCT */
    00152         /* 2048/1920 */
    00153         b_scale = 1;
    00154         scale = COEF_CONST(1.0666666666666667);
    00155     }
    00156 #endif
    00157 #endif
    00158
    00159     /* pre-IFFT complex multiplication */
    00160     for (k = 0; k < N4; k++)
    00161     {
    00162         ComplexMult(&IM(Z1[k]), &RE(Z1[k]),
    00163             X_in[2*k], X_in[N2 - 1 - 2*k], RE(sincos[k]), IM(sincos[k]));
    00164     }
    00165
    00166 #ifdef PROFILE
    00167     count1 = faad_get_ts();
    00168 #endif
    00169
    00170     /* complex IFFT, any non-scaling FFT can be used here */
    00171     cfftb(mdct->cfft, Z1);
    00172
    00173 #ifdef PROFILE
    00174     count1 = faad_get_ts() - count1;
    00175 #endif
    00176
    00177     /* post-IFFT complex multiplication */
    00178     for (k = 0; k < N4; k++)
    00179     {
    00180         RE(x) = RE(Z1[k]);
    00181         IM(x) = IM(Z1[k]);
    00182         ComplexMult(&IM(Z1[k]), &RE(Z1[k]),
    00183             IM(x), RE(x), RE(sincos[k]), IM(sincos[k]));
    00184
    00185 #ifdef ALLOW_SMALL_FRAMELENGTH
    00186 #ifdef FIXED_POINT
    00187         /* non-power of 2 MDCT scaling */
    00188         if (b_scale)
    00189         {
    00190             RE(Z1[k]) = MUL_C(RE(Z1[k]), scale);
    00191             IM(Z1[k]) = MUL_C(IM(Z1[k]), scale);
    00192         }
    00193 #endif
    00194 #endif
    00195     }
    00196
    00197     /* reordering */
    00198     for (k = 0; k < N8; k+=2)
    00199     {
    00200         X_out[              2*k] =  IM(Z1[N8 +     k]);
    00201         X_out[          2 + 2*k] =  IM(Z1[N8 + 1 + k]);
    00202
    00203         X_out[          1 + 2*k] = -RE(Z1[N8 - 1 - k]);
    00204         X_out[          3 + 2*k] = -RE(Z1[N8 - 2 - k]);
    00205
    00206         X_out[N4 +          2*k] =  RE(Z1[         k]);
    00207         X_out[N4 +    + 2 + 2*k] =  RE(Z1[     1 + k]);
    00208
    00209         X_out[N4 +      1 + 2*k] = -IM(Z1[N4 - 1 - k]);
    00210         X_out[N4 +      3 + 2*k] = -IM(Z1[N4 - 2 - k]);
    00211
    00212         X_out[N2 +          2*k] =  RE(Z1[N8 +     k]);
    00213         X_out[N2 +    + 2 + 2*k] =  RE(Z1[N8 + 1 + k]);
    00214
    00215         X_out[N2 +      1 + 2*k] = -IM(Z1[N8 - 1 - k]);
    00216         X_out[N2 +      3 + 2*k] = -IM(Z1[N8 - 2 - k]);
    00217
    00218         X_out[N2 + N4 +     2*k] = -IM(Z1[         k]);
    00219         X_out[N2 + N4 + 2 + 2*k] = -IM(Z1[     1 + k]);
    00220
    00221         X_out[N2 + N4 + 1 + 2*k] =  RE(Z1[N4 - 1 - k]);
    00222         X_out[N2 + N4 + 3 + 2*k] =  RE(Z1[N4 - 2 - k]);
    00223     }
    00224
    00225 #ifdef PROFILE
    00226     count2 = faad_get_ts() - count2;
    00227     mdct->fft_cycles += count1;
    00228     mdct->cycles += (count2 - count1);
    00229 #endif
    00230 }
    00231
    00232 #ifdef LTP_DEC
    00233 void faad_mdct(mdct_info *mdct, real_t *X_in, real_t *X_out)
    00234 {
    00235     uint16_t k;
    00236
    00237     complex_t x;
    00238     ALIGN complex_t Z1[512];
    00239     complex_t *sincos = mdct->sincos;
    00240
    00241     uint16_t N  = mdct->N;
    00242     uint16_t N2 = N >> 1;
    00243     uint16_t N4 = N >> 2;
    00244     uint16_t N8 = N >> 3;
    00245
    00246 #ifndef FIXED_POINT
    00247         real_t scale = REAL_CONST(N);
    00248 #else
    00249         real_t scale = REAL_CONST(4.0/N);
    00250 #endif
    00251
    00252 #ifdef ALLOW_SMALL_FRAMELENGTH
    00253 #ifdef FIXED_POINT
    00254     /* detect non-power of 2 */
    00255     if (N & (N-1))
    00256     {
    00257         /* adjust scale for non-power of 2 MDCT */
    00258         /* *= sqrt(2048/1920) */
    00259         scale = MUL_C(scale, COEF_CONST(1.0327955589886444));
    00260     }
    00261 #endif
    00262 #endif
    00263
    00264     /* pre-FFT complex multiplication */
    00265     for (k = 0; k < N8; k++)
    00266     {
    00267         uint16_t n = k << 1;
    00268         RE(x) = X_in[N - N4 - 1 - n] + X_in[N - N4 +     n];
    00269         IM(x) = X_in[    N4 +     n] - X_in[    N4 - 1 - n];
    00270
    00271         ComplexMult(&RE(Z1[k]), &IM(Z1[k]),
    00272             RE(x), IM(x), RE(sincos[k]), IM(sincos[k]));
    00273
    00274         RE(Z1[k]) = MUL_R(RE(Z1[k]), scale);
    00275         IM(Z1[k]) = MUL_R(IM(Z1[k]), scale);
    00276
    00277         RE(x) =  X_in[N2 - 1 - n] - X_in[        n];
    00278         IM(x) =  X_in[N2 +     n] + X_in[N - 1 - n];
    00279
    00280         ComplexMult(&RE(Z1[k + N8]), &IM(Z1[k + N8]),
    00281             RE(x), IM(x), RE(sincos[k + N8]), IM(sincos[k + N8]));
    00282
    00283         RE(Z1[k + N8]) = MUL_R(RE(Z1[k + N8]), scale);
    00284         IM(Z1[k + N8]) = MUL_R(IM(Z1[k + N8]), scale);
    00285     }
    00286
    00287     /* complex FFT, any non-scaling FFT can be used here  */
    00288     cfftf(mdct->cfft, Z1);
    00289
    00290     /* post-FFT complex multiplication */
    00291     for (k = 0; k < N4; k++)
    00292     {
    00293         uint16_t n = k << 1;
    00294         ComplexMult(&RE(x), &IM(x),
    00295             RE(Z1[k]), IM(Z1[k]), RE(sincos[k]), IM(sincos[k]));
    00296
    00297         X_out[         n] = -RE(x);
    00298         X_out[N2 - 1 - n] =  IM(x);
    00299         X_out[N2 +     n] = -IM(x);
    00300         X_out[N  - 1 - n] =  RE(x);
    00301     }
    00302 }
    00303 #endif