lpc.cpp

在vc下的lpc源代码,此文件包含组成 lpc 的所有应用程序

// lpc.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include "iostream"
#include "math.h"
#include "stdlib.h"
#include "stdio.h"
#include "float.h"
#include "string.h"
using namespace std;

#define LPC_SAMPLES_PER_FRAME	    240		//每桢样本数160
#define LPC_FILTORDER			10		//预测系数10个
#define BUFLEN			((LPC_SAMPLES_PER_FRAME * 3) / 2)		//缓存,每次取出240个样本进行处理
struct lpc_encoder_state {
	float   w_s[BUFLEN], w_y[BUFLEN], w_h[BUFLEN], w_w[BUFLEN];
    float   fa[6], u, u1, yp1, yp2;
    int     vuv;
};

#define M_PI 3.14159265358979323846
#define FS				8000.0f					/* 采样率 */
#define DOWN			5   					/* 下采样,5个样本用一个代替*/
#define PITCHORDER		4						/* 5-1 */
#define FC				600.0f					/* Pitch analyzer filter cutoff 为最大基音频率的2倍，采样率*/
#define MINPIT			50.0f					/* Minimum pitch 频率*/
#define MAXPIT			300.0f					/* Maximum pitch 频率*/

#define MINPER			(int)(FS/(DOWN*MAXPIT)+.5f)		/* Minimum period 样本个数*/
#define MAXPER			(int)(FS/(DOWN*MINPIT)+.5f)		/* Maximum period 样本个数*/

#define WSCALE			1.5863f	/* Energy loss due to windowing，加窗后，能量损失，采用该系数进行能量补偿 */
#define GAIN_ADJUST		0.98f




static void lpc_auto_correl1(float *w, int n, float *r)//最长周期所占样本数,逆滤波器用
{
    int i, k;
    
    for (k=0; k <= MAXPER; k++, n--) 
	{
        r[k] = 0.0;
        for (i=0; i < n; i++)
            r[k] += (w[i] *  w[i+k]);
    }
}



static void lpc_auto_correl2(float *w, int n, float *r)//求预测系数用
{
    int i, k;
    
    for (k=0; k <= LPC_FILTORDER; k++, n--) 
	{
        r[k] = 0.0;
        for (i=0; i < n; i++) 
		{
            r[k] += (w[i] *  w[i+k]);
        }
    }
}



static void lpc_auto_correl3(float *w, int n, float *r)//求逆滤波器用
{
    int i, k;
    
    for (k=0; k <= PITCHORDER; k++, n--) 
	{
        r[k] = 0.0;
        for (i=0; i < n; i++)
		{
            r[k] += (w[i] *  w[i+k]);
        }
    }
}



static void lpc_durbin(float *r, int p, float *k, /*@out@*/float *g)//相关系数转成k，输出误差，作为解码的激励源
{
    int i, j;
    float a[LPC_FILTORDER + 1], at[LPC_FILTORDER + 1], e;
    
    for (i = 0; i <= p; i++)
        a[i] = at[i] = 0.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] -= a[j] * r[i - j];
        }
        if (fabs(e) < FLT_EPSILON) 
		{
            e = 0.0f;
            break;
        }
        k[i] /= e;
        a[i] = k[i];
        for (j = 1; j < i; j++)
            a[j] = at[j] + k[i] * at[i - j];
        e *= 1.0f - k[i] * k[i];
    }
    
	//for(int h=1;h<10;h++)
	//	cout<<a[h]<<endl;
	//cout<<k[10]<<endl;
    if (e < FLT_EPSILON) 
	{
        e = 0.0f;
    }
    *g = (float) sqrt(e);
}

static void lpc_inverse_filter(float *w, float *k)//求残差，其结果用来进行基因周期检测
{
    int i, j;
    float b[PITCHORDER + 1], bp[PITCHORDER + 1], f[PITCHORDER + 1];
    
    for (i = 0; i <= PITCHORDER; i++)
        b[i] = f[i] = bp[i] = 0.0;
    
    for (i = 0; i < BUFLEN / DOWN; i++) 
	{
        f[0] = b[0] = w[i];
        for (j = 1; j <= PITCHORDER; j++) 
		{
            f[j] = f[j - 1] + k[j] * bp[j - 1];
            b[j] = k[j] * f[j - 1] + bp[j - 1];
            bp[j - 1] = b[j - 1];
        }
        w[i] = f[PITCHORDER];
    }
}



static void lpc_calc_pitch(float *w, /*@out@*/float *per, int *vuv)
{
    int i, j, rpos;
    float d[BUFLEN / DOWN]; /* 下采样后的结果 */
    float k[PITCHORDER + 1];
    float r[MAXPER + 1];
    float g, rmax;
    float rval, rm, rp;
    float a, b, c, x, y;
   
    for (i = 0, j = 0; i < BUFLEN; j++) {
        d[j] = 0;
        for (rpos = 0; rpos < DOWN; rpos++) {
            d[j] += w[i++];
        }
   	
    }
    lpc_auto_correl3(d, BUFLEN / DOWN, r);//求四阶相关系数
    lpc_durbin(r, PITCHORDER, k, &g);//转化系数
    lpc_inverse_filter(d, k);//由预测系数对下采样后的样本进行预测，求出残差
    lpc_auto_correl1(d, BUFLEN / DOWN, r);//对残差求相关系数
    rpos = 0;
    rmax = 0.0;
    for (i = MINPER; i <= MAXPER; i++) {//求出相关系数最大的点
        if (r[i] > rmax) {
            rmax = r[i];
            rpos = i;
        }
    }
    //由于下采样了，得出的基因周期误差比较大，所以必须进行插值拟合，找出真正的基音周期
    rm = r[rpos - 1];
    rp = r[rpos + 1];
    
    a = 0.5f * rm - rmax + 0.5f * rp;
    b = -0.5f * rm * (2.0f * rpos + 1.0f) + 
        2.0f * rpos * rmax + 0.5f * rp * (1.0f - 2.0f * rpos);
    c = 0.5f * rm * (rpos * rpos + rpos) +
        rmax * (1.0f - rpos * rpos) + 0.5f * rp * (rpos * rpos - rpos);
    
    x = -b / (2.0f * a);
    y = a * x * x + b * x + c;
    x *= DOWN;//恢复到原来的样本点
    
    rmax = y;
    if (r[0] == 0.0f)//若完全不相关，则没有周期性
	{
        rval = 1.0;//表示不是元音
    } 
	else 
	{
        rval = rmax / r[0];//可能为元音
    }


    if ((rval >= 0.4 || (*vuv == 3 && rval >= 0.3)) && (x > 0.0f))
	{
        *per = x;
        *vuv = (*vuv & 1) * 2 + 1;
    } 
	else
	{
        *per = 0.0;
        *vuv = (*vuv & 1) * 2;
    }
}

//lpc_encoder_state *create_lpc_encoder_state()
//{
//    lpc_encoder_state *state;
//    
//    state = (lpc_encoder_state *)malloc(sizeof(lpc_encoder_state));
//    
//    return state;
//}

void init_lpc_encoder_state(lpc_encoder_state *st)//主要工作是低通滤波
{
    int i;
    float r, v, w, wcT;
    
    for (i = 0; i < BUFLEN; i++)
	{
        st->w_s[i] = 0.0;//样本点值
        st->w_h[i] = (float) (WSCALE * (0.54f - 0.46f *(float)cos(2.0 * M_PI * i / (BUFLEN - 1.0))));//hanming 窗
    }


    wcT = (float) (2 * M_PI * FC / FS);
    r = (float) (0.36891079 * wcT);
    v = (float) (0.18445539 * wcT);
    w = (float) (0.92307712 * wcT);
    st->fa[1] = (float) (-exp(-r));
    st->fa[2] = 1.0f + st->fa[1];
    st->fa[3] = (float) (-2.0 * exp(-v) * cos(w));
    st->fa[4] = (float) exp(-2.0 * v);
    st->fa[5] = 1.0f + st->fa[3] + st->fa[4];
    
    st->u1 = 0.0;
    st->yp1 = 0.0;
    st->yp2 = 0.0;
    st->vuv = 0;
}

int lpc_encode(const short *in, unsigned char *out, lpc_encoder_state *st)
{
    int 	i, j;
    float	r[LPC_FILTORDER + 1];
    float	per, G, k[LPC_FILTORDER + 1];
    
    for (j = 0;j< BUFLEN;j++/*i = 0, j = BUFLEN - LPC_SAMPLES_PER_FRAME; i < LPC_SAMPLES_PER_FRAME; i++, j++*/) 
	{
        st->w_s[j] = (float) (GAIN_ADJUST * ((*in++) / 32768.0f));//换到1范围内处理
        
        st->u = st->fa[2] * st->w_s[j] - st->fa[1] * st->u1;
        st->w_y[j] = st->fa[5] * st->u1 - st->fa[3] * st->yp1 - st->fa[4] * st->yp2;
        st->u1 = st->u;
        st->yp2 = st->yp1;
        st->yp1 = st->w_y[j];
    }
    
    lpc_calc_pitch(st->w_y, &per, &st->vuv);
    
    for (i = 0; i < BUFLEN; i++)
        st->w_w[i] = st->w_s[i] * st->w_h[i];
    lpc_auto_correl2(st->w_w, BUFLEN, r);
    lpc_durbin(r, LPC_FILTORDER, k, &G);
	//cout<<G<<" "<<per<<endl;
	
    //per = per  * 2;
    if(per > 255.0f)
    {
        per = 255.0f;
    }
    out[0] = (unsigned char)(per);//输出
    
    i = (int)(G * 256);
    if(i > 255) i = 255; //此时一般是辅音
    out[1] = (char)i;//输出

    for (i = 0; i < LPC_FILTORDER; i++) {
        float u = k[i + 1];
        
        if (u < -0.9999f) //这种情况出现比较少
		{
            u = -0.9999f;
        }
		else if (u > 0.9999f)
		{
            u = 0.9999f;
        }
        out[i + 2] =(unsigned char) ((char)(128.0f * u));//输出
    }
    
    return true;
}

void destroy_lpc_encoder_state(lpc_encoder_state *st)
{
    if(st != NULL)
    {
        free(st);
        st = NULL;
    }
}










struct lpc_decoder_state{
		float Oldper, OldG, Oldk[LPC_FILTORDER], bp[LPC_FILTORDER + 1];
		int pitchctr;
        float exc;
		int tri;
		int tri_1;
		int tri_2;
};


//lpc_decoder_state *create_lpc_decoder_state()
//{
//    lpc_decoder_state *state;
//    
//    state = (lpc_decoder_state *)malloc(sizeof(lpc_decoder_state));
//    
//    return state;
//}

void init_lpc_decoder_state(lpc_decoder_state *st)//解压初始化工作
{
    int i;
    
    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;
	st->tri=0;
	st->tri_1=0;
	st->tri_2=0;

}


__inline int PNrandom (void)//产生白噪声
{
    static int the_random = 2385;
    
    if((the_random & 1) == 1)
    {
        the_random = ((the_random >> 1) ^ 0xfff6) | 0x8000;
    }
    else
    { 
        the_random >>= 1;
    }
    
    return(the_random);
}

int lpc_decode(unsigned char *in,short *out, lpc_decoder_state *st)//解压方法:一般是接受本次参数,对上次参数进行合成
{
    int i;
    register float u, f, per, G, NewG, Ginc, Newper, perinc;
    float k[LPC_FILTORDER], Newk[LPC_FILTORDER];
    float bp0, bp1, bp2, bp3, bp4, bp5, bp6, bp7, bp8, bp9, bp10;
    float b, kj;
    
    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];