Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 1 | |
Jean-Marc Valin | d69c1cb | 2009-12-28 00:34:29 -0500 | [diff] [blame] | 2 | #ifndef NEW_PLC |
| 3 | #define NEW_PLC |
| 4 | #endif |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 5 | |
Jean-Marc Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 6 | #ifdef FIXED_POINT |
Jean-Marc Valin | bd82ca8 | 2010-06-17 00:11:54 -0400 | [diff] [blame] | 7 | static celt_word32 frac_div32(celt_word32 a, celt_word32 b) |
| 8 | { |
Jean-Marc Valin | 3a9699e | 2010-06-17 00:35:26 -0400 | [diff] [blame^] | 9 | 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 Valin | bd82ca8 | 2010-06-17 00:11:54 -0400 | [diff] [blame] | 17 | 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 Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 22 | #else |
Jean-Marc Valin | bd82ca8 | 2010-06-17 00:11:54 -0400 | [diff] [blame] | 23 | #define frac_div32(a,b) ((float)(a)/(b)) |
Jean-Marc Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 24 | #endif |
| 25 | |
| 26 | |
Jean-Marc Valin | d69c1cb | 2009-12-28 00:34:29 -0500 | [diff] [blame] | 27 | float _celt_lpc( |
Jean-Marc Valin | 74128be | 2010-01-01 09:33:17 -0500 | [diff] [blame] | 28 | celt_word16 *_lpc, /* out: [0...p-1] LPC coefficients */ |
Jean-Marc Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 29 | const celt_word32 *ac, /* in: [0...p] autocorrelation values */ |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 30 | int p |
| 31 | ) |
| 32 | { |
| 33 | int i, j; |
Jean-Marc Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 34 | celt_word32 r; |
| 35 | celt_word32 error = ac[0]; |
Jean-Marc Valin | 74128be | 2010-01-01 09:33:17 -0500 | [diff] [blame] | 36 | #ifdef FIXED_POINT |
Jean-Marc Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 37 | celt_word32 lpc[LPC_ORDER]; |
Jean-Marc Valin | 74128be | 2010-01-01 09:33:17 -0500 | [diff] [blame] | 38 | #else |
| 39 | float *lpc = _lpc; |
| 40 | #endif |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 41 | |
Jean-Marc Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 42 | for (i = 0; i < p; i++) |
| 43 | lpc[i] = 0; |
| 44 | if (ac[0] != 0) |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 45 | { |
Jean-Marc Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 46 | 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 Valin | bd82ca8 | 2010-06-17 00:11:54 -0400 | [diff] [blame] | 52 | r = -frac_div32(SHL32(rr,3), error); |
Jean-Marc Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 53 | /* 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 Valin | bd82ca8 | 2010-06-17 00:11:54 -0400 | [diff] [blame] | 65 | /* Bail out once we get 30 dB gain */ |
| 66 | #ifdef FIXED_POINT |
| 67 | if (error<SHR32(ac[0],10)) |
Jean-Marc Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 68 | break; |
Jean-Marc Valin | bd82ca8 | 2010-06-17 00:11:54 -0400 | [diff] [blame] | 69 | #else |
| 70 | if (error<.001*ac[0]) |
| 71 | break; |
| 72 | #endif |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 73 | } |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 74 | } |
Jean-Marc Valin | 303b3b6 | 2009-12-30 22:40:24 -0500 | [diff] [blame] | 75 | #ifdef FIXED_POINT |
| 76 | for (i=0;i<p;i++) |
Jean-Marc Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 77 | _lpc[i] = ROUND16(lpc[i],16); |
Jean-Marc Valin | 303b3b6 | 2009-12-30 22:40:24 -0500 | [diff] [blame] | 78 | #endif |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 79 | return error; |
| 80 | } |
| 81 | |
Jean-Marc Valin | 303b3b6 | 2009-12-30 22:40:24 -0500 | [diff] [blame] | 82 | void fir(const celt_word16 *x, |
Jean-Marc Valin | 74128be | 2010-01-01 09:33:17 -0500 | [diff] [blame] | 83 | const celt_word16 *num, |
Jean-Marc Valin | 303b3b6 | 2009-12-30 22:40:24 -0500 | [diff] [blame] | 84 | celt_word16 *y, |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 85 | int N, |
| 86 | int ord, |
Jean-Marc Valin | 74128be | 2010-01-01 09:33:17 -0500 | [diff] [blame] | 87 | celt_word16 *mem) |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 88 | { |
| 89 | int i,j; |
| 90 | |
| 91 | for (i=0;i<N;i++) |
| 92 | { |
Jean-Marc Valin | 74128be | 2010-01-01 09:33:17 -0500 | [diff] [blame] | 93 | celt_word32 sum = SHL32(EXTEND32(x[i]), SIG_SHIFT); |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 94 | for (j=0;j<ord;j++) |
| 95 | { |
Jean-Marc Valin | 74128be | 2010-01-01 09:33:17 -0500 | [diff] [blame] | 96 | sum += MULT16_16(num[j],mem[j]); |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 97 | } |
| 98 | for (j=ord-1;j>=1;j--) |
| 99 | { |
| 100 | mem[j]=mem[j-1]; |
| 101 | } |
| 102 | mem[0] = x[i]; |
Jean-Marc Valin | 74128be | 2010-01-01 09:33:17 -0500 | [diff] [blame] | 103 | y[i] = ROUND16(sum, SIG_SHIFT); |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 104 | } |
| 105 | } |
| 106 | |
Jean-Marc Valin | d69c1cb | 2009-12-28 00:34:29 -0500 | [diff] [blame] | 107 | void iir(const celt_word32 *x, |
Jean-Marc Valin | 74128be | 2010-01-01 09:33:17 -0500 | [diff] [blame] | 108 | const celt_word16 *den, |
Jean-Marc Valin | d69c1cb | 2009-12-28 00:34:29 -0500 | [diff] [blame] | 109 | celt_word32 *y, |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 110 | int N, |
| 111 | int ord, |
Jean-Marc Valin | 74128be | 2010-01-01 09:33:17 -0500 | [diff] [blame] | 112 | celt_word16 *mem) |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 113 | { |
| 114 | int i,j; |
| 115 | for (i=0;i<N;i++) |
| 116 | { |
Jean-Marc Valin | 74128be | 2010-01-01 09:33:17 -0500 | [diff] [blame] | 117 | celt_word32 sum = x[i]; |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 118 | for (j=0;j<ord;j++) |
| 119 | { |
Jean-Marc Valin | 74128be | 2010-01-01 09:33:17 -0500 | [diff] [blame] | 120 | sum -= MULT16_16(den[j],mem[j]); |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 121 | } |
| 122 | for (j=ord-1;j>=1;j--) |
| 123 | { |
| 124 | mem[j]=mem[j-1]; |
| 125 | } |
Jean-Marc Valin | 74128be | 2010-01-01 09:33:17 -0500 | [diff] [blame] | 126 | mem[0] = ROUND16(sum,SIG_SHIFT); |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 127 | y[i] = sum; |
| 128 | } |
| 129 | } |
| 130 | |
| 131 | void _celt_autocorr( |
Jean-Marc Valin | 303b3b6 | 2009-12-30 22:40:24 -0500 | [diff] [blame] | 132 | const celt_word16 *x, /* in: [0...n-1] samples x */ |
Jean-Marc Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 133 | celt_word32 *ac, /* out: [0...lag-1] ac values */ |
Jean-Marc Valin | d69c1cb | 2009-12-28 00:34:29 -0500 | [diff] [blame] | 134 | const celt_word16 *window, |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 135 | int overlap, |
| 136 | int lag, |
| 137 | int n |
| 138 | ) |
| 139 | { |
| 140 | float d; |
| 141 | int i; |
Jean-Marc Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 142 | float scale=1; |
Thorvald Natvig | b8002a0 | 2009-12-11 13:19:09 +0100 | [diff] [blame] | 143 | VARDECL(float, xx); |
| 144 | SAVE_STACK; |
| 145 | ALLOC(xx, n, float); |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 146 | for (i=0;i<n;i++) |
| 147 | xx[i] = x[i]; |
| 148 | for (i=0;i<overlap;i++) |
| 149 | { |
Jean-Marc Valin | d69c1cb | 2009-12-28 00:34:29 -0500 | [diff] [blame] | 150 | xx[i] *= (1./Q15ONE)*window[i]; |
| 151 | xx[n-i-1] *= (1./Q15ONE)*window[i]; |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 152 | } |
Jean-Marc Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 153 | #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 Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 162 | while (lag>=0) |
| 163 | { |
| 164 | for (i = lag, d = 0; i < n; i++) |
| 165 | d += x[i] * x[i-lag]; |
Jean-Marc Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 166 | ac[lag] = d*scale; |
| 167 | /*printf ("%f ", ac[lag]);*/ |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 168 | lag--; |
| 169 | } |
Jean-Marc Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 170 | /*printf ("\n");*/ |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 171 | ac[0] += 10; |
Jean-Marc Valin | 456eab2 | 2010-06-16 22:38:57 -0400 | [diff] [blame] | 172 | |
Thorvald Natvig | b8002a0 | 2009-12-11 13:19:09 +0100 | [diff] [blame] | 173 | RESTORE_STACK; |
Jean-Marc Valin | e14fe90 | 2009-12-11 00:07:31 -0500 | [diff] [blame] | 174 | } |