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