ntt_tf_proc_spectrum.c

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

/*****************************************************************************/
/* This software module was originally developed by                          */
/*   Naoki Iwakami (NTT)                                                     */
/* and edited by                                                             */
/*   Naoki Iwakami (NTT) on 1997-04-18,                                      */
/*   Akio Jin (NTT) 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.                                                        */
/*****************************************************************************/


/* 25-aug-1997  NI  added bandwidth control */
/* 01-oct-1997  HP  "iscount++" bug fix in line 627 & 631 */
/*                  (as reported by T. Moriya, NTT) */
/* 14-oct-1997  TM  mods by Takehiro Moriya (NTT) */


#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_bwp_enc.h"         /* structs */
#include "nok_ltp_common.h"      /* structs */
#include "tf_mainStruct.h"       /* structs */
#include "tns.h"                 /* structs */

#include "ntt_conf.h"
#include "ntt_encode.h"
#include "ntt_tools.h"
#include  "ntt_relsp.h"
#include "ntt_enc_para.h"

void ntt_tf_proc_spectrum(double spectrum[],       /* spectrum */
			  ntt_INDEX  *index,       /* code indices*/
			  /* Output */
			  double lpc_spectrum[],   /* LPC spectrum */
			  double bark_env[],       /* Bark-scale envelope*/
			  double pitch_sequence[],/*periodic peak components*/
			  double gain[])          /* gain factor */
{
  double tc[ntt_T_FR_MAX], pwt[ntt_T_FR_MAX];
  int    current_block, i_ch, i_ma, top;
  static int InitFlag = 1;
  
  /*--- Parameters ---*/
  double band_lower, band_upper;
  int    nfr, nsf, n_ch;
  int    n_crb, *bark_tbl;
  int    fw_ndiv, fw_length, fw_len_max, fw_cb_size;
  double fw_prcpt;
  double *fw_codev;
  int    n_pr;
  int    ifr, isf, subtop, ismp, block_size;
  double *lsp_code;
  double *lsp_fgcode;
  int    lsp_csize[2], lsp_cdim[2];
  int    flag_blim_l;
  int    p_w_type;
  /*--- Memories ---*/
  /* Bark-scale envelope */
  /* LPC spectrum */
  static double prev_buf[MAX_TIME_CHANNELS][ntt_MA_NP][ntt_N_PR_MAX+1];
  static int counter;
  /*--- 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;
    bark_tbl = index->nttDataBase->ntt_crb_tbl;
    /* quantization (Bark-scale envelope) */
    fw_ndiv = ntt_FW_N_DIV;
    fw_cb_size = ntt_FW_CB_SIZE;
    fw_length = ntt_FW_CB_LEN;
    fw_codev = index->nttDataBase->ntt_fwcodev;
    fw_len_max = ntt_FW_CB_LEN;
    fw_prcpt = 0.4;
    /* envelope calculation (Bark-scale envelope) */
    current_block = ntt_BLOCK_LONG;
    /* bandwidth control */
    flag_blim_l = 0;
    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) */
    fw_cb_size = 1;
    current_block = ntt_BLOCK_SHORT;
    /* bandwidth control */
    flag_blim_l = 1;
    break;
  default:
        fprintf( stderr,"Fatal error! %d: no such window type.\n", index->w_type );
        exit(1);
  }
  block_size = nfr*nsf;

  /* LPC spectrum */
  n_pr = ntt_N_PR;
  lsp_code   = index->nttDataBase->lsp_code_base;
  lsp_fgcode = index->nttDataBase->lsp_fgcode_base;
  lsp_csize[0] = ntt_NC0;
  lsp_csize[1] = ntt_NC1;
  lsp_cdim[0] = ntt_N_PR;
  lsp_cdim[1] = ntt_N_PR;

    p_w_type = index->nttDataBase->prev_w_type;
  
  /*--- Encoding tools ---*/
  /* LPC spectrum encoding */
  /** new version moriya informed **/
       ntt_scalebase_enc_lpc_spectrum(
                            nfr, nsf, 
			    index->block_size_samples,
			    n_ch, n_pr, lsp_code, lsp_fgcode,
                            lsp_csize, lsp_cdim, 
			    index->nttDataBase->prev_lsp_code,
			    /* prev_buf, */
			    index->nttDataBase->ma_np, 
                            spectrum, index->w_type, current_block,
                            index->lsp, lpc_spectrum,
                            index->nttDataBase->bandLower, 
			    index->nttDataBase->bandUpper,
			    index->nttDataBase->cos_TT);

 /* bandwidth control */
  flag_blim_l =0;
  for (i_ch=0; i_ch<n_ch; i_ch++){
    top = i_ch * nsf * nfr;
    if ( index->bandlimit == 1){
      index->blim_h[i_ch] = ntt_BLIM_STEP_H/2;
      band_upper = 
	index->nttDataBase->bandUpper * 
	  (1. - (1.-ntt_CUT_M_H) *
	   (double)index->blim_h[i_ch]/(double)ntt_BLIM_STEP_H);
    }
    else{
      index->blim_h[i_ch] = 0;
      band_upper = index->nttDataBase->bandUpper;
    }
    if (( index->bandlimit == 1)) index->blim_l[i_ch] = 1;
    else                           index->blim_l[i_ch] = 0;
    band_lower = index->nttDataBase->bandLower 
	       + ntt_CUT_M_L * index->blim_l[i_ch];
    ntt_freq_domain_pre_process(nfr, nsf, band_lower, band_upper,
				spectrum+top, lpc_spectrum+top,
				spectrum+top, lpc_spectrum+top);
  }
  /* periodic peak components encoding */
  index->pitch_q =0;
  ntt_zerod(index->block_size_samples*index->numChannel, pitch_sequence);

  if(tvq_debug_level>5)
   fprintf(stderr, "UUUUU ntt_tf_proc_spectrum tvq_ppc_enable %5d \n",
   index->tvq_ppc_enable);

if(index->tvq_ppc_enable){
    if((index->nttDataBase->ntt_VQTOOL_BITS > ntt_PIT_TBITperCH*index->numChannel*3)){
    ntt_enc_pitch(spectrum, lpc_spectrum, index->w_type,
		  index->numChannel,
		  index->block_size_samples,
		  index->isampf,
		  index->nttDataBase->bandUpper,
		  index->nttDataBase->ntt_codevp0,
		  index->nttDataBase->pleave0,
		  index->nttDataBase->pleave1,
		  index->pit, index->pls, index->pgain, pitch_sequence);
{ /*pitch_q*/
 double tmp[ntt_N_FR_MAX*MAX_TIME_CHANNELS];
 double dis1, dis2;
 for (i_ch=0; i_ch<n_ch; i_ch++){
    top = i_ch * block_size;
    for(isf=0; isf<nsf; isf++){
      subtop = top + isf*nfr;
      /* make input data */
      for (ismp=0; ismp<nfr; ismp++){
       tmp[ismp+subtop] =
         spectrum[ismp+subtop] * lpc_spectrum[ismp+subtop]
           - pitch_sequence[ismp+subtop];
     }
    }
 }
      dis1 = 0.0;
      dis2 = 0.0;
 for (i_ch=0; i_ch<n_ch; i_ch++){
    top = i_ch * block_size;
    for(isf=0; isf<nsf; isf++){
      subtop = top + isf*nfr;
      for (ismp=0; ismp<nfr; ismp++){
       dis1 += tmp[ismp+subtop]*tmp[ismp+subtop];
       dis2 += pitch_sequence[ismp+subtop]*pitch_sequence[ismp+subtop];
      } 
    }
 }
 if(((dis1+dis2)/(index->nttDataBase->ntt_VQTOOL_BITS)
    <  dis2/ (ntt_PIT_TBITperCH*index->numChannel)) && 
    (index->nttDataBase->ntt_VQTOOL_BITS > ntt_PIT_TBITperCH*index->numChannel*3) ) index->pitch_q = 1;
 else {
    index->pitch_q =0;
    ntt_zerod(index->block_size_samples*index->numChannel, pitch_sequence);
 }

if(tvq_debug_level>5)
  fprintf(stderr, "%5d %5d %5d %12.3f %12.5f %5d YYYYY\n", index->pitch_q, ntt_PIT_TBITperCH*index->numChannel, index->nttDataBase->ntt_VQTOOL_BITS, 
dis1+dis2, dis2, index->w_type);

}  /* pitch_q*/ 
/**UUVV    index->pitch_q =0; */ 
/**UUVV    ntt_zerod(index->block_size_samples*index->numChannel, pitch_sequence);*/ 
  }
  else{
    ntt_zerod(index->block_size_samples*index->numChannel, pitch_sequence);
  }
} /* if flag_ppc */

/* Bark-scale envelope encoding */
 if(fw_cb_size==1){
    for(ifr=0; ifr<nfr*nsf*n_ch; ifr++){
       bark_env[ifr]=1.0;
    }
    for (i_ch=0; i_ch<n_ch; i_ch++){
    top = i_ch * block_size;
    for(isf=0; isf<nsf; isf++){
      subtop = top + isf*nfr;
      /* make input data */
      for (ismp=0; ismp<nfr; ismp++){
       tc[ismp+subtop] =
         spectrum[ismp+subtop] * lpc_spectrum[ismp+subtop]
           - pitch_sequence[ismp+subtop];
     }
    }
  }
 }
 else{
    ntt_enc_bark_env(nfr, nsf, n_ch, n_crb, bark_tbl,
                   index->nttDataBase->bandUpper, 
		   fw_ndiv, fw_cb_size, fw_length,
                   fw_codev, fw_len_max, fw_prcpt,
                   index->nttDataBase->prev_fw_code,
		   /*fw_env_memory, fw_pdot, */
		   &p_w_type,
                   spectrum, lpc_spectrum, pitch_sequence, current_block,
                   tc, pwt, index->fw, index->fw_alf, bark_env);
    if(index->w_type == EIGHT_SHORT_SEQUENCE)
	index->nttDataBase->prev_w_type = ntt_BLOCK_SHORT;
    else
	index->nttDataBase->prev_w_type = ntt_BLOCK_LONG;
  }
  /* gain encoding */
  ntt_enc_gain_t(nfr, nsf, n_ch, bark_env, tc, index->pow, gain);
}


void ntt_enc_bark_env(/* Parameters */
		      int    nfr,               /* block length */
		      int    nsf,               /* number of sub frames */
		      int    n_ch,
		      int    n_crb,
		      int    *bark_tbl,
		      double bandUpper,
		      int    ndiv,              /* number of interleave division */
		      int    cb_size,           /* codebook size */
		      int    length,
		      double *codev,            /* code book */
		      int    len_max,           /* codevector memory length */
		      double prcpt,
		      int    prev_fw_code[],
		      /*
		      double *env_memory,
		      double *pdot,
		      */
		      int    *p_w_type,
		      /* Input */
		      double spectrum[],        /* spectrum */
		      double lpc_spectrum[],    /* LPC spectrum */
		      double pitch_sequence[],  /* periodic peak components */
		      int    current_block,     /* block length type */
		      /* Output */
		      double tc[ntt_T_FR_MAX],  /* flattened spectrum */
		      double pwt[ntt_T_FR_MAX], /* perceptual weight */
		      int    index_fw[],        /* indices for Bark-scale enveope coding */
		      int    ind_fw_alf[],
		      double bark_env[])        /* Bark-scale envelope */
{
  /*--- Variables ---*/
  int    ismp, i_ch, top, isf, subtop, pred_sw, block_size;
  /*--- Initialization ---*/
  block_size = nfr * nsf;
  /*--- Encoding process ---*/
  for (i_ch=0; i_ch<n_ch; i_ch++){
    top = i_ch * block_size;
    /* make weight */
    for (ismp=0; ismp<block_size; ismp++){
      pwt[ismp+top] = 1. / pow(lpc_spectrum[ismp+top], 0.5);
    }
    /* set MA prediction switch */
    if (*p_w_type == current_block) pred_sw = 1;
    else                           pred_sw = 0;
    
    for(isf=0; isf<nsf; isf++){
      subtop = top + isf*nfr;
      /* make input data */
      for (ismp=0; ismp<nfr; ismp++){
	tc[ismp+subtop] =
	  spectrum[ismp+subtop] * lpc_spectrum[ismp+subtop]
	    - pitch_sequence[ismp+subtop];
      }
      /* envelope coding */
      ntt_fwpred(nfr, n_crb, bark_tbl, bandUpper, ndiv,
		 cb_size, length, codev, len_max, prcpt,
		 prev_fw_code+(i_ch*nsf+isf)*ndiv,
		 /* env_memory, pdot, */
		 pwt+subtop, bark_env+subtop, tc+subtop,
		 &index_fw[(i_ch*nsf+isf)*ndiv],
		 &ind_fw_alf[isf+i_ch*nsf], i_ch, pred_sw);
      
      pred_sw = 1;
    }
  }
  
  *p_w_type = current_block;
}

void ntt_enc_gaim(/* Input */
		  double power,      /* Power to be encoded */
		  /* Output */
		  int    *index_pow, /* Power code index */
		  double *gain,      /* Gain */
		  double *norm)      /* Power normalize factor */
{
  double mu_power, dmu_power, dec_power;
  
  /* encoding */
  mu_power = ntt_mulaw(power, ntt_SUB_AMP_MAX, ntt_MU);  /* mu-Law power */
  *index_pow = (int)floor((mu_power)/ ntt_SUB_STEP); /* quantization */

  /* local decoding */
  dmu_power = *index_pow * ntt_SUB_STEP +ntt_SUB_STEP/2.;
  dec_power = ntt_mulawinv(dmu_power, ntt_SUB_AMP_MAX, ntt_MU);
  *gain= dec_power/ ntt_SUB_AMP_NM;                       /* calculate gain */
  *norm = ntt_SUB_AMP_NM/(dec_power+0.1);                 /* calculate norm. */
}

void ntt_enc_gain(/* Input */
		  double power,      /* Power to be encoded */
		  /* Output */
		  int    *index_pow, /* Power code index */
		  double *gain,      /* Gain */
		  double *norm)      /* Power normalize factor */
{
    /*--- Variables --*/
    double mu_power, dmu_power, dec_power;

    mu_power = ntt_mulaw(power, ntt_AMP_MAX, ntt_MU);           /* mu-Law power */

    *index_pow = (int)floor((mu_power)/ ntt_STEP);  /* quantization */

    dmu_power = *index_pow * ntt_STEP + ntt_STEP/2.;        /* decoded mu-Law power */
    dec_power = ntt_mulawinv(dmu_power, ntt_AMP_MAX, ntt_MU);   /* decoded power */

    *gain= dec_power * ntt_I_AMP_NM;                    /* calculate gain */
    *norm = ntt_AMP_NM/(dec_power+0.1);                 /* calculate norm. */
}

void ntt_enc_gain_t(/* Parameters */
		    int    nfr,          /* block length */
		    int    nsf,          /* number of sub frames */
		    int    n_ch,
		    /* Input */
		    double bark_env[],   /* Bark-scale envelope */
                    /* In/Out */
		    double tc[],         /* flattened spectrum */
		    /* Output */ 
		    int    index_pow[],  /* gain code indices */
		    double gain[])       /* gain factor */
{
  int    ismp, top, subtop, isf, iptop, gtop, i_ch, block_size;
  double norm, pwr;
  double sub_pwr_tmp, sub_pwr[ntt_N_SHRT_MAX];
  double g_gain, g_norm;

  block_size = nfr * nsf;
  
  /*--- Encoding process ---*/
  for (i_ch=0; i_ch<n_ch; i_ch++){
    top = i_ch * block_size;
    /*--- Gain calculation ---*/
    pwr = 0.1;
    for(isf=0; isf<nsf; isf++){
      subtop = top + isf*nfr;
      for ( ismp=0; ismp<nfr; ismp++ ){
	tc[ismp+subtop] /= bark_env[ismp+subtop];
      }
      sub_pwr_tmp
	=(ntt_dotdd(nfr, tc+subtop, tc+subtop)+0.1)/(double)nfr;
      pwr += sub_pwr_tmp;
      sub_pwr[isf] =sqrt( sub_pwr_tmp );
    }
    pwr = sqrt(pwr/(double)nsf);

    /*--- Quantization ---*/
    if (isf == 1){
      /* global gain */
      ntt_enc_gain(pwr,&index_pow[i_ch],&gain[i_ch],&g_norm);
    }
    else{
      /* global gain */
      iptop = (nsf+1)*i_ch;
      gtop  = nsf*i_ch;
      ntt_enc_gain(pwr,&index_pow[iptop],&g_gain,&g_norm);
      gain[gtop] = g_gain;

      /* subband gain */
      for(isf=0; isf<nsf; isf++){
	subtop = top + isf*nfr;
	ntt_enc_gaim(sub_pwr[isf]*g_norm, &index_pow[iptop+isf+1],
		     &gain[gtop+isf], &norm);
	gain[gtop+isf] *= g_gain;
	norm *= g_norm;
      }
    }
  }
}

void ntt_enc_gair_p(int  *index_pgain, /* Output: Power code index */
	      double powG,        /* Input : power of the weighted target */
	      double gain2,       /* Input : squared gain of the target */
	      double *g_opt,      /* In/Out: optimum gain */                
	      double nume,/*Input: cross power or the orig. and quant. resid.*/
	      double denom)       /* Input : power of the orig. regid */
{
    double norm;
    double acc2_2, g_opt_2, acc2;
    int index_tmp;
    double pgain_step;
    
    pgain_step   = ntt_PGAIN_MAX / ntt_PGAIN_NSTEP;

    /*--- gain re-quantization ---*/
    ntt_enc_pgain(*g_opt, &index_tmp, g_opt, &norm);
    *index_pgain = index_tmp;    
    acc2= powG*gain2 -2.* *g_opt*nume + *g_opt* *g_opt*denom;
    if(index_tmp != 0){
	g_opt_2= ntt_mulawinv((index_tmp-1)*pgain_step+pgain_step/2., ntt_PGAIN_MAX, ntt_PGAIN_MU);
	acc2_2= powG*gain2 -2.*g_opt_2*nume+g_opt_2*g_opt_2*denom;
	if(acc2_2 < acc2){
	    acc2= acc2_2;
	    *index_pgain = ntt_max(index_tmp -1, 0);
	    *g_opt= g_opt_2;
	    return; 
	}
    }
}