melp_ana.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].

*/
/*
    Name: melp_ana.c
    Description: MELP analysis
    Inputs:
      speech[] - input speech signal
    Outputs: 
      *par - MELP parameter structure
    Returns: void
*/

/* compiler include files */

#include <stdio.h>
#include <math.h>
#include "mathhalf.h"
#include "melp.h"
#include "spbstd.h"
#include "lpc.h"
#include "mat.h"
#include "vq.h"
#include "fs.h"
#include "pit.h"
#include "math_lib.h"
#include "constant.h"
#include "wmops.h"

/* compiler constants */

#define BEGIN 0
#define END 1

#define HF_CORR_Q15 4       /* 0.0001220703125 in Q15 */
#define LAG_WIN 0.005  /* not used */

#define PDECAY_Q15 31129            /* (0.95*(1<<15)) */
#define PDECAY_PITCH_Q7 320         /* ((0.05*DEFAULT_PITCH)*(1<<7)) */
#define PEAK_THRESH_Q12 5488        /* (1.34*(1<<12)) */
#define PEAK_THR2_Q12 6553          /* (1.6*(1<<12)) */
#define SILENCE_DB_Q8 7680          /* (30.0*(1<<8)) */
#define MAX_ORD LPF_ORD
#define FRAME_BEG (PITCHMAX-(FRAME/2))
#define FRAME_END (FRAME_BEG+FRAME)
#define PITCH_BEG (FRAME_END-PITCHMAX)
#define PITCH_FR ((2*PITCHMAX)+1)
#define IN_BEG (PITCH_BEG+PITCH_FR-FRAME)
#define SIG_LENGTH (LPF_ORD+PITCH_FR)


/* external memory references */
 
extern Shortword win_cof[];

extern Shortword lpf_num[];
extern Shortword lpf_den[];

extern Shortword msvq_cb[];
extern Shortword msvq_cb_mean[];
extern Shortword fsvq_cb[];
extern Shortword fsvq_weighted;

extern FILE *fp_indat,*fp_outdat;

/* memory definitions */
 
static Shortword sigbuf[SIG_LENGTH];
static Shortword speech[IN_BEG+FRAME];
static Shortword dcdelin[DC_ORD];
static Shortword dcdelout_hi[DC_ORD];
static Shortword dcdelout_lo[DC_ORD];
static Shortword lpfsp_delin[LPF_ORD];
static Shortword lpfsp_delout[LPF_ORD];
static Shortword pitch_avg;
static Shortword fpitch[2];
static struct msvq_param vq_par;  /* MSVQ parameters */
static struct msvq_param fs_vq_par;  /* Fourier series VQ parameters */
static Shortword w_fs[NUM_HARM];

Shortword frames=0;


void melp_ana(Shortword sp_in[],struct melp_param *par)
{

    Shortword i;
    Shortword begin;
    Shortword sub_pitch;
    Shortword temp,pcorr,bpthresh;
    Shortword r[LPC_ORD+1],refc[LPC_ORD+1],lpc[LPC_ORD+1];
    Shortword weights[LPC_ORD];
    Shortword section;
    Shortword temp_delin[LPF_ORD];
    Shortword temp_delout[LPF_ORD];
    /* char in_num[20];*/

    /* Remove DC from input speech */
    dc_rmv(sp_in,&speech[IN_BEG],dcdelin,dcdelout_hi,dcdelout_lo,FRAME);

#if (COMPLEXITY_COUNT)
complexity_count();
#endif

    /* Copy input speech to pitch window and lowpass filter */
    v_equ(&sigbuf[LPF_ORD],&speech[PITCH_BEG],PITCH_FR);
    for (section=0; section<LPF_ORD/2; section++) {
      iir_2nd_s(&sigbuf[LPF_ORD],&lpf_den[section*3],&lpf_num[section*3],
                &sigbuf[LPF_ORD],&lpfsp_delin[section*2],
                &lpfsp_delout[section*2],FRAME);
      /* save delay buffers for the next overlapping frame */
      for (i = section*2; i < section*2+2; i++) {
        temp_delin[i] = lpfsp_delin[i];      data_move();
        temp_delout[i] = lpfsp_delout[i];      data_move();
      }
      iir_2nd_s(&sigbuf[LPF_ORD+FRAME],&lpf_den[section*3],
                &lpf_num[section*3],&sigbuf[LPF_ORD+FRAME],
                &lpfsp_delin[section*2],&lpfsp_delout[section*2],
                PITCH_FR-FRAME);
      /* restore delay buffers for the next overlapping frame */
      for (i = section*2; i < section*2+2; i++) {
        lpfsp_delin[i] = temp_delin[i];      data_move();
        lpfsp_delout[i] = temp_delout[i];      data_move();
      }
    }
#if (COMPLEXITY_COUNT)
complexity_count();
#endif

    /* Perform global pitch search at frame end on lowpass speech signal */
    /* Note: avoid short pitches due to formant tracking */
    f_pitch_scale(&sigbuf[LPF_ORD],&sigbuf[LPF_ORD],PITCH_FR);
    fpitch[END] = find_pitch(&sigbuf[LPF_ORD+(PITCH_FR/2)],&temp,
			     (2*PITCHMIN),PITCHMAX,PITCHMAX);
    fpitch[END] = shl(fpitch[END],7);     /* fpitch in Q7 */

#if (COMPLEXITY_COUNT)
complexity_count();
#endif

    /* Perform bandpass voicing analysis for end of frame */
    bpvc_ana(&speech[FRAME_END], fpitch, &par->bpvc[0], &sub_pitch);

#if (COMPLEXITY_COUNT)
complexity_count();
#endif

    /* Force jitter if lowest band voicing strength is weak */    
    if (par->bpvc[0] < VJIT_Q14)      /* par->bpvc in Q14 */
	par->jitter = (Shortword)MAX_JITTER_Q15;   /* par->jitter in Q15 */
    else
	par->jitter = (Shortword)0;

    /* Calculate LPC for end of frame */
    lpc_autocorr(&speech[(FRAME_END-(LPC_FRAME/2))],win_cof,r,
                 HF_CORR_Q15,(Shortword)LAG_WIN,LPC_ORD,LPC_FRAME);

    lpc[0] = ONE_Q12;        /* 1 in Q12 */
    lpc_schur(r,lpc,refc,LPC_ORD);     /* lpc in Q12, refc in Q15 */

#if (COMPLEXITY_COUNT)
complexity_count();
#endif

    /* Calculate LPC residual */
    zerflt(&speech[PITCH_BEG],lpc,&sigbuf[LPF_ORD],LPC_ORD,PITCH_FR);

#if (COMPLEXITY_COUNT)
complexity_count();
#endif

    /* Check peakiness of residual signal */
    begin = (LPF_ORD+(PITCHMAX/2));
    temp = peakiness(&sigbuf[begin],PITCHMAX);    /* temp in Q12 */

#if (COMPLEXITY_COUNT)
complexity_count();
#endif

    /* Peakiness: force lowest band to be voiced  */
    if (temp > PEAK_THRESH_Q12) {
	par->bpvc[0] = 1<<14;
    }

    /* Extreme peakiness: force second and third bands to be voiced */
    if (temp > PEAK_THR2_Q12) {
	par->bpvc[1] = 1<<14;
	par->bpvc[2] = 1<<14;
    }

    /* Calculate overall frame pitch using lowpass filtered residual */
    par->pitch = pitch_ana(&speech[FRAME_END], &sigbuf[LPF_ORD+PITCHMAX], 
			   sub_pitch,pitch_avg,&pcorr);
#if (COMPLEXITY_COUNT)
complexity_count();
#endif

    bpthresh = (Shortword)BPTHRESH_Q14;

    /* Calculate gain of input speech for each gain subframe */
    for (i = 0; i < NUM_GAINFR; i++) {
	if (par->bpvc[0] > bpthresh) {

	    /* voiced mode: pitch synchronous window length */
	    temp = sub_pitch;

	    par->gain[i] = gain_ana(&speech[FRAME_BEG+(i+1)*GAINFR],
				    temp,MIN_GAINFR,PITCHMAX_X2);

	}
	else {
	    temp = (Shortword)GAIN_PITCH_Q7;

	    par->gain[i] = gain_ana(&speech[FRAME_BEG+(i+1)*GAINFR],
				    temp,0,PITCHMAX_X2);

	}
    }
#if (COMPLEXITY_COUNT)
complexity_count();
#endif

    /* Update average pitch value */
    if (par->gain[NUM_GAINFR-1] > SILENCE_DB_Q8)    /* par->gain in Q8 */
      temp = pcorr;
    else
      temp = 0;
    /* pcorr in Q14 */
    pitch_avg = p_avg_update(par->pitch,temp,(Shortword)VMIN_Q14);

    /* Calculate Line Spectral Frequencies */
    lpc_pred2lsp(lpc,par->lsf,LPC_ORD);

#if (COMPLEXITY_COUNT)
complexity_count();
#endif

    /* Force minimum LSF bandwidth (separation) */
    lpc_clamp(par->lsf,BWMIN_Q15,LPC_ORD);

    /* Quantize MELP parameters to 2400 bps and generate bitstream */
    /* Quantize LSF's with MSVQ */
    vq_lspw(weights, &par->lsf[1], lpc, LPC_ORD);

    msvq_enc(&par->lsf[1], weights, &par->lsf[1], vq_par);

#if (COMPLEXITY_COUNT)
complexity_count();
#endif

frames++;

    par->msvq_index = vq_par.indices;

    /* Force minimum LSF bandwidth (separation) */
    lpc_clamp(par->lsf,BWMIN_Q15,LPC_ORD);

    /* Quantize logarithmic pitch period */
    /* Reserve all zero code for completely unvoiced */
    par->pitch = log10_fxp(par->pitch,7);      /* par->pitch in Q12 */
    quant_u(&par->pitch,&par->pitch_index,(Shortword)PIT_QLO_Q12,
	    (Shortword)PIT_QUP_Q12,PIT_QLEV_M1,PIT_QLEV_M1_Q8,1,7);
    /* convert pitch back to linear in Q7 */
    par->pitch = pow10_fxp(par->pitch,7);

    /* Quantize gain terms with uniform log quantizer	*/
    q_gain(par->gain, par->gain_index,(Shortword)GN_QLO_Q8,
	   (Shortword)GN_QUP_Q8,GN_QLEV_M1,GN_QLEV_M1_Q10,0,5);
    
    /* Quantize jitter */
    /* quant_u(&par->jitter,&par->jit_index,0,MAX_JITTER_Q15,2); */
    if (par->jitter < shr((Shortword)MAX_JITTER_Q15,1)) {
      par->jitter = 0;
      par->jit_index = 0;
    }
    else {
      par->jitter = (Shortword)MAX_JITTER_Q15;
      par->jit_index = 1;
    }

    /* Quantize bandpass voicing */
    par->uv_flag = q_bpvc(&par->bpvc[0],&par->bpvc_index,bpthresh,
			  NUM_BANDS);
#if (COMPLEXITY_COUNT)
complexity_count();
#endif

    /*	Calculate Fourier coefficients of residual signal from quantized LPC */
    fill(par->fs_mag,ONE_Q13,NUM_HARM);
    if (par->bpvc[0] > bpthresh) {
	lpc_lsp2pred(par->lsf,lpc,LPC_ORD);

	zerflt(&speech[(FRAME_END-(LPC_FRAME/2))],lpc,sigbuf,
	       LPC_ORD,LPC_FRAME);
	window(sigbuf,win_cof,sigbuf,LPC_FRAME);

	find_harm(sigbuf, par->fs_mag, par->pitch, NUM_HARM, LPC_FRAME);

    }
#if (COMPLEXITY_COUNT)
complexity_count();
#endif
    
    /* quantize Fourier coefficients */
    /* pre-weight vector, then use Euclidean distance */
    window_Q(&par->fs_mag[0],w_fs,&par->fs_mag[0],NUM_HARM,14);
    fsvq_enc(&par->fs_mag[0], &par->fs_mag[0], fs_vq_par);

#if (COMPLEXITY_COUNT)
complexity_count();
#endif

    /* Set MELP indeces to point to same array */
    par->fsvq_index = fs_vq_par.indices;

    /* Update MSVQ information */
    par->msvq_stages = vq_par.num_stages;
    par->msvq_bits = vq_par.num_bits;

    /* Write channel bitstream */
    melp_chn_write(par);

    /* Update delay buffers for next frame */
    v_equ(&speech[0],&speech[FRAME],IN_BEG);
    fpitch[BEGIN] = fpitch[END];

#if (COMPLEXITY_COUNT)
complexity_count();
#endif

} /* melp_ana */



/* 
 * melp_ana_init: perform initialization 
 */


void melp_ana_init()
{

    Shortword j;

    bpvc_ana_init(FRAME,PITCHMIN,PITCHMAX,NUM_BANDS,2,MINLENGTH);
    pitch_ana_init(PITCHMIN,PITCHMAX,FRAME,LPF_ORD,MINLENGTH);
    p_avg_init((Shortword)PDECAY_Q15,DEFAULT_PITCH_Q7,
	       (Shortword)PDECAY_PITCH_Q7,3);

    v_zap(speech,IN_BEG+FRAME);
    pitch_avg=DEFAULT_PITCH_Q7;
    fill(fpitch,DEFAULT_PITCH_Q7,2);
    v_zap(lpfsp_delin,LPF_ORD);
    v_zap(lpfsp_delout,LPF_ORD);

    /* Initialize multi-stage vector quantization (read codebook) */
	
    vq_par.num_best = MSVQ_M;
    vq_par.num_stages = 4;
    vq_par.dimension = 10;

    /* 
     * Allocate memory for number of levels per stage and indices
     * and for number of bits per stage 
     */
 
    MEM_ALLOC(MALLOC,vq_par.num_levels,vq_par.num_stages,Shortword);
    MEM_ALLOC(MALLOC,vq_par.indices,vq_par.num_stages,Shortword);
    MEM_ALLOC(MALLOC,vq_par.num_bits,vq_par.num_stages,Shortword);
	
    vq_par.num_levels[0] = 128;
    vq_par.num_levels[1] = 64;
    vq_par.num_levels[2] = 64;
    vq_par.num_levels[3] = 64;
	
    vq_par.num_bits[0] = 7;
    vq_par.num_bits[1] = 6;
    vq_par.num_bits[2] = 6;
    vq_par.num_bits[3] = 6;
	
    vq_par.cb = msvq_cb;
    vq_par.cb_mean = msvq_cb_mean;

    /* Initialize Fourier magnitude vector quantization (read codebook) */
	
    fs_vq_par.num_best = 1;
    fs_vq_par.num_stages = 1;
    fs_vq_par.dimension = NUM_HARM;

    /* 
     * Allocate memory for number of levels per stage and indices
     * and for number of bits per stage 
     */
 
    MEM_ALLOC(MALLOC,fs_vq_par.num_levels,fs_vq_par.num_stages,Shortword);
    MEM_ALLOC(MALLOC,fs_vq_par.indices,fs_vq_par.num_stages,Shortword);
    MEM_ALLOC(MALLOC,fs_vq_par.num_bits,fs_vq_par.num_stages,Shortword);

    fs_vq_par.num_levels[0] = FS_LEVELS;
    fs_vq_par.num_bits[0] = FS_BITS;
    fs_vq_par.cb = fsvq_cb;
	
    /* Initialize fixed MSE weighting and inverse of weighting */
	
    vq_fsw(w_fs, NUM_HARM, X60_Q9);
	
    /* Pre-weight codebook (assume single stage only) */
	
    if (fsvq_weighted == 0)
      {
	  fsvq_weighted = 1;
	  for (j = 0; j < fs_vq_par.num_levels[0]; j++)
            window_Q(&fs_vq_par.cb[j*NUM_HARM],w_fs,
		     &fs_vq_par.cb[j*NUM_HARM],NUM_HARM,14);
      }

    /* Initialize wr_array and wi_array */
    fs_init();
} /* melp_ana_init */