线性预测与Levinson-Durbin算法实现

　　在学习信号处理的时候，线性预测是一个比较难理解的知识点，为了加快很多朋友的理解，这里给出Levinson-Durbin算法的线性预测实现和一个测试Demo，Demo中很明确的把输入信号、预测信号、预测误差打印了出来，这样就能以最直观的方式，把线性预测的实现与作用展示出来。话不多说，直接上代码！

　　

typedef float　OsFlt32;
typedef int　　OsInt32;

OsFlt32 lpc(const OsFlt32 *r,OsInt32 p,OsFlt32 *a)
{
    OsInt32 i,j;
    OsFlt32 err;

    if(0 == r[0])
    {
        for(i = 0; i < p; i++) a[i] = 0;
        return 0;
    }
    a[0] = 1.0;
    err = r[0];
    for(i = 0; i < p; i++)
    {
        OsFlt32 lambda = 0.0;
        for(j = 0; j <= i; j++)
            lambda -= a[j]*r[i+1-j];
        lambda /= err;
        // Update LPC coefficients and total error
        for(j = 0; j <= (i+1)/2; j++)
        {
            OsFlt32 temp = a[i+1-j] + lambda * a[j];
            a[j] = a[j] + lambda * a[i+1-j];
            a[i+1-j] = temp;
        }
        err *= (1.0 - lambda*lambda);
    }
    return err;
}

void autocorr(const OsFlt32 *x,OsInt32 n,OsFlt32 *r,OsInt32 k)
{
    OsFlt32 d;
    OsInt32 i,p;

    for(p = 0; p <= k; p++)
    {
        for(i = 0,d = 0.0; i < n-p; i++)
            d += x[i] * x[i+p];
        r[p] = d;
    }
}

#include "lpc.h"

int main(int argc,char **argv)
{
    OsInt32 nLen = 128;
    OsFlt32 *pOriginal,*pPredicted;
    OsInt32 i,j;
    const OsInt32 order = 4;
    OsFlt32 R[order+1] = {0.0};
    OsFlt32 A[order+1] = {0.0};
    OsFlt32 error;

    pOriginal   = (OsFlt32*)calloc(nLen,sizeof(OsFlt32));
    pPredicted  = (OsFlt32*)calloc(nLen,sizeof(OsFlt32));

    for(i = 0; i < nLen; i++)
        pOriginal[i] = sin(i*0.01) + 0.75 * sin(i*0.03) + 0.5 * sin(i*0.05) + 0.25 * sin(i*0.11);

    autocorr(pOriginal,nLen,R,order);
    lpc(R,order,A);
    
    for(i = 1; i <= order; i++)
        A[i-1] = A[i];

    for(i = order; i < nLen; i++)
    {
        pPredicted[i] = 0.0;
        for(j = 0; j < order; j++)
            pPredicted[i] -= A[j] * pOriginal[i-1-j];
    }
    
    error = 0;
    for(i = order; i < nLen; i++)
    {
        double delta = pPredicted[i] - pOriginal[i];
        printf( "Index: %.2d / Original: %.6f / Predicted: %.6f
",i,pOriginal[i],pPredicted[i]);
        error += delta * delta;
    }
    printf("Forward Linear Prediction Approximation Error: %f
",error);

    free(pPredicted);
    free(pOriginal);
    return 0;
}