diff options
author | erdgeist@erdgeist.org <erdgeist@bauklotz.fritz.box> | 2019-07-04 23:26:09 +0200 |
---|---|---|
committer | erdgeist@erdgeist.org <erdgeist@bauklotz.fritz.box> | 2019-07-04 23:26:09 +0200 |
commit | f02dfce6e6c34b3d8a7b8a0e784b506178e331fa (patch) | |
tree | 45556e6104242d4702689760433d7321ae74ec17 /lsp.c |
stripdown of version 0.9
Diffstat (limited to 'lsp.c')
-rw-r--r-- | lsp.c | 321 |
1 files changed, 321 insertions, 0 deletions
@@ -0,0 +1,321 @@ | |||
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 | ||
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 | |||
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 | |||
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 | |||