ntt_tf_proc_spectrum.c

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

  /*--- Variables ---*/
  int    j, vqword;
  double wegt[ntt_N_PR_MAX+1];
  double flsp[ntt_N_PR_MAX+1];
  double target_buf[ntt_N_MODE_MAX][ntt_N_PR_MAX];
  double fg_sum[ntt_N_MODE_MAX][ntt_N_PR_MAX];
  double out_vec[ntt_N_PR_MAX];
  int    i_mode, i_ma;
  int ntt_isp[ntt_NSP_MAX+1];
  double  pred_vec[ntt_N_PR_MAX], buf_prev[ntt_N_PR_MAX+1];
  double  lspq[ntt_N_PR_MAX+1];

/*--- Initialization ---*/
  /* set the band-split table */
  ntt_set_isp(nsp, n_pr, ntt_isp);

  /*--- Make the codebook with the prediction ---*/
  for(i_mode=0; i_mode<ntt_N_MODE; i_mode++){
    for(j=0; j<n_pr; j++){
      fg_sum[i_mode][j] =1.0;
      if(ma_np==1) fg_sum[i_mode][j] -= fgcode[i_mode][0][j];
    }
  }
  for(j=0; j<n_pr; j++){ pred_vec[j] = 0.0; }

  ntt_redec(n_pr, prev_lsp_code, csize,  nsp, code,
		      fg_sum[0], pred_vec, lspq, buf_prev );

  ntt_movdd( n_pr, freq+1, flsp+1 );
  ntt_get_wegt(n_pr, flsp, wegt );
  for(j=1; j<nsp; j++) {
    wegt[ntt_isp[j]] *= 1.2;
    wegt[ntt_isp[j]+1] *= 1.2;
}

  for (i_mode=0; i_mode< ntt_N_MODE; i_mode++){
    for(j=0; j<n_pr; j++) target_buf[i_mode][j] = flsp[j+1];
    if(ma_np==1){
      for(j=0; j<n_pr; j++){
        target_buf[i_mode][j] -= fgcode[i_mode][0][j]*buf_prev[j];
      }
    }
    for(j=0; j<n_pr; j++)
      target_buf[i_mode][j] /= fg_sum[i_mode][j];
  }
  ntt_relspwed(n_pr, flsp, wegt, freqout, ntt_LSP_NSTAGE,
               csize, code, fg_sum,
               target_buf, out_vec, &vqword , index, nsp );
}


void ntt_lsptowt(/* Parameters */
		 int nfr,        /* block length */
		 int n_pr,
		 int block_size_samples,
		 /* Input */
		 double lsp[],
		 /* Output */
		 double wt[],
		 double *cos_TT)
{
    register double  a,b,c,d,cosz2;
    double sa, sb,sc,sd;
    double coslsp[ntt_N_PR_MAX];
    int i, mag;
    double *lspp, *wp1, *wp2, *cp;

    mag =2*block_size_samples/nfr;
    for(i=0; i<n_pr; i++) coslsp[i] =  2.* cos(lsp[i+1]);

    wp1= wt;
    wp2= wt+nfr-1;
    cp = cos_TT+mag;
    switch(n_pr){
    case 20 :
	 for(i=0;i<nfr/2;i++) { 
           lspp = coslsp;
	   cosz2=2. * *cp;
	   a = (cosz2 - *lspp);
	   c = (cosz2 + *(lspp++));
	   b = (cosz2 - *lspp);
	   d = (cosz2 + *(lspp++));
	 { a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++)); }
	 { a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++)); }
	 { a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++)); }
	 { a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++)); }
	 { a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++)); }
	 { a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++)); }
	 { a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++)); }
	 { a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++)); }
	 { a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++)); }
	 sa=a*a; sb=b*b; sc=c*c; sd=d*d;
	 *(wp1++) = sa+sb +  (sa-sb)* *cp;
	 *(wp2--) = sc+sd + (-sc+sd)* *cp;
	 cp += mag*2;
      }
    break;
    case 16 :
	 for(i=0;i<nfr/2;i++) { 
           lspp = coslsp;
	   cosz2=2. * *cp;
	   a = (cosz2 - *lspp);
	   c = (cosz2 + *(lspp++));
	   b = (cosz2 - *lspp);
	   d = (cosz2 + *(lspp++));
	 { a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++)); }
	 { a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++)); }
	 { a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++)); }
	 { a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++)); }
	 { a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++)); }
	 { a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++)); }
	 { a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++)); }
	 sa=a*a; sb=b*b; sc=c*c; sd=d*d;
	 *(wp1++) = sa+sb +  (sa-sb)* *cp;
	 *(wp2--) = sc+sd + (-sc+sd)* *cp;
	 cp += mag*2;
      }
    break;
    default:
	 for(i=0;i<nfr/2;i++) { 
           lspp = coslsp;
	   cosz2=2. * *cp;
	   a = (cosz2 - *lspp);
	   c = (cosz2 + *(lspp++));
	   b = (cosz2 - *lspp);
	   d = (cosz2 + *(lspp++));
         while(lspp<coslsp+n_pr)     /* n_pr should be even */ 
	 {
	   a *= (cosz2 - *lspp);
	   c *= (cosz2 + *(lspp++));
	   b *= (cosz2 - *lspp);
	   d *= (cosz2 + *(lspp++));
	 }
	 sa=a*a; sb=b*b; sc=c*c; sd=d*d;
	 *(wp1++) = sa+sb +  (sa-sb)* *cp;
	 *(wp2--) = sc+sd + (-sc+sd)* *cp;
	 cp += mag*2;
     }
   }
}



/***************************************/
/* quantization in ntt_mulaw-scale         */
/* coded by K. Mano    28/Jul/1988     */
/***************************************/

#define DISP_OVF 0

double ntt_mulaw(double x, double xmax,double mu) /* Input */
{
   double  y;
   if (x>xmax){
#if (DISP_OVF)
       fprintf(stderr, "ntt_mulaw(): %lf, %lf: Input overflow.\n", x, xmax);
#endif
       x=xmax;
   }
#if 0
   else if (x< -xmax){
#if (DISP_OVF)
       fprintf(stderr, "ntt_mulaw(): %lf, %lf: Input underflow.\n", x, -xmax);
#endif
       x= -xmax;
   }
#endif
   else if ( x <=0.0){
       x = 0.0; 
   }
   y = xmax*log10(1.+mu*fabs(x)/xmax)/log10(1.+mu); 
   if (x<0.) y= -y;
   return(y);
}


void ntt_opt_gain_p(
	      int    i_div,    /* Input : division number */
	      int    cb_len,   /* Param. : codebook length */
	      int    numChannel,
	      double gt[],     /* Input : gain */
	      double d_targetv[], /* Input : target vector */
	      double sig_l[],  /* Input : decoded residual */
	      double d_wt[],   /* Input : weight */
	      double nume[],   /* Output: cross power (total)*/
	      double gain[],   /* Input : gain of the residual */
	      double denom[],  /* Output: decoded resid. power without gain */
	      double nume0[],  /* Output: cross power (divided) */
	      double denom0[], /* Output: decoded resid. power without gain */
              short *pleave0, /*ntt_SMP_ACCp*/
              short *pleave1) /*ntt_BIT_REVp*/

{
    double fact[MAX_TIME_CHANNELS], acc;
    int i_smp, point, points, serial;
     
    /*serial = ntt_SMP_ACCp[i_div];  */
      serial = *(pleave0+i_div); 
   if( numChannel ==1){
	points = 0;
	for(i_smp=0; i_smp<cb_len; i_smp++){
	    acc =  gain[points]*d_targetv[i_smp]*sig_l[i_smp] ;
	    nume[points] +=  acc;
	    nume0[i_div] +=  acc;
	    acc = sig_l[i_smp]*sig_l[i_smp]*d_wt[i_smp];
	    denom[points] += acc;
	    denom0[i_div] += acc;
	}
    }
    else{
	for(points=0; points<numChannel; points++)
	    fact[points] = 1./(gt[points])/(gt[points]);
	for(i_smp=0; i_smp<cb_len; i_smp++){
	    /* point = ntt_BIT_REVp[serial++]; */
	    point = pleave1[serial++];
	    points = (point)/ntt_N_FR_P;
	    acc =  gain[points]*d_targetv[i_smp]*sig_l[i_smp]*fact[points] ;
	    nume[points] +=  acc;
	    nume0[i_div] +=  acc;
	    acc = sig_l[i_smp]*sig_l[i_smp]*d_wt[i_smp]*fact[points];
	    denom[points] += acc;
	    denom0[i_div] += acc;
	}
    }
}


void ntt_pre_dot_p(int    cb_len,      /* Param. : codebook length */
	       int    i_div,       /* Input : division number */
	       double targetv[],   /* Input : target vector */
	       double d_targetv[], /* Output: target subvector (weighted) */
	       double wt[],        /* Input : weighting vector */
	       double mwt[],       /* Input : perceptual controlled weight */
	       double gt[],        /* Input : gain */
	       double powG[],      /* Output: power of the weighted target */
	       double d_wt[],      /* Output: weighting subvector */
               int   numChannel,
	       short  *pleave0,
               short  *pleave1)

{
    int   i_smp, points, point, serial;
    register double acc, accb;
    
    /* Make divided target and weighting vector */
      serial = *(pleave0+i_div); 
    if(numChannel==1){
	points = 0;
	for ( i_smp=0; i_smp<cb_len; i_smp++ ) {
	    point = pleave1[serial++]; 
	    accb = mwt[point];
	    acc = targetv[point]* (wt[point] * accb);
	    accb *= gt[points];
	    d_targetv[i_smp] = acc *gt[points]*accb;
	    d_wt[i_smp] = accb*accb;
	    powG[points] += acc* acc;
	}
    }
    else{
	for ( i_smp=0; i_smp<cb_len; i_smp++ ) {
	    /*point = ntt_BIT_REVp[serial++]; */
	    point = pleave1[serial++]; 
	    points = (point) / ntt_N_FR_P;
	    accb = mwt[point];
	    acc = targetv[point]* (wt[point] * accb);
	    accb *= gt[points];
	    d_targetv[i_smp] = acc *gt[points]*accb;
	    d_wt[i_smp] = accb*accb;
	    powG[points] += acc* acc;
	}
    }
}


void ntt_relspwed(int    n_pr,
                  double lsp[ntt_N_PR_MAX+1],
                  double wegt[ntt_N_PR_MAX+1],
                  double lspq[ntt_N_PR_MAX+1],
                  int    nstage,
                  int    csize[],
                  double code[][ntt_N_PR_MAX],
                  double fg_sum[ntt_N_MODE_MAX][ntt_N_PR_MAX],
                  double target_buf[ntt_N_MODE_MAX][ntt_N_PR_MAX],
                  double out_vec[ntt_N_PR_MAX],
                  int    *vqword,
                  int    index_lsp[],
                  int    nsp)
{
  int      i, j, k, off_code, off_dim, off, code1, code2;
  double   dist, dmin, diff;
  int   kst, kend, index[ntt_N_MODE_MAX][NCD], index0;
  int   nbit0, nbit, k1;
  double dmin1, dmin2, tdist[NCD*ntt_N_MODE_MAX];
  int  tindex[NCD*ntt_N_MODE_MAX], mode_index, mode;
  double bufsave, rbuf[ntt_N_MODE_MAX][ntt_N_PR_MAX];
  double buf[ntt_N_PR_MAX], d[ntt_NC0_MAX];
  double lspbuf[ntt_N_PR_MAX+1];
  int i_mode, idim;
  int index_sp[ntt_NSP_MAX];
  register double acc, acc1;
  double *p_code, *p_rbuf;
/*
  int ntt_isp[ntt_NSP_MAX];
*/
  int ntt_isp[ntt_NSP_MAX+1];

  ntt_set_isp(nsp, n_pr, ntt_isp);
  for(i_mode=0; i_mode<ntt_N_MODE; i_mode++)
    ntt_movdd(n_pr, target_buf[i_mode], rbuf[i_mode]);
  off_code = 0;
  kst = 0;
  kend = nstage;
  off_dim = 0;
  off = csize[kst];
  nbit0 = (int)floor(log((double)csize[kst])/log(2.0));
  nbit  = (int)floor(log((double)csize[kend-1])/log(2.0) + 0.00001);
  for(mode = 0; mode<ntt_N_MODE; mode++){
    for(i=off_code; i<csize[kst]+off_code; i++){
      acc =0.0;
      acc1 =0.0;
      p_code = code[i];
      p_rbuf = rbuf[mode];
      for(j=off_dim; j<n_pr;j++){
        acc += *(p_rbuf++) * *(p_code);
        acc1+= *(p_code) * *(p_code++);
    }
      d[i-off_code] =  -acc*2.0 + acc1;

  }
    for(j=0; j<NCD; j++){
      for(i=1, index[mode][j]=0, dmin= d[0]; i<csize[kst]; i++){
        if(d[i]<dmin){
          dmin=d[i];
          index[mode][j]=i;
      }
    }
      d[index[mode][j]]=1.e38;
  }
    for(k=0; k<NCD; k++){
      code1=off_code+index[mode][k];
      for(dmin1=1.e38, k1 = 0; k1<csize[kend-1]; k1++){
        code2=off_code+k1+off;
        for(j=ntt_isp[0]; j<ntt_isp[1]; j++){
          buf[j]=(code[code1][j]+code[code2][j]);
      }
        ntt_dist_lsp(ntt_isp[1], rbuf[mode], buf, wegt, &dist);
        if(dist<dmin1) { dmin1 = dist; index_sp[0] = k1; }
        for(j=ntt_isp[0]; j<ntt_isp[1]; j++){
          buf[j]=(code[code1][j]+code[off_code+index_sp[0]+off][j]);
      }
    }
      bufsave = buf[ntt_isp[1]-1];
      for(idim=1; idim<nsp-1; idim++){

        for(dmin2=1.e38, k1 = 0; k1<csize[kend-1]; k1++){
          code2=off_code+k1+off;
          buf[ntt_isp[idim]-1] =bufsave;

          for(j=ntt_isp[idim]; j<ntt_isp[idim+1]; j++){
            buf[j]=(code[code1][j]+code[code2][j]);
	}
          ntt_dist_lsp(ntt_isp[idim+1]-ntt_isp[idim],rbuf[mode]+ntt_isp[idim],
                       buf+ntt_isp[idim],wegt+ntt_isp[idim],&dist);
          if(dist<dmin2) { dmin2 = dist; index_sp[idim] = k1; }
      }
        for(j=ntt_isp[idim]; j<ntt_isp[idim+1]; j++){
          buf[j]=(code[code1][j]+code[off_code+index_sp[idim]+off][j]);
      }
        bufsave = buf[ntt_isp[idim+1]-1];
    }

      for(dmin2=1.e38, k1 = 0; k1<csize[kend-1]; k1++){
        code2=off_code+k1+off;
        buf[ntt_isp[nsp-1]-1] =bufsave;
        for(j=ntt_isp[nsp-1]; j<ntt_isp[nsp]; j++){
          buf[j]=(code[code1][j]+code[code2][j]);
      }
        ntt_dist_lsp(ntt_isp[nsp]-ntt_isp[nsp-1],rbuf[mode]+ntt_isp[nsp-1],
                     buf+ntt_isp[nsp-1],wegt+ntt_isp[nsp-1],&dist);
        if(dist<dmin2) { dmin2 = dist; index_sp[nsp-1] = k1; }
    }
      for(j=ntt_isp[nsp-1]; j<ntt_isp[nsp]; j++){
        buf[j]=(code[code1][j]+code[off_code+index_sp[nsp-1]+off][j]);
    }
      if(buf[0] < L_LIMIT) {
        diff = L_LIMIT - buf[0];
        buf[0] += diff*1.2;
    }
      if(buf[n_pr-1] > M_LIMIT ){
        diff = buf[n_pr-1]-M_LIMIT;
        buf[n_pr-1] -= diff*1.2;
    }

      ntt_check_lsp(n_pr-1, buf+1, MIN_GAP);

      lspbuf[0] = 0.0;
      for(j=0; j<n_pr; j++){
        lspbuf[j+1] = lsp[j+1] +
          (buf[j] - rbuf[mode][j] ) * fg_sum[mode][j] ;
    }
      if(lspbuf[1] < L_LIMIT) {
        diff = L_LIMIT - lspbuf[1];
        lspbuf[1] += diff*1.2;
    }
      if(lspbuf[n_pr] > M_LIMIT ){
        diff = lspbuf[n_pr]-M_LIMIT;
        lspbuf[n_pr] -= diff*1.2;
    }
      ntt_check_lsp(n_pr-1, lspbuf+2, MIN_GAP);
      ntt_check_lsp(n_pr-1, lspbuf+2, MIN_GAP*0.95);
      ntt_check_lsp_sort(n_pr-1, lspbuf+2);
      ntt_dist_lsp(n_pr, lsp+1, lspbuf+1, wegt, &dist);

      tdist[k+mode*NCD] = dist;
      tindex[k+mode*NCD] =  0;
      for(idim=0; idim<nsp; idim++){
        tindex[k+mode*NCD] +=(index_sp[idim]<<(nbit*(nsp-idim-1)));
    }
  }
}
  for(dmin= 1.e38, mode=0; mode<ntt_N_MODE ; mode ++){
    for(k=0; k<NCD; k++){
      if(tdist[k+mode*NCD]<dmin) {
        dmin = tdist[k+mode*NCD];
        index0 = k;
        mode_index = mode;
    }
  }
}
  *vqword = (mode_index<<(nbit*nsp+nbit0))
    | (index[mode_index][index0]<<(nbit*nsp))
      | tindex[index0+mode_index*NCD];
  index_lsp[0] = mode_index;
  index_lsp[1] = index[mode_index][index0];