lpc.c 5.27 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
/*
  Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann,
  Technische Universitaet Berlin

  Any use of this software is permitted provided that this notice is not
  removed and that neither the authors nor the Technische Universitaet Berlin
  are deemed to have made any representations as to the suitability of this
  software for any purpose nor are held responsible for any defects of
  this software.  THERE IS ABSOLUTELY NO WARRANTY FOR THIS SOFTWARE.

  As a matter of courtesy, the authors request to be informed about uses
  this software has found, about bugs in this software, and about any
  improvements that may be of general interest.

  Berlin, 28.11.1994
  Jutta Degener
  Carsten Bormann
jmvalin's avatar
jmvalin committed
18 19


20
   Code modified by Jean-Marc Valin
jmvalin's avatar
jmvalin committed
21 22

   Speex License:
jmvalin's avatar
jmvalin committed
23

jm's avatar
jm committed
24 25 26
   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions
   are met:
Tristan Matthews's avatar
Tristan Matthews committed
27

jm's avatar
jm committed
28 29
   - Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.
Tristan Matthews's avatar
Tristan Matthews committed
30

jm's avatar
jm committed
31 32 33
   - Redistributions in binary form must reproduce the above copyright
   notice, this list of conditions and the following disclaimer in the
   documentation and/or other materials provided with the distribution.
Tristan Matthews's avatar
Tristan Matthews committed
34

jm's avatar
jm committed
35 36 37
   - Neither the name of the Xiph.org Foundation nor the names of its
   contributors may be used to endorse or promote products derived from
   this software without specific prior written permission.
Tristan Matthews's avatar
Tristan Matthews committed
38

jm's avatar
jm committed
39 40 41 42 43 44 45 46 47 48 49
   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
jmvalin's avatar
jmvalin committed
50 51
*/

52 53 54 55
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

56 57
#ifndef DISABLE_ENCODER

58
#include "lpc.h"
59

60 61 62
#ifdef BFIN_ASM
#include "lpc_bfin.h"
#endif
63

64
/* LPC analysis
65 66 67 68 69 70 71 72 73
 *
 * The next two functions calculate linear prediction coefficients
 * and/or the related reflection coefficients from the first P_MAX+1
 * values of the autocorrelation function.
 */

/* Invented by N. Levinson in 1947, modified by J. Durbin in 1959.
 */

jm's avatar
jm committed
74
/* returns minimum mean square error    */
jm's avatar
jm committed
75
spx_word32_t _spx_lpc(
jm's avatar
jm committed
76
spx_coef_t       *lpc, /* out: [0...p-1] LPC coefficients      */
jm's avatar
jm committed
77
const spx_word16_t *ac,  /* in:  [0...p] autocorrelation values  */
jm's avatar
jm committed
78 79 80
int          p
)
{
Tristan Matthews's avatar
Tristan Matthews committed
81
   int i, j;
jm's avatar
jm committed
82
   spx_word16_t r;
jm's avatar
jm committed
83
   spx_word16_t error = ac[0];
jm's avatar
jm committed
84 85 86 87

   for (i = 0; i < p; i++) {

      /* Sum up this iteration's reflection coefficient */
jm's avatar
jm committed
88
      spx_word32_t rr = NEG32(SHL32(EXTEND32(ac[i + 1]),13));
Tristan Matthews's avatar
Tristan Matthews committed
89
      for (j = 0; j < i; j++)
90
         rr = SUB32(rr,MULT16_16(lpc[j],ac[i - j]));
91
#ifdef FIXED_POINT
92
      r = DIV32_16(rr+PSHR32(error,1),ADD16(error,8));
93 94 95
#else
      r = rr/(error+.003*ac[0]);
#endif
jm's avatar
jm committed
96
      /*  Update LPC coefficients and total error */
jm's avatar
jm committed
97
      lpc[i] = r;
Tristan Matthews's avatar
Tristan Matthews committed
98
      for (j = 0; j < (i+1)>>1; j++)
jm's avatar
jm committed
99
      {
100 101 102 103 104 105
         spx_word16_t tmp1, tmp2;
         /* It could be that j == i-1-j, in which case, we're updating the same value twice, which is OK */
         tmp1 = lpc[j];
         tmp2 = lpc[i-1-j];
         lpc[j]     = MAC16_16_P13(tmp1,r,tmp2);
         lpc[i-1-j] = MAC16_16_P13(tmp2,r,tmp1);
jm's avatar
jm committed
106 107
      }

108
      error = SUB16(error,MULT16_16_Q13(r,MULT16_16_Q13(error,r)));
jm's avatar
jm committed
109 110 111 112 113
   }
   return error;
}


114 115
#ifdef FIXED_POINT

jm's avatar
jm committed
116 117 118 119 120 121
/* Compute the autocorrelation
 *                      ,--,
 *              ac(i) = >  x(n) * x(n-i)  for all n
 *                      `--'
 * for lags between 0 and lag-1, and x == 0 outside 0...n-1
 */
122

123
#ifndef OVERRIDE_SPEEX_AUTOCORR
jm's avatar
jm committed
124
void _spx_autocorr(
125
const spx_word16_t *x,   /*  in: [0...n-1] samples x   */
jm's avatar
jm committed
126
spx_word16_t       *ac,  /* out: [0...lag-1] ac values */
Tristan Matthews's avatar
Tristan Matthews committed
127
int          lag,
jm's avatar
jm committed
128 129 130
int          n
)
{
jm's avatar
jm committed
131
   spx_word32_t d;
jm's avatar
jm committed
132
   int i, j;
jm's avatar
jm committed
133 134
   spx_word32_t ac0=1;
   int shift, ac_shift;
Tristan Matthews's avatar
Tristan Matthews committed
135

jm's avatar
jm committed
136
   for (j=0;j<n;j++)
jm's avatar
jm committed
137
      ac0 = ADD32(ac0,SHR32(MULT16_16(x[j],x[j]),8));
138
   ac0 = ADD32(ac0,n);
jm's avatar
jm committed
139 140 141 142 143 144 145 146 147 148 149 150
   shift = 8;
   while (shift && ac0<0x40000000)
   {
      shift--;
      ac0 <<= 1;
   }
   ac_shift = 18;
   while (ac_shift && ac0<0x40000000)
   {
      ac_shift--;
      ac0 <<= 1;
   }
Tristan Matthews's avatar
Tristan Matthews committed
151 152


jm's avatar
jm committed
153 154 155 156 157
   for (i=0;i<lag;i++)
   {
      d=0;
      for (j=i;j<n;j++)
      {
jm's avatar
jm committed
158
         d = ADD32(d,SHR32(MULT16_16(x[j],x[j-i]), shift));
jm's avatar
jm committed
159
      }
Tristan Matthews's avatar
Tristan Matthews committed
160

jm's avatar
jm committed
161
      ac[i] = SHR32(d, ac_shift);
jm's avatar
jm committed
162 163
   }
}
164
#endif
jm's avatar
jm committed
165 166 167 168 169 170


#else



171 172 173 174 175 176
/* Compute the autocorrelation
 *                      ,--,
 *              ac(i) = >  x(n) * x(n-i)  for all n
 *                      `--'
 * for lags between 0 and lag-1, and x == 0 outside 0...n-1
 */
177
void _spx_autocorr(
178
const spx_word16_t *x,   /*  in: [0...n-1] samples x   */
179
float       *ac,  /* out: [0...lag-1] ac values */
Tristan Matthews's avatar
Tristan Matthews committed
180
int          lag,
181 182
int          n
)
183
{
184 185
   float d;
   int i;
Tristan Matthews's avatar
Tristan Matthews committed
186
   while (lag--)
187
   {
Tristan Matthews's avatar
Tristan Matthews committed
188
      for (i = lag, d = 0; i < n; i++)
189
         d += x[i] * x[i-lag];
190 191
      ac[lag] = d;
   }
jm's avatar
jm committed
192
   ac[0] += 10;
193
}
jm's avatar
jm committed
194 195 196 197

#endif


198
#endif /* DISABLE_ENCODER */