ntt_tf_proc_spectrum.c

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

void ntt_enc_pgain(/* Input */
		   double power,      /* Power to be encoded */
		   /* Output */
		   int    *index, /* Power code index */
		   double *pgain,
		   double *norm)      /* Gain */
{
    /*--- Variables --*/
    double mu_power;
    int itmp;
    double pgain_step;
    
    pgain_step   = ntt_PGAIN_MAX / ntt_PGAIN_NSTEP;

    mu_power = ntt_mulaw(power, ntt_PGAIN_MAX, ntt_PGAIN_MU);           /* mu-Law power */
    itmp = (int)floor((mu_power)/ pgain_step);  /* quantization */
    *index = itmp;  /* quantization */
    
    ntt_dec_pgain( *index, pgain);
    *norm = 1./(*pgain+1.e-3);

}


void ntt_enc_pitch(/* Input */
	       double spectrum[],         /* spectrum */
	       double lpc_spectrum[],     /* LPC spectrum */
	       int    w_type,             /* block type */
	       int    numChannel,
	       int     block_size_samples,
	       int     isampf,
	       double bandUpper,
	       double *pcode,
	       short *pleave0,
	       short *pleave1,
	     /* Output */
	       int    index_pit[],        /* indices for periodic peak components coding */
	       int    index_pls[],        
	       int    index_pgain[],
	       double pitch_sequence[])   /* periodic peak components */


{
    /*--- Variables ---*/
    double tc[ntt_T_FR_MAX], pwt[ntt_T_FR_MAX];
    double pit_seq[ntt_T_FR_MAX];
    double pit_pack[ntt_N_FR_P_MAX*MAX_TIME_CHANNELS];
    double tc_pack[ntt_N_FR_P_MAX*MAX_TIME_CHANNELS];
    double wtbuf[ntt_T_FR_MAX];
    double pwtbuf[ntt_T_FR_MAX];
    double norm[MAX_TIME_CHANNELS], dtmp, dpow; 
    int    ismp, i_sup, top, ptop;
    double pgain[MAX_TIME_CHANNELS];
    double pgainq[MAX_TIME_CHANNELS];

    switch(w_type){
    case ONLY_LONG_SEQUENCE:
    case LONG_START_SEQUENCE:
    case LONG_STOP_SEQUENCE:
	/* LPC flattening */
	for (i_sup=0; i_sup<numChannel; i_sup++){
	    top = i_sup * block_size_samples;
	    for (ismp=0; ismp<block_size_samples; ismp++){
		tc[ismp+top] = spectrum[ismp+top] * lpc_spectrum[ismp+top];
	    }
	}
	/* pitch components extraction */
	for (i_sup=0; i_sup<numChannel; i_sup++){
	    top = i_sup * block_size_samples;
	    ptop = i_sup * ntt_N_FR_P;
            ntt_sear_pitch(tc+top, lpc_spectrum+top,
		    block_size_samples, isampf, bandUpper, &index_pit[i_sup]);
            ntt_extract_pitch(index_pit[i_sup], tc+top, block_size_samples,
		    isampf, bandUpper, tc_pack+ptop);
            ntt_extract_pitch(index_pit[i_sup], lpc_spectrum+top, 
		    block_size_samples, isampf, bandUpper, wtbuf+ptop);
	    dtmp = 1.e-3+ntt_dotdd(ntt_N_FR_P, tc_pack+ptop, tc_pack+ptop);
	    dpow = sqrt(dtmp/ntt_N_FR_P);
	    ntt_enc_pgain( dpow, &index_pgain[i_sup], &pgain[i_sup], &norm[i_sup]);
	    ntt_dec_pgain(index_pgain[i_sup], &pgainq[i_sup]);
	    for (ismp=0; ismp<block_size_samples; ismp++){
		pwt[ismp+top] = 1. / pow(lpc_spectrum[ismp+top], 0.5);
            }
	    ntt_extract_pitch(index_pit[i_sup], pwt+top, block_size_samples,
	         isampf, bandUpper, pwtbuf+ptop);
            ntt_extract_pitch(index_pit[i_sup],spectrum+top,block_size_samples, 
		 isampf, bandUpper, pit_pack+ptop);
	    dtmp = 1./(pgainq[i_sup]+1.e-5);
	    for (ismp=0; ismp<ntt_N_FR_P; ismp++){
	      pit_pack[ptop+ismp] *= dtmp;
	    }
            pgain[i_sup] = pow(pgain[i_sup] , ntt_GAMMA_CH);	
	}

	/*--- Quantization of the pitch value ---*/
	ntt_wvq_pitch(pit_pack, wtbuf, pwtbuf, pgain,
		  pcode, pleave0, pleave1, numChannel,
		  index_pls, pgainq, index_pgain);
	/*--- Making the pitch sequence ---*/
	ntt_dec_pit_seq(index_pit, index_pls, numChannel,
	       block_size_samples, isampf, bandUpper,
	       pcode, pleave1, pit_seq);
	for (i_sup=0; i_sup<numChannel; i_sup++){
	    top = i_sup * block_size_samples;
	    for(ismp=0; ismp<block_size_samples; ismp++){
		pitch_sequence[ismp+top] = pgainq[i_sup] * pit_seq[ismp+top];
	    }
	}
	break;
    case EIGHT_SHORT_SEQUENCE:
	ntt_zerod(block_size_samples*numChannel, pitch_sequence);
	break;
    }


}

void     ntt_extract_pitch(
	      /* Input */
	      int     index_pit, 
	      double  pit_seq[],
	      int     block_size_samples,
	      int     isampf,
	      double  bandUpper,
	      /* Output */
	      double  pit_pack[]) 
{
    /*--- Variables ---*/
    int	   i_smp;
    static int InitFlag = 1;
    float dtmp, pitch;
    int npcount, iscount;
    int ii,jj;
    float fcmin, fcmax, fctmp0, fctmp1;
    float tmpnp0, tmpnp1;
    float bandwidth;
    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;
    InitFlag = 0;
    bandwidth=4.0;
    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*/	
        dtmp = (float)index_pit / (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;
        iscount=0;

        for (jj=0; jj<ntt_N_FR_P; jj++){
              pit_pack[jj] = 0.0;
        }

        for (jj=0; jj<npcount/2; jj++){
              pit_pack[iscount++] = pit_seq[jj];
        }

        for (ii=0; (ii<ntt_N_FR_P) && (iscount<ntt_N_FR_P); ii++){
           dtmp = pitch*(float)(ii+1);
           dtmp += 0.5;
           i_smp = (int)dtmp;
           if(i_smp+(npcount-1)/2+1 < block_size_samples) {
             for (jj=-npcount/2; (jj<(npcount-1)/2+1) && (iscount < ntt_N_FR_P);
 jj++, iscount++){
                pit_pack[iscount] = pit_seq[jj+i_smp];
             }
           }
        }

}

void ntt_freq_domain_pre_process(/* Parameters */
				 int    nfr,              /* block length */
				 int    nsf,              /* number of sub frames */
				 double band_lower,
				 double band_upper,
				 /* Input */
				 double spectrum[],       /* spectrum */
				 double lpc_spectrum[],   /* LPC spectrum */
				 /* Output */
				 double spectrum_out[],   /* spectrum */
				 double lpc_spectrum_out[]) /* LPC spectrum */
{
  /*--- Variables ---*/
  int isf, subtop, start, stop, ismp, width;
  float ftmp;

  /*--- Band limit ---*/
  ftmp = (float)band_lower;
  ftmp *= (float)nfr;
  start =  (int)ftmp;
  ftmp = (float)band_upper-(float)band_lower;
  ftmp *= (float)nfr;
  width  = (int)ftmp;
  stop  = start+width;
  for (isf=0; isf<nsf; isf++){
    subtop = isf * nfr;
    /* lower cut */
    for (ismp=0; ismp<start; ismp++){
      spectrum_out[ismp+subtop] = 0.;
      lpc_spectrum_out[ismp+subtop] = 1.e5;
    }
    /* signal band */
    for (ismp=start; ismp<stop; ismp++){
      spectrum_out[ismp+subtop] = spectrum[ismp+subtop];
      lpc_spectrum_out[ismp+subtop] = lpc_spectrum[ismp+subtop];
    }
    /* upper cut */
    for (ismp=stop; ismp<nfr; ismp++){
      spectrum_out[ismp+subtop] = 0.;
      lpc_spectrum_out[ismp+subtop] = 1.e5;
    }
  }
}


void ntt_fw_div(int    ndiv,    /* Param. : number of interleave division */
		int    length,  /* Param. : codebook length */
		double prcpt,   /* Param. : perceptual weight */
		double env[],   /* Input  : envelope vector */
		double denv[],  /* Output : envelope subvector */
		double wt[],    /* Input  : weighting factor vector */
		double dwt[],   /* Output : weighting factor subvector */
		int    idiv)    /* Param. : division number */
{
  int ismp, ismp2;

  for ( ismp=0; ismp<length; ismp++ ){
    ismp2 = idiv + ismp * ndiv;
    denv[ismp] = env[ismp2];
    dwt[ismp] = pow(wt[ismp2], prcpt);
  }
}


void ntt_fw_sear(int    cb_size, /* Param  : codebook size */
		 int    length,  /* Param  : codevector length */
		 double *codev,  /* Param  : codebook */
		 int    len_max, /* Param  : codevector memory length */
		 double denv[],  /* Input  : envelope subvector */
		 double dwt[],   /* Input  : weight subvector */
		 int    *index,  /* Output : quantization index */
		 double *dmin)   /* Output : quant. weighted distortion */
{
  /*--- Variables ---*/
  double dtemp, ldmin;
  int    ismp, icb;
  double dist_buf;
  
  ldmin = 1.e50;
  for (icb=0; icb<cb_size; icb++){
    /*--- Calculate the weighted distortion measure ---*/
    dist_buf = 0.;
    for (ismp=0; ismp<length; ismp++){
      dtemp = denv[ismp] - codev[icb*len_max+ismp];
      dist_buf += dwt[ismp] * dtemp * dtemp;
    }
    
    /*--- Search for the optimum code vector ---*/
    if ( dist_buf < (ldmin) ){
      ldmin = dist_buf;
      *index = icb;
    }
  }
  *dmin = ldmin;
}


void ntt_fwalfenc(double alf,    /* Input  : AR prediction coefficients */
	      double *alfq,  /* Output : quantized AR pred. coefficients */
	      int *ind)      /* Output : the index */
{
    /*--- Variables ---*/
    double expand;
    
    /*--- Main operation ---*/
    expand = alf / (double) ntt_FW_ALF_STEP+0.5;
    *ind = ntt_min( (int)expand, ntt_FW_ARQ_NSTEP-1 );
    *ind = ntt_max(*ind, 0);
    *alfq = (double)*ind*ntt_FW_ALF_STEP;
}

void ntt_fwpred(int    nfr,  /* Param:  block size */
		int    n_crb,       /* Param:  number of Bark-scale band */
		int    *bark_tbl,   /* Param:  Bark-scale table */
		double  bandUpper,
		int    ndiv,        /* Param:  number of interleave division */
		int    cb_size,     /* Param:  codebook size */
		int    length,      /* Param:  codevector length */
		double *codev,      /* Param:  codebook */
		int    len_max,     /* Param:  codevector memory length */
		double prcpt,
		int     prev_fw_code[],
		/*double *penv_tmp, *//*Param:  memory for Bark-scale envelope*/
		/*double *pdot_tmp, *//*Param:  memory for env. calculation */
		double iwt[],       /* Input:  LPC envelope */
		double pred[],      /* Output : Prediction factor */
		double tc[],        /* Input:  Residual signal */
		int    index_fw[],  /* Output: Envelope VQ indices */
		int    *ind_fw_alf, /* Output: AR prediction coefficient index */
		int    i_sup,       /* Input:  Channel number */
		int    pred_sw)     /* Input:  prediction switch */
{
    /*--- Variables ---*/
    double        bwt[ntt_N_CRB_MAX], env[ntt_N_CRB_MAX];
    double        fwgain, alf, alfq;
    double        penv[ntt_N_CRB_MAX];
    double        envsave[ntt_N_CRB_MAX];
    register double acc, cdot, acc2, col0, col1, dtmp;
    int           top, top2, ismp, iblow, ibhigh, ib;
    int  n_crb_eff;
    float ftmp;

    /*--- Initialization ---*/
    top  = i_sup*nfr;
    top2 = i_sup*n_crb;
    ntt_fwex( prev_fw_code, ndiv, length, codev, len_max, penv );
    /*--- Calculate envelope ---*/
    iblow=0; fwgain=0.; ib=0;
     n_crb_eff=0;
     for (ib=0; ib<n_crb; ib++){
	ibhigh = bark_tbl[ib];
	/* reset the accumulations */
	acc  = 0.0;
	acc2 = 0.0;
	/* bandle the input coefficient lines into the bark-scale band */
	ftmp = (float)bandUpper;
	ftmp *= (float)nfr;
	for(ismp=iblow; (ismp<ibhigh) && (ismp<(int)ftmp); ismp++){
	    cdot =  tc[ismp]*tc[ismp];/* calc. power */
	    acc +=  cdot; /* acc. the current and previous power */
	    acc2 = acc2 + iwt[ismp];  /* accumulate the weighting factor */
	}
	env[ib] = sqrt(acc/((double)(ismp-iblow)+0.001)) + 0.1; /* envelope in bark*/
	bwt[ib] = acc2/((double)(ismp-iblow)+0.001); 	/* weighting factor in bark scale */
	if(env[ib]>0.1){
	  n_crb_eff++;
	  fwgain += env[ib];	/* accumulate the envelope power */
	}
	iblow = ibhigh;
    }
    /* gain to normalize the bark-scale envelope */
    if(n_crb_eff != 0){
	   fwgain = (double) n_crb_eff / fwgain;
    }
    else{
	  fwgain = 0.;
    }

    /*--- Normalize the envelope, and remove the offset ---*/
    for (ib=0; ib<n_crb_eff; ib++){
	env[ib] = env[ib]*fwgain - 1.0;
	envsave[ib]= env[ib];
    }
    for (ib=n_crb_eff; ib<n_crb; ib++){
	env[ib] = 0.0;
	envsave[ib]= env[ib];
    }

    /*--- Calculate the AR prediction coefficient and quantize it ---*/
    col0 = 0.1;
    col1 = 0.;
    if (pred_sw == 1){
      /* calculate auto correlation and cross correlation */
      for (ib=0; ib<n_crb; ib++){
	col0 += penv[ib]*penv[ib]*bwt[ib];
	col1 += penv[ib]*env[ib]*bwt[ib];
      }
      alf = col1/col0;   /* normalize the cross correlation coefficient */
      ntt_fwalfenc( alf, &alfq, ind_fw_alf ); /* quantize the cor. coefficient */
      /* AR prediction (cancellation using cor. coeff.) */
      for ( ib=0; ib<n_crb; ib++ ){
	env[ib] -= alfq*penv[ib];
      }
    }
    else{ /* in stop frame */
      alf = 0.;
      ntt_fwalfenc( alf, &alfq, ind_fw_alf );
    }
    /*--- Quantization of the envelope ---*/
    ntt_fwvq(ndiv, cb_size, length, codev, len_max, prcpt,
	     env, bwt, index_fw);

    /*--- Local decoding ---*/
    ntt_fwex( index_fw, ndiv, length, codev, len_max, env );
    for (ib=ndiv*length; ib<n_crb; ib++) env[ib] = 0.;
    for (ib=0, ismp=0; ib<n_crb; ib++){
	ibhigh = bark_tbl[ib];
	dtmp = ntt_max(env[ib]+alfq*penv[ib]+1., ntt_FW_LLIM);
	while(ismp<ibhigh){
	    pred[ismp] = dtmp;
	    ismp++;
	}
    }
    /*--- Store the current envelope ---
    ntt_movdd(n_crb, env, penv); */
}


void ntt_fwvq(int    ndiv,     /* Param  : number of interleave division */
	      int    cb_size,  /* Param  : codebook size */
	      int    length,   /* Param  : codevector length */
	      double *codev,   /* Param  : codebook */
	      int    len_max,  /* Param  : codevector memory length */
	      double prcpt,    /* Param  : perceptual control factor */
	      double env[],    /* Input  : the envelope */
	      double wt[],     /* Input  : the weighting factor */
	      int    index[])  /* Output : indexes */
{
    /*--- Variables ---*/
    int    idiv;
    double denv[ntt_FW_CB_LEN_MX], dwt[ntt_FW_CB_LEN_MX], dmin;

    /*--- Initialization ---*/

    /*--- Main operation ---*/
    for ( idiv=0; idiv<ndiv; idiv++ ){
	/*--- Divide the input frame into subframe ---*/
	ntt_fw_div(ndiv, length, prcpt,
		   env, denv, wt, dwt, idiv);
	/*--- search ---*/
	ntt_fw_sear(cb_size, length,
		    codev, len_max,
		    denv, dwt, &index[idiv], &dmin );

    }
}

void ntt_get_wegt(/* Parameters */
		  int n_pr,
                  /* Input */
		  double flsp[],
		  /* Output */
		  double wegt[])
{
    int    i;
    double tmp;
    double llimit, ulimit;
    
    llimit = 0.4/(double)n_pr;
    ulimit = 1.-0.4/(double)n_pr;

    tmp     = ( ( flsp[2]-PI*llimit )-1.0);
    wegt[1] = tmp * tmp * 10.+1.0;
    if (tmp > 0.0 ) wegt[1] = 1.0;

    for ( i=2; i<=n_pr-1; i++ ) {
	tmp     = -flsp[i-1]+flsp[i+1]- 1.0;
        wegt[i] = tmp * tmp * 10.+1.0;
        if(tmp > 0.0 ) wegt[i] = 1.0;
    }
    tmp = -flsp[n_pr-1]+PI*ulimit- 1.0;
    wegt[n_pr] = tmp* tmp * 10.+1.0;
    if (tmp > 0.0 ) wegt[n_pr] = 1.0;
}



void ntt_lsp_encw(int    n_pr,
                  double code[][ntt_N_PR_MAX],
                  double fgcode[][ntt_MA_NP][ntt_N_PR_MAX],
                  int    *csize,
                  int    prev_lsp_code[],
		  /*
		  double buf_prev[ntt_MA_NP][ntt_N_PR_MAX+1],
                  */
                  int    ma_np,
		  double freq[],     /* Input  : LSP parameters */
                  double freqout[],  /* Output : quantized LSP parameters */
                  int index[],       /* Output : the optimum LSP code index */
                  int  nsp)
{