1 /*---------------------------------------------------------------------------*\
2
3 FILE........: lsp.c
4 AUTHOR......: David Rowe
5 DATE CREATED: 24/2/93
6
7
8 This file contains functions for LPC to LSP conversion and LSP to
9 LPC conversion. Note that the LSP coefficients are not in radians
10 format but in the x domain of the unit circle.
11
12 \*---------------------------------------------------------------------------*/
13
14 /*
15 Copyright (C) 2009 David Rowe
16
17 All rights reserved.
18
19 This program is free software; you can redistribute it and/or modify
20 it under the terms of the GNU Lesser General Public License version 2.1, as
21 published by the Free Software Foundation. This program is
22 distributed in the hope that it will be useful, but WITHOUT ANY
23 WARRANTY; without even the implied warranty of MERCHANTABILITY or
24 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
25 License for more details.
26
27 You should have received a copy of the GNU Lesser General Public License
28 along with this program; if not, see <http://www.gnu.org/licenses/>.
29 */
30
31 #include "defines.h"
32 #include "lsp.h"
33 #include <math.h>
34 #include <stdio.h>
35 #include <stdlib.h>
36
37 /*---------------------------------------------------------------------------*\
38
39 Introduction to Line Spectrum Pairs (LSPs)
40 ------------------------------------------
41
42 LSPs are used to encode the LPC filter coefficients {ak} for
43 transmission over the channel. LSPs have several properties (like
44 less sensitivity to quantisation noise) that make them superior to
45 direct quantisation of {ak}.
46
47 A(z) is a polynomial of order lpcrdr with {ak} as the coefficients.
48
49 A(z) is transformed to P(z) and Q(z) (using a substitution and some
50 algebra), to obtain something like:
51
52 A(z) = 0.5[P(z)(z+z^-1) + Q(z)(z-z^-1)] (1)
53
54 As you can imagine A(z) has complex zeros all over the z-plane. P(z)
55 and Q(z) have the very neat property of only having zeros _on_ the
56 unit circle. So to find them we take a test point z=exp(jw) and
57 evaluate P (exp(jw)) and Q(exp(jw)) using a grid of points between 0
58 and pi.
59
60 The zeros (roots) of P(z) also happen to alternate, which is why we
61 swap coefficients as we find roots. So the process of finding the
62 LSP frequencies is basically finding the roots of 5th order
63 polynomials.
64
65 The root so P(z) and Q(z) occur in symmetrical pairs at +/-w, hence
66 the name Line Spectrum Pairs (LSPs).
67
68 To convert back to ak we just evaluate (1), "clocking" an impulse
69 thru it lpcrdr times gives us the impulse response of A(z) which is
70 {ak}.
71
72 \*---------------------------------------------------------------------------*/
73
74 /*---------------------------------------------------------------------------*\
75
76 FUNCTION....: cheb_poly_eva()
77 AUTHOR......: David Rowe
78 DATE CREATED: 24/2/93
79
80 This function evalutes a series of chebyshev polynomials
81
82 FIXME: performing memory allocation at run time is very inefficient,
83 replace with stack variables of MAX_P size.
84
85 \*---------------------------------------------------------------------------*/
86
87
88 static float
cheb_poly_eva(float * coef,float x,int order)89 cheb_poly_eva(float *coef,float x,int order)
90 /* float coef[] coefficients of the polynomial to be evaluated */
91 /* float x the point where polynomial is to be evaluated */
92 /* int order order of the polynomial */
93 {
94 int i;
95 float *t,*u,*v,sum;
96 float T[(order / 2) + 1];
97
98 /* Initialise pointers */
99
100 t = T; /* T[i-2] */
101 *t++ = 1.0;
102 u = t--; /* T[i-1] */
103 *u++ = x;
104 v = u--; /* T[i] */
105
106 /* Evaluate chebyshev series formulation using iterative approach */
107
108 for(i=2;i<=order/2;i++)
109 *v++ = (2*x)*(*u++) - *t++; /* T[i] = 2*x*T[i-1] - T[i-2] */
110
111 sum=0.0; /* initialise sum to zero */
112 t = T; /* reset pointer */
113
114 /* Evaluate polynomial and return value also free memory space */
115
116 for(i=0;i<=order/2;i++)
117 sum+=coef[(order/2)-i]**t++;
118
119 return sum;
120 }
121
122
123 /*---------------------------------------------------------------------------*\
124
125 FUNCTION....: lpc_to_lsp()
126 AUTHOR......: David Rowe
127 DATE CREATED: 24/2/93
128
129 This function converts LPC coefficients to LSP coefficients.
130
131 \*---------------------------------------------------------------------------*/
132
lpc_to_lsp(float * a,int order,float * freq,int nb,float delta)133 int lpc_to_lsp (float *a, int order, float *freq, int nb, float delta)
134 /* float *a lpc coefficients */
135 /* int order order of LPC coefficients (10) */
136 /* float *freq LSP frequencies in radians */
137 /* int nb number of sub-intervals (4) */
138 /* float delta grid spacing interval (0.02) */
139 {
140 float psuml,psumr,psumm,temp_xr,xl,xr,xm = 0;
141 float temp_psumr;
142 int i,j,m,flag,k;
143 float *px; /* ptrs of respective P'(z) & Q'(z) */
144 float *qx;
145 float *p;
146 float *q;
147 float *pt; /* ptr used for cheb_poly_eval()
148 whether P' or Q' */
149 int roots=0; /* number of roots found */
150 float Q[order + 1];
151 float P[order + 1];
152
153 flag = 1;
154 m = order/2; /* order of P'(z) & Q'(z) polynimials */
155
156 /* Allocate memory space for polynomials */
157
158 /* determine P'(z)'s and Q'(z)'s coefficients where
159 P'(z) = P(z)/(1 + z^(-1)) and Q'(z) = Q(z)/(1-z^(-1)) */
160
161 px = P; /* initilaise ptrs */
162 qx = Q;
163 p = px;
164 q = qx;
165 *px++ = 1.0;
166 *qx++ = 1.0;
167 for(i=1;i<=m;i++){
168 *px++ = a[i]+a[order+1-i]-*p++;
169 *qx++ = a[i]-a[order+1-i]+*q++;
170 }
171 px = P;
172 qx = Q;
173 for(i=0;i<m;i++){
174 *px = 2**px;
175 *qx = 2**qx;
176 px++;
177 qx++;
178 }
179 px = P; /* re-initialise ptrs */
180 qx = Q;
181
182 /* Search for a zero in P'(z) polynomial first and then alternate to Q'(z).
183 Keep alternating between the two polynomials as each zero is found */
184
185 xr = 0; /* initialise xr to zero */
186 xl = 1.0; /* start at point xl = 1 */
187
188
189 for(j=0;j<order;j++){
190 if(j%2) /* determines whether P' or Q' is eval. */
191 pt = qx;
192 else
193 pt = px;
194
195 psuml = cheb_poly_eva(pt,xl,order); /* evals poly. at xl */
196 flag = 1;
197 while(flag && (xr >= -1.0)){
198 xr = xl - delta ; /* interval spacing */
199 psumr = cheb_poly_eva(pt,xr,order);/* poly(xl-delta_x) */
200 temp_psumr = psumr;
201 temp_xr = xr;
202
203 /* if no sign change increment xr and re-evaluate
204 poly(xr). Repeat til sign change. if a sign change has
205 occurred the interval is bisected and then checked again
206 for a sign change which determines in which interval the
207 zero lies in. If there is no sign change between poly(xm)
208 and poly(xl) set interval between xm and xr else set
209 interval between xl and xr and repeat till root is located
210 within the specified limits */
211
212 if(((psumr*psuml)<0.0) || (psumr == 0.0)){
213 roots++;
214
215 psumm=psuml;
216 for(k=0;k<=nb;k++){
217 xm = (xl+xr)/2; /* bisect the interval */
218 psumm=cheb_poly_eva(pt,xm,order);
219 if(psumm*psuml>0.){
220 psuml=psumm;
221 xl=xm;
222 }
223 else{
224 psumr=psumm;
225 xr=xm;
226 }
227 }
228
229 /* once zero is found, reset initial interval to xr */
230 freq[j] = (xm);
231 xl = xm;
232 flag = 0; /* reset flag for next search */
233 }
234 else{
235 psuml=temp_psumr;
236 xl=temp_xr;
237 }
238 }
239 }
240
241 /* convert from x domain to radians */
242
243 for(i=0; i<order; i++) {
244 freq[i] = acosf(freq[i]);
245 }
246
247 return(roots);
248 }
249
250 /*---------------------------------------------------------------------------*\
251
252 FUNCTION....: lsp_to_lpc()
253 AUTHOR......: David Rowe
254 DATE CREATED: 24/2/93
255
256 This function converts LSP coefficients to LPC coefficients. In the
257 Speex code we worked out a way to simplify this significantly.
258
259 \*---------------------------------------------------------------------------*/
260
lsp_to_lpc(float * lsp,float * ak,int order)261 void lsp_to_lpc(float *lsp, float *ak, int order)
262 /* float *freq array of LSP frequencies in radians */
263 /* float *ak array of LPC coefficients */
264 /* int order order of LPC coefficients */
265
266
267 {
268 int i,j;
269 float xout1,xout2,xin1,xin2;
270 float *pw,*n1,*n2,*n3,*n4 = 0;
271 float freq[order];
272 float Wp[(order * 4) + 2];
273
274 /* convert from radians to the x=cos(w) domain */
275
276 for(i=0; i<order; i++)
277 freq[i] = cosf(lsp[i]);
278
279 pw = Wp;
280
281 /* initialise contents of array */
282
283 for(i=0;i<=4*(order/2)+1;i++){ /* set contents of buffer to 0 */
284 *pw++ = 0.0;
285 }
286
287 /* Set pointers up */
288
289 pw = Wp;
290 xin1 = 1.0;
291 xin2 = 1.0;
292
293 /* reconstruct P(z) and Q(z) by cascading second order polynomials
294 in form 1 - 2xz(-1) +z(-2), where x is the LSP coefficient */
295
296 for(j=0;j<=order;j++){
297 for(i=0;i<(order/2);i++){
298 n1 = pw+(i*4);
299 n2 = n1 + 1;
300 n3 = n2 + 1;
301 n4 = n3 + 1;
302 xout1 = xin1 - 2*(freq[2*i]) * *n1 + *n2;
303 xout2 = xin2 - 2*(freq[2*i+1]) * *n3 + *n4;
304 *n2 = *n1;
305 *n4 = *n3;
306 *n1 = xin1;
307 *n3 = xin2;
308 xin1 = xout1;
309 xin2 = xout2;
310 }
311 xout1 = xin1 + *(n4+1);
312 xout2 = xin2 - *(n4+2);
313 ak[j] = (xout1 + xout2)*0.5;
314 *(n4+1) = xin1;
315 *(n4+2) = xin2;
316
317 xin1 = 0.0;
318 xin2 = 0.0;
319 }
320 }
321
322