1 /* Copyright (c) 2009-2010 Xiph.Org Foundation 2 Written by Jean-Marc Valin */ 3 /* 4 Redistribution and use in source and binary forms, with or without 5 modification, are permitted provided that the following conditions 6 are met: 7 8 - Redistributions of source code must retain the above copyright 9 notice, this list of conditions and the following disclaimer. 10 11 - Redistributions in binary form must reproduce the above copyright 12 notice, this list of conditions and the following disclaimer in the 13 documentation and/or other materials provided with the distribution. 14 15 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 16 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 17 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 18 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER 19 OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 20 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 21 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 22 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 23 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 24 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 25 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 26 */ 27 28 #ifdef HAVE_CONFIG_H 29 #include "config.h" 30 #endif 31 32 #include "celt_lpc.h" 33 #include "stack_alloc.h" 34 #include "mathops.h" 35 36 void _celt_lpc( 37 opus_val16 *_lpc, /* out: [0...p-1] LPC coefficients */ 38 const opus_val32 *ac, /* in: [0...p] autocorrelation values */ 39 int p 40 ) 41 { 42 int i, j; 43 opus_val32 r; 44 opus_val32 error = ac[0]; 45 #ifdef FIXED_POINT 46 opus_val32 lpc[LPC_ORDER]; 47 #else 48 float *lpc = _lpc; 49 #endif 50 51 for (i = 0; i < p; i++) 52 lpc[i] = 0; 53 if (ac[0] != 0) 54 { 55 for (i = 0; i < p; i++) { 56 /* Sum up this iteration's reflection coefficient */ 57 opus_val32 rr = 0; 58 for (j = 0; j < i; j++) 59 rr += MULT32_32_Q31(lpc[j],ac[i - j]); 60 rr += SHR32(ac[i + 1],3); 61 r = -frac_div32(SHL32(rr,3), error); 62 /* Update LPC coefficients and total error */ 63 lpc[i] = SHR32(r,3); 64 for (j = 0; j < (i+1)>>1; j++) 65 { 66 opus_val32 tmp1, tmp2; 67 tmp1 = lpc[j]; 68 tmp2 = lpc[i-1-j]; 69 lpc[j] = tmp1 + MULT32_32_Q31(r,tmp2); 70 lpc[i-1-j] = tmp2 + MULT32_32_Q31(r,tmp1); 71 } 72 73 error = error - MULT32_32_Q31(MULT32_32_Q31(r,r),error); 74 /* Bail out once we get 30 dB gain */ 75 #ifdef FIXED_POINT 76 if (error<SHR32(ac[0],10)) 77 break; 78 #else 79 if (error<.001f*ac[0]) 80 break; 81 #endif 82 } 83 } 84 #ifdef FIXED_POINT 85 for (i=0;i<p;i++) 86 _lpc[i] = ROUND16(lpc[i],16); 87 #endif 88 } 89 90 void celt_fir(const opus_val16 *x, 91 const opus_val16 *num, 92 opus_val16 *y, 93 int N, 94 int ord, 95 opus_val16 *mem) 96 { 97 int i,j; 98 99 for (i=0;i<N;i++) 100 { 101 opus_val32 sum = SHL32(EXTEND32(x[i]), SIG_SHIFT); 102 for (j=0;j<ord;j++) 103 { 104 sum += MULT16_16(num[j],mem[j]); 105 } 106 for (j=ord-1;j>=1;j--) 107 { 108 mem[j]=mem[j-1]; 109 } 110 mem[0] = x[i]; 111 y[i] = ROUND16(sum, SIG_SHIFT); 112 } 113 } 114 115 void celt_iir(const opus_val32 *x, 116 const opus_val16 *den, 117 opus_val32 *y, 118 int N, 119 int ord, 120 opus_val16 *mem) 121 { 122 int i,j; 123 for (i=0;i<N;i++) 124 { 125 opus_val32 sum = x[i]; 126 for (j=0;j<ord;j++) 127 { 128 sum -= MULT16_16(den[j],mem[j]); 129 } 130 for (j=ord-1;j>=1;j--) 131 { 132 mem[j]=mem[j-1]; 133 } 134 mem[0] = ROUND16(sum,SIG_SHIFT); 135 y[i] = sum; 136 } 137 } 138 139 void _celt_autocorr( 140 const opus_val16 *x, /* in: [0...n-1] samples x */ 141 opus_val32 *ac, /* out: [0...lag-1] ac values */ 142 const opus_val16 *window, 143 int overlap, 144 int lag, 145 int n 146 ) 147 { 148 opus_val32 d; 149 int i; 150 VARDECL(opus_val16, xx); 151 SAVE_STACK; 152 ALLOC(xx, n, opus_val16); 153 celt_assert(n>0); 154 celt_assert(overlap>=0); 155 for (i=0;i<n;i++) 156 xx[i] = x[i]; 157 for (i=0;i<overlap;i++) 158 { 159 xx[i] = MULT16_16_Q15(x[i],window[i]); 160 xx[n-i-1] = MULT16_16_Q15(x[n-i-1],window[i]); 161 } 162 #ifdef FIXED_POINT 163 { 164 opus_val32 ac0=0; 165 int shift; 166 for(i=0;i<n;i++) 167 ac0 += SHR32(MULT16_16(xx[i],xx[i]),9); 168 ac0 += 1+n; 169 170 shift = celt_ilog2(ac0)-30+10; 171 shift = (shift+1)/2; 172 for(i=0;i<n;i++) 173 xx[i] = VSHR32(xx[i], shift); 174 } 175 #endif 176 while (lag>=0) 177 { 178 for (i = lag, d = 0; i < n; i++) 179 d += xx[i] * xx[i-lag]; 180 ac[lag] = d; 181 /*printf ("%f ", ac[lag]);*/ 182 lag--; 183 } 184 /*printf ("\n");*/ 185 ac[0] += 10; 186 187 RESTORE_STACK; 188 } 189
