ntt_tf_proc_spectrum.c

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

  for(idim=0; idim<nsp; idim++){
    index_lsp[idim+2] =
      (tindex[index0+mode_index*NCD]>>(nbit*(nsp-idim-1)))&(csize[kend-1]-1) ;
}

  code1=off_code+index[mode_index][index0];
  for(idim=0; idim<nsp; idim++){
    for(j=ntt_isp[idim]; j<ntt_isp[idim+1]; j++){
      lspq[j+1] =
        (code[code1][j] +code[off_code+off+index_lsp[2+idim]][j]);
  }
}
  if(lspq[1] < L_LIMIT) {
    diff = L_LIMIT - lspq[1];
    lspq[1] += diff*1.2;
}
  if(lspq[n_pr] > M_LIMIT ){
    diff = lspq[n_pr]-M_LIMIT;
    lspq[n_pr] -= diff*1.2;
}

  ntt_check_lsp(n_pr-1, lspq+2, MIN_GAP);

  for(j=off_dim; j<n_pr; j++) buf[j] = lspq[j+1];
  for(j=off_dim; j<n_pr; j++){
    lspq[j+1] = lsp[j+1] +
      (lspq[j+1] - rbuf[mode_index][j] ) *
        fg_sum[mode_index][j] ;
    out_vec[j] = buf[j];
}
  if(lspq[1] < L_LIMIT) {
    diff = L_LIMIT - lspq[1];
    lspq[1] += diff*1.2;
}
  if(lspq[n_pr] > M_LIMIT ){
    diff = lspq[n_pr]-M_LIMIT;
    lspq[n_pr] -= diff*1.2;
}
  for(i=1; i<n_pr; i++){
    if(lspq[i] > lspq[i+1] - MIN_GAP){
      diff = lspq[i] - lspq[i+1];
      lspq[i]   -= (diff+MIN_GAP)/2.;
      lspq[i+1] += (diff+MIN_GAP)/2.;
  }
}
  for(i=1; i<n_pr; i++){
    if(lspq[i] > lspq[i+1] - MIN_GAP*0.95){
      diff = lspq[i] - lspq[i+1];
      /**/        fprintf(stderr,"%12.5f %5d %12.5f PERMUTATION \n",
                          diff, i, lsp[i+1]-lsp[i]);
      /**/
      lspq[i]   -= (diff+MIN_GAP*0.95)/2.;
      lspq[i+1] += (diff+MIN_GAP*0.95)/2.;
  }
}
  ntt_check_lsp_sort(n_pr-1, lspq+2);
}






void ntt_sear_ch_p(/* Parameters */
             int    cb_len,      /* codebook length */
	     int    pol_bits0,   /* number of polarity bits */
	     int    cb_size0,    /* codebook size */
	     int    n_can,       /* number of candidates */
             /* Input */
	     double d_targetv[], /* target subvector (weighted) */
	     double *codevv, /* codebook */
	     double d_wt[],      /* divided weight */
	     /* Output */
	     double can_code_targ[ntt_POLBITS_P][ntt_N_CAN_MAX],
	     double can_code_sign[ntt_POLBITS_P][ntt_N_CAN_MAX],
	     int    can_ind[])
{
    /*--- Variables ---*/
    int i_cb, i_pol, i_can, j_can, ismp, tb_top, spl_pnt, t_ind;
    double cros_max_t[ntt_N_CAN_MAX], dtmp;
    double code_targ[ntt_POLBITS][ntt_PIT_CB_SIZE_MAX], 
    code_targ_tmp[ntt_PIT_CB_SIZE_MAX];
    register double accb, accc;



    for (i_cb=0; i_cb<cb_size0; i_cb++){
       accc = 0.;
       accb = 0.;
      {
       double *p_code, *p_targ, *p_wt;
       p_code= codevv+i_cb*ntt_CB_LEN_P_MAX;
       p_targ=d_targetv;
       p_wt = d_wt;
       for (ismp=0; ismp<cb_len; ismp++){
           accb += *(p_targ++) * (*p_code);
           accc += *(p_wt++)  * *(p_code) * *(p_code++);
       }
      }
      code_targ[0][i_cb] = accb;
      if((pol_bits0 == 1)&&(accb<=0.0)){
           code_targ_tmp[i_cb] = -4.0*accb-accc;
      }
      else{
           code_targ_tmp[i_cb] =  4.0*accb-accc;
      }
    }

    /*--- pre-selection search ---*/
    if (n_can == ntt_PIT_CB_SIZE){ /* full search */
	i_pol=0;
	for( i_can=0; i_can<n_can; i_can++){
	    can_ind[i_can] = i_can;
	    can_code_targ[i_pol][i_can] = code_targ_tmp[i_can];
	    can_code_sign[i_pol][i_can] = 
		code_targ[i_pol][i_can] > 0. ? 1.0 : -1.0;
	}/* end 1st-ch preselection */
    }
    else{
	cros_max_t[0] = -1.e30;
	can_ind[0] = 0;
	tb_top=0;
	/*--- search for the insertion point ---*/
	for (i_cb=0; i_cb<cb_size0; i_cb++){
	    dtmp = code_targ_tmp[i_cb];
	    i_can=tb_top; j_can=0;
	    if (dtmp>cros_max_t[tb_top]){
		i_can=tb_top; j_can=0;
		while(i_can>j_can){
		    spl_pnt = (i_can-j_can)/2 + j_can;
		    if (cros_max_t[spl_pnt]>dtmp){
			j_can = spl_pnt+1;
		    }
		    else{
			i_can = spl_pnt;
		    }
		}
		tb_top = ntt_min(tb_top+1, n_can-1);
		for (j_can=tb_top; j_can>i_can; j_can--){
		    cros_max_t[j_can] = cros_max_t[j_can-1];
		    can_ind[j_can] = can_ind[j_can-1];
		}
		cros_max_t[i_can] = dtmp;
		can_ind[i_can] = i_cb;
	    }
	} /* end 1st-ch preselection */

	
	/*--- Make output ---*/
	i_pol = 0;
	for (i_can=0; i_can<n_can; i_can++){
	    t_ind = can_ind[i_can];
	    can_code_targ[i_pol][i_can] = code_targ_tmp[t_ind];
	    can_code_sign[i_pol][i_can] = 
		code_targ[i_pol][t_ind] > 0. ? 1.0 : -1.0;
	}
    }
}


#define ntt_GAMMA_P  0.05 

void ntt_sear_pitch(
		    /* Input */
		    double tc[],
		    double lpc_spectrum[],  /* LPC spectrum */
		    int   block_size_samples,
		    int   isampf,
		    double bandUpper,
		    /* Output */
		    int    *index_pit)
{
  /*--- Variables ---*/
  int ismp, ii, jj;
  int npcount, ifc;
  int  iscount;
  float pmax,  pitch;
  float dtmp, fcmax, fcmin;
  double ddtmp;
  double tcbuf[ntt_N_FR_MAX];
  float bandwidth;
  float tmpnp0, tmpnp1;
  float fctmp0, fctmp1;

  /*--- Calculate the cepstrum ---*/
  bandwidth =4.;
  if(isampf==8) bandwidth=1.5;
  if(isampf>=11) bandwidth=2.;
  if(isampf>=22) bandwidth=4.;
  if(bandwidth < (float)(1./bandUpper)) bandwidth = (float)(1./bandUpper);/*pitch_q*/	
    fctmp1 = (float)block_size_samples/(float)isampf;
    fctmp0 = 1./(float)log(2.);
    fcmin = (float)log(fctmp1*0.2)*fctmp0;
    fcmax = (float)log(fctmp1*2.4)*fctmp0;
  for(ii=0; ii<block_size_samples; ii++) 
    tcbuf[ii] = tc[ii]*tc[ii]*pow((lpc_spectrum[ii]+0.001), -ntt_GAMMA_P*2.0);

  pmax= -1.0e30;
  for (ifc=0; ifc<(1<<ntt_BASF_BIT); ifc++){
        dtmp = (float)ifc / (float) ntt_BASF_STEP;
        fctmp0 =  fcmax-fcmin;
        dtmp = dtmp * fctmp0;
        dtmp = dtmp  + fcmin;
        pitch = (float)pow(2., dtmp);
        tmpnp0 = pitch*bandwidth;
        tmpnp1 = tmpnp0/(float)block_size_samples;
        tmpnp0 = (float)ntt_N_FR_P*tmpnp1;
        npcount = (int)tmpnp0;
    ddtmp=0.;
    iscount=0;
    for (jj=0; jj<npcount/2; jj++){
      iscount ++;
      ddtmp += tcbuf[jj];
    }
    for (ii=0; ii<(ntt_N_FR_P) && ( iscount<ntt_N_FR_P) ; ii++){
      dtmp = pitch*(float)(ii+1);
      dtmp += 0.5;
      ismp = (int)dtmp;
      if(ismp+(npcount-1)/2+1<block_size_samples){
	for (jj=-npcount/2; (jj<(npcount-1)/2+1) && (iscount<ntt_N_FR_P);
	     jj++, iscount++){
	  ddtmp += tcbuf[ismp+jj];
	}
      }
    }	
    
    if(pmax<ddtmp){
      pmax=ddtmp;
      *index_pit = ifc;
    }
  }

}


void ntt_sear_x_p(/* Parameters */
	    int    i_div,      /* division number */
	    int    cb_len,     /* codebook length */
	    int    pol_bits0,  /* polarity bits (0 ch) */
	    int    pol_bits1,  /* polarity bits (1 ch) */
	    int    can_ind0[], /* candidate indexes (0 ch) */
	    int    n_can0,     /* number of candidates (0 ch) */
	    int    can_ind1[], /* candidate indexes (1 ch) */
	    int    n_can1,     /* number of candidates (1 ch) */
            int    numChannel,
	    /* Output */
	    int    index[],    /* quantization index */
            /* Input */
	    double can_code_targ0[ntt_POLBITS_P][ntt_N_CAN_MAX], 
	    double can_code_targ1[ntt_POLBITS_P][ntt_N_CAN_MAX], 
	    double can_code_sign0[ntt_POLBITS_P][ntt_N_CAN_MAX], 
	    double can_code_sign1[ntt_POLBITS_P][ntt_N_CAN_MAX], 
	    double d_wt[],     /* weight subvector */
            double *pcode,
	    /* Output */
	    double sig_l[],    /* decoded residual */
	    double *dist_min)  /* quantization distortion */
{    
    /*--- Variables ---*/
    int i_pol, i_can, j_can, i_smp, index0, index1, i_can0, i_can1;
    double pol0[ntt_POLBITS_P], pol1[ntt_POLBITS_P];
    double cross, dist;
    double *pt_i, *pt_j;
    double weightbuf[ntt_CB_LEN_P_MAX];
    int n_can0_act, n_can1_act;
    double ttmp;

		i_pol = 0;
    if (pol_bits1 == 1){
	*dist_min = 1.e50;
        n_can0_act = n_can0;
        n_can1_act = n_can1;
	for(i_can=0; i_can<n_can0_act; i_can++){
            ttmp = -can_code_targ0[i_pol][i_can];
	    pt_i = pcode + can_ind0[i_can]*ntt_CB_LEN_P_MAX;
	    for(i_smp=0; i_smp<cb_len; i_smp++){
	        weightbuf[i_smp] = 2. * *(pt_i++) * d_wt[i_smp]
		                 * can_code_sign0[i_pol][i_can];
	    }
	    for(j_can=0; j_can<n_can1_act; j_can++){
		/*--- Calculate the distortion measure ---*/
	        pt_j = pcode +(can_ind1[j_can]+ntt_PIT_CB_SIZE)
			       *ntt_CB_LEN_P_MAX;
		cross =0.0;
	        for(i_smp=0; i_smp<cb_len; i_smp++){
	            cross += weightbuf[i_smp]* *(pt_j++);
	        }
		cross *= can_code_sign1[i_pol][j_can];

        	dist = cross + ttmp - can_code_targ1[i_pol][j_can];

		/*--- Compare the distortion with the minimum ---*/
		if ( dist < *dist_min ){
		    i_can0 = i_can;
		    i_can1 = j_can;
		    *dist_min = dist;
		}
	    } /* j_can */
	}  /* i_can */
	
	/*--- Make output indexes ---*/
	index[i_div] = index0 = can_ind0[i_can0];
	index[i_div+ntt_N_DIV_PperCH*numChannel] = index1 = can_ind1[i_can1];
	if(can_code_sign0[0][i_can0] <= 0) 
	    index[i_div] += ((0x1)<<(ntt_MAXBIT_SHAPE_P));
	if(can_code_sign1[0][i_can1] <= 0) 
	    index[i_div+ntt_N_DIV_PperCH*numChannel] +=  
	     ((0x1)<<(ntt_MAXBIT_SHAPE_P));
	
    }
    /**/
    else {
	*dist_min = 1.e50;
        n_can0_act = n_can0;
        n_can1_act = n_can1;
        if(n_can0> 8){
	  i_can=0;
          while( (can_code_targ0[0][n_can0-1] >
                  can_code_targ1[0][n_can1-1-i_can]) &&(n_can1_act>8) ){
            i_can++;
            n_can1_act --;
          }
          while( (can_code_targ0[0][n_can0-1-i_can] <
                  can_code_targ1[0][n_can1-1]) && (n_can0_act>8) ){
            i_can++;
            n_can0_act --;
          }
        }
        for(i_can=0; i_can<n_can0_act; i_can++){
            ttmp = -can_code_targ0[i_pol][i_can];
	    pt_i = pcode +can_ind0[i_can]*ntt_CB_LEN_P_MAX;
	    for(i_smp=0; i_smp<cb_len; i_smp++){
	        weightbuf[i_smp] = 2. * *(pt_i++) * d_wt[i_smp]
		                 * can_code_sign0[i_pol][i_can];
	    }
	    for(j_can=0; j_can<n_can1_act; j_can++){
		
		/*--- Calculate the distortion measure ---*/
	        pt_j = pcode + (can_ind1[j_can]+ntt_PIT_CB_SIZE)
			       *ntt_CB_LEN_P_MAX;
		cross =0.0;
	        for(i_smp=0; i_smp<cb_len; i_smp++){
	            cross += weightbuf[i_smp]* *(pt_j++);
	        }
        	dist = cross + ttmp - can_code_targ1[i_pol][j_can];

		/*--- Compare the distortion with the minimum ---*/
		if ( dist < *dist_min ){
		    i_can0 = i_can;
		    i_can1 = j_can;
		    *dist_min = dist;
		}
	    } /* j_can */
	}  /* i_can */
	/*--- Make output indexes ---*/
	index[i_div] = index0 = can_ind0[i_can0];
	index[i_div+ntt_N_DIV_PperCH*numChannel] = index1 = can_ind1[i_can1];
	if(pol_bits0 ==1){
	    if(can_code_sign0[0][i_can0] <= 0) 
	    index[i_div] += ((0x1)<<(ntt_MAXBIT_SHAPE_P));
        }
    }
    
    
    /*--- Make the polarity indexes ---*/
    for(i_pol=0; i_pol<ntt_POLBITS_P; i_pol++){
	pol0[i_pol] =1- 2*((index[i_div] >> (i_pol+ntt_MAXBIT_SHAPE_P)) & 0x1);
	pol1[i_pol] =1- 2*((index[i_div+ntt_N_DIV_PperCH*numChannel]>>(i_pol+ntt_MAXBIT_SHAPE_P)) & 0x1);
    }

    /*--- Local decoding ---*/
    for(i_pol=0; i_pol<ntt_POLBITS_P; i_pol++){
	for ( i_smp=i_pol; i_smp<cb_len; i_smp +=ntt_POLBITS_P ){
	    sig_l[i_smp] = 
		(pol0[i_pol]*pcode[index0*ntt_CB_LEN_P_MAX+i_smp]
		 +pol1[i_pol]*
		  pcode[(index1+ntt_PIT_CB_SIZE)*ntt_CB_LEN_P_MAX+i_smp])*0.5;
	}
    }
}


#define EPS  1.e-10

void ntt_wvq_pitch(
	   double targetv[],     /* Input : target vector */
	   double wt[],          /* Input : weighting vector */
	   double pwt[],         /* Input : perceptual controlled weight */
	   double pgain[],       /* Input : gain */
	   double *pcode,
	   short *pleave0,
	   short *pleave1,
	   int    numChannel,
	   int    index_pls[],   /* Output: quantization indexes */
	   double pgainq[],      /* In/Out : gain */
	   int    index_pgain[]) /* Output: gain index */
{
    /*--- Variables ---*/
    int		i_div;
    double	d_targetv[ntt_CB_LEN_P_MAX], d_wt[ntt_CB_LEN_P_MAX];
    double	dist_min;
    int         can_ind0[ntt_N_CAN_MAX], can_ind1[ntt_N_CAN_MAX];
    double      can_code_targ0[ntt_POLBITS_P][ntt_N_CAN_MAX], 
		can_code_targ1[ntt_POLBITS_P][ntt_N_CAN_MAX];
    double      can_code_sign0[ntt_POLBITS_P][ntt_N_CAN_MAX];
    double      can_code_sign1[ntt_POLBITS_P][ntt_N_CAN_MAX];
    static double      gain2[MAX_TIME_CHANNELS];
    static int bits_p[ntt_N_DIV_P_MAX], length_p[ntt_N_DIV_P_MAX];
    int         cb_len, cb_size0, cb_size1;
    double   sig_l[ntt_CB_LEN_P_MAX];
    double   nume[MAX_TIME_CHANNELS], denom[MAX_TIME_CHANNELS];
    double   g_opt[MAX_TIME_CHANNELS];
    int      pol_bits0, pol_bits1;
    double   nume0[ntt_N_DIV_P_MAX], denom0[ntt_N_DIV_P_MAX];
    double   powG[MAX_TIME_CHANNELS];
    int      i_sup;
    double   acc_tmp;
    static double ntt_codevp1[ntt_PIT_CB_SIZE_MAX][ntt_CB_LEN_P_MAX];
    double   *ptr;
    int index_pgain0;

    ntt_vec_lenp(numChannel, bits_p, length_p);
    for(i_sup=0; i_sup<numChannel; i_sup++){ 
	nume[i_sup] = denom[i_sup] = 0.;
	gain2[i_sup] = pgain[i_sup]*pgain[i_sup];
	powG[i_sup] = 0.0;
    }

    /*--- Main operation ---*/
    acc_tmp=0.0;
    for ( i_div=0; i_div<ntt_N_DIV_PperCH*numChannel; i_div++ ){
	nume0[i_div]= denom0[i_div] =0.0;
	/*--- set codebook lengths and sizes ---*/
	cb_len = length_p[i_div];
	cb_size0 = 0x1 << ntt_MAXBIT_SHAPE_P;
	cb_size1 = 0x1 << ntt_MAXBIT_SHAPE_P;
	pol_bits0 = ((bits_p[i_div]+1)/2)-ntt_MAXBIT_SHAPE_P;
	pol_bits1 = ((bits_p[i_div])/2)-ntt_MAXBIT_SHAPE_P;

	/*--- divide vectors into subvectors ---*/
	ntt_pre_dot_p(cb_len, i_div, targetv, d_targetv,
		wt, pwt, pgain, powG, d_wt, 
		numChannel, pleave0, pleave1);

	/*--- Best codebook search ---*/
	/* pre-selection 0 */
        ntt_sear_ch_p(cb_len,pol_bits0, cb_size0, ntt_N_CAN_P, d_targetv, 
	       pcode, d_wt,
	       can_code_targ0, can_code_sign0,
	       can_ind0);
	/* pre-selection 1 */
        ntt_sear_ch_p(cb_len,pol_bits1, cb_size1, ntt_N_CAN_P, d_targetv, 
	       pcode+ntt_PIT_CB_SIZE*ntt_CB_LEN_P_MAX, d_wt,
	       can_code_targ1, can_code_sign1,
	       can_ind1);
	/* main selection */
	ntt_sear_x_p(i_div, cb_len, pol_bits0, pol_bits1,
	       can_ind0, ntt_N_CAN_P, can_ind1, ntt_N_CAN_P, 
	       numChannel, index_pls,
	       can_code_targ0, can_code_targ1,
	       can_code_sign0, can_code_sign1, 
	       d_wt, pcode, sig_l, &dist_min);
	acc_tmp += dist_min;
	/* Prepair for calculating the optimum gain */
        ntt_opt_gain_p(i_div,cb_len, numChannel, 
		 pgain,d_targetv,sig_l,d_wt,nume,
                 pgainq, denom, nume0, denom0, pleave0, pleave1);
    }

    /*--- gain re-quantization ---*/
    for(i_sup=0; i_sup<numChannel; i_sup++){
	g_opt[i_sup] = nume[i_sup]/(denom[i_sup]+EPS); /* calc. opt. gain */
	/* re-quantization */
	ntt_enc_gair_p(&index_pgain0,
		 powG[i_sup], gain2[i_sup], &g_opt[i_sup],
		 nume[i_sup], denom[i_sup]);
	index_pgain[i_sup] = index_pgain0;
        ntt_dec_pgain(index_pgain[i_sup], &pgainq[i_sup]); 
    } /* end i_sup */
}