fixed-point: starting conversion of the new PLC
[opus.git] / libcelt / plc.c
1
2 #ifndef NEW_PLC
3 #define NEW_PLC
4 #endif
5
6 float _celt_lpc(
7       float       *lpc, /* out: [0...p-1] LPC coefficients      */
8 const float *ac,  /* in:  [0...p] autocorrelation values  */
9 int          p
10 )
11 {
12    int i, j;  
13    float r;
14    float error = ac[0];
15
16    if (ac[0] == 0)
17    {
18       for (i = 0; i < p; i++)
19          lpc[i] = 0;
20       return 0;
21    }
22    
23    for (i = 0; i < p; i++) {
24       
25       /* Sum up this iteration's reflection coefficient */
26       float rr = -ac[i + 1];
27       for (j = 0; j < i; j++) 
28          rr = rr - lpc[j]*ac[i - j];
29       r = rr/(error+1e-15);
30       /*  Update LPC coefficients and total error */
31       lpc[i] = r;
32       for (j = 0; j < i>>1; j++) 
33       {
34          float tmp  = lpc[j];
35          lpc[j]     = lpc[j    ] + r*lpc[i-1-j];
36          lpc[i-1-j] = lpc[i-1-j] + r*tmp;
37       }
38       if (i & 1) 
39          lpc[j] = lpc[j] + lpc[j]*r;
40       
41       error = error - r*r*error;
42       if (error<.00001*ac[0])
43          break;
44    }
45 #ifdef FIXED_POINT
46    for (i=0;i<p;i++)
47       lpc[i] = (1./4096)*floor(.5+4096*lpc[i]);
48 #endif
49    return error;
50 }
51
52 void fir(const celt_word16 *x,
53          const float *num,
54          celt_word16 *y,
55          int N,
56          int ord,
57          float *mem)
58 {
59    int i,j;
60
61    for (i=0;i<N;i++)
62    {
63       float sum = x[i];
64       for (j=0;j<ord;j++)
65       {
66          sum += num[j]*mem[j];
67       }
68       for (j=ord-1;j>=1;j--)
69       {
70          mem[j]=mem[j-1];
71       }
72       mem[0] = x[i];
73       y[i] = sum;
74    }
75 }
76
77 void iir(const celt_word32 *x,
78          const float *den,
79          celt_word32 *y,
80          int N,
81          int ord,
82          float *mem)
83 {
84    int i,j;
85    for (i=0;i<N;i++)
86    {
87       float sum = x[i];
88       for (j=0;j<ord;j++)
89       {
90          sum -= den[j]*mem[j];
91       }
92       for (j=ord-1;j>=1;j--)
93       {
94          mem[j]=mem[j-1];
95       }
96       mem[0] = sum;
97       y[i] = sum;
98    }
99 }
100
101 void _celt_autocorr(
102                    const celt_word16 *x,   /*  in: [0...n-1] samples x   */
103                    float       *ac,  /* out: [0...lag-1] ac values */
104                    const celt_word16       *window,
105                    int          overlap,
106                    int          lag, 
107                    int          n
108                   )
109 {
110    float d;
111    int i;
112    VARDECL(float, xx);
113    SAVE_STACK;
114    ALLOC(xx, n, float);
115    for (i=0;i<n;i++)
116       xx[i] = x[i];
117    for (i=0;i<overlap;i++)
118    {
119       xx[i] *= (1./Q15ONE)*window[i];
120       xx[n-i-1] *= (1./Q15ONE)*window[i];
121    }
122    while (lag>=0)
123    {
124       for (i = lag, d = 0; i < n; i++) 
125          d += x[i] * x[i-lag];
126       ac[lag] = d;
127       lag--;
128    }
129    ac[0] += 10;
130    RESTORE_STACK;
131 }