openlpc.fixed.c

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

        xv40 = xv41; 	
#ifdef FAST_FILTERS        
        xv41 = ((u * 3) >> 1) +((u * 43) >> 9);		/* e(n)	, add 4 dB to compensate attenuation */
        yv40 = yv41; 
        yv41 = ((yv40 * 11) >> 4) + ((yv40 * 7) >> 10) 	 /* u(n) */
            + ((xv41 * 23) >> 5) + ((xv41 * 7) >> 11)
            - ((xv40 * 11) >> 4) - ((xv40 * 7) >> 10);
#else
        xv41 = fixmul32(u, ftofix32(1.584));
        yv40 = yv41; 
        yv41 = fixmul32(ftofix32(TAU/(1.0f+RHO+TAU)), yv40) 
            + fixmul32(ftofix32((RHO+TAU)/(1.0f+RHO+TAU)), xv41)
            - fixmul32(ftofix32(TAU/(1.0f+RHO+TAU)), xv40);
#endif
        u = yv41;
        
        st->s[j] = u;
    }
    st->xv2[0] = xv20;
    st->xv2[1] = xv21;
    st->yv2[0] = yv20;
    st->yv2[1] = yv21;
    st->xv4[0] = xv40;
    st->xv4[1] = xv41;
    st->yv4[0] = yv40;
    st->yv4[1] = yv41;
#endif
    
    /* operate windowing s[] -> w[] */
    
    for (i=0; i < st->buflen; i++)
        st->w[i] = fixmul32(st->s[i], st->h[i]);
    
    /* compute LPC coeff. from autocorrelation (first 11 values) of windowed data */
    auto_correl2(st->w, st->buflen, st->r);
    durbin(st->r, LPC_FILTORDER, k, &gain);
    
    /* calculate pitch */
    calc_pitch(st->y, st->framelen, &per1);                 /* first 2/3 of buffer */ 
    calc_pitch(st->y + st->buflen - st->framelen, st->framelen, &per2); /* last 2/3 of buffer */
    if(per1 > 0 && per2 > 0)
        per = (per1+per2) / 2;
    else if(per1 > 0)
        per = per1;
    else if(per2 > 0)
        per = per2;
    else
        per = 0;

    /* logarithmic q.: 0 = MINPER, 256 = MAXPER */
    parm[0] = (unsigned char)(per == 0? 0 : (unsigned char)fixtoi32(fixdiv32(fixlog32(fixdiv32(per, itofix32(REAL_MINPER))), logmaxminper) * 256));

#ifdef LINEAR_G_Q
    i = fixtoi32(gain * 128);
    if(i > 255)
        i = 255;
#else
    i = fixtoi32(256 * fixlog32(itofix32(1) + fixmul32(ftofix32((2.718-1.f)/10.f), gain))); /* deriv = 5.82 allowing to reserve 2 bits */
    if(i > 255) i = 255;	 /* reached when gain = 10 */
    i = (i+2) & 0xfc;
#endif
    
    parm[1] = (unsigned char)i;
    
    if(per1 > 0)
        parm[1] |= 1;
    if(per2 > 0)
        parm[1] |= 2;
    
    for(j=2; j < sizeofparm; j++)
        parm[j] = 0;
    
    for (i=0; i < LPC_FILTORDER; i++) {
        int bitamount = parambits[i];
        int bitc8 = 8-bitamount;
        int q = (1 << bitc8);  /* quantum: 1, 2, 4... */
        fixed32 u = k[i+1];
        int iu;
        
#ifdef ARCSIN_Q
        if(i < 2) u = fixmul32(fixasin32(u), ftofix32(2.f/M_PI));
#endif
        u *= 127;
        if(u < 0) 
            u += ftofix32(0.6) * q;
        else
            u += ftofix32(0.4) * q; /* highly empirical! */
        
        iu = fixtoi32(u);
        iu = iu & 0xff; /* keep only 8 bits */
        
        /* make room at the left of parm array shifting left */
        for(j=sizeofparm-1; j >= 3; j--) {
            parm[j] = (unsigned char)((parm[j] << bitamount) | (parm[j-1] >> bitc8));
        }
        parm[2] = (unsigned char)((parm[2] << bitamount) | (iu >> bitc8)); /* parm[2] */
    }
    bcopy(st->s + st->framelen, st->s, (st->buflen - st->framelen)*sizeof(st->s[0]));
    bcopy(st->y + st->framelen, st->y, (st->buflen - st->framelen)*sizeof(st->y[0]));
    
    return sizeofparm;
}

openlpc_decoder_state *create_openlpc_decoder_state(void)
{
    openlpc_decoder_state *state;
    
    state = (openlpc_decoder_state *)malloc(sizeof(openlpc_decoder_state));
    
    return state;
}

void init_openlpc_decoder_state(openlpc_decoder_state *st, int framelen)
{
    int i, j;
    
    st->Oldper = 0;
    st->OldG = 0;
    for (i = 0; i <= LPC_FILTORDER; i++) {
        st->Oldk[i] = 0;
        st->bp[i] = 0;
    }
    st->pitchctr = 0;
    st->exc = 0;
    logmaxminper = fixlog32(fixdiv32(itofix32(MAXPER), itofix32(MINPER)));
    
    for(i=0, j=0; i<sizeof(parambits) / sizeof(parambits[0]); i++) {
        j += parambits[i];
    }
    sizeofparm = (j + 7) / 8 + 2;
    
    /* test for a valid frame len? */
    st->framelen = framelen;
    st->buflen = framelen * 3 / 2;
    st->gainadj = fixsqrt32(itofix32(3) / st->buflen);
}

#define MIDTAP 1
#define MAXTAP 4
static short y[MAXTAP+1]={-21161, -8478, 30892,-10216, 16950};
static int j=MIDTAP, k=MAXTAP;

__inline int random16 (void)
{
    int the_random;
    
    y[k] = (short)(y[k] + y[j]);
    
    the_random = y[k];
    k--;
    if (k < 0) k = MAXTAP;
    j--;
    if (j < 0) j = MAXTAP;
    
    return(the_random);
}

/* LPC Synthesis (decoding) */

int openlpc_decode(unsigned char *parm, short *buf, openlpc_decoder_state *st)
{
    int i, j, flen=st->framelen;
    fixed32 per, gain, k[LPC_FILTORDER+1];
    fixed32 u, NewG, Ginc, Newper, perinc;
    fixed32 Newk[LPC_FILTORDER+1], kinc[LPC_FILTORDER+1];
    fixed32 gainadj;
    int hframe;
    fixed32 hper[2];
    int ii;
    fixed32 bp0, bp1, bp2, bp3, bp4, bp5, bp6, bp7, bp8, bp9, bp10;
    fixed32 stgain;
    
    bp0 = st->bp[0];
    bp1 = st->bp[1];
    bp2 = st->bp[2];
    bp3 = st->bp[3];
    bp4 = st->bp[4];
    bp5 = st->bp[5];
    bp6 = st->bp[6];
    bp7 = st->bp[7];
    bp8 = st->bp[8];
    bp9 = st->bp[9];
    bp10 = st->bp[10];
    stgain = st->gainadj;
    
    per = itofix32(parm[0]);
    
    per = (fixed32)(per == 0? 0: REAL_MINPER * fixexp32(fixmul32(per/256, logmaxminper)));
    
    hper[0] = hper[1] = per;
    
    if((parm[1] & 0x1) == 0) hper[0] = 0;
    if((parm[1] & 0x2) == 0) hper[1] = 0;
    
#ifdef LINEAR_G_Q
    gain = itofix32(parm[1]) / 128;
#else
    gain = itofix32(parm[1]) / 256;
    gain = fixdiv32((fixexp32(gain) - itofix32(1)), ftofix32((2.718-1.f)/10));
#endif
    
    k[0] = 0;
    
    for (i=LPC_FILTORDER-1; i >= 0; i--) {
        int bitamount = parambits[i];
        int bitc8 = 8-bitamount;
        /* casting to char should set the sign properly */
        char c = (char)(parm[2] << bitc8);
        
        for(j=2; j<sizeofparm; j++)
            parm[j] = (unsigned char)((parm[j] >> bitamount) | (parm[j+1] << bitc8)); 
        
        k[i+1] = itofix32(c) / 128;
#ifdef ARCSIN_Q
        if(i<2) k[i+1] = fixsin32(fixmul32(ftofix32(M_PI/2), k[i+1]));
#endif
    }
    
    /* k[] are the same in the two subframes */
    for (i=1; i <= LPC_FILTORDER; i++) {
        Newk[i] = st->Oldk[i];
        kinc[i] = (k[i] - st->Oldk[i]) / flen;
    }
    
    /* Loop on two half frames */
    
    for(hframe=0, ii=0; hframe<2; hframe++) {
        
        Newper = st->Oldper;
        NewG = st->OldG;
        
        Ginc = (gain - st->OldG) / (flen / 2);
        per = hper[hframe];
        
        if (per == 0) {			 /* if unvoiced */
            gainadj = stgain;
        } else {
            gainadj = fixsqrt32(per / st->buflen);
        }
        
        /* Interpolate period ONLY if both old and new subframes are voiced, gain and K always */ 
        
        if (st->Oldper != 0 && per != 0) {
            perinc = (per - st->Oldper) / (flen / 2);
        } else {
            perinc = 0; 
            Newper = per;
        }
        
        if (Newper == 0) st->pitchctr = 0;
        
        for (i=0; i < flen / 2; i++, ii++) {
            fixed32 kj;

            if (Newper == 0) {
                u = fixmul32((random16() << (PRECISION - 15 - 1)), fixmul32(NewG, gainadj)); 
            } else {			/* voiced: send a delta every per samples */
                /* triangular excitation */
                if (st->pitchctr == 0) {
                    st->exc = fixmul32(NewG, gainadj >> 2);
                    st->pitchctr = fixtoi32(Newper);
                } else {
                    st->exc -= fixmul32(fixdiv32(NewG, Newper), gainadj >> 1);
                    st->pitchctr--;
                }
                u = st->exc;
            }
            /* excitation */
            kj = Newk[10];
            u -= fixmul32(kj, bp9);
            bp10 = bp9 + fixmul32(kj, u);
            
            kj = Newk[9];
            u -= fixmul32(kj, bp8);
            bp9 = bp8 + fixmul32(kj, u);
            
            kj = Newk[8];
            u -= fixmul32(kj, bp7);
            bp8 = bp7 + fixmul32(kj, u);
            
            kj = Newk[7];
            u -= fixmul32(kj, bp6);
            bp7 = bp6 + fixmul32(kj, u);
            
            kj = Newk[6];
            u -= fixmul32(kj, bp5);
            bp6 = bp5 + fixmul32(kj, u);
            
            kj = Newk[5];
            u -= fixmul32(kj, bp4);
            bp5 = bp4 + fixmul32(kj, u);
            
            kj = Newk[4];
            u -= fixmul32(kj, bp3);
            bp4 = bp3 + fixmul32(kj, u);
            
            kj = Newk[3];
            u -= fixmul32(kj, bp2);
            bp3 = bp2 + fixmul32(kj, u);
            
            kj = Newk[2];
            u -= fixmul32(kj, bp1);
            bp2 = bp1 + fixmul32(kj, u);
            
            kj = Newk[1];
            u -= fixmul32(kj, bp0);
            bp1 = bp0 + fixmul32(kj, u);
            
            bp0 = u;
            
            if (u  < ftofix32(-0.9999)) {
                u = ftofix32(-0.9999);
            } else if (u > ftofix32(0.9999)) {
                u = ftofix32(0.9999);
            }
            buf[ii] = (short)(u >> (PRECISION - 15));
            
            Newper += perinc;
            NewG += Ginc;
            for (j=1; j <= LPC_FILTORDER; j++) Newk[j] += kinc[j];
            
        }
        
        st->Oldper = per;
        st->OldG = gain;
    }
    st->bp[0] = bp0;
    st->bp[1] = bp1;
    st->bp[2] = bp2;
    st->bp[3] = bp3;
    st->bp[4] = bp4;
    st->bp[5] = bp5;
    st->bp[6] = bp6;
    st->bp[7] = bp7;
    st->bp[8] = bp8;
    st->bp[9] = bp9;
    st->bp[10] = bp10;
    
    for (j=1; j <= LPC_FILTORDER; j++) st->Oldk[j] = k[j];
    
    return flen;
}

void destroy_openlpc_decoder_state(openlpc_decoder_state *st)
{
    if(st != NULL)
    {
        free(st);
        st = NULL;
    }
}