imdct.c

语音压缩算法

/**********************************************************************
  This software module was originally developed by
  Bernhard Grill (University of Erlangen)
  and edited by
  Huseyin Kemal Cakmak (Fraunhofer Gesellschaft IIS)
  and edited by
  Bernhard Grill (University of Erlangen),
  Takashi Koike (Sony Corporation)
  Naoki Iwakami (Nippon Telegraph and Telephone)

  FN3 LN3 (CN3),
  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. 
  **********************************************************************/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "block.h"               /* handler, defines, enums */
#include "buffersHandle.h"       /* handler, defines, enums */
#include "concealmentHandle.h"   /* handler, defines, enums */
#include "mod_bufHandle.h"       /* handler, defines, enums */
#include "reorderspecHandle.h"   /* handler, defines, enums */
#include "resilienceHandle.h"    /* handler, defines, enums */
#include "tf_mainHandle.h"       /* handler, defines, enums */

#include "all.h"                 /* structs */
#include "nok_ltp_common.h"      /* structs */
#include "obj_descr.h"           /* structs */
#include "tf_mainStruct.h"       /* structs */

#include "allVariables.h"        /* variables */

#include "common_m4a.h"
/*#include "tf_main.h"*/
#include "transfo.h"
#include "dolby_win.h"
#include "dolby_win_ssr.h"

#define DPI	3.14159265358979323846264338327950288
#define USE_FAST_IMDCT 1

static double zero = 0;

/* This is actually not an fft but implements just the dft formula,
   and should therefore be replaced */

void fft( double in[], double out[], int len )
{
  float tmp1, tmp2;
  int n, m;
  tmp1 = 2*DPI/len;
  for( m=0; m<len; m++ ) {
    tmp2 = tmp1*m;
    out[2*m] = out[2*m+1] = 0;
    for( n=0; n<len; n++ ) {
      out[2*m]   += in[2*n] * cos( tmp2*n ) - in[2*n+1] * sin( tmp2*n );
      out[2*m+1] -= in[2*n] * sin( tmp2*n ) + in[2*n+1] * cos( tmp2*n );
    }
  }
}

static void vcopy( double src[], double dest[], int inc_src, int inc_dest, int vlen )
{
  int i;

  for( i=0; i<vlen-1; i++ ) {
    *dest = *src;
    dest += inc_dest;
    src  += inc_src;
  }
  if (vlen) /* just for bounds-checkers sake */
    *dest = *src;

}

static void vmult( double src1[], double src2[], double dest[], 
            int inc_src1, int inc_src2, int inc_dest, int vlen )
{
  int i;

  for( i=0; i<vlen-1; i++ ) {
    *dest = *src1 * *src2;
    dest += inc_dest;
    src1 += inc_src1;
    src2 += inc_src2;
  }
  if (i<vlen)
    *dest = *src1 * *src2;

}

static void vadd( double src1[], double src2[], double dest[], 
            int inc_src1, int inc_src2, int inc_dest, int vlen )
{
  int i;

  for( i=0; i<vlen; i++ ) {
    *dest = *src1 + *src2;
    dest += inc_dest;
    src1 += inc_src1;
    src2 += inc_src2;
  }
}
 
/*****************************************************************************

    functionname: Izero
    description:  calculates the modified Bessel function of the first kind
    returns:      value of Bessel function
    input:        argument
    output:       

*****************************************************************************/
static double Izero(double x)
{
  const double IzeroEPSILON = 1E-41;  /* Max error acceptable in Izero */
  double sum, u, halfx, temp;
  int n;
  
  sum = u = n = 1;
  halfx = x/2.0;
  do {
    temp = halfx/(double)n;
    n += 1;
    temp *= temp;
    u *= temp;
    sum += u;
  } while (u >= IzeroEPSILON*sum);

  return(sum);
}



/*****************************************************************************

    functionname: CalculateKBDWindow
    description:  calculates the window coefficients for the Kaiser-Bessel
                  derived window
    returns:        
    input:        window length, alpha
    output:       window coefficients

*****************************************************************************/
static void CalculateKBDWindow(double* win, double alpha, int length)
{
#ifndef M_PI
#define M_PI 3.14159265358979323846264338327950288
#endif
  int i;
  double IBeta;
  double tmp;
  double sum = 0.0;

  alpha *= M_PI;
  IBeta = 1.0/Izero(alpha);
 
  /* calculate lower half of Kaiser Bessel window */
  for(i=0; i<(length>>1); i++) {
    tmp = 4.0*(double)i/(double)length - 1.0;
    win[i] = Izero(alpha*sqrt(1.0-tmp*tmp))*IBeta;
    sum += win[i];
  }

  sum = 1.0/sum;
  tmp = 0.0;

  /* calculate lower half of window */
  for(i=0; i<(length>>1); i++) {
    tmp += win[i];
    win[i] = sqrt(tmp*sum);
  }
}

/* Calculate window */
static void calc_window( double window[], int len, WINDOW_SHAPE wfun_select)

{
  int i;
  extern double dolby_win_1024[];
  extern double dolby_win_960[];
  extern double dolby_win_128[];
  extern double dolby_win_120[];

  switch(wfun_select) {
  case WS_FHG: 
    for( i=0; i<len; i++ ) 
      window[i] = sin( ((i+1)-0.5) * DPI / (2*len) );
    break;
    
  case WS_DOLBY: 
    switch(len)
      {
      case BLOCK_LEN_SHORT_S:
        for( i=0; i<len; i++ ) 
          window[i] = dolby_win_120[i]; 
        break;
      case BLOCK_LEN_SHORT:
        for( i=0; i<len; i++ ) 
          window[i] = dolby_win_128[i]; 
        break;
      case BLOCK_LEN_SHORT_SSR:
        for( i=0; i<len; i++ ) 
          window[i] = dolby_win_32[i]; 
        break;
      case BLOCK_LEN_MEDIUM_S:
      case BLOCK_LEN_MEDIUM:
        for( i=0; i<len; i++ ) 
          window[i] = sin( ((double)(i+1)-0.5) * DPI / (double)(2*len) ); /* Using WS_FHG */
        break;
      case BLOCK_LEN_LONG_S:
        for( i=0; i<len; i++ ) 
          window[i] = dolby_win_960[i]; 
        break;
      case BLOCK_LEN_LONG:
        for( i=0; i<len; i++ ) 
          window[i] = dolby_win_1024[i]; 
        break;
      case BLOCK_LEN_LONG_SSR:
        for( i=0; i<len; i++ ) 
          window[i] = dolby_win_256[i]; 
        break;
      case 4:
        memcpy(window, dolby_win_4, 4*sizeof(double));
        break;
      case 16:
        memcpy(window, dolby_win_16, 16*sizeof(double));
        break;
      default:
        CalculateKBDWindow(window, 6.0, 2*len);
        /* 	  CommonExit(1,"Unsupported window size: %d", len); */
        break;
      }
    break;

  default:
    CommonExit(1,"Unsupported window shape: %d", wfun_select);
    break;
  }
}


/* %%%%%%%%%%%%%%%%%% IMDCT - STUFF %%%%%%%%%%%%%%%%*/

#define MAX_SHIFT_LEN_LONG 4096

void imdct(double in_data[], double out_data[], int len)
{
  vcopy(in_data, out_data, 1, 1, len/2);
  MDCT(out_data, len, len/2, -1);
}

void freq2buffer(
  double           p_in_data[], 
  double           p_out_data[],
  double           p_overlap[],
  WINDOW_SEQUENCE  windowSequence,
  int              nlong,            /* shift length for long windows   */
  int              nmed,             /* shift length for medium windows */
  int              nshort,           /* shift length for short windows  */
  WINDOW_SHAPE     wfun_select,      /* offers the possibility to select different window functions */
  WINDOW_SHAPE     wfun_select_prev, /* YB : 971113 */
  Imdct_out	   overlap_select,   /* select imdct output *TK*	*/
				     /* switch (overlap_select) {	*/
				     /* case OVERLAPPED_MODE:		*/
				     /*   p_out_data[]			*/
				     /*   = overlapped and added signal */
				     /*		(bufferlength: nlong)	*/
				     /* case NON_OVERLAPPED_MODE:		*/
				     /*   p_out_data[]			*/
				     /*   = non overlapped signal	*/
				     /*		(bufferlength: 2*nlong)	*/
  int              num_short_win     /* number of short windows to      */
                                     /* transform                       */
  )
{
  double           *o_buf, transf_buf[ 2*MAX_SHIFT_LEN_LONG ];
  double           overlap_buf[ 2*MAX_SHIFT_LEN_LONG ]; 
 
  double           window_long[MAX_SHIFT_LEN_LONG]; 
  double           window_long_prev[MAX_SHIFT_LEN_LONG]; 
  double           window_med[MAX_SHIFT_LEN_LONG]; 
  double           window_med_prev[MAX_SHIFT_LEN_LONG]; 
  double           window_short[MAX_SHIFT_LEN_LONG]; 
  double           window_short_prev[MAX_SHIFT_LEN_LONG]; 
  double           *window_short_prev_ptr;
 
  double  *fp; 
  int      k; 
 
  int nflat_ls    = (nlong-nshort)/ 2; 
  int transfak_ls =  nlong/nshort; 
  int transfak_lm =  nlong/nmed; 
#if 0 /* not used anymore */
  int nflat_lm    = (nlong-nmed)  / 2; 
  int nflat_ms    = (nmed-nshort) / 2; 
#endif
 
  window_short_prev_ptr=window_short_prev ; 


  if( (nlong%nshort) || (nlong > MAX_SHIFT_LEN_LONG) || (nshort > MAX_SHIFT_LEN_LONG/2) ) { 
    CommonExit( 1, "mdct_synthesis: Problem with window length" ); 
  } 
  if( (nlong%nmed) || (nmed%nshort) || (nmed > MAX_SHIFT_LEN_LONG/2) ) { 
    CommonExit( 1, "mdct_synthesis: Problem with window length" ); 
  } 
  if( transfak_lm%2 ) { 
    CommonExit( 1, "mdct_synthesis: Problem with window length" ); 
  } 
  if( windowSequence==EIGHT_SHORT_SEQUENCE
      && ( (num_short_win <= 0) || (num_short_win > transfak_ls) ) ) {
    CommonExit( 1, "mdct_synthesis: Problem with number of short windows" ); 
  } 

  calc_window( window_long,      nlong, wfun_select ); 
  calc_window( window_long_prev, nlong, wfun_select_prev ); 
  calc_window( window_med, nmed, wfun_select ); 
  calc_window( window_med_prev, nmed, wfun_select_prev );
  calc_window( window_short,      nshort, wfun_select ); 
  calc_window( window_short_prev_ptr, nshort, wfun_select_prev ); 

  /* Assemble overlap buffer */ 
  vcopy( p_overlap, overlap_buf, 1, 1, nlong ); 
  o_buf = overlap_buf; 
 
#if 0
  printf("\n input imdct: \n%d:",0);
  for (i=0;i<1024;i++){
#if 0
    p_in_data[i] = 0;
    if (i==50)
      p_in_data[i] = 500;    
#endif
    printf(" %d",(int)p_in_data[i]);
    if (i%50 == 0)
      printf("\n%d:",i);
  }
#endif

  /* Separate action for each Block Type */ 
   switch( windowSequence ) { 
   case ONLY_LONG_SEQUENCE : 
    imdct( p_in_data, transf_buf, 2*nlong ); 
    vmult( transf_buf, window_long_prev, transf_buf, 1, 1, 1, nlong ); 
    if (overlap_select != NON_OVERLAPPED_MODE) {
      vadd( transf_buf, o_buf, o_buf, 1, 1, 1, nlong ); 
      vmult( transf_buf+nlong, window_long+nlong-1, o_buf+nlong, 1, -1, 1, nlong ); 
    }
    else { /* overlap_select == NON_OVERLAPPED_MODE */