ref: e79f5009f833fbbe65b213c59e2e5c27e1ee1ce1
dir: /libfaad/mdct.c/
/* ** FAAD - Freeware Advanced Audio Decoder ** Copyright (C) 2002 M. Bakker ** ** 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. ** ** $Id: mdct.c,v 1.9 2002/04/23 21:08:26 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. * Thanks to the FFTW library this (I)MDCT can now be used for any data * size n, where n is divisible by 8. * */ #include "common.h" #include <stdlib.h> #include <assert.h> /* uses fftw (http://www.fftw.org) for very fast arbitrary-n FFT and IFFT */ #include <fftw.h> #include "mdct.h" void faad_mdct_init(mdct_info *mdct, uint16_t N) { uint16_t k; assert(N % 8 == 0); mdct->N = N; mdct->sincos = (faad_sincos*)malloc(N/4*sizeof(faad_sincos)); mdct->Z1 = (fftw_complex*)malloc(N/4*sizeof(fftw_complex)); mdct->Z2 = (fftw_complex*)malloc(N/4*sizeof(fftw_complex)); for (k = 0; k < N/4; k++) { real_t angle = 2.0 * M_PI * (k + 1.0/8.0)/(real_t)N; mdct->sincos[k].sin = -sin(angle); mdct->sincos[k].cos = -cos(angle); } mdct->plan_backward = fftw_create_plan(N/4, FFTW_BACKWARD, FFTW_ESTIMATE); #ifdef LTP_DEC mdct->plan_forward = fftw_create_plan(N/4, FFTW_FORWARD, FFTW_ESTIMATE); #endif } void faad_mdct_end(mdct_info *mdct) { fftw_destroy_plan(mdct->plan_backward); #ifdef LTP_DEC fftw_destroy_plan(mdct->plan_forward); #endif if (mdct->Z2) free(mdct->Z2); if (mdct->Z1) free(mdct->Z1); if (mdct->sincos) free(mdct->sincos); } void faad_imdct(mdct_info *mdct, real_t *X_in, real_t *X_out) { uint16_t k; fftw_complex *Z1 = mdct->Z1; fftw_complex *Z2 = mdct->Z2; faad_sincos *sincos = mdct->sincos; uint16_t N = mdct->N; uint16_t N2 = N >> 1; uint16_t N4 = N >> 2; uint16_t N8 = N >> 3; real_t fac = 2.0/(real_t)N; /* pre-IFFT complex multiplication */ for (k = 0; k < N4; k++) { uint16_t n = k << 1; real_t x0 = X_in[ n]; real_t x1 = X_in[N2 - 1 - n]; Z1[k].re = x1 * sincos[k].cos - x0 * sincos[k].sin; Z1[k].im = x0 * sincos[k].cos + x1 * sincos[k].sin; } /* complex IFFT */ fftw_one(mdct->plan_backward, Z1, Z2); /* post-IFFT complex multiplication */ for (k = 0; k < N4; k++) { real_t zr = Z2[k].re; real_t zi = Z2[k].im; Z2[k].re = fac * (zr * sincos[k].cos - zi * sincos[k].sin); Z2[k].im = fac * (zi * sincos[k].cos + zr * sincos[k].sin); } /* reordering */ for (k = 0; k < N8; k++) { uint16_t n = k << 1; X_out[ n] = -Z2[N8 + k].im; X_out[ 1 + n] = Z2[N8 - 1 - k].re; X_out[N4 + n] = -Z2[ k].re; X_out[N4 + 1 + n] = Z2[N4 - 1 - k].im; X_out[N2 + n] = -Z2[N8 + k].re; X_out[N2 + 1 + n] = Z2[N8 - 1 - k].im; X_out[N2 + N4 + n] = Z2[ k].im; X_out[N2 + N4 + 1 + n] = -Z2[N4 - 1 - k].re; } } #ifdef LTP_DEC void faad_mdct(mdct_info *mdct, real_t *X_in, real_t *X_out) { uint16_t k; fftw_complex *Z1 = mdct->Z1; fftw_complex *Z2 = mdct->Z2; faad_sincos *sincos = mdct->sincos; uint16_t N = mdct->N; uint16_t N2 = N >> 1; uint16_t N4 = N >> 2; uint16_t N8 = N >> 3; /* pre-FFT complex multiplication */ for (k = 0; k < N8; k++) { uint16_t n = k << 1; real_t zr = X_in[N - N4 - 1 - n] + X_in[N - N4 + n]; real_t zi = X_in[ N4 + n] - X_in[ N4 - 1 - n]; Z1[k ].re = -zr * sincos[k ].cos - zi * sincos[k ].sin; Z1[k ].im = -zi * sincos[k ].cos + zr * sincos[k ].sin; zr = X_in[ N2 - 1 - n] - X_in[ n]; zi = X_in[ N2 + n] + X_in[N - 1 - n]; Z1[k + N8].re = -zr * sincos[k + N8].cos - zi * sincos[k + N8].sin; Z1[k + N8].im = -zi * sincos[k + N8].cos + zr * sincos[k + N8].sin; } /* complex FFT */ fftw_one(mdct->plan_forward, Z1, Z2); /* post-FFT complex multiplication */ for (k = 0; k < N4; k++) { uint16_t n = k << 1; real_t zr = 2.0 * (Z2[k].re * sincos[k].cos + Z2[k].im * sincos[k].sin); real_t zi = 2.0 * (Z2[k].im * sincos[k].cos - Z2[k].re * sincos[k].sin); X_out[ n] = -zr; X_out[N2 - 1 - n] = zi; X_out[N2 + n] = -zi; X_out[N - 1 - n] = zr; } } #endif