openlpc.c

HawkVoiceDI Project created and maintained by Phil Frisbie, Jr. <[e-mail]> HawkVoiceD

        yv12 = (float)((xv10 + xv12) - 2 * xv11
            + ( -0.8948742499f * yv10) + ( 1.8890389823f * yv11));
        
        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 */
        
        /*st->xv3[0] = (float)(u / 2.127814584e+001);*/ /* GAIN */
        xv30 = (float)(u * 0.04699658f); /* GAIN */
        yv30 = yv31;
        yv31 = yv32; 
        yv32 = xv30 + (float)(( -0.7166152306f * yv30) + (1.6696186545f * yv31));
        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++) {
        float 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) */
        xv21 = (float)(u * 1.584f);		/* e(n)	, add 4 dB to compensate attenuation */
        yv20 = yv21;
        yv21 = (float)(TAU/(1.f+RHO+TAU) * yv20 	 /* u(n) */
            + (RHO+TAU)/(1.f+RHO+TAU) * xv21
            - TAU/(1.f+RHO+TAU) * xv20);
        u = yv21;
        
        /* cascaded copy of handcoded filter: 1 zero at 640 Hz, 1 pole at 3200 */
        xv40 = xv41;
        xv41 = (float)(u * 1.584f);
        yv40 = yv41;
        yv41 = (float)(TAU/(1.f+RHO+TAU) * yv40
            + (RHO+TAU)/(1.f+RHO+TAU) * xv41
            - TAU/(1.f+RHO+TAU) * xv40);
        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] = 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)(log(per/(REAL_MINPER)) / logmaxminper * (1<<8)));
    
#ifdef LINEAR_G_Q
    i = gain * (1<<7);
    if(i > 255) 	/* bug fix by EM */
        i = 255;
#else
    i = (int)(float)(256.0f * log(1 + (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... */
        float u = k[i+1];
        int iu;
#ifdef ARCSIN_Q
        if(i < 2) u = (float)(asin(u)*2.f/M_PI);
#endif
        u *= 127;
        if(u < 0)
            u += (0.6f * q);
        else
            u += (0.4f * q); /* highly empirical! */
        
        iu = lrintf(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.0f;
    st->OldG = 0.0f;
    for (i = 0; i <= LPC_FILTORDER; i++) {
        st->Oldk[i] = 0.0f;
        st->bp[i] = 0.0f;
    }
    st->pitchctr = 0;
    st->exc = 0.0f;
    logmaxminper = (float)log((float)MAXPER/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;
}

/* LPC Synthesis (decoding) */

int openlpc_decode(unsigned char *parm, short *buf, openlpc_decoder_state *st)
{
    int i, j, flen=st->framelen;
    float per, gain, k[LPC_FILTORDER+1];
    float f, u, newgain, Ginc, Newper, perinc;
    float Newk[LPC_FILTORDER+1], kinc[LPC_FILTORDER+1];
    float gainadj;
    int hframe;
    float hper[2];
    int ii;
    float bp0, bp1, bp2, bp3, bp4, bp5, bp6, bp7, bp8, bp9, bp10;
    
    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];
    
    per = (float)(parm[0]);
    
    per = (float)(per == 0? 0: REAL_MINPER * exp(per/(1<<8) * 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 = (float)parm[1] / (1<<7);
#else
    gain = (float)parm[1] / 256.f;
    gain = (float)((exp(gain) - 1)/((2.718-1.f)/10));
#endif
    
    k[0] = 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] = ((float)c / (1<<7));
#ifdef ARCSIN_Q
        if(i<2) k[i+1] = (float)sin(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;
        newgain = st->OldG;
        
        Ginc = (gain - st->OldG) / (flen/2);
        per = hper[hframe];
        
        if (per == 0.0) {			 /* if unvoiced */
            gainadj = /* 1.5874 * */ (float)sqrt(3.0f/st->buflen);
        } else {
            gainadj = (float)sqrt(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.0f; 
            Newper = per;
        }
        
        if (Newper == 0.f) st->pitchctr = 0;
        
        for (i=0; i < flen/2; i++, ii++) {
            float b, kj;
            if (Newper == 0.f) {
                u = (float)(((rand()*(1/(1.0f+RAND_MAX))) - 0.5f ) * newgain * gainadj); 
            } else {			/* voiced: send a delta every per samples */
                /* triangular excitation */
                if (st->pitchctr == 0) {
                    st->exc = newgain * 0.25f * gainadj;
                    st->pitchctr = (int) Newper;
                } else {
                    st->exc -= newgain/Newper * 0.5f * gainadj;
                    st->pitchctr--;
                }
                u = st->exc;
            }
            f = u;
	        /* excitation */
            b = bp9;
            kj = Newk[10];
            f -= kj * bp9;
            bp10 = bp9 + kj * f;

            kj = Newk[9];
            f -= kj * bp8;
            bp9 = bp8 + kj * f;

            kj = Newk[8];
            f -= kj * bp7;
            bp8 = bp7 + kj * f;

            kj = Newk[7];
            f -= kj * bp6;
            bp7 = bp6 + kj * f;

            kj = Newk[6];
            f -= kj * bp5;
            bp6 = bp5 + kj * f;

            kj = Newk[5];
            f -= kj * bp4;
            bp5 = bp4 + kj * f;

            kj = Newk[4];
            f -= kj * bp3;
            bp4 = bp3 + kj * f;

            kj = Newk[3];
            f -= kj * bp2;
            bp3 = bp2 + kj * f;

            kj = Newk[2];
            f -= kj * bp1;
            bp2 = bp1 + kj * f;

            kj = Newk[1];
            f -= kj * bp0;
            bp1 = bp0 + kj * f;

            bp0 = f;
            u = f;

            if (u  < -0.9999f) {
                u = -0.9999f;
            } else if (u > 0.9999f) {
                u = 0.9999f;
            }

            buf[ii] = (short)lrintf(u * 32767.0f);

            Newper += perinc;
            newgain += 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;
    }
}