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