transfo.c

MPEG-4编解码的实现(包括MPEG4视音频编解码)

/************************* MPEG-2 NBC Audio Decoder **************************
 *                                                                           *
"This software module was originally developed by
AT&T, Dolby Laboratories, Fraunhofer Gesellschaft IIS in the course of
development of the MPEG-2 NBC/MPEG-4 Audio standard ISO/IEC 13818-7,
14496-1,2 and 3. This software module is an implementation of a part of one or more
MPEG-2 NBC/MPEG-4 Audio tools as specified by the MPEG-2 NBC/MPEG-4
Audio standard. ISO/IEC  gives users of the MPEG-2 NBC/MPEG-4 Audio
standards free license to this software module or modifications thereof for use in
hardware or software products claiming conformance to the MPEG-2 NBC/MPEG-4
Audio  standards. Those intending to use this software module in hardware or
software products are advised that this use may infringe existing patents.
The original developer of this software module and his/her company, the subsequent
editors and their companies, and ISO/IEC have no liability for use of this software
module or modifications thereof in an implementation. Copyright is not released for
non MPEG-2 NBC/MPEG-4 Audio conforming products.The original developer
retains full right to use the code for his/her  own purpose, assign or donate the
code to a third party and to inhibit third party from using the code for non
MPEG-2 NBC/MPEG-4 Audio conforming products. This copyright notice must
be included in all copies or derivative works."
Copyright(c)1996.
 *                                                                           *
 ****************************************************************************/
/*
 * $Id: transfo.c,v 1.5 2002/01/09 22:25:41 wmay Exp $
 */

#include <math.h>
#include "all.h"
#include "transfo.h"
#include "fastfft.h"

/*
#define INTEL_SPL
*/

#ifdef INTEL_SPL
#define nsp_UsesAll
#include <nsp.h>
#include <nspwin.h>

/* choose the library that best fits your processor */
#pragma comment (lib, "nsppxl.lib")
#endif


#ifndef M_PI
#define M_PI        3.14159265358979323846
#endif

#ifndef M_PI_2
#define M_PI_2      1.57079632679489661923
#endif

void MDCT_Long(faacDecHandle hDecoder, fftw_real *data)
{
#ifdef INTEL_SPL
    SCplx FFT_data[512];
#else
    fftw_complex FFTarray[512];    /* the array for in-place FFT */
#endif
    fftw_real tempr, tempi, c, s, cold, cfreq, sfreq; /* temps for pre and post twiddle */
/*  fftw_real freq = 0.0030679616611450911f; */
    fftw_real fac,cosfreq8,sinfreq8;
    int i, n;
    int b = 2048 >> 1;
    int N4 = 2048 >> 2;
    int N2 = 2048 >> 1;
    int a = 2048 - b;
    int a2 = a >> 1;
    int a4 = a >> 2;
    int b4 = b >> 2;
    int unscambled;


    /* Choosing to allocate 2/N factor to Inverse Xform! */
    fac = 2.; /* 2 from MDCT inverse  to forward */

    /* prepare for recurrence relation in pre-twiddle */
    cfreq = 0.99999529123306274f;
    sfreq = 0.0030679567717015743f;
    cosfreq8 = 0.99999994039535522f;
    sinfreq8 = 0.00038349519824312089f;

    c = cosfreq8;
    s = sinfreq8;

    for (i = 0; i < N4; i++)
    {
        /* calculate real and imaginary parts of g(n) or G(p) */

    n = 2048 / 2 - 1 - 2 * i;
        if (i < b4) {
      tempr = data [a2 + n] + data [2048 + a2 - 1 - n]; /* use second form of e(n) for n = N / 2 - 1 - 2i */
        } else {
            tempr = data [a2 + n] - data [a2 - 1 - n]; /* use first form of e(n) for n = N / 2 - 1 - 2i */
        }
        n = 2 * i;
        if (i < a4) {
            tempi = data [a2 + n] - data [a2 - 1 - n]; /* use first form of e(n) for n=2i */
        } else {
      tempi = data [a2 + n] + data [2048 + a2 - 1 - n]; /* use second form of e(n) for n=2i*/
        }

        /* calculate pre-twiddled FFT input */
#ifdef INTEL_SPL
    FFT_data[i].re = tempr * c + tempi * s;
    FFT_data[i].im = tempi * c - tempr * s;
#else
        FFTarray [i].re = tempr * c + tempi * s;
        FFTarray [i].im = tempi * c - tempr * s;
#endif

        /* use recurrence to prepare cosine and sine for next value of i */
        cold = c;
        c = c * cfreq - s * sfreq;
        s = s * cfreq + cold * sfreq;
    }

    /* Perform in-place complex FFT of length N/4 */
#ifdef INTEL_SPL
    nspcFft(FFT_data, 9, NSP_Forw);
#else
    pfftw_512(FFTarray);
#endif


    /* prepare for recurrence relations in post-twiddle */
    c = cosfreq8;
    s = sinfreq8;

    /* post-twiddle FFT output and then get output data */
    for (i = 0; i < N4; i++)
    {

      /* get post-twiddled FFT output  */
      /* Note: fac allocates 4/N factor from IFFT to forward and inverse */
#ifdef INTEL_SPL
      tempr = fac * (FFT_data[i].re * c + FFT_data[i].im * s);
      tempi = fac * (FFT_data[i].im * c - FFT_data[i].re * s);
#else
      unscambled = hDecoder->unscambled512[i];

      tempr = fac * (FFTarray [unscambled].re * c + FFTarray [unscambled].im * s);
      tempi = fac * (FFTarray [unscambled].im * c - FFTarray [unscambled].re * s);
#endif

      /* fill in output values */
      data [2 * i] = -tempr;   /* first half even */
      data [N2 - 1 - 2 * i] = tempi;  /* first half odd */
      data [N2 + 2 * i] = -tempi;  /* second half even */
      data [2048 - 1 - 2 * i] = tempr;  /* second half odd */

      /* use recurrence to prepare cosine and sine for next value of i */
      cold = c;
      c = c * cfreq - s * sfreq;
      s = s * cfreq + cold * sfreq;
    }
}

void MDCT_Short(faacDecHandle hDecoder, fftw_real *data)
{
#ifdef INTEL_SPL
    SCplx FFT_data[64];
#else
    fftw_complex FFTarray[64];    /* the array for in-place FFT */
#endif
    fftw_real tempr, tempi, c, s, cold, cfreq, sfreq; /* temps for pre and post twiddle */
/*  fftw_real freq = 0.024543693289160728f; */
    fftw_real fac,cosfreq8,sinfreq8;
    int i, n;
    int b = 256 >> 1;
    int N4 = 256 >> 2;
    int N2 = 256 >> 1;
    int a = 256 - b;
    int a2 = a >> 1;
    int a4 = a >> 2;
    int b4 = b >> 2;
    int unscambled;


    /* Choosing to allocate 2/N factor to Inverse Xform! */
    fac = 2.; /* 2 from MDCT inverse  to forward */

    /* prepare for recurrence relation in pre-twiddle */
    cfreq = 0.99969881772994995f;
    sfreq = 0.024541229009628296f;
    cosfreq8 = 0.99999529123306274f;
    sinfreq8 = 0.0030679568483393833f;

    c = cosfreq8;
    s = sinfreq8;

    for (i = 0; i < N4; i++)
    {
        /* calculate real and imaginary parts of g(n) or G(p) */

    n = 256 / 2 - 1 - 2 * i;
        if (i < b4) {
      tempr = data [a2 + n] + data [256 + a2 - 1 - n]; /* use second form of e(n) for n = N / 2 - 1 - 2i */
        } else {
            tempr = data [a2 + n] - data [a2 - 1 - n]; /* use first form of e(n) for n = N / 2 - 1 - 2i */
        }
        n = 2 * i;
        if (i < a4) {
            tempi = data [a2 + n] - data [a2 - 1 - n]; /* use first form of e(n) for n=2i */
        } else {
      tempi = data [a2 + n] + data [256 + a2 - 1 - n]; /* use second form of e(n) for n=2i*/
        }

        /* calculate pre-twiddled FFT input */
#ifdef INTEL_SPL
    FFT_data[i].re = tempr * c + tempi * s;
    FFT_data[i].im = tempi * c - tempr * s;
#else
        FFTarray [i].re = tempr * c + tempi * s;
        FFTarray [i].im = tempi * c - tempr * s;
#endif

        /* use recurrence to prepare cosine and sine for next value of i */
        cold = c;
        c = c * cfreq - s * sfreq;
        s = s * cfreq + cold * sfreq;
    }

    /* Perform in-place complex FFT of length N/4 */
#ifdef INTEL_SPL
    nspcFft(FFT_data, 6, NSP_Forw);
#else
    pfftw_64(FFTarray);
#endif

    /* prepare for recurrence relations in post-twiddle */
    c = cosfreq8;
    s = sinfreq8;

    /* post-twiddle FFT output and then get output data */
    for (i = 0; i < N4; i++) {

        /* get post-twiddled FFT output  */
        /* Note: fac allocates 4/N factor from IFFT to forward and inverse */
#ifdef INTEL_SPL
    tempr = fac * (FFT_data[i].re * c + FFT_data[i].im * s);
    tempi = fac * (FFT_data[i].im * c - FFT_data[i].re * s);
#else
        unscambled = hDecoder->unscambled64[i];

    tempr = fac * (FFTarray [unscambled].re * c + FFTarray [unscambled].im * s);
    tempi = fac * (FFTarray [unscambled].im * c - FFTarray [unscambled].re * s);
#endif

        /* fill in output values */
        data [2 * i] = -tempr;   /* first half even */
        data [N2 - 1 - 2 * i] = tempi;  /* first half odd */
        data [N2 + 2 * i] = -tempi;  /* second half even */
    data [256 - 1 - 2 * i] = tempr;  /* second half odd */

        /* use recurrence to prepare cosine and sine for next value of i */
        cold = c;
        c = c * cfreq - s * sfreq;
        s = s * cfreq + cold * sfreq;
    }
}


void IMDCT_Long(faacDecHandle hDecoder, fftw_real *data)
{
#ifdef INTEL_SPL
  SCplx FFT_data[512];    /* the array for in-place FFT */
#else
    fftw_complex FFTarray[512];    /* the array for in-place FFT */
#endif
    fftw_real tempr, tempi, c, s, cold, cfreq, sfreq; /* temps for pre and post twiddle */

/*  fftw_real freq = 0.0030679616611450911f; */
    fftw_real fac, cosfreq8, sinfreq8;
    int i;
  int Nd2 = 2048 >> 1;
  int Nd4 = 2048 >> 2;
  int Nd8 = 2048 >> 3;
    int unscambled;

    /* Choosing to allocate 2/N factor to Inverse Xform! */
  fac = 0.0009765625f;

    /* prepare for recurrence relation in pre-twiddle */
  cfreq = 0.99999529123306274f;
  sfreq = 0.0030679567717015743f;
  cosfreq8 = 0.99999994039535522f;
  sinfreq8 = 0.00038349519824312089f;

    c = cosfreq8;
    s = sinfreq8;


    for (i = 0; i < Nd4; i++) {

        /* calculate real and imaginary parts of g(n) or G(p) */
        tempr = -data [2 * i];
        tempi = data [Nd2 - 1 - 2 * i];

        /* calculate pre-twiddled FFT input */
#ifdef INTEL_SPL
    FFT_data[i].re = tempr * c - tempi * s;
    FFT_data[i].im = tempi * c + tempr * s;
#else
        unscambled = hDecoder->unscambled512[i];

        FFTarray [unscambled].re = tempr * c - tempi * s;
        FFTarray [unscambled].im = tempi * c + tempr * s;
#endif

        /* use recurrence to prepare cosine and sine for next value of i */
        cold = c;
        c = c * cfreq - s * sfreq;
        s = s * cfreq + cold * sfreq;
    }

    /* Perform in-place complex IFFT of length N/4 */
#ifdef INTEL_SPL
  nspcFft(FFT_data, 9, NSP_Inv | NSP_NoScale);
#else
  pfftwi_512(FFTarray);
#endif

    /* prepare for recurrence relations in post-twiddle */
    c = cosfreq8;
    s = sinfreq8;

    /* post-twiddle FFT output and then get output data */
    for (i = 0; i < Nd4; i++) {

        /* get post-twiddled FFT output  */
#ifdef INTEL_SPL
  tempr = fac * (FFT_data[i].re * c - FFT_data[i].im * s);
      tempi = fac * (FFT_data[i].im * c + FFT_data[i].re * s);
#else
    tempr = fac * (FFTarray[i].re * c - FFTarray[i].im * s);
        tempi = fac * (FFTarray[i].im * c + FFTarray[i].re * s);
#endif

        /* fill in output values */
        data [Nd2 + Nd4 - 1 - 2 * i] = tempr;
        if (i < Nd8) {
            data [Nd2 + Nd4 + 2 * i] = tempr;
        } else {
            data [2 * i - Nd4] = -tempr;
        }
        data [Nd4 + 2 * i] = tempi;
        if (i < Nd8) {
            data [Nd4 - 1 - 2 * i] = -tempi;
        } else {
      data [Nd4 + 2048 - 1 - 2*i] = tempi;
        }

        /* use recurrence to prepare cosine and sine for next value of i */
        cold = c;
        c = c * cfreq - s * sfreq;
        s = s * cfreq + cold * sfreq;
  }
}


void IMDCT_Short(faacDecHandle hDecoder, fftw_real *data)
{
#ifdef INTEL_SPL
  SCplx FFT_data[64];    /* the array for in-place FFT */
#else
    fftw_complex FFTarray[64];    /* the array for in-place FFT */
#endif
    fftw_real tempr, tempi, c, s, cold, cfreq, sfreq; /* temps for pre and post twiddle */
/*  fftw_real freq = 0.024543693289160728f; */
    fftw_real fac, cosfreq8, sinfreq8;
    int i;
    int Nd2 = 256 >> 1;
    int Nd4 = 256 >> 2;
    int Nd8 = 256 >> 3;
    int unscambled;

    /* Choosing to allocate 2/N factor to Inverse Xform! */
  fac = 0.0078125f; /* remaining 2/N from 4/N IFFT factor */

    /* prepare for recurrence relation in pre-twiddle */
  cfreq = 0.99969881772994995f;
  sfreq = 0.024541229009628296f;
  cosfreq8 = 0.99999529123306274f;
  sinfreq8 = 0.0030679568483393833f;

    c = cosfreq8;
    s = sinfreq8;

    for (i = 0; i < Nd4; i++)
  {

        /* calculate real and imaginary parts of g(n) or G(p) */
        tempr = -data [2 * i];
        tempi = data [Nd2 - 1 - 2 * i];

        /* calculate pre-twiddled FFT input */
#ifdef INTEL_SPL
    FFT_data[i].re = tempr * c - tempi * s;
    FFT_data[i].im = tempi * c + tempr * s;
#else
        unscambled = hDecoder->unscambled64[i];

        FFTarray [unscambled].re = tempr * c - tempi * s;
        FFTarray [unscambled].im = tempi * c + tempr * s;
#endif

        /* use recurrence to prepare cosine and sine for next value of i */
        cold = c;
        c = c * cfreq - s * sfreq;
        s = s * cfreq + cold * sfreq;
    }

    /* Perform in-place complex IFFT of length N/4 */

#ifdef INTEL_SPL
  nspcFft(FFT_data, 6, NSP_Inv | NSP_NoScale);
#else
    pfftwi_64(FFTarray);
#endif

    /* prepare for recurrence relations in post-twiddle */
    c = cosfreq8;
    s = sinfreq8;

    /* post-twiddle FFT output and then get output data */
    for (i = 0; i < Nd4; i++) {

        /* get post-twiddled FFT output  */
#ifdef INTEL_SPL
      tempr = fac * (FFT_data[i].re * c - FFT_data[i].im * s);
      tempi = fac * (FFT_data[i].im * c + FFT_data[i].re * s);
#else
    tempr = fac * (FFTarray[i].re * c - FFTarray[i].im * s);
        tempi = fac * (FFTarray[i].im * c + FFTarray[i].re * s);
#endif

        /* fill in output values */
        data [Nd2 + Nd4 - 1 - 2 * i] = tempr;
        if (i < Nd8) {
            data [Nd2 + Nd4 + 2 * i] = tempr;
        } else {
            data [2 * i - Nd4] = -tempr;
        }
        data [Nd4 + 2 * i] = tempi;
        if (i < Nd8) {
            data [Nd4 - 1 - 2 * i] = -tempi;
        } else {
      data [Nd4 + 256 - 1 - 2*i] = tempi;
        }

        /* use recurrence to prepare cosine and sine for next value of i */
        cold = c;
        c = c * cfreq - s * sfreq;
        s = s * cfreq + cold * sfreq;
    }
}