melp_sub.c

It is source code for Melp2.4kps vocoder using dsp tms320vc55x of ti

/*

*2.4 kbps MELP Proposed Federal Standard speech coder
*
*TMS320C5x assembly 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 following code was hand optimized for the Texas Instuments
* TMS320C5x DSP by DSPCon, Inc.  For information, please contact DSPCon
* at:
* 
*                       DSPCon, Inc.
*                       380 Foothill Road
*                       Bridgewater, New Jersey 08807
*                       [phone]
*                       [e-mail]
*                       www.dspcon.com
*
*************************************************************************/




/*

  melp_sub.c: MELP-specific subroutines

*/


#include <math.h>
#include "spbstd.h"
#include "mathhalf.h"
#include "math_lib.h"
#include "mat.h"
#include "dsp_sub.h"
#include "melp_sub.h"
#include "pit.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 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

/* Constants */
#define X01_Q14    (0.1*(1<<14))

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 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)
{
    iir_2nd_d(sigin, dc_den, dc_num, sigout,
	       delin,  delout_hi, delout_lo,
	      frame);

    iir_2nd_d(sigout,&dc_den[3],&dc_num[3],sigout,
	      &delin[2],&delout_hi[2],&delout_lo[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 X05_Q6         (0.5 * (1<<6))
#define TEN_Q11        (10 * (1<<11))
#define N2_Q11         (-2*(1<<11))
#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 *tmp_sig;
    Shortword S_temp1, S_temp2;
    Longword L_temp;

    static  Shortword   temp_buf[350];   /* Must be at least 250 */

    /* initialize scaled pitch and minlength */
    pitch_Q6 = shr(pitch,1);      /* Q7->Q6 */    

    tmp_minlength = shl(minlength,6);    

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

      flength = add(flength,pitch_Q6);      
    }

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

    /* increment complexity for if statement */
    
    if (length > maxlength) {
      /* divide by 2 */
      length = shr(length,1);    
    }

    /* 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));


#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);

    
    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 */
    
    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);    
    }
    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);    
    }

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


    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 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);
}

/*