Making new PLC code work in fixed-point even though it's still using float
[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    return error;
46 }
47
48 void fir(const float *x,
49          const float *num,
50          float *y,
51          int N,
52          int ord,
53          float *mem)
54 {
55    int i,j;
56
57    for (i=0;i<N;i++)
58    {
59       float sum = x[i];
60       for (j=0;j<ord;j++)
61       {
62          sum += num[j]*mem[j];
63       }
64       for (j=ord-1;j>=1;j--)
65       {
66          mem[j]=mem[j-1];
67       }
68       mem[0] = x[i];
69       y[i] = sum;
70    }
71 }
72
73 void iir(const celt_word32 *x,
74          const float *den,
75          celt_word32 *y,
76          int N,
77          int ord,
78          float *mem)
79 {
80    int i,j;
81    for (i=0;i<N;i++)
82    {
83       float sum = x[i];
84       for (j=0;j<ord;j++)
85       {
86          sum -= den[j]*mem[j];
87       }
88       for (j=ord-1;j>=1;j--)
89       {
90          mem[j]=mem[j-1];
91       }
92       mem[0] = sum;
93       y[i] = sum;
94    }
95 }
96
97 void _celt_autocorr(
98                    const float *x,   /*  in: [0...n-1] samples x   */
99                    float       *ac,  /* out: [0...lag-1] ac values */
100                    const celt_word16       *window,
101                    int          overlap,
102                    int          lag, 
103                    int          n
104                   )
105 {
106    float d;
107    int i;
108    VARDECL(float, xx);
109    SAVE_STACK;
110    ALLOC(xx, n, float);
111    for (i=0;i<n;i++)
112       xx[i] = x[i];
113    for (i=0;i<overlap;i++)
114    {
115       xx[i] *= (1./Q15ONE)*window[i];
116       xx[n-i-1] *= (1./Q15ONE)*window[i];
117    }
118    while (lag>=0)
119    {
120       for (i = lag, d = 0; i < n; i++) 
121          d += x[i] * x[i-lag];
122       ac[lag] = d;
123       lag--;
124    }
125    ac[0] += 10;
126    RESTORE_STACK;
127 }