首页 > 学院 > 开发设计 > 正文

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

2019-11-06 09:57:45
字体:
来源:转载
供稿:网友

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

[html] view plain copy 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;      }  }  [html] view plain copy #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/n",i,pOriginal[i],pPredicted[i]);          error += delta * delta;      }      printf("Forward Linear Prediction Approximation Error: %f/n",error);        free(pPredicted);      free(pOriginal);      return 0;  } 
发表评论 共有条评论
用户名: 密码:
验证码: 匿名发表