1 /*---------------------------------------------------------------------------*\
4 AUTHOR......: David Rowe
7 Heavily modified by Jean-Marc Valin (c) 2002-2006 (fixed-point,
8 optimizations, additional functions, ...)
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
17 Redistribution and use in source and binary forms, with or without
18 modification, are permitted provided that the following conditions
21 - Redistributions of source code must retain the above copyright
22 notice, this list of conditions and the following disclaimer.
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.
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.
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.
45 /*---------------------------------------------------------------------------*\
47 Introduction to Line Spectrum Pairs (LSPs)
48 ------------------------------------------
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}.
55 A(z) is a polynomial of order lpcrdr with {ak} as the coefficients.
57 A(z) is transformed to P(z) and Q(z) (using a substitution and some
58 algebra), to obtain something like:
60 A(z) = 0.5[P(z)(z+z^-1) + Q(z)(z-z^-1)] (1)
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
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
73 The root so P(z) and Q(z) occur in symmetrical pairs at +/-w, hence
74 the name Line Spectrum Pairs (LSPs).
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
80 \*---------------------------------------------------------------------------*/
83 #include "config-speex.h"
88 #include "stack_alloc.h"
89 #include "math_approx.h"
92 #define M_PI 3.14159265358979323846 /* pi */
101 #define FREQ_SCALE 16384
103 /*#define ANGLE2X(a) (32768*cos(((a)/8192.)))*/
104 #define ANGLE2X(a) (SHL16(spx_cos(a),2))
106 /*#define X2ANGLE(x) (acos(.00006103515625*(x))*LSP_SCALING)*/
107 #define X2ANGLE(x) (spx_acos(x))
110 #include "lsp_bfin.h"
115 /*#define C1 0.99940307
116 #define C2 -0.49558072
117 #define C3 0.03679168*/
119 #define FREQ_SCALE 1.
120 #define ANGLE2X(a) (spx_cos(a))
121 #define X2ANGLE(x) (acos(x))
126 /*---------------------------------------------------------------------------*\
128 FUNCTION....: cheb_poly_eva()
130 AUTHOR......: David Rowe
131 DATE CREATED: 24/2/93
133 This function evaluates a series of Chebyshev polynomials
135 \*---------------------------------------------------------------------------*/
137 #ifndef SPEEX_DISABLE_ENCODER
141 #ifndef OVERRIDE_CHEB_POLY_EVA
142 static inline spx_word32_t
cheb_poly_eva(
143 spx_word16_t
*coef
, /* P or Q coefs in Q13 format */
144 spx_word16_t x
, /* cos of freq (-1.0 to 1.0) in Q14 format */
145 int m
, /* LPC order/2 */
153 /*Prevents overflows*/
159 /* Initialise values */
163 /* Evaluate Chebyshev series formulation usin g iterative approach */
164 sum
= ADD32(EXTEND32(coef
[m
]), EXTEND32(MULT16_16_P14(coef
[m
-1],x
)));
168 b0
= SUB16(MULT16_16_Q13(x
,b0
), b1
);
170 sum
= ADD32(sum
, EXTEND32(MULT16_16_P14(coef
[m
-i
],b0
)));
179 static float cheb_poly_eva(spx_word32_t
*coef
, spx_word16_t x
, int m
, char *stack
)
184 /* Initial conditions */
190 /* Calculate the b_(k) */
193 tmp
=b0
; /* tmp holds the previous value of b0 */
194 b0
=x
*b0
-b1
+coef
[m
-k
]; /* b0 holds its new value based on b0 and b1 */
195 b1
=tmp
; /* b1 holds the previous value of b0 */
198 return(-b1
+.5*x
*b0
+coef
[m
]);
202 /*---------------------------------------------------------------------------*\
204 FUNCTION....: lpc_to_lsp()
206 AUTHOR......: David Rowe
207 DATE CREATED: 24/2/93
209 This function converts LPC coefficients to LSP
212 \*---------------------------------------------------------------------------*/
215 #define SIGN_CHANGE(a,b) (((a)&0x70000000)^((b)&0x70000000)||(b==0))
217 #define SIGN_CHANGE(a,b) (((a)*(b))<0.0)
221 int lpc_to_lsp (spx_coef_t
*a
,int lpcrdr
,spx_lsp_t
*freq
,int nb
,spx_word16_t delta
, char *stack
)
222 /* float *a lpc coefficients */
223 /* int lpcrdr order of LPC coefficients (10) */
224 /* float *freq LSP frequencies in the x domain */
225 /* int nb number of sub-intervals (4) */
226 /* float delta grid spacing interval (0.02) */
230 spx_word16_t temp_xr
,xl
,xr
,xm
=0;
231 spx_word32_t psuml
,psumr
,psumm
,temp_psumr
/*,temp_qsumr*/;
233 VARDECL(spx_word32_t
*Q
); /* ptrs for memory allocation */
234 VARDECL(spx_word32_t
*P
);
235 VARDECL(spx_word16_t
*Q16
); /* ptrs for memory allocation */
236 VARDECL(spx_word16_t
*P16
);
237 spx_word32_t
*px
; /* ptrs of respective P'(z) & Q'(z) */
241 spx_word16_t
*pt
; /* ptr used for cheb_poly_eval()
243 int roots
=0; /* DR 8/2/94: number of roots found */
244 flag
= 1; /* program is searching for a root when,
245 1 else has found one */
246 m
= lpcrdr
/2; /* order of P'(z) & Q'(z) polynomials */
248 /* Allocate memory space for polynomials */
249 ALLOC(Q
, (m
+1), spx_word32_t
);
250 ALLOC(P
, (m
+1), spx_word32_t
);
252 /* determine P'(z)'s and Q'(z)'s coefficients where
253 P'(z) = P(z)/(1 + z^(-1)) and Q'(z) = Q(z)/(1-z^(-1)) */
255 px
= P
; /* initialise ptrs */
264 *px
++ = SUB32(ADD32(EXTEND32(a
[i
]),EXTEND32(a
[lpcrdr
-i
-1])), *p
++);
265 *qx
++ = ADD32(SUB32(EXTEND32(a
[i
]),EXTEND32(a
[lpcrdr
-i
-1])), *q
++);
271 /*if (fabs(*px)>=32768)
272 speex_warning_int("px", *px);
273 if (fabs(*qx)>=32768)
274 speex_warning_int("qx", *qx);*/
280 /* The reason for this lies in the way cheb_poly_eva() is implemented for fixed-point */
281 P
[m
] = PSHR32(P
[m
],3);
282 Q
[m
] = PSHR32(Q
[m
],3);
287 *px
++ = (a
[i
]+a
[lpcrdr
-1-i
]) - *p
++;
288 *qx
++ = (a
[i
]-a
[lpcrdr
-1-i
]) + *q
++;
300 px
= P
; /* re-initialise ptrs */
303 /* now that we have computed P and Q convert to 16 bits to
304 speed up cheb_poly_eval */
306 ALLOC(P16
, m
+1, spx_word16_t
);
307 ALLOC(Q16
, m
+1, spx_word16_t
);
315 /* Search for a zero in P'(z) polynomial first and then alternate to Q'(z).
316 Keep alternating between the two polynomials as each zero is found */
318 xr
= 0; /* initialise xr to zero */
319 xl
= FREQ_SCALE
; /* start at point xl = 1 */
321 for(j
=0;j
<lpcrdr
;j
++){
322 if(j
&1) /* determines whether P' or Q' is eval. */
327 psuml
= cheb_poly_eva(pt
,xl
,m
,stack
); /* evals poly. at xl */
329 while(flag
&& (xr
>= -FREQ_SCALE
)){
331 /* Modified by JMV to provide smaller steps around x=+-1 */
333 dd
= MULT16_16_Q15(delta
,SUB16(FREQ_SCALE
, MULT16_16_Q14(MULT16_16_Q14(xl
,xl
),14000)));
334 if (psuml
<512 && psuml
>-512)
337 dd
=delta
*(1-.9*xl
*xl
);
341 xr
= SUB16(xl
, dd
); /* interval spacing */
342 psumr
= cheb_poly_eva(pt
,xr
,m
,stack
);/* poly(xl-delta_x) */
346 /* if no sign change increment xr and re-evaluate poly(xr). Repeat til
348 if a sign change has occurred the interval is bisected and then
349 checked again for a sign change which determines in which
350 interval the zero lies in.
351 If there is no sign change between poly(xm) and poly(xl) set interval
352 between xm and xr else set interval between xl and xr and repeat till
353 root is located within the specified limits */
355 if(SIGN_CHANGE(psumr
,psuml
))
362 xm
= ADD16(PSHR16(xl
,1),PSHR16(xr
,1)); /* bisect the interval */
364 xm
= .5*(xl
+xr
); /* bisect the interval */
366 psumm
=cheb_poly_eva(pt
,xm
,m
,stack
);
367 /*if(psumm*psuml>0.)*/
368 if(!SIGN_CHANGE(psumm
,psuml
))
378 /* once zero is found, reset initial interval to xr */
379 freq
[j
] = X2ANGLE(xm
);
381 flag
= 0; /* reset flag for next search */
392 #endif /* SPEEX_DISABLE_ENCODER */
394 /*---------------------------------------------------------------------------*\
396 FUNCTION....: lsp_to_lpc()
398 AUTHOR......: David Rowe
399 DATE CREATED: 24/2/93
401 Converts LSP coefficients to LPC coefficients.
403 \*---------------------------------------------------------------------------*/
407 void lsp_to_lpc(spx_lsp_t
*freq
,spx_coef_t
*ak
,int lpcrdr
, char *stack
)
408 /* float *freq array of LSP frequencies in the x domain */
409 /* float *ak array of LPC coefficients */
410 /* int lpcrdr order of LPC coefficients */
413 spx_word32_t xout1
,xout2
,xin
;
414 spx_word32_t mult
, a
;
415 VARDECL(spx_word16_t
*freqn
);
416 VARDECL(spx_word32_t
**xp
);
417 VARDECL(spx_word32_t
*xpmem
);
418 VARDECL(spx_word32_t
**xq
);
419 VARDECL(spx_word32_t
*xqmem
);
424 Reconstruct P(z) and Q(z) by cascading second order polynomials
425 in form 1 - 2cos(w)z(-1) + z(-2), where w is the LSP frequency.
426 In the time domain this is:
428 y(n) = x(n) - 2cos(w)x(n-1) + x(n-2)
430 This is what the ALLOCS below are trying to do:
432 int xp[m+1][lpcrdr+1+2]; // P matrix in QIMP
433 int xq[m+1][lpcrdr+1+2]; // Q matrix in QIMP
435 These matrices store the output of each stage on each row. The
436 final (m-th) row has the output of the final (m-th) cascaded
437 2nd order filter. The first row is the impulse input to the
438 system (not written as it is known).
440 The version below takes advantage of the fact that a lot of the
441 outputs are zero or known, for example if we put an inpulse
442 into the first section the "clock" it 10 times only the first 3
443 outputs samples are non-zero (it's an FIR filter).
446 ALLOC(xp
, (m
+1), spx_word32_t
*);
447 ALLOC(xpmem
, (m
+1)*(lpcrdr
+1+2), spx_word32_t
);
449 ALLOC(xq
, (m
+1), spx_word32_t
*);
450 ALLOC(xqmem
, (m
+1)*(lpcrdr
+1+2), spx_word32_t
);
452 for(i
=0; i
<=m
; i
++) {
453 xp
[i
] = xpmem
+ i
*(lpcrdr
+1+2);
454 xq
[i
] = xqmem
+ i
*(lpcrdr
+1+2);
457 /* work out 2cos terms in Q14 */
459 ALLOC(freqn
, lpcrdr
, spx_word16_t
);
460 for (i
=0;i
<lpcrdr
;i
++)
461 freqn
[i
] = ANGLE2X(freq
[i
]);
463 #define QIMP 21 /* scaling for impulse */
465 xin
= SHL32(EXTEND32(1), (QIMP
-1)); /* 0.5 in QIMP format */
467 /* first col and last non-zero values of each row are trivial */
478 /* 2nd row (first output row) is trivial */
480 xp
[1][3] = -MULT16_32_Q14(freqn
[0],xp
[0][2]);
481 xq
[1][3] = -MULT16_32_Q14(freqn
[1],xq
[0][2]);
485 /* now generate remaining rows */
489 for(j
=1;j
<2*(i
+1)-1;j
++) {
490 mult
= MULT16_32_Q14(freqn
[2*i
],xp
[i
][j
+1]);
491 xp
[i
+1][j
+2] = ADD32(SUB32(xp
[i
][j
+2], mult
), xp
[i
][j
]);
492 mult
= MULT16_32_Q14(freqn
[2*i
+1],xq
[i
][j
+1]);
493 xq
[i
+1][j
+2] = ADD32(SUB32(xq
[i
][j
+2], mult
), xq
[i
][j
]);
496 /* for last col xp[i][j+2] = xq[i][j+2] = 0 */
498 mult
= MULT16_32_Q14(freqn
[2*i
],xp
[i
][j
+1]);
499 xp
[i
+1][j
+2] = SUB32(xp
[i
][j
], mult
);
500 mult
= MULT16_32_Q14(freqn
[2*i
+1],xq
[i
][j
+1]);
501 xq
[i
+1][j
+2] = SUB32(xq
[i
][j
], mult
);
504 /* process last row to extra a{k} */
506 for(j
=1;j
<=lpcrdr
;j
++) {
509 /* final filter sections */
510 a
= PSHR32(xp
[m
][j
+2] + xout1
+ xq
[m
][j
+2] - xout2
, shift
);
514 /* hard limit ak's to +/- 32767 */
516 if (a
< -32767) a
= -32767;
517 if (a
> 32767) a
= 32767;
526 void lsp_to_lpc(spx_lsp_t
*freq
,spx_coef_t
*ak
,int lpcrdr
, char *stack
)
527 /* float *freq array of LSP frequencies in the x domain */
528 /* float *ak array of LPC coefficients */
529 /* int lpcrdr order of LPC coefficients */
534 float xout1
,xout2
,xin1
,xin2
;
536 float *pw
,*n1
,*n2
,*n3
,*n4
=NULL
;
537 VARDECL(float *x_freq
);
540 ALLOC(Wp
, 4*m
+2, float);
543 /* initialise contents of array */
545 for(i
=0;i
<=4*m
+1;i
++){ /* set contents of buffer to 0 */
549 /* Set pointers up */
555 ALLOC(x_freq
, lpcrdr
, float);
556 for (i
=0;i
<lpcrdr
;i
++)
557 x_freq
[i
] = ANGLE2X(freq
[i
]);
559 /* reconstruct P(z) and Q(z) by cascading second order
560 polynomials in form 1 - 2xz(-1) +z(-2), where x is the
563 for(j
=0;j
<=lpcrdr
;j
++){
565 for(i
=0;i
<m
;i
++,i2
+=2){
570 xout1
= xin1
- 2.f
*x_freq
[i2
] * *n1
+ *n2
;
571 xout2
= xin2
- 2.f
*x_freq
[i2
+1] * *n3
+ *n4
;
579 xout1
= xin1
+ *(n4
+1);
580 xout2
= xin2
- *(n4
+2);
582 ak
[j
-1] = (xout1
+ xout2
)*0.5f
;
596 /*Makes sure the LSPs are stable*/
597 void lsp_enforce_margin(spx_lsp_t
*lsp
, int len
, spx_word16_t margin
)
600 spx_word16_t m
= margin
;
601 spx_word16_t m2
= 25736-margin
;
607 for (i
=1;i
<len
-1;i
++)
609 if (lsp
[i
]<lsp
[i
-1]+m
)
612 if (lsp
[i
]>lsp
[i
+1]-m
)
613 lsp
[i
]= SHR16(lsp
[i
],1) + SHR16(lsp
[i
+1]-m
,1);
618 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
)
621 spx_word16_t tmp
= DIV32_16(SHL32(EXTEND32(1 + subframe
),14),nb_subframes
);
622 spx_word16_t tmp2
= 16384-tmp
;
625 interp_lsp
[i
] = MULT16_16_P14(tmp2
,old_lsp
[i
]) + MULT16_16_P14(tmp
,new_lsp
[i
]);
631 /*Makes sure the LSPs are stable*/
632 void lsp_enforce_margin(spx_lsp_t
*lsp
, int len
, spx_word16_t margin
)
635 if (lsp
[0]<LSP_SCALING
*margin
)
636 lsp
[0]=LSP_SCALING
*margin
;
637 if (lsp
[len
-1]>LSP_SCALING
*(M_PI
-margin
))
638 lsp
[len
-1]=LSP_SCALING
*(M_PI
-margin
);
639 for (i
=1;i
<len
-1;i
++)
641 if (lsp
[i
]<lsp
[i
-1]+LSP_SCALING
*margin
)
642 lsp
[i
]=lsp
[i
-1]+LSP_SCALING
*margin
;
644 if (lsp
[i
]>lsp
[i
+1]-LSP_SCALING
*margin
)
645 lsp
[i
]= .5f
* (lsp
[i
] + lsp
[i
+1]-LSP_SCALING
*margin
);
650 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
)
653 float tmp
= (1.0f
+ subframe
)/nb_subframes
;
656 interp_lsp
[i
] = (1-tmp
)*old_lsp
[i
] + tmp
*new_lsp
[i
];