1 /*---------------------------------------------------------------------------*\
2 Original copyright
3 FILE........: lsp.c
4 AUTHOR......: David Rowe
5 DATE CREATED: 24/2/93
6
7 Heavily modified by Jean-Marc Valin (c) 2002-2006 (fixed-point,
8 optimizations, additional functions, ...)
9
10 This file contains functions for converting Linear Prediction
11 Coefficients (LPC) to Line Spectral Pair (LSP) and back. Note that the
12 LSP coefficients are not in radians format but in the x domain of the
13 unit circle.
14
15 Speex License:
16
17 Redistribution and use in source and binary forms, with or without
18 modification, are permitted provided that the following conditions
19 are met:
20
21 - Redistributions of source code must retain the above copyright
22 notice, this list of conditions and the following disclaimer.
23
24 - Redistributions in binary form must reproduce the above copyright
25 notice, this list of conditions and the following disclaimer in the
26 documentation and/or other materials provided with the distribution.
27
28 - Neither the name of the Xiph.org Foundation nor the names of its
29 contributors may be used to endorse or promote products derived from
30 this software without specific prior written permission.
31
32 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
33 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
34 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
35 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR
36 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
37 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
38 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
39 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
40 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
41 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
42 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
43 */
44
45 /*---------------------------------------------------------------------------*\
46
47 Introduction to Line Spectrum Pairs (LSPs)
48 ------------------------------------------
49
50 LSPs are used to encode the LPC filter coefficients {ak} for
51 transmission over the channel. LSPs have several properties (like
52 less sensitivity to quantisation noise) that make them superior to
53 direct quantisation of {ak}.
54
55 A(z) is a polynomial of order lpcrdr with {ak} as the coefficients.
56
57 A(z) is transformed to P(z) and Q(z) (using a substitution and some
58 algebra), to obtain something like:
59
60 A(z) = 0.5[P(z)(z+z^-1) + Q(z)(z-z^-1)] (1)
61
62 As you can imagine A(z) has complex zeros all over the z-plane. P(z)
63 and Q(z) have the very neat property of only having zeros _on_ the
64 unit circle. So to find them we take a test point z=exp(jw) and
65 evaluate P (exp(jw)) and Q(exp(jw)) using a grid of points between 0
66 and pi.
67
68 The zeros (roots) of P(z) also happen to alternate, which is why we
69 swap coefficients as we find roots. So the process of finding the
70 LSP frequencies is basically finding the roots of 5th order
71 polynomials.
72
73 The root so P(z) and Q(z) occur in symmetrical pairs at +/-w, hence
74 the name Line Spectrum Pairs (LSPs).
75
76 To convert back to ak we just evaluate (1), "clocking" an impulse
77 thru it lpcrdr times gives us the impulse response of A(z) which is
78 {ak}.
79
80 \*---------------------------------------------------------------------------*/
81
82 #ifdef HAVE_CONFIG_H
83 #include "config.h"
84 #endif
85
86 #include <math.h>
87 #include "lsp.h"
88 #include "stack_alloc.h"
89 #include "math_approx.h"
90
91 #ifndef M_PI
92 #define M_PI 3.14159265358979323846 /* pi */
93 #endif
94
95 #ifndef NULL
96 #define NULL 0
97 #endif
98
99 #ifdef FIXED_POINT
100
101 #define FREQ_SCALE 16384
102
103 /*#define ANGLE2X(a) (32768*cos(((a)/8192.)))*/
104 #define ANGLE2X(a) (SHL16(spx_cos(a),2))
105
106 /*#define X2ANGLE(x) (acos(.00006103515625*(x))*LSP_SCALING)*/
107 #define X2ANGLE(x) (spx_acos(x))
108
109 #ifdef BFIN_ASM
110 #include "lsp_bfin.h"
111 #endif
112
113 #else
114
115 /*#define C1 0.99940307
116 #define C2 -0.49558072
117 #define C3 0.03679168*/
118
119 #define FREQ_SCALE 1.
120 #define ANGLE2X(a) (spx_cos(a))
121 #define X2ANGLE(x) (acos(x))
122
123 #endif
124
125
126 /*---------------------------------------------------------------------------*\
127
128 FUNCTION....: cheb_poly_eva()
129
130 AUTHOR......: David Rowe
131 DATE CREATED: 24/2/93
132
133 This function evaluates a series of Chebyshev polynomials
134
135 \*---------------------------------------------------------------------------*/
136
137 #ifdef FIXED_POINT
138
139 #ifndef OVERRIDE_CHEB_POLY_EVA
cheb_poly_eva(spx_word16_t * coef,spx_word16_t x,int m,char * stack)140 static inline spx_word32_t cheb_poly_eva(
141 spx_word16_t *coef, /* P or Q coefs in Q13 format */
142 spx_word16_t x, /* cos of freq (-1.0 to 1.0) in Q14 format */
143 int m, /* LPC order/2 */
144 char *stack
145 )
146 {
147 int i;
148 spx_word16_t b0, b1;
149 spx_word32_t sum;
150
151 /*Prevents overflows*/
152 if (x>16383)
153 x = 16383;
154 if (x<-16383)
155 x = -16383;
156
157 /* Initialise values */
158 b1=16384;
159 b0=x;
160
161 /* Evaluate Chebyshev series formulation usin g iterative approach */
162 sum = ADD32(EXTEND32(coef[m]), EXTEND32(MULT16_16_P14(coef[m-1],x)));
163 for(i=2;i<=m;i++)
164 {
165 spx_word16_t tmp=b0;
166 b0 = SUB16(MULT16_16_Q13(x,b0), b1);
167 b1 = tmp;
168 sum = ADD32(sum, EXTEND32(MULT16_16_P14(coef[m-i],b0)));
169 }
170
171 return sum;
172 }
173 #endif
174
175 #else
176
cheb_poly_eva(spx_word32_t * coef,spx_word16_t x,int m,char * stack)177 static float cheb_poly_eva(spx_word32_t *coef, spx_word16_t x, int m, char *stack)
178 {
179 int k;
180 float b0, b1, tmp;
181
182 /* Initial conditions */
183 b0=0; /* b_(m+1) */
184 b1=0; /* b_(m+2) */
185
186 x*=2;
187
188 /* Calculate the b_(k) */
189 for(k=m;k>0;k--)
190 {
191 tmp=b0; /* tmp holds the previous value of b0 */
192 b0=x*b0-b1+coef[m-k]; /* b0 holds its new value based on b0 and b1 */
193 b1=tmp; /* b1 holds the previous value of b0 */
194 }
195
196 return(-b1+.5*x*b0+coef[m]);
197 }
198 #endif
199
200 /*---------------------------------------------------------------------------*\
201
202 FUNCTION....: lpc_to_lsp()
203
204 AUTHOR......: David Rowe
205 DATE CREATED: 24/2/93
206
207 This function converts LPC coefficients to LSP
208 coefficients.
209
210 \*---------------------------------------------------------------------------*/
211
212 #ifdef FIXED_POINT
213 #define SIGN_CHANGE(a,b) (((a)&0x70000000)^((b)&0x70000000)||(b==0))
214 #else
215 #define SIGN_CHANGE(a,b) (((a)*(b))<0.0)
216 #endif
217
218
lpc_to_lsp(spx_coef_t * a,int lpcrdr,spx_lsp_t * freq,int nb,spx_word16_t delta,char * stack)219 int lpc_to_lsp (spx_coef_t *a,int lpcrdr,spx_lsp_t *freq,int nb,spx_word16_t delta, char *stack)
220 /* float *a lpc coefficients */
221 /* int lpcrdr order of LPC coefficients (10) */
222 /* float *freq LSP frequencies in the x domain */
223 /* int nb number of sub-intervals (4) */
224 /* float delta grid spacing interval (0.02) */
225
226
227 {
228 spx_word16_t temp_xr,xl,xr,xm=0;
229 spx_word32_t psuml,psumr,psumm,temp_psumr/*,temp_qsumr*/;
230 int i,j,m,flag,k;
231 VARDECL(spx_word32_t *Q); /* ptrs for memory allocation */
232 VARDECL(spx_word32_t *P);
233 VARDECL(spx_word16_t *Q16); /* ptrs for memory allocation */
234 VARDECL(spx_word16_t *P16);
235 spx_word32_t *px; /* ptrs of respective P'(z) & Q'(z) */
236 spx_word32_t *qx;
237 spx_word32_t *p;
238 spx_word32_t *q;
239 spx_word16_t *pt; /* ptr used for cheb_poly_eval()
240 whether P' or Q' */
241 int roots=0; /* DR 8/2/94: number of roots found */
242 flag = 1; /* program is searching for a root when,
243 1 else has found one */
244 m = lpcrdr/2; /* order of P'(z) & Q'(z) polynomials */
245
246 /* Allocate memory space for polynomials */
247 ALLOC(Q, (m+1), spx_word32_t);
248 ALLOC(P, (m+1), spx_word32_t);
249
250 /* determine P'(z)'s and Q'(z)'s coefficients where
251 P'(z) = P(z)/(1 + z^(-1)) and Q'(z) = Q(z)/(1-z^(-1)) */
252
253 px = P; /* initialise ptrs */
254 qx = Q;
255 p = px;
256 q = qx;
257
258 #ifdef FIXED_POINT
259 *px++ = LPC_SCALING;
260 *qx++ = LPC_SCALING;
261 for(i=0;i<m;i++){
262 *px++ = SUB32(ADD32(EXTEND32(a[i]),EXTEND32(a[lpcrdr-i-1])), *p++);
263 *qx++ = ADD32(SUB32(EXTEND32(a[i]),EXTEND32(a[lpcrdr-i-1])), *q++);
264 }
265 px = P;
266 qx = Q;
267 for(i=0;i<m;i++)
268 {
269 /*if (fabs(*px)>=32768)
270 speex_warning_int("px", *px);
271 if (fabs(*qx)>=32768)
272 speex_warning_int("qx", *qx);*/
273 *px = PSHR32(*px,2);
274 *qx = PSHR32(*qx,2);
275 px++;
276 qx++;
277 }
278 /* The reason for this lies in the way cheb_poly_eva() is implemented for fixed-point */
279 P[m] = PSHR32(P[m],3);
280 Q[m] = PSHR32(Q[m],3);
281 #else
282 *px++ = LPC_SCALING;
283 *qx++ = LPC_SCALING;
284 for(i=0;i<m;i++){
285 *px++ = (a[i]+a[lpcrdr-1-i]) - *p++;
286 *qx++ = (a[i]-a[lpcrdr-1-i]) + *q++;
287 }
288 px = P;
289 qx = Q;
290 for(i=0;i<m;i++){
291 *px = 2**px;
292 *qx = 2**qx;
293 px++;
294 qx++;
295 }
296 #endif
297
298 px = P; /* re-initialise ptrs */
299 qx = Q;
300
301 /* now that we have computed P and Q convert to 16 bits to
302 speed up cheb_poly_eval */
303
304 ALLOC(P16, m+1, spx_word16_t);
305 ALLOC(Q16, m+1, spx_word16_t);
306
307 for (i=0;i<m+1;i++)
308 {
309 P16[i] = P[i];
310 Q16[i] = Q[i];
311 }
312
313 /* Search for a zero in P'(z) polynomial first and then alternate to Q'(z).
314 Keep alternating between the two polynomials as each zero is found */
315
316 xr = 0; /* initialise xr to zero */
317 xl = FREQ_SCALE; /* start at point xl = 1 */
318
319 for(j=0;j<lpcrdr;j++){
320 if(j&1) /* determines whether P' or Q' is eval. */
321 pt = Q16;
322 else
323 pt = P16;
324
325 psuml = cheb_poly_eva(pt,xl,m,stack); /* evals poly. at xl */
326 flag = 1;
327 while(flag && (xr >= -FREQ_SCALE)){
328 spx_word16_t dd;
329 /* Modified by JMV to provide smaller steps around x=+-1 */
330 #ifdef FIXED_POINT
331 dd = MULT16_16_Q15(delta,SUB16(FREQ_SCALE, MULT16_16_Q14(MULT16_16_Q14(xl,xl),14000)));
332 if (psuml<512 && psuml>-512)
333 dd = PSHR16(dd,1);
334 #else
335 dd=delta*(1-.9*xl*xl);
336 if (fabs(psuml)<.2)
337 dd *= .5;
338 #endif
339 xr = SUB16(xl, dd); /* interval spacing */
340 psumr = cheb_poly_eva(pt,xr,m,stack);/* poly(xl-delta_x) */
341 temp_psumr = psumr;
342 temp_xr = xr;
343
344 /* if no sign change increment xr and re-evaluate poly(xr). Repeat til
345 sign change.
346 if a sign change has occurred the interval is bisected and then
347 checked again for a sign change which determines in which
348 interval the zero lies in.
349 If there is no sign change between poly(xm) and poly(xl) set interval
350 between xm and xr else set interval between xl and xr and repeat till
351 root is located within the specified limits */
352
353 if(SIGN_CHANGE(psumr,psuml))
354 {
355 roots++;
356
357 psumm=psuml;
358 for(k=0;k<=nb;k++){
359 #ifdef FIXED_POINT
360 xm = ADD16(PSHR16(xl,1),PSHR16(xr,1)); /* bisect the interval */
361 #else
362 xm = .5*(xl+xr); /* bisect the interval */
363 #endif
364 psumm=cheb_poly_eva(pt,xm,m,stack);
365 /*if(psumm*psuml>0.)*/
366 if(!SIGN_CHANGE(psumm,psuml))
367 {
368 psuml=psumm;
369 xl=xm;
370 } else {
371 psumr=psumm;
372 xr=xm;
373 }
374 }
375
376 /* once zero is found, reset initial interval to xr */
377 freq[j] = X2ANGLE(xm);
378 xl = xm;
379 flag = 0; /* reset flag for next search */
380 }
381 else{
382 psuml=temp_psumr;
383 xl=temp_xr;
384 }
385 }
386 }
387 return(roots);
388 }
389
390 /*---------------------------------------------------------------------------*\
391
392 FUNCTION....: lsp_to_lpc()
393
394 AUTHOR......: David Rowe
395 DATE CREATED: 24/2/93
396
397 Converts LSP coefficients to LPC coefficients.
398
399 \*---------------------------------------------------------------------------*/
400
401 #ifdef FIXED_POINT
402
lsp_to_lpc(spx_lsp_t * freq,spx_coef_t * ak,int lpcrdr,char * stack)403 void lsp_to_lpc(spx_lsp_t *freq,spx_coef_t *ak,int lpcrdr, char *stack)
404 /* float *freq array of LSP frequencies in the x domain */
405 /* float *ak array of LPC coefficients */
406 /* int lpcrdr order of LPC coefficients */
407 {
408 int i,j;
409 spx_word32_t xout1,xout2,xin;
410 spx_word32_t mult, a;
411 VARDECL(spx_word16_t *freqn);
412 VARDECL(spx_word32_t **xp);
413 VARDECL(spx_word32_t *xpmem);
414 VARDECL(spx_word32_t **xq);
415 VARDECL(spx_word32_t *xqmem);
416 int m = lpcrdr>>1;
417
418 /*
419
420 Reconstruct P(z) and Q(z) by cascading second order polynomials
421 in form 1 - 2cos(w)z(-1) + z(-2), where w is the LSP frequency.
422 In the time domain this is:
423
424 y(n) = x(n) - 2cos(w)x(n-1) + x(n-2)
425
426 This is what the ALLOCS below are trying to do:
427
428 int xp[m+1][lpcrdr+1+2]; // P matrix in QIMP
429 int xq[m+1][lpcrdr+1+2]; // Q matrix in QIMP
430
431 These matrices store the output of each stage on each row. The
432 final (m-th) row has the output of the final (m-th) cascaded
433 2nd order filter. The first row is the impulse input to the
434 system (not written as it is known).
435
436 The version below takes advantage of the fact that a lot of the
437 outputs are zero or known, for example if we put an inpulse
438 into the first section the "clock" it 10 times only the first 3
439 outputs samples are non-zero (it's an FIR filter).
440 */
441
442 ALLOC(xp, (m+1), spx_word32_t*);
443 ALLOC(xpmem, (m+1)*(lpcrdr+1+2), spx_word32_t);
444
445 ALLOC(xq, (m+1), spx_word32_t*);
446 ALLOC(xqmem, (m+1)*(lpcrdr+1+2), spx_word32_t);
447
448 for(i=0; i<=m; i++) {
449 xp[i] = xpmem + i*(lpcrdr+1+2);
450 xq[i] = xqmem + i*(lpcrdr+1+2);
451 }
452
453 /* work out 2cos terms in Q14 */
454
455 ALLOC(freqn, lpcrdr, spx_word16_t);
456 for (i=0;i<lpcrdr;i++)
457 freqn[i] = ANGLE2X(freq[i]);
458
459 #define QIMP 21 /* scaling for impulse */
460
461 xin = SHL32(EXTEND32(1), (QIMP-1)); /* 0.5 in QIMP format */
462
463 /* first col and last non-zero values of each row are trivial */
464
465 for(i=0;i<=m;i++) {
466 xp[i][1] = 0;
467 xp[i][2] = xin;
468 xp[i][2+2*i] = xin;
469 xq[i][1] = 0;
470 xq[i][2] = xin;
471 xq[i][2+2*i] = xin;
472 }
473
474 /* 2nd row (first output row) is trivial */
475
476 xp[1][3] = -MULT16_32_Q14(freqn[0],xp[0][2]);
477 xq[1][3] = -MULT16_32_Q14(freqn[1],xq[0][2]);
478
479 xout1 = xout2 = 0;
480
481 /* now generate remaining rows */
482
483 for(i=1;i<m;i++) {
484
485 for(j=1;j<2*(i+1)-1;j++) {
486 mult = MULT16_32_Q14(freqn[2*i],xp[i][j+1]);
487 xp[i+1][j+2] = ADD32(SUB32(xp[i][j+2], mult), xp[i][j]);
488 mult = MULT16_32_Q14(freqn[2*i+1],xq[i][j+1]);
489 xq[i+1][j+2] = ADD32(SUB32(xq[i][j+2], mult), xq[i][j]);
490 }
491
492 /* for last col xp[i][j+2] = xq[i][j+2] = 0 */
493
494 mult = MULT16_32_Q14(freqn[2*i],xp[i][j+1]);
495 xp[i+1][j+2] = SUB32(xp[i][j], mult);
496 mult = MULT16_32_Q14(freqn[2*i+1],xq[i][j+1]);
497 xq[i+1][j+2] = SUB32(xq[i][j], mult);
498 }
499
500 /* process last row to extra a{k} */
501
502 for(j=1;j<=lpcrdr;j++) {
503 int shift = QIMP-13;
504
505 /* final filter sections */
506 a = PSHR32(xp[m][j+2] + xout1 + xq[m][j+2] - xout2, shift);
507 xout1 = xp[m][j+2];
508 xout2 = xq[m][j+2];
509
510 /* hard limit ak's to +/- 32767 */
511
512 if (a < -32767) a = -32767;
513 if (a > 32767) a = 32767;
514 ak[j-1] = (short)a;
515
516 }
517
518 }
519
520 #else
521
lsp_to_lpc(spx_lsp_t * freq,spx_coef_t * ak,int lpcrdr,char * stack)522 void lsp_to_lpc(spx_lsp_t *freq,spx_coef_t *ak,int lpcrdr, char *stack)
523 /* float *freq array of LSP frequencies in the x domain */
524 /* float *ak array of LPC coefficients */
525 /* int lpcrdr order of LPC coefficients */
526
527
528 {
529 int i,j;
530 float xout1,xout2,xin1,xin2;
531 VARDECL(float *Wp);
532 float *pw,*n1,*n2,*n3,*n4=NULL;
533 VARDECL(float *x_freq);
534 int m = lpcrdr>>1;
535
536 ALLOC(Wp, 4*m+2, float);
537 pw = Wp;
538
539 /* initialise contents of array */
540
541 for(i=0;i<=4*m+1;i++){ /* set contents of buffer to 0 */
542 *pw++ = 0.0;
543 }
544
545 /* Set pointers up */
546
547 pw = Wp;
548 xin1 = 1.0;
549 xin2 = 1.0;
550
551 ALLOC(x_freq, lpcrdr, float);
552 for (i=0;i<lpcrdr;i++)
553 x_freq[i] = ANGLE2X(freq[i]);
554
555 /* reconstruct P(z) and Q(z) by cascading second order
556 polynomials in form 1 - 2xz(-1) +z(-2), where x is the
557 LSP coefficient */
558
559 for(j=0;j<=lpcrdr;j++){
560 int i2=0;
561 for(i=0;i<m;i++,i2+=2){
562 n1 = pw+(i*4);
563 n2 = n1 + 1;
564 n3 = n2 + 1;
565 n4 = n3 + 1;
566 xout1 = xin1 - 2.f*x_freq[i2] * *n1 + *n2;
567 xout2 = xin2 - 2.f*x_freq[i2+1] * *n3 + *n4;
568 *n2 = *n1;
569 *n4 = *n3;
570 *n1 = xin1;
571 *n3 = xin2;
572 xin1 = xout1;
573 xin2 = xout2;
574 }
575 xout1 = xin1 + *(n4+1);
576 xout2 = xin2 - *(n4+2);
577 if (j>0)
578 ak[j-1] = (xout1 + xout2)*0.5f;
579 *(n4+1) = xin1;
580 *(n4+2) = xin2;
581
582 xin1 = 0.0;
583 xin2 = 0.0;
584 }
585
586 }
587 #endif
588
589
590 #ifdef FIXED_POINT
591
592 /*Makes sure the LSPs are stable*/
lsp_enforce_margin(spx_lsp_t * lsp,int len,spx_word16_t margin)593 void lsp_enforce_margin(spx_lsp_t *lsp, int len, spx_word16_t margin)
594 {
595 int i;
596 spx_word16_t m = margin;
597 spx_word16_t m2 = 25736-margin;
598
599 if (lsp[0]<m)
600 lsp[0]=m;
601 if (lsp[len-1]>m2)
602 lsp[len-1]=m2;
603 for (i=1;i<len-1;i++)
604 {
605 if (lsp[i]<lsp[i-1]+m)
606 lsp[i]=lsp[i-1]+m;
607
608 if (lsp[i]>lsp[i+1]-m)
609 lsp[i]= SHR16(lsp[i],1) + SHR16(lsp[i+1]-m,1);
610 }
611 }
612
613
lsp_interpolate(spx_lsp_t * old_lsp,spx_lsp_t * new_lsp,spx_lsp_t * interp_lsp,int len,int subframe,int nb_subframes)614 void lsp_interpolate(spx_lsp_t *old_lsp, spx_lsp_t *new_lsp, spx_lsp_t *interp_lsp, int len, int subframe, int nb_subframes)
615 {
616 int i;
617 spx_word16_t tmp = DIV32_16(SHL32(EXTEND32(1 + subframe),14),nb_subframes);
618 spx_word16_t tmp2 = 16384-tmp;
619 for (i=0;i<len;i++)
620 {
621 interp_lsp[i] = MULT16_16_P14(tmp2,old_lsp[i]) + MULT16_16_P14(tmp,new_lsp[i]);
622 }
623 }
624
625 #else
626
627 /*Makes sure the LSPs are stable*/
lsp_enforce_margin(spx_lsp_t * lsp,int len,spx_word16_t margin)628 void lsp_enforce_margin(spx_lsp_t *lsp, int len, spx_word16_t margin)
629 {
630 int i;
631 if (lsp[0]<LSP_SCALING*margin)
632 lsp[0]=LSP_SCALING*margin;
633 if (lsp[len-1]>LSP_SCALING*(M_PI-margin))
634 lsp[len-1]=LSP_SCALING*(M_PI-margin);
635 for (i=1;i<len-1;i++)
636 {
637 if (lsp[i]<lsp[i-1]+LSP_SCALING*margin)
638 lsp[i]=lsp[i-1]+LSP_SCALING*margin;
639
640 if (lsp[i]>lsp[i+1]-LSP_SCALING*margin)
641 lsp[i]= .5f* (lsp[i] + lsp[i+1]-LSP_SCALING*margin);
642 }
643 }
644
645
lsp_interpolate(spx_lsp_t * old_lsp,spx_lsp_t * new_lsp,spx_lsp_t * interp_lsp,int len,int subframe,int nb_subframes)646 void lsp_interpolate(spx_lsp_t *old_lsp, spx_lsp_t *new_lsp, spx_lsp_t *interp_lsp, int len, int subframe, int nb_subframes)
647 {
648 int i;
649 float tmp = (1.0f + subframe)/nb_subframes;
650 for (i=0;i<len;i++)
651 {
652 interp_lsp[i] = (1-tmp)*old_lsp[i] + tmp*new_lsp[i];
653 }
654 }
655
656 #endif
657