ntt_scale_coder.c

MPEG2/MPEG4编解码参考程序(实现了MPEG4的部分功能)

/*****************************************************************************/
/* This software module was originally developed by                          */
/*   Akio Jin      (NTT),                                                    */
/*   Takeshi Norimatsu,                                                      */
/*   Mineo Tsushima,                                                         */
/*   and Tomokazu Ishikawa (Matsushita Electric Industrial Co Ltd.)          */
/* and edited by                                                             */
/*   Naoki Iwakami (NTT) on 1997-07-17,                                      */
/*   Akio Jin (NTT),                                                         */
/*   Mineo Tsushima, (Matsushita Electric Industrial Co Ltd.)                */
/*   and Tomokazu Ishikawa (Matsushita Electric Industrial Co Ltd.)          */
/*   on 1997-10-23,                                                          */
/* 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 <stdio.h>
#include <math.h>

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

#include "nok_ltp_common.h"      /* structs */
#include "tf_mainStruct.h"       /* structs */
#include "tns.h"                 /* structs */

#include "nok_bwp_enc.h"
#include "ntt_conf.h"
#include "ntt_scale_conf.h"
#include "ntt_scale_enc_para.h"
#include "ntt_encode.h"
#include "ntt_scale_encode.h"
#include "ntt_tools.h"
#include "mat_def_ntt.h"


void ntt_scale_enc_lpc_spectrum
  (/* Parameters */
   int    nfr,
   int    nsf,
   int    block_size_samples,
   int    n_ch,
   int    n_pr,
   double *lsp_code,
   double *lsp_fgcode,
   int    prev_lsp_code[],
   /*
   double prev_buf[MAX_TIME_CHANNELS][ntt_MA_NP][ntt_N_PR_MAX+1],
   */
   int    ma_np,
   /* Input */
   double spectrum[],
   int    w_type,
   int    current_block,
   /* Output */
   int    index_lsp[MAX_TIME_CHANNELS][ntt_LSP_NIDX_MAX],
   double lpc_spectrum[],
   /* Scalable layer ID */
   int iscl,
   double band_lower[],
   double band_upper[],
   double *cos_TT
   )
{
  /*--- Variables ---*/
  double resid, alf[MAX_TIME_CHANNELS][ntt_N_PR_MAX+1],
  lsp[ntt_N_PR_MAX+1], lspq[MAX_TIME_CHANNELS*(ntt_N_PR_MAX+1)];
  double spectrum_wide[ntt_T_FR_MAX],
         lpc_spectrum_re[ntt_T_FR_MAX];
  int    i_sup, top, subtop, w_type_t, isf, lsptop, ismp;
  int    block_size, nfr_l, nfr_u, nfr_lu;
  int lsp_csize[2];
  float ftmp ;

  block_size = nfr * nsf;

  lsp_csize[0] = ntt_NC0_SCL;
  lsp_csize[1] = ntt_NC1_SCL;
  /*--- encode LSP coefficients ---*/
  if (current_block == ntt_BLOCK_LONG) w_type_t = w_type;
  else                                 w_type_t = ONLY_LONG_SEQUENCE;
  for ( i_sup=0; i_sup<n_ch; i_sup++ ){
       ftmp = (float)band_upper[i_sup];
       ftmp *= (float)nfr;
       nfr_u  = (int)ftmp;
       ftmp = (float)band_lower[i_sup];
       ftmp *= (float)nfr;
       nfr_l  = (int)ftmp;
       ftmp = (float)band_upper[i_sup] - (float)band_lower[i_sup];
       ftmp *= (float)nfr;
       nfr_lu = (int)ftmp;
       top    = i_sup*block_size;
       ntt_zerod(block_size, spectrum_wide+top);
       ntt_setd(block_size,1.e5, lpc_spectrum+top);

       ftmp = (float)band_upper[i_sup] - (float)band_lower[i_sup];
       ftmp = 1./ftmp;
       /*--- band extension for nallow band lpc (by A.Jin 1997.6.9) ---*/
       for(isf=0; isf<nsf; isf++){
	    float fftmp;
	    subtop = top + isf*nfr;
	    for(ismp=nfr_l; ismp<nfr_l+nfr_lu; ismp++){
		 fftmp = ftmp;
		 fftmp *= (float)(ismp-nfr_l);
		 spectrum_wide[(int)fftmp+subtop] = spectrum[ismp+subtop];
	    }
       }

       /* calculate LSP coefficients */
       ntt_scale_fgetalf(nfr, nsf, block_size_samples,
			 spectrum_wide+top,alf[i_sup],&resid,
			 ntt_N_PR_SCL, ntt_BAND_EXPN_SCL, cos_TT);
       ntt_alflsp(ntt_N_PR_SCL, alf[i_sup], lsp);
       /* quantize LSP coefficients */
       ntt_lsp_encw(ntt_N_PR_SCL, (double (*)[ntt_N_PR_MAX])lsp_code,
		    (double (*)[ntt_MA_NP][ntt_N_PR_MAX])lsp_fgcode,
		    lsp_csize,
		    prev_lsp_code+i_sup*ntt_LSP_NIDX_MAX,
		    /* prev_buf[i_sup], */
		    ma_np, 
		    lsp, lspq+i_sup*ntt_N_PR_SCL,
		    index_lsp[i_sup], ntt_LSP_SPLIT_SCL);
  }

  /*--- reconstruct LPC spectrum ---*/
  for (i_sup=0; i_sup<n_ch; i_sup++){
       ftmp = (float)band_lower[i_sup];
       ftmp *= (float)nfr;
       nfr_l  =(int)ftmp;
       ftmp = (float)band_upper[i_sup]-(float)band_lower[i_sup];
       ftmp *= (float)nfr;
       nfr_lu = (int)ftmp;
       nfr_u  = nfr_l+nfr_lu;
       top = i_sup * block_size;
       lsptop = i_sup * n_pr;
       ntt_lsptowt(nfr, n_pr, block_size_samples,
		   lspq+lsptop, lpc_spectrum+top, cos_TT);

       /*--- band re-extension for nallow band lpc (by A.Jin 1997.6.9) ---*/
       ftmp = (float)band_upper[i_sup] - (float)band_lower[i_sup];
       ftmp = 1./ftmp;
       for(ismp=0; ismp<nfr_lu; ismp++){
	    float fftmp;
	    fftmp = ftmp*(float)ismp;
	    lpc_spectrum_re[ismp+top+nfr_l] = lpc_spectrum[(int)fftmp+top];
       }
       ntt_movdd(nfr_lu,lpc_spectrum_re+top+nfr_l, lpc_spectrum+top+nfr_l);
       ntt_setd(nfr_l, 999., lpc_spectrum+top);
       ntt_setd(nfr- nfr_u, 999., lpc_spectrum+top+nfr_u);

       for (isf=1; isf<nsf; isf++){
	    subtop = top + isf*nfr;
	    ntt_movdd(nfr, lpc_spectrum+top, lpc_spectrum+subtop);
       }
  }
}


#define ntt_GAMMA  0.95 

void ntt_scale_fgetalf(/* Parameters */
		 int    nfr,    /* subframe size */
		 int    nsf,    /* number of subframes */
		 int    block_size_samples,
		 /* Input */
		 double *spectrum,  /* MDCT spectrum */
		 /* Output */
		 double *alf,   /* LPC coefficients */
		 double *resid, /* residual power */
		 /* Input */
		 int n_pr,  
		 double band_exp,
		 double *ntt_cos_TT) 
{
    double wlag[ntt_N_PR_MAX + 1];
    double powspec[ntt_N_FR_MAX];
    int    unstable;
    int    ip;
    double ref[ntt_N_PR_MAX+1];         
    double cep[(ntt_N_PR_MAX+1)*2];        
    double refb[ntt_N_PR_MAX+1];      
    double alfb[ntt_N_PR_MAX+1];     
    double gamma; 			   
    double          
	cor[ntt_N_PR_MAX+1],
        pow_counter;
    double 
      xr[ntt_N_FR_MAX*2];
    int    ismp, isf;

    /*--- Initialization ---*/
    /*--- set lag window ---*/
    ntt_lagwdw( wlag, n_pr, band_exp );

    /*--- Calculate power spectrum ---*/
    pow_counter = 0.;
    for (ismp=0; ismp<nfr; ismp++) powspec[ismp] = 0.;
    for (isf=0; isf<nsf; isf++){
      for (ismp=0; ismp<nfr; ismp++){
	powspec[ismp] += spectrum[ismp+isf*nfr]*spectrum[ismp+isf*nfr];
	pow_counter   += powspec[ismp]*1.e-6;
      }
    }
    /*--- convert low power signals to zero(for certain LPC) ---*/
    if(pow_counter<1.e-2){
      ntt_zerod(nfr, powspec);
    }
    /*--- Calculate auto correlation coefficients ---*/
 {
 double tmp, theta;
 int jsmp, pointer, pointer0;
     for(jsmp=0; jsmp<n_pr+1; jsmp++){
       tmp=0.0;
       pointer =0; pointer0=0; theta=1.0;
       for(ismp=0; ismp<nfr; ismp++){
        tmp += powspec[ismp]*theta;
        pointer += jsmp*block_size_samples*4/nfr;
        if(pointer >= block_size_samples*8) pointer -= block_size_samples*8;
        if(0<= pointer && pointer <block_size_samples*2 ) {
          pointer0 = pointer;
	  theta = ntt_cos_TT[pointer0];
        }
        else if(block_size_samples*2 <= pointer && pointer <block_size_samples*4 ){
            pointer0 = block_size_samples*4-pointer;
            if(pointer0==block_size_samples*2) theta =0.0;
	    else  theta = -ntt_cos_TT[pointer0];
        }
        else if(block_size_samples*4<= pointer && pointer <block_size_samples*6 ){
            pointer0 = pointer-block_size_samples*4;
            theta = -ntt_cos_TT[pointer0];
        }
        else if(block_size_samples*6<= pointer && pointer <block_size_samples*8 ){
            pointer0 = block_size_samples*8-pointer;
            if(pointer0==block_size_samples*2) theta =0.0;
	    else  theta = ntt_cos_TT[pointer0];
        }
       }
        xr[jsmp]= tmp;
     }
 }
    ntt_movdd(n_pr+1, xr, cor);
    for (ip=1; ip<=n_pr; ip++){
      cor[ip] /= (cor[0] + 1.e-5)*1.02;
    }
    cor[0] = 1.0;

    /*--- LPC analysis --*/
    ntt_mulddd(n_pr + 1, cor, wlag, cor);
    ntt_corref(n_pr, cor, alfb ,refb, resid); alfb[0] = 1.;

    /*--- LPC parameters to cepstrum parameters ---*/
    ntt_alfcep(n_pr, alfb, cep, (n_pr)*2);
   
    /*--- Divide cepstrum parameters by 2 (sqrt in the spectrum domain) ---*/
    ntt_mulcdd((n_pr)*2, 0.5, cep+1, cep+1);

    /*--- Cepstrum parameters to LPC alpha parameters ---*/
    ntt_cep2alf((n_pr)*2, n_pr, cep, alf);

    /*--- Stability check ---*/
    alf[0]=1.0;
    do{
       ntt_alfref(n_pr, alf, ref, resid);
       unstable =0;
       for(ip=1; ip<=n_pr; ip++){
          if(ref[ip] > 0.999 || ref[ip] < -0.999) unstable =1;
          if(ref[ip] > 0.999 || ref[ip] < -0.999)
	       fprintf(stderr,
		       "ntt_scale_fgetalf(): unstable ref. pram. order:%5d, value: %7.4f\n",
		       ip, ref[ip]);
       }
       if(unstable) {
         gamma= ntt_GAMMA;
         for(ip=1; ip<=n_pr; ip++){
           alf[ip] *= gamma;
           gamma *= ntt_GAMMA;
         }
       }
    }while(unstable);

}

void ntt_scale_tf_perceptual_model(/* Input */
				   double lpc_spectrum[], /* LPC spectrum */
				   double bark_env[],     /* Bark-scale envelope */
				   double gain[],         /* gain factor */
				   int       w_type,      /* block type */
				   double spectrum[],     /* specturm */
				   double pitch_sequence[], /* pitch peak components */
				   int    numChannel,
				   int    block_size_samples,
				   /* Output */
				   double perceptual_weight[]) /* perceptual weight */
{
    /*--- Variables ---*/
    int    ismp, i_sup, top, isf, subtop;
    double ratio;
    double ratio0;

    switch(w_type){
    case ONLY_LONG_SEQUENCE:
    case LONG_START_SEQUENCE:
    case LONG_STOP_SEQUENCE:
	for (i_sup=0; i_sup<numChannel; i_sup++){
	    top = i_sup * block_size_samples;
	    for (ismp=0; ismp<block_size_samples; ismp++){
		perceptual_weight[ismp+top] = (1./lpc_spectrum[ismp+top]);
	    }
            ntt_prcptw(perceptual_weight+top, bark_env+top,
	    ntt_GAMMA_W2, ntt_GAMMA_W_MIC2, block_size_samples);
            ratio = pow((gain[i_sup]+0.0001), ntt_GAMMA_CH)/(gain[i_sup]+0.0001);
	    for (ismp=0; ismp<block_size_samples; ismp++){
		perceptual_weight[ismp+top] *=
		    lpc_spectrum[ismp+top]/bark_env[ismp+top]*ratio;
	    }

	}
	break;
    case EIGHT_SHORT_SEQUENCE:
	for (i_sup=0; i_sup<numChannel; i_sup++){
	    top = i_sup * block_size_samples;
	    for (isf=0; isf<ntt_N_SHRT; isf++){
		subtop = top + isf * block_size_samples/8;
		for (ismp=0; ismp<block_size_samples/8; ismp++){
		    perceptual_weight[ismp+subtop] =
		      (bark_env[ismp+subtop] / lpc_spectrum[ismp+subtop]);
		}
	    }
	}
        ntt_prcptw_s(perceptual_weight, gain, perceptual_weight,
	      ntt_GAMMA_W_S2, ntt_GAMMA_W_S_T2,
	      numChannel, block_size_samples);
	for (i_sup=0; i_sup<numChannel; i_sup++){
	    top = i_sup * block_size_samples;
	    for (isf=0; isf<ntt_N_SHRT; isf++){
		subtop = top + isf * block_size_samples/8;
		for (ismp=0; ismp<block_size_samples/8; ismp++){
		    perceptual_weight[ismp+subtop] *=
			lpc_spectrum[ismp+subtop] / bark_env[ismp+subtop];
		}
	    }
	}
	break;
    }
}


void ntt_scale_tf_proc_spectrum(/* In/Out */
                                double spectrum[],
                                /* In/Out */
                                ntt_INDEX  *index,
                                double   *cos_TT,
				/* Output */
                                double lpc_spectrum[],
                                double bark_env[],
                                double pitch_sequence[],
                                double gain[],
                                /* In/Out */
                                int iscl,
                                double lpc_spectrum_org[],
                                double lpc_spectrum_mod[])
{
  double tc[ntt_T_FR_MAX], pwt[ntt_T_FR_MAX];
  int    current_block, i_ch, i_ma, ii;
  double band_l[MAX_TIME_CHANNELS], band_u[MAX_TIME_CHANNELS];
  static int InitFlag = 1;

  /*--- Parameters ---*/
  double band_lower, band_upper;
  int    nfr, nsf, n_ch;
  int    n_crb, bark_tbl[128];
  int    fw_ndiv, fw_length, fw_len_max, fw_cb_size;
  double fw_prcpt;
  double *fw_codev;
  /* double *fw_env_memory, *fw_pdot; */
  int    n_pr;
  double *lsp_code;
  double *lsp_fgcode;
  int    lsp_csize[2], lsp_cdim[2];
  /*--- Memories ---*/
  /* Bark-scale envelope */
  int   p_w_type[ntt_NSclLay_MAX];
  /* LPC spectrum */
  static double
    prev_buf[ntt_NSclLay_MAX][MAX_TIME_CHANNELS][ntt_MA_NP][ntt_N_PR_MAX+1];
  int fw_nbit;
  int ifr;
  int isf, top, subtop, ismp, block_size;

  /*--- Set parameters ---*/
  
  switch(index->w_type){
  case ONLY_LONG_SEQUENCE: case LONG_START_SEQUENCE: case LONG_STOP_SEQUENCE:
    /* subframe */
    nfr = index->block_size_samples;
    nsf = 1;
    n_ch = index->numChannel;
    /* Bark-scale table (Bark-scale envelope) */
    n_crb    = ntt_N_CRB_SCL;
    /* bark_tbl = ntt_crb_tbl_scl[iscl]; */
    
    mat_scale_set_shift_para2(iscl, bark_tbl, index);

    /* quantization (Bark-scale envelope) */
    fw_ndiv = ntt_FW_N_DIV_SCL;
    fw_cb_size = ntt_FW_CB_SIZE_SCL;
    fw_length  = ntt_FW_CB_LEN_SCL;
    fw_len_max = ntt_FW_CB_LEN_SCL;
    fw_codev = index->nttDataScl->ntt_fwcodev_scl;
    fw_prcpt = 0.4;
    /*fw_env_memory = fw_env_memory_long[iscl]; */
    /* envelope calculation (Bark-scale envelope) */
    current_block = ntt_BLOCK_LONG;
    fw_nbit = ntt_FW_N_BIT_SCL;
    break;
 case EIGHT_SHORT_SEQUENCE:
    /* subframe */
    nfr = index->block_size_samples/8;
    nsf = ntt_N_SHRT;
    n_ch = index->numChannel;
    /* Bark-scale table (Bark-scale envelope) */