melp_sub.c

音频编解码 2400 bps MELP语音编解码器 定点算法

/*

2.4 kbps MELP Proposed Federal Standard speech coder

Fixed-point C code, version 1.0

Copyright (c) 1998, Texas Instruments, Inc.  

Texas Instruments has intellectual property rights on the MELP
algorithm.  The Texas Instruments contact for licensing issues for
commercial and non-government use is William Gordon, Director,
Government Contracts, Texas Instruments Incorporated, Semiconductor
Group (phone [phone]).

The fixed-point version of the voice codec Mixed Excitation Linear
Prediction (MELP) is based on specifications on the C-language software
simulation contained in GSM 06.06 which is protected by copyright and
is the property of the European Telecommunications Standards Institute
(ETSI). This standard is available from the ETSI publication office
tel. [phone]. ETSI has granted a license to United States
Department of Defense to use the C-language software simulation contained
in GSM 06.06 for the purposes of the development of a fixed-point
version of the voice codec Mixed Excitation Linear Prediction (MELP).
Requests for authorization to make other use of the GSM 06.06 or
otherwise distribute or modify them need to be addressed to the ETSI
Secretariat fax: [phone].

*/
/*

  melp_sub.c: MELP-specific subroutines

*/

#include <stdio.h>
#include <math.h>
#include "spbstd.h"
#include "mathhalf.h"
#include "wmops.h"
#include "math_lib.h"
#include "mat.h"
#include "dsp_sub.h"
#include "melp_sub.h"
#include "pit.h"
#include "constant.h"

/*
    Name: bpvc_ana.c
    Description: Bandpass voicing analysis
    Inputs:
      speech[] - input speech signal
      fpitch[] - initial (floating point) pitch estimates
    Outputs: 
      bpvc[] - bandpass voicing decisions
      pitch[] - frame pitch estimates
    Returns: void 

    Copyright (c) 1995,1997 by Texas Instruments, Inc.  All rights reserved.
*/

/* Q values
   --------
   speech - Q0
   fpitch - Q7
   bpvc - Q14
   pitch - Q7

*/

/* Filter orders */
#define ENV_ORD 2
#define BPF_ORD 6


static Shortword PITCHMAX;
static Shortword PITCHMIN;
static Shortword PITCHMAX_Q7;
static Shortword PITCHMIN_Q7;
static Shortword FRAME;
static Shortword NUM_BANDS;
static Shortword PIT_BEG;
static Shortword PIT_P_FR;
static Shortword PIT_FR_BEG;
static Shortword FIRST_CNTR;
static Shortword BPF_BEG;
static Shortword NUM_PITCHES;
static Shortword LMIN;

/* Static memory */
static Shortword **bpfsp;
static Shortword **bpfdelin;
static Shortword **bpfdelout;
static Shortword **envdel;
static Shortword *envdel2;
static Shortword *sigbuf;

/* External variables */
extern Shortword bpf_num[], bpf_den[];

void bpvc_ana(Shortword speech[], Shortword fpitch[], Shortword bpvc[], 
	      Shortword pitch[])
{

    Shortword pcorr, temp, scale;
    Shortword j,section;
    Shortword filt_index;
	
    /* Filter lowest band and estimate pitch */
    v_equ(&sigbuf[PIT_BEG],&bpfsp[0][0],(Shortword)(PIT_P_FR-FRAME));
    v_equ(&sigbuf[PIT_BEG+PIT_P_FR-FRAME],
	  &speech[PIT_FR_BEG+PIT_P_FR-FRAME],FRAME);

    for (section=0; section<BPF_ORD/2; section++) {
      iir_2nd_s(&sigbuf[PIT_BEG+PIT_P_FR-FRAME],
		&bpf_den[section*3],&bpf_num[section*3],
		&sigbuf[PIT_BEG+PIT_P_FR-FRAME],&bpfdelin[0][section*2],
		&bpfdelout[0][section*2],FRAME);
    }
    v_equ(&bpfsp[0][0],&sigbuf[PIT_BEG+FRAME],(Shortword)(PIT_P_FR-FRAME));

    /* Scale signal for pitch correlations */
    f_pitch_scale(&sigbuf[PIT_BEG],&sigbuf[PIT_BEG],PIT_P_FR);

    *pitch = frac_pch(&sigbuf[FIRST_CNTR],&bpvc[0],fpitch[0],5,
		      PITCHMIN,PITCHMAX,PITCHMIN_Q7,PITCHMAX_Q7,LMIN);

    for (j = 1; j < NUM_PITCHES; j++) {
      temp = frac_pch(&sigbuf[FIRST_CNTR],&pcorr,fpitch[j],5,
		      PITCHMIN,PITCHMAX,PITCHMIN_Q7,PITCHMAX_Q7,LMIN);
	
      /* choose largest correlation value */
      if (pcorr > bpvc[0]) {
	*pitch = temp;
	bpvc[0] = pcorr; 
      }
    }
    
    /* Calculate bandpass voicing for frames */
    
    for (j = 1; j < NUM_BANDS; j++) {

      /* Bandpass filter input speech */
      v_equ(&sigbuf[PIT_BEG],&bpfsp[j][0],(Shortword)(PIT_P_FR-FRAME));
      v_equ(&sigbuf[PIT_BEG+PIT_P_FR-FRAME],
	    &speech[PIT_FR_BEG+PIT_P_FR-FRAME],FRAME);

      for (section=0; section<BPF_ORD/2; section++) {
	  filt_index = j*(BPF_ORD/2)*3 + section*3;
	  iir_2nd_s(&sigbuf[PIT_BEG+PIT_P_FR-FRAME],
		    &bpf_den[filt_index],&bpf_num[filt_index],
		    &sigbuf[PIT_BEG+PIT_P_FR-FRAME],&bpfdelin[j][section*2],
		    &bpfdelout[j][section*2],FRAME);
      }
      v_equ(&bpfsp[j][0],&sigbuf[PIT_BEG+FRAME],(Shortword)(PIT_P_FR-FRAME));

      /* Scale signal for pitch correlations */
      scale = f_pitch_scale(&sigbuf[PIT_BEG],&sigbuf[PIT_BEG],PIT_P_FR);

      /* Check correlations for each frame */
      temp = frac_pch(&sigbuf[FIRST_CNTR],&bpvc[j],*pitch,0,
		      PITCHMIN,PITCHMAX,PITCHMIN_Q7,PITCHMAX_Q7,LMIN);

      /* Calculate envelope of bandpass filtered input speech */
      /* update delay buffers without scaling */
      temp = shr(envdel2[j],scale);
      envdel2[j] = shr(sigbuf[PIT_BEG+FRAME-1],(Shortword)(-scale));
      v_equ_shr(&sigbuf[PIT_BEG-ENV_ORD],&envdel[j][0],scale,ENV_ORD);
      envelope(&sigbuf[PIT_BEG],temp,&sigbuf[PIT_BEG],PIT_P_FR);
      v_equ_shr(&envdel[j][0],&sigbuf[PIT_BEG+FRAME-ENV_ORD],
		(Shortword)(-scale),ENV_ORD);

      /* Scale signal for pitch correlations */
      f_pitch_scale(&sigbuf[PIT_BEG],&sigbuf[PIT_BEG],PIT_P_FR);

      /* Check correlations for each frame */
      temp = frac_pch(&sigbuf[FIRST_CNTR],&pcorr,*pitch,0,
		      PITCHMIN,PITCHMAX,PITCHMIN_Q7,PITCHMAX_Q7,LMIN);

      /* reduce envelope correlation */
      pcorr = sub(pcorr,(Shortword)X01_Q14);

      if (pcorr > bpvc[j])
	bpvc[j] = pcorr;
    }

} /* bpvc_ana */

void bpvc_ana_init(Shortword fr, Shortword pmin, Shortword pmax, 
		   Shortword nbands, Shortword num_p, Shortword lmin)
{

    /* Initialize constants */
    FRAME = fr;
    PITCHMIN = pmin;
    PITCHMAX = pmax;
    PITCHMIN_Q7 = shl(pmin,7);
    PITCHMAX_Q7 = shl(pmax,7);
    NUM_BANDS = nbands;
    NUM_PITCHES = num_p;
    LMIN = lmin;
    PIT_BEG = BPF_ORD;
    PIT_P_FR = ((2*PITCHMAX)+1);
    PIT_FR_BEG = (-PITCHMAX);
    FIRST_CNTR = (PIT_BEG+PITCHMAX);
    BPF_BEG = (PIT_BEG+PIT_P_FR-FRAME);

    /* Allocate memory */
    MEM_2ALLOC(malloc,bpfdelin,NUM_BANDS,BPF_ORD,Shortword);
    v_zap(&bpfdelin[0][0],(Shortword)(NUM_BANDS*BPF_ORD));
    MEM_2ALLOC(malloc,bpfdelout,NUM_BANDS,BPF_ORD,Shortword);
    v_zap(&bpfdelout[0][0],(Shortword)(NUM_BANDS*BPF_ORD));
    MEM_2ALLOC(malloc,bpfsp,NUM_BANDS,PIT_P_FR-FRAME,Shortword);
    v_zap(&bpfsp[0][0],(Shortword)(NUM_BANDS*(PIT_P_FR-FRAME)));

    MEM_2ALLOC(malloc,envdel,NUM_BANDS,ENV_ORD,Shortword);
    v_zap(&envdel[0][0],(Shortword)(NUM_BANDS*ENV_ORD));
    MEM_ALLOC(malloc,envdel2,NUM_BANDS,Shortword);
    v_zap(envdel2,NUM_BANDS);

    /* Allocate scratch buffer */
    MEM_ALLOC(malloc,sigbuf,PIT_BEG+PIT_P_FR,Shortword);
}


/*
    Name: dc_rmv.c
    Description: remove DC from input signal
    Inputs: 
      sigin[] - input signal
      dcdel[] - filter delay history (size DC_ORD)
      frame - number of samples to filter
    Outputs: 
      sigout[] - output signal
      dcdel[] - updated filter delay history
    Returns: void 
    See_Also:

    Copyright (c) 1995,1997 by Texas Instruments, Inc.  All rights reserved.
*/

/* Filter order */
#define DC_ORD 4

/* DC removal filter */
/* 4th order Chebychev Type II 60 Hz removal filter */
/* cutoff=60 Hz, stop=-30 dB */
/* matlab commands: [z,p,k]=cheby2(4,30,0.015,'high');  */
/*                      sos=zp2sos(z,p,k);        */
/* order of sections swapped */
/* Q13 */
static Shortword dc_num[(DC_ORD/2)*3] = {
           7769, -15534, 7769,
           8007, -15999, 8007};

/* sign of coefficients for dc_den is reversed */
static Shortword dc_den[(DC_ORD/2)*3] = {
           -8192, 15521, -7358,
	   -8192, 15986, -7836};

void dc_rmv(Shortword sigin[],Shortword sigout[],Shortword delin[],
	    Shortword delout_hi[],Shortword delout_lo[],
	    Shortword frame)
{
  Shortword section;

  /* Remove DC from input speech */
  v_equ(sigout,sigin,frame);
  for (section=0; section<DC_ORD/2; section++) {
    iir_2nd_d(sigout,&dc_den[section*3],&dc_num[section*3],sigout,
	      &delin[section*2],&delout_hi[section*2],&delout_lo[section*2],
	      frame);
  }
}


/*
    Name: gain_ana.c
    Description: analyze gain level for input signal
    Inputs: 
      sigin[] - input signal
      pitch - pitch value (for pitch synchronous window)
      minlength - minimum window length
      maxlength - maximum window length
    Outputs: 
    Returns: log gain in dB
    See_Also:

    Copyright (c) 1997 by Texas Instruments, Inc.  All rights reserved.

    Q values
    --------
    pitch - Q7
    sigin (speech) - Q0

*/

#define CONSTANT_SCALE 0
#define MINGAIN        0
#define LOG10OF2X2     1233         /* 2*log10(2) in Q11 */

Shortword gain_ana(Shortword sigin[], Shortword pitch, 
		   Shortword minlength, Shortword maxlength)
{
    Shortword length;    /* Q0 */
    Shortword flength;   /* Q6 */
    Shortword gain;      /* Q8 */
    Shortword pitch_Q6;  /* Q6 */
    Shortword scale;
    Shortword tmp_minlength;
    Shortword sigbegin;
    Shortword *temp_buf, *tmp_sig;
    Shortword S_temp1, S_temp2;
    Longword L_temp;
    
    /* initialize scaled pitch and minlength */
    pitch_Q6 = shr(pitch,1);      /* Q7->Q6 */    data_move();

    tmp_minlength = shl(minlength,6);    data_move();

    /* Find shortest pitch multiple window length (floating point) */
    flength = pitch_Q6;    data_move();
    while (flength < tmp_minlength) {
      /* increment complexity for while statement */
      compare_nonzero();

      flength = add(flength,pitch_Q6);      data_move();
    }

    /* Convert window length to integer and check against maximum */
    length = shr(add(flength, (Shortword)X05_Q6), 6);     data_move();

    /* increment complexity for if statement */
    compare_nonzero();
    if (length > maxlength) {
      /* divide by 2 */
      length = shr(length,1);    data_move();
    }
    
    /* Calculate RMS gain in dB */
    /*gain = 10.0*log10(0.01+(v_magsq(&sigin[-(length/2)],length)/length));*/
    
    /* use this adaptive scaling if more precision needed at low 
       levels.  Seemed like it wasn't worth the mips */

    sigbegin = negate(shr(length,1));

    /* Right shift input signal and put in temp buffer */
    MEM_ALLOC(MALLOC,temp_buf,length,Shortword);

#if(CONSTANT_SCALE)
    scale=3;
#else
    scale = 3;
    v_equ_shr(temp_buf,&sigin[sigbegin],scale,length);
    L_temp = L_v_magsq(temp_buf,length,0,1);

    compare_zero();
    if (L_temp) {
      S_temp1 = norm_l(L_temp);
      scale = sub(scale,shr(S_temp1,1));
      if (scale < 0)
	scale = 0;
    }
    else
      scale = 0;

#endif

    if (scale)
      v_equ_shr(temp_buf,&sigin[sigbegin],scale,length);
    else
      v_equ(temp_buf,&sigin[sigbegin],length);
    tmp_sig = temp_buf - sigbegin;


    /* increment complexity for if statement */
    compare_zero();
    if (scale) {
      L_temp = L_v_magsq(&tmp_sig[sigbegin],length,0,0);
      S_temp1 = L_log10_fxp(L_temp,0);
      S_temp2 = L_log10_fxp(L_deposit_l(length),0);
      S_temp1 = sub(S_temp1,S_temp2);
      /* shift right to counter-act for the shift in L_mult */
      S_temp2 = extract_l(L_shr(L_mult(scale,LOG10OF2X2),1));
      S_temp1 = add(S_temp1,S_temp2);
      S_temp1 = mult(TEN_Q11,S_temp1);
      gain = shl(S_temp1,1);    data_move();
    }
    else {
      L_temp = L_v_magsq(&tmp_sig[sigbegin],length,0,0);
      /* Add one to avoid log of a zero value */
      if (L_temp == 0)
	S_temp1 = N2_Q11;
      else
	S_temp1 = L_log10_fxp(L_temp,0);
      S_temp2 = L_log10_fxp(L_deposit_l(length),0);
      S_temp1 = sub(S_temp1,S_temp2);
      S_temp1 = mult(TEN_Q11,S_temp1);
      gain = shl(S_temp1,1);    data_move();
    }

    /* increment complexity for if statement */
    compare_nonzero();
    if (gain < MINGAIN) {
      gain = MINGAIN;    data_move();
    }

    MEM_FREE(FREE,temp_buf);

    return(gain);
} /* gain_ana */

/*
    Name: lin_int_bnd.c
    Description: Linear interpolation within bounds
    Inputs:
      x - input X value
      xmin - minimum X value
      xmax - maximum X value
      ymin - minimum Y value
      ymax - maximum Y value
    Returns: y - interpolated and bounded y value

    Copyright (c) 1995,1997 by Texas Instruments, Inc.  All rights reserved.

    Q values:
    x,xmin,xmax - Q8
    y,ymin,ymax - Q15
*/

Shortword lin_int_bnd(Shortword x,Shortword xmin,Shortword xmax,
		      Shortword ymin,Shortword ymax)

{
    Shortword y;
    Shortword temp1,temp2;
			
    if (x <= xmin)
      y = ymin;

    else if (x >= xmax)
      y = ymax;

    else {
      /* y = ymin + (x-xmin)*(ymax-ymin)/(xmax-xmin); */
      temp1 = sub(x,xmin);
      temp2 = sub(ymax,ymin);
      temp1 = mult(temp1,temp2);     /* temp1 in Q8 */
      temp2 = sub(xmax,xmin);
      /* temp1 always smaller than temp2 */
      temp1 = divide_s(temp1,temp2);    /* temp1 in Q15 */
      y = add(ymin,temp1);
    }

    return(y);
}

/*
    Name: noise_est.c
    Description: Estimate long-term noise floor
    Inputs:
      gain - input gain (in dB)
      noise_gain - current noise gain estimate
      up - maximum up stepsize
      down - maximum down stepsize
      min - minimum allowed gain
      max - maximum allowed gain
    Outputs:
      noise_gain - updated noise gain estimate
    Returns: void

    Copyright (c) 1995,1997 by Texas Instruments, Inc.  All rights reserved.

    Q values: