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

 1 typedef float　OsFlt32;
 2 typedef int　　OsInt32;
 3 
 4 OsFlt32 lpc(const OsFlt32 *r,OsInt32 p,OsFlt32 *a)
 5 {
 6     OsInt32 i,j;
 7     OsFlt32 err;
 8 
 9     if(0 == r[0])
10     {
11         for(i = 0; i < p; i++) a[i] = 0;
12         return 0;
13     }
14     a[0] = 1.0;
15     err = r[0];
16     for(i = 0; i < p; i++)
17     {
18         OsFlt32 lambda = 0.0;
19         for(j = 0; j <= i; j++)
20             lambda -= a[j]*r[i+1-j];
21         lambda /= err;
22         // Update LPC coefficients and total error
23         for(j = 0; j <= (i+1)/2; j++)
24         {
25             OsFlt32 temp = a[i+1-j] + lambda * a[j];
26             a[j] = a[j] + lambda * a[i+1-j];
27             a[i+1-j] = temp;
28         }
29         err *= (1.0 - lambda*lambda);
30     }
31     return err;
32 }
33 
34 void autocorr(const OsFlt32 *x,OsInt32 n,OsFlt32 *r,OsInt32 k)
35 {
36     OsFlt32 d;
37     OsInt32 i,p;
38 
39     for(p = 0; p <= k; p++)
40     {
41         for(i = 0,d = 0.0; i < n-p; i++)
42             d += x[i] * x[i+p];
43         r[p] = d;
44     }
45 }

 1 #include "lpc.h"
 2 
 3 int main(int argc,char **argv)
 4 {
 5     OsInt32 nLen = 128;
 6     OsFlt32 *pOriginal,*pPredicted;
 7     OsInt32 i,j;
 8     const OsInt32 order = 4;
 9     OsFlt32 R[order+1] = {0.0};
10     OsFlt32 A[order+1] = {0.0};
11     OsFlt32 error;
12 
13     pOriginal   = (OsFlt32*)calloc(nLen,sizeof(OsFlt32));
14     pPredicted  = (OsFlt32*)calloc(nLen,sizeof(OsFlt32));
15 
16     for(i = 0; i < nLen; i++)
17         pOriginal[i] = sin(i*0.01) + 0.75 * sin(i*0.03) + 0.5 * sin(i*0.05) + 0.25 * sin(i*0.11);
18 
19     autocorr(pOriginal,nLen,R,order);
20     lpc(R,order,A);
21     
22     for(i = 1; i <= order; i++)
23         A[i-1] = A[i];
24 
25     for(i = order; i < nLen; i++)
26     {
27         pPredicted[i] = 0.0;
28         for(j = 0; j < order; j++)
29             pPredicted[i] -= A[j] * pOriginal[i-1-j];
30     }
31     
32     error = 0;
33     for(i = order; i < nLen; i++)
34     {
35         double delta = pPredicted[i] - pOriginal[i];
36         printf( "Index: %.2d / Original: %.6f / Predicted: %.6f\n",i,pOriginal[i],pPredicted[i]);
37         error += delta * delta;
38     }
39     printf("Forward Linear Prediction Approximation Error: %f\n",error);
40 
41     free(pPredicted);
42     free(pOriginal);
43     return 0;
44 }