blob: 334da8fb18366c3cc383352c1f78c83212747c46 [file] [log] [blame]
Jean-Marc Valine14fe902009-12-11 00:07:31 -05001
Jean-Marc Valind69c1cb2009-12-28 00:34:29 -05002#ifndef NEW_PLC
3#define NEW_PLC
4#endif
Jean-Marc Valine14fe902009-12-11 00:07:31 -05005
Jean-Marc Valin456eab22010-06-16 22:38:57 -04006#ifdef FIXED_POINT
Jean-Marc Valinbd82ca82010-06-17 00:11:54 -04007static celt_word32 frac_div32(celt_word32 a, celt_word32 b)
8{
Jean-Marc Valin3a9699e2010-06-17 00:35:26 -04009 celt_word16 rcp;
10 celt_word32 result, rem;
11 int shift = 30-celt_ilog2(b);
12 a = SHL32(a,shift);
13 b = SHL32(b,shift);
14
15 /* 16-bit reciprocal */
16 rcp = ROUND16(celt_rcp(ROUND16(b,16)),2);
Jean-Marc Valinbd82ca82010-06-17 00:11:54 -040017 result = SHL32(MULT16_32_Q15(rcp, a),1);
18 rem = a-MULT32_32_Q31(result, b);
19 result += SHL32(MULT16_32_Q15(rcp, rem),1);
20 return result;
21}
Jean-Marc Valin456eab22010-06-16 22:38:57 -040022#else
Jean-Marc Valinbd82ca82010-06-17 00:11:54 -040023#define frac_div32(a,b) ((float)(a)/(b))
Jean-Marc Valin456eab22010-06-16 22:38:57 -040024#endif
25
26
Jean-Marc Valind69c1cb2009-12-28 00:34:29 -050027float _celt_lpc(
Jean-Marc Valin74128be2010-01-01 09:33:17 -050028 celt_word16 *_lpc, /* out: [0...p-1] LPC coefficients */
Jean-Marc Valin456eab22010-06-16 22:38:57 -040029const celt_word32 *ac, /* in: [0...p] autocorrelation values */
Jean-Marc Valine14fe902009-12-11 00:07:31 -050030int p
31)
32{
33 int i, j;
Jean-Marc Valin456eab22010-06-16 22:38:57 -040034 celt_word32 r;
35 celt_word32 error = ac[0];
Jean-Marc Valin74128be2010-01-01 09:33:17 -050036#ifdef FIXED_POINT
Jean-Marc Valin456eab22010-06-16 22:38:57 -040037 celt_word32 lpc[LPC_ORDER];
Jean-Marc Valin74128be2010-01-01 09:33:17 -050038#else
39 float *lpc = _lpc;
40#endif
Jean-Marc Valine14fe902009-12-11 00:07:31 -050041
Jean-Marc Valin456eab22010-06-16 22:38:57 -040042 for (i = 0; i < p; i++)
43 lpc[i] = 0;
44 if (ac[0] != 0)
Jean-Marc Valine14fe902009-12-11 00:07:31 -050045 {
Jean-Marc Valin456eab22010-06-16 22:38:57 -040046 for (i = 0; i < p; i++) {
47 /* Sum up this iteration's reflection coefficient */
48 celt_word32 rr = 0;
49 for (j = 0; j < i; j++)
50 rr += MULT32_32_Q31(lpc[j],ac[i - j]);
51 rr += SHR32(ac[i + 1],3);
Jean-Marc Valinbd82ca82010-06-17 00:11:54 -040052 r = -frac_div32(SHL32(rr,3), error);
Jean-Marc Valin456eab22010-06-16 22:38:57 -040053 /* Update LPC coefficients and total error */
54 lpc[i] = SHR32(r,3);
55 for (j = 0; j < (i+1)>>1; j++)
56 {
57 celt_word32 tmp1, tmp2;
58 tmp1 = lpc[j];
59 tmp2 = lpc[i-1-j];
60 lpc[j] = tmp1 + MULT32_32_Q31(r,tmp2);
61 lpc[i-1-j] = tmp2 + MULT32_32_Q31(r,tmp1);
62 }
63
64 error = error - MULT32_32_Q31(MULT32_32_Q31(r,r),error);
Jean-Marc Valinbd82ca82010-06-17 00:11:54 -040065 /* Bail out once we get 30 dB gain */
66#ifdef FIXED_POINT
67 if (error<SHR32(ac[0],10))
Jean-Marc Valin456eab22010-06-16 22:38:57 -040068 break;
Jean-Marc Valinbd82ca82010-06-17 00:11:54 -040069#else
70 if (error<.001*ac[0])
71 break;
72#endif
Jean-Marc Valine14fe902009-12-11 00:07:31 -050073 }
Jean-Marc Valine14fe902009-12-11 00:07:31 -050074 }
Jean-Marc Valin303b3b62009-12-30 22:40:24 -050075#ifdef FIXED_POINT
76 for (i=0;i<p;i++)
Jean-Marc Valin456eab22010-06-16 22:38:57 -040077 _lpc[i] = ROUND16(lpc[i],16);
Jean-Marc Valin303b3b62009-12-30 22:40:24 -050078#endif
Jean-Marc Valine14fe902009-12-11 00:07:31 -050079 return error;
80}
81
Jean-Marc Valin303b3b62009-12-30 22:40:24 -050082void fir(const celt_word16 *x,
Jean-Marc Valin74128be2010-01-01 09:33:17 -050083 const celt_word16 *num,
Jean-Marc Valin303b3b62009-12-30 22:40:24 -050084 celt_word16 *y,
Jean-Marc Valine14fe902009-12-11 00:07:31 -050085 int N,
86 int ord,
Jean-Marc Valin74128be2010-01-01 09:33:17 -050087 celt_word16 *mem)
Jean-Marc Valine14fe902009-12-11 00:07:31 -050088{
89 int i,j;
90
91 for (i=0;i<N;i++)
92 {
Jean-Marc Valin74128be2010-01-01 09:33:17 -050093 celt_word32 sum = SHL32(EXTEND32(x[i]), SIG_SHIFT);
Jean-Marc Valine14fe902009-12-11 00:07:31 -050094 for (j=0;j<ord;j++)
95 {
Jean-Marc Valin74128be2010-01-01 09:33:17 -050096 sum += MULT16_16(num[j],mem[j]);
Jean-Marc Valine14fe902009-12-11 00:07:31 -050097 }
98 for (j=ord-1;j>=1;j--)
99 {
100 mem[j]=mem[j-1];
101 }
102 mem[0] = x[i];
Jean-Marc Valin74128be2010-01-01 09:33:17 -0500103 y[i] = ROUND16(sum, SIG_SHIFT);
Jean-Marc Valine14fe902009-12-11 00:07:31 -0500104 }
105}
106
Jean-Marc Valind69c1cb2009-12-28 00:34:29 -0500107void iir(const celt_word32 *x,
Jean-Marc Valin74128be2010-01-01 09:33:17 -0500108 const celt_word16 *den,
Jean-Marc Valind69c1cb2009-12-28 00:34:29 -0500109 celt_word32 *y,
Jean-Marc Valine14fe902009-12-11 00:07:31 -0500110 int N,
111 int ord,
Jean-Marc Valin74128be2010-01-01 09:33:17 -0500112 celt_word16 *mem)
Jean-Marc Valine14fe902009-12-11 00:07:31 -0500113{
114 int i,j;
115 for (i=0;i<N;i++)
116 {
Jean-Marc Valin74128be2010-01-01 09:33:17 -0500117 celt_word32 sum = x[i];
Jean-Marc Valine14fe902009-12-11 00:07:31 -0500118 for (j=0;j<ord;j++)
119 {
Jean-Marc Valin74128be2010-01-01 09:33:17 -0500120 sum -= MULT16_16(den[j],mem[j]);
Jean-Marc Valine14fe902009-12-11 00:07:31 -0500121 }
122 for (j=ord-1;j>=1;j--)
123 {
124 mem[j]=mem[j-1];
125 }
Jean-Marc Valin74128be2010-01-01 09:33:17 -0500126 mem[0] = ROUND16(sum,SIG_SHIFT);
Jean-Marc Valine14fe902009-12-11 00:07:31 -0500127 y[i] = sum;
128 }
129}
130
131void _celt_autocorr(
Jean-Marc Valin303b3b62009-12-30 22:40:24 -0500132 const celt_word16 *x, /* in: [0...n-1] samples x */
Jean-Marc Valin456eab22010-06-16 22:38:57 -0400133 celt_word32 *ac, /* out: [0...lag-1] ac values */
Jean-Marc Valind69c1cb2009-12-28 00:34:29 -0500134 const celt_word16 *window,
Jean-Marc Valine14fe902009-12-11 00:07:31 -0500135 int overlap,
136 int lag,
137 int n
138 )
139{
140 float d;
141 int i;
Jean-Marc Valin456eab22010-06-16 22:38:57 -0400142 float scale=1;
Thorvald Natvigb8002a02009-12-11 13:19:09 +0100143 VARDECL(float, xx);
144 SAVE_STACK;
145 ALLOC(xx, n, float);
Jean-Marc Valine14fe902009-12-11 00:07:31 -0500146 for (i=0;i<n;i++)
147 xx[i] = x[i];
148 for (i=0;i<overlap;i++)
149 {
Jean-Marc Valind69c1cb2009-12-28 00:34:29 -0500150 xx[i] *= (1./Q15ONE)*window[i];
151 xx[n-i-1] *= (1./Q15ONE)*window[i];
Jean-Marc Valine14fe902009-12-11 00:07:31 -0500152 }
Jean-Marc Valin456eab22010-06-16 22:38:57 -0400153#ifdef FIXED_POINT
154 {
155 float ac0=0;
156 for(i=0;i<n;i++)
157 ac0 += x[i]*x[i];
158 ac0+=10;
159 scale = 2000000000/ac0;
160 }
161#endif
Jean-Marc Valine14fe902009-12-11 00:07:31 -0500162 while (lag>=0)
163 {
164 for (i = lag, d = 0; i < n; i++)
165 d += x[i] * x[i-lag];
Jean-Marc Valin456eab22010-06-16 22:38:57 -0400166 ac[lag] = d*scale;
167 /*printf ("%f ", ac[lag]);*/
Jean-Marc Valine14fe902009-12-11 00:07:31 -0500168 lag--;
169 }
Jean-Marc Valin456eab22010-06-16 22:38:57 -0400170 /*printf ("\n");*/
Jean-Marc Valine14fe902009-12-11 00:07:31 -0500171 ac[0] += 10;
Jean-Marc Valin456eab22010-06-16 22:38:57 -0400172
Thorvald Natvigb8002a02009-12-11 13:19:09 +0100173 RESTORE_STACK;
Jean-Marc Valine14fe902009-12-11 00:07:31 -0500174}