线性预测与Levinson-Durbin算法实现
2016-06-10 18:02
543 查看
在学习信号处理的时候,线性预测是一个比较难理解的知识点,为了加快很多朋友的理解,这里给出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;
}
}
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\n",i,pOriginal[i],pPredicted[i]); error += delta * delta; } printf("Forward Linear Prediction Approximation Error: %f\n",error); free(pPredicted); free(pOriginal); return 0; }
相关文章推荐
- LPC入门教程
- MUDOS详解
- LPC基础讨论--数据类型--object
- lpc1343 UART
- LPC 1343 ADC
- 语音的线性预测系数(Linear Prediction Coefficient,LPC)
- LPCXpresso54114基于MDK的初试
- Node.js 实用工具
- Android intent flag启动模式
- “no public key available” while upgrading using update-manager
- POJ 2762 Going from u to v or from v to u?【强连通Kosaraju+拓扑排序】
- 基于QtQuick2.0应用程序运行于XP系统的诸多问题
- rhel7文件的归档与压缩
- centos 忘记用户名密码该怎么修改用户密码
- LMS、NLMS最优步长理论分析与Speex回声消除可能的改进想法
- mysql修改表名,字段,增加字段,删除字段
- WPF学习- AllowDrop 用户控件启用拖放功能
- symfony could not load type 'datetime'
- character-RNN模型介绍以及代码解析
- 第二章 对象、消息、运行期