openlpc.fixed.c

LPC10 2.4Kpbs 语音压缩定点运算C语言源程序OPENLPC

    temp >>= PRECISION;
    temp = (temp - ftofix32(.0170881256f)) * x;
    temp >>= PRECISION;
    temp = (temp + ftofix32(.0308918810f)) * x;
    temp >>= PRECISION;
    temp = (temp - ftofix32(.0501743046f)) * x;
    temp >>= PRECISION;
    temp = (temp + ftofix32(.0889789874f)) * x;
    temp >>= PRECISION;
    temp = (temp - ftofix32(.2145988016f)) * x;
    temp >>= PRECISION;
    temp = (temp + ftofix32(1.570796305f)) * fixsqrt32(itofix32(1) - x);
    temp >>= PRECISION;

    return sign * (ftofix32(M_PI/2) - (fixed32)temp);
}

#define PREEMPH

#define bcopy(a, b, n)	  memmove(b, a, n)

#define LPC_FILTORDER   10
#define FS              8000.0	/* Sampling rate */
#define MAXWINDOW       1000	/* Max analysis window length */


typedef struct openlpc_e_state{
    int     framelen, buflen;
    fixed32 s[MAXWINDOW], y[MAXWINDOW], h[MAXWINDOW];
    fixed32 xv1[3], yv1[3],
            xv2[2], yv2[2],
            xv3[1], yv3[3],
            xv4[2], yv4[2];
    fixed32 w[MAXWINDOW], r[LPC_FILTORDER+1];
} openlpc_e_state_t;

typedef struct openlpc_d_state{
    fixed32 Oldper, OldG, Oldk[LPC_FILTORDER + 1];
    fixed32 bp[LPC_FILTORDER+1];
    fixed32 exc;
    fixed32 gainadj;
    int pitchctr, framelen, buflen;
} openlpc_d_state_t;

#define FC		200.0	/* Pitch analyzer filter cutoff */
#define DOWN		5	/* Decimation for pitch analyzer */
#define MINPIT		40.0	/* Minimum pitch (observed: 74) */
#define MAXPIT		320.0	/* Maximum pitch (observed: 250) */

#define MINPER		(int)(FS / (DOWN * MAXPIT) + .5)	/* Minimum period  */
#define MAXPER		(int)(FS / (DOWN * MINPIT) + .5)	/* Maximum period  */

#define REAL_MINPER	 (DOWN * MINPER) /* converted to samples units */

#define WSCALE		1.5863	/* Energy loss due to windowing */

#define BITS_FOR_LPC 38

#define ARCSIN_Q /* provides better quantization of first two k[] at low bitrates */

#if BITS_FOR_LPC == 38
/* (38 bit LPC-10, 2.7 Kbit/s @ 20ms, 2.4 Kbit/s @ 22.5 ms */
static int parambits[LPC_FILTORDER] = {6,5,5,4,4,3,3,3,3,2};
#elif BITS_FOR_LPC == 32
/* (32 bit LPC-10, 2.4 Kbit/s, not so good */
static int parambits[LPC_FILTORDER] = {5,5,5,4,3,3,2,2,2,1};
#else /* BITS_FOR_LPC == 80	*/
/* 80-bit LPC10, 4.8 Kbit/s */
static int parambits[LPC_FILTORDER] = {8,8,8,8,8,8,8,8,8,8};
#endif

static fixed32 logmaxminper;
static int sizeofparm;	/* computed by lpc_init */

static void auto_correl1(fixed32 *w, int n, fixed32 *r)
{
    int i, k;
    fixed64 temp, temp2;
    
    for (k=0; k <= MAXPER; k++, n--) {
        temp = 0;
        for (i=0; i < n; i++) {
            temp2 = w[i];
            temp += temp2 * w[i+k];
        }
        r[k] = (fixed32)(temp >> PRECISION);
    }
}

static void auto_correl2(fixed32 *w, int n, fixed32 *r)
{
    int i, k;
    fixed64 temp, temp2;
    
    for (k=0; k <= LPC_FILTORDER; k++, n--) {
        temp = 0;
        for (i=0; i < n; i++) {
            temp2 = w[i];
            temp += temp2 * w[i+k];
        }
        r[k] = (fixed32)(temp >> PRECISION);
    }
}

static void durbin(fixed32 r[], int p, fixed32 k[], fixed32 *g)
{
    int i, j;
    fixed32 a[LPC_FILTORDER+1], at[LPC_FILTORDER+1], e;
    
    for (i=0; i <= p; i++)
        a[i] = at[i] = 0;
    
    e = r[0];
    for (i=1; i <= p; i++) {
        k[i] = -r[i];
        for (j=1; j < i; j++) {
            at[j] = a[j];
            k[i] -= fixmul32(a[j], r[i-j]);
        }
        if (e == 0) {  /* fix by John Walker */
            *g = 0;
            return;
        }
        k[i] = fixdiv32(k[i], e);
        a[i] = k[i];
        for (j=1; j < i; j++)
            a[j] = at[j] + fixmul32(k[i], at[i-j]);
        e = fixmul32(e, (itofix32(1) - fixmul32(k[i], k[i])));
    }
    if (e < 0) {
        e = 0; /* fix by John Walker */
    }
    *g = fixsqrt32(e);
}

static void calc_pitch(fixed32 w[], int len, fixed32 *per)
{
    int i, j, rpos;
    fixed32 d[MAXWINDOW / DOWN], r[MAXPER + 1], rmax;
    fixed32 rval, rm, rp;
    fixed32 x, y;
    fixed32 vthresh;
    
    /* decimation */
    for (i=0, j=0; i < len; i+=DOWN)
        d[j++] = w[i];
    
    auto_correl1(d, len / DOWN, r);
    
    /* find peak between MINPER and MAXPER */
    x = itofix32(1);
    rpos = 0;
    rmax = 0;
    
    for (i = 1; i <= MAXPER; i++) {
        rm = r[i-1];
        rp = r[i+1];
        y = rm+r[i]+rp; /* find max of integral from i-1 to i+1 */
        if (y > rmax && r[i] > rm && r[i] > rp &&  i > MINPER) {
            rmax = y;
            rpos = i;
        }
    }
    
    /* consider adjacent values */
    rm = r[rpos-1];
    rp = r[rpos+1];
    
    if(rpos > 0) {
        x = fixdiv32(((rpos-1) * rm + rpos * r[rpos] + (rpos+1) * rp), (rm+r[rpos]+rp)); 
    }
    /* normalize, so that 0. < rval < 1. */ 
    rval = (r[0] == 0 ? 0 : fixdiv32(r[rpos], r[0]));
    
    
    /* periods near the low boundary and at low volumes
    are usually spurious and 
    manifest themselves as annoying mosquito buzzes */
    
    *per = 0;	/* default: unvoiced */
    if ( x > itofix32(MINPER) &&  /* x could be < MINPER or even < 0 if rpos == MINPER */
        x < itofix32(MAXPER + 1) /* same story */
        ) {
        
        vthresh = ftofix32(0.6); 
        if(r[0] > ftofix32(0.002))	   /* at low volumes (< 0.002), prefer unvoiced */ 
            vthresh = ftofix32(0.25);    /* drop threshold at high volumes */
        
        if(rval > vthresh)
            *per = x * DOWN;
    }
}

/* Initialization of various parameters */

openlpc_encoder_state *create_openlpc_encoder_state(void)
{
    openlpc_encoder_state *state;
    
    state = (openlpc_encoder_state *)malloc(sizeof(openlpc_encoder_state));
    
    return state;
}


void init_openlpc_encoder_state(openlpc_encoder_state *st, int framelen)
{
    int i, j;
    
    st->framelen = framelen;
    memset(st->y, 0, sizeof(st->y));
    st->buflen = framelen * 3 / 2;
    /*  (st->buflen > MAXWINDOW) return -1;*/
    
    for(i=0, j=0; i<sizeof(parambits)/sizeof(parambits[0]); i++) {
        j += parambits[i];
    }
    sizeofparm = (j + 7) / 8 + 2;
    for (i = 0; i < st->buflen; i++) {
        st->s[i] = 0;
        /* this is only calculated once, but used each frame, */
        /* so we will use floating point for accuracy */
        st->h[i] = ftofix32(WSCALE*(0.54 - 0.46 * cos(2 * M_PI * i / (st->buflen-1.0))));
    }
    /* init the filters */
    st->xv1[0] = st->xv1[1] = st->xv1[2] = st->yv1[0] = st->yv1[1] = st->yv1[2] = 0;
    st->xv2[0] = st->xv2[1] = st->yv2[0] = st->yv2[1] = 0;
    st->xv3[0] = st->yv3[0] = st->yv3[1] = st->yv3[2] = 0;
    st->xv4[0] = st->xv4[1] = st->yv4[0] = st->yv4[1] = 0;

    logmaxminper = fixlog32(fixdiv32(itofix32(MAXPER), itofix32(MINPER)));
}

void destroy_openlpc_encoder_state(openlpc_encoder_state *st)
{
    if(st != NULL)
    {
        free(st);
        st = NULL;
    }
}

/* LPC Analysis (compression) */

int openlpc_encode(const short *buf, unsigned char *parm, openlpc_encoder_state *st)
{
    int i, j;
    fixed32 per, gain, k[LPC_FILTORDER+1];
    fixed32 per1, per2;
    fixed32 xv10, xv11, xv12, yv10, yv11, yv12, xv30, yv30, yv31, yv32;
#ifdef PREEMPH
    fixed32 xv20, xv21, yv20, yv21, xv40, xv41, yv40, yv41;
#endif
    
    xv10 = st->xv1[0];
    xv11 = st->xv1[1];
    xv12 = st->xv1[2];
    yv10 = st->yv1[0];
    yv11 = st->yv1[1];
    yv12 = st->yv1[2];
    xv30 = st->xv3[0];
    yv30 = st->yv3[0];
    yv31 = st->yv3[1];
    yv32 = st->yv3[2];
    /* convert short data in buf[] to signed lin. data in s[] and prefilter */
    for (i=0, j=st->buflen - st->framelen; i < st->framelen; i++, j++) {
        
        /* special handling here for the intitial conversion */
        fixed32 u = (fixed32)(buf[i] << (PRECISION - 15));
        
        /* Anti-hum 2nd order Butterworth high-pass, 100 Hz corner frequency */
        /* Digital filter designed by mkfilter/mkshape/gencode   A.J. Fisher
        mkfilter -Bu -Hp -o 2 -a 0.0125 -l -z */
        
        xv10 = xv11;
        xv11 = xv12; 
#ifdef FAST_FILTERS        
        xv12 = ((u * 15) >> 4) + (u >> 7) + ((u * 11) >> 14); /* /GAIN */
        yv10 = yv11;
        yv11 = yv12;
        yv12 = (fixed32)((xv10 + xv12) - (xv11 + xv11)
            - ((yv10 * 7) >> 3) - ((yv10 * 5) >> 8)
            + ((yv11 * 15) >> 3) + (yv11 >> 6) );
#else
        xv12 = fixmul32(u, ftofix32(0.94597831)); /* /GAIN */
        yv10 = yv11;
        yv11 = yv12;
        yv12 = (fixed32)((xv10 + xv12) - (xv11 + xv11)
            + fixmul32(ftofix32(-0.8948742499), yv10) + fixmul32(ftofix32(1.8890389823), yv11));
#endif
        u = st->s[j] = yv12;	/* also affects input of next stage, to the LPC filter synth */
        
        /* low-pass filter s[] -> y[] before computing pitch */
        /* second-order Butterworth low-pass filter, corner at 300 Hz */
        /* Digital filter designed by mkfilter/mkshape/gencode   A.J. Fisher
        MKFILTER.EXE -Bu -Lp -o 2 -a 0.0375 -l -z */
#ifdef FAST_FILTERS        
        xv30 = ((u * 3) >> 6) + (u >> 13); /* GAIN */
        yv30 = yv31;
        yv31 = yv32;
        yv32 = xv30 - ((yv30 * 23) >> 5) + (yv30 >> 9)
            + ((yv31 * 107) >> 6) - (yv31 >> 9);
#else
        xv30 = fixmul32(u, ftofix32(0.04699658)); /* GAIN */
        yv30 = yv31;
        yv31 = yv32; 
        yv32 = xv30 + fixmul32(ftofix32(-0.7166152306), yv30) + fixmul32(ftofix32(1.6696186545), yv31);
#endif
        st->y[j] = yv32;
    }
    st->xv1[0] = xv10;
    st->xv1[1] = xv11;
    st->xv1[2] = xv12;
    st->yv1[0] = yv10;
    st->yv1[1] = yv11;
    st->yv1[2] = yv12;
    st->xv3[0] = xv30;
    st->yv3[0] = yv30;
    st->yv3[1] = yv31;
    st->yv3[2] = yv32;
#ifdef PREEMPH
    /* operate optional preemphasis s[] -> s[] on the newly arrived frame */
    xv20 = st->xv2[0];
    xv21 = st->xv2[1];
    yv20 = st->yv2[0];
    yv21 = st->yv2[1];
    xv40 = st->xv4[0];
    xv41 = st->xv4[1];
    yv40 = st->yv4[0];
    yv41 = st->yv4[1];
    for (j=st->buflen - st->framelen; j < st->buflen; j++) {
        fixed32 u = st->s[j];
        
        /* handcoded filter: 1 zero at 640 Hz, 1 pole at 3200 */
#define TAU (FS / 3200.f)
#define RHO (0.1f)
        xv20 = xv21; 	/* e(n-1) */
#ifdef FAST_FILTERS        
        xv21 = ((u * 3) >> 1) +((u * 43) >> 9);		/* e(n)	, add 4 dB to compensate attenuation */
        yv20 = yv21; 
        yv21 = ((yv20 * 11) >> 4) + ((yv20 * 7) >> 10) 	 /* u(n) */
            + ((xv21 * 23) >> 5) + ((xv21 * 7) >> 11)
            - ((xv20 * 11) >> 4) - ((xv20 * 7) >> 10);
#else
        xv21 = fixmul32(u, ftofix32(1.584));		/* e(n)	, add 4 dB to compensate attenuation */
        yv20 = yv21; 
        yv21 = fixmul32(ftofix32(TAU/(1.0f+RHO+TAU)), yv20) 	 /* u(n) */
            + fixmul32(ftofix32((RHO+TAU)/(1.0f+RHO+TAU)), xv21)
            - fixmul32(ftofix32(TAU/(1.0f+RHO+TAU)), xv20);
#endif
        u = yv21;
        
        /* cascaded copy of handcoded filter: 1 zero at 640 Hz, 1 pole at 3200 */