diff options
Diffstat (limited to 'lsp.c')
-rw-r--r-- | lsp.c | 354 |
1 files changed, 173 insertions, 181 deletions
@@ -28,12 +28,14 @@ | |||
28 | along with this program; if not, see <http://www.gnu.org/licenses/>. | 28 | along with this program; if not, see <http://www.gnu.org/licenses/>. |
29 | */ | 29 | */ |
30 | 30 | ||
31 | #include "defines.h" | ||
32 | #include "lsp.h" | 31 | #include "lsp.h" |
32 | |||
33 | #include <math.h> | 33 | #include <math.h> |
34 | #include <stdio.h> | 34 | #include <stdio.h> |
35 | #include <stdlib.h> | 35 | #include <stdlib.h> |
36 | 36 | ||
37 | #include "defines.h" | ||
38 | |||
37 | /*---------------------------------------------------------------------------*\ | 39 | /*---------------------------------------------------------------------------*\ |
38 | 40 | ||
39 | Introduction to Line Spectrum Pairs (LSPs) | 41 | Introduction to Line Spectrum Pairs (LSPs) |
@@ -77,49 +79,45 @@ | |||
77 | AUTHOR......: David Rowe | 79 | AUTHOR......: David Rowe |
78 | DATE CREATED: 24/2/93 | 80 | DATE CREATED: 24/2/93 |
79 | 81 | ||
80 | This function evalutes a series of chebyshev polynomials | 82 | This function evaluates a series of chebyshev polynomials |
81 | 83 | ||
82 | FIXME: performing memory allocation at run time is very inefficient, | 84 | FIXME: performing memory allocation at run time is very inefficient, |
83 | replace with stack variables of MAX_P size. | 85 | replace with stack variables of MAX_P size. |
84 | 86 | ||
85 | \*---------------------------------------------------------------------------*/ | 87 | \*---------------------------------------------------------------------------*/ |
86 | 88 | ||
87 | 89 | static float cheb_poly_eva(float *coef, float x, int order) | |
88 | static float | ||
89 | cheb_poly_eva(float *coef,float x,int order) | ||
90 | /* float coef[] coefficients of the polynomial to be evaluated */ | 90 | /* float coef[] coefficients of the polynomial to be evaluated */ |
91 | /* float x the point where polynomial is to be evaluated */ | 91 | /* float x the point where polynomial is to be evaluated */ |
92 | /* int order order of the polynomial */ | 92 | /* int order order of the polynomial */ |
93 | { | 93 | { |
94 | int i; | 94 | int i; |
95 | float *t,*u,*v,sum; | 95 | float *t, *u, *v, sum; |
96 | float T[(order / 2) + 1]; | 96 | float T[(order / 2) + 1]; |
97 | 97 | ||
98 | /* Initialise pointers */ | 98 | /* Initialise pointers */ |
99 | 99 | ||
100 | t = T; /* T[i-2] */ | 100 | t = T; /* T[i-2] */ |
101 | *t++ = 1.0; | 101 | *t++ = 1.0; |
102 | u = t--; /* T[i-1] */ | 102 | u = t--; /* T[i-1] */ |
103 | *u++ = x; | 103 | *u++ = x; |
104 | v = u--; /* T[i] */ | 104 | v = u--; /* T[i] */ |
105 | 105 | ||
106 | /* Evaluate chebyshev series formulation using iterative approach */ | 106 | /* Evaluate chebyshev series formulation using iterative approach */ |
107 | 107 | ||
108 | for(i=2;i<=order/2;i++) | 108 | for (i = 2; i <= order / 2; i++) |
109 | *v++ = (2*x)*(*u++) - *t++; /* T[i] = 2*x*T[i-1] - T[i-2] */ | 109 | *v++ = (2 * x) * (*u++) - *t++; /* T[i] = 2*x*T[i-1] - T[i-2] */ |
110 | 110 | ||
111 | sum=0.0; /* initialise sum to zero */ | 111 | sum = 0.0; /* initialise sum to zero */ |
112 | t = T; /* reset pointer */ | 112 | t = T; /* reset pointer */ |
113 | 113 | ||
114 | /* Evaluate polynomial and return value also free memory space */ | 114 | /* Evaluate polynomial and return value also free memory space */ |
115 | 115 | ||
116 | for(i=0;i<=order/2;i++) | 116 | for (i = 0; i <= order / 2; i++) sum += coef[(order / 2) - i] * *t++; |
117 | sum+=coef[(order/2)-i]**t++; | ||
118 | 117 | ||
119 | return sum; | 118 | return sum; |
120 | } | 119 | } |
121 | 120 | ||
122 | |||
123 | /*---------------------------------------------------------------------------*\ | 121 | /*---------------------------------------------------------------------------*\ |
124 | 122 | ||
125 | FUNCTION....: lpc_to_lsp() | 123 | FUNCTION....: lpc_to_lsp() |
@@ -130,121 +128,118 @@ cheb_poly_eva(float *coef,float x,int order) | |||
130 | 128 | ||
131 | \*---------------------------------------------------------------------------*/ | 129 | \*---------------------------------------------------------------------------*/ |
132 | 130 | ||
133 | int lpc_to_lsp (float *a, int order, float *freq, int nb, float delta) | 131 | int lpc_to_lsp(float *a, int order, float *freq, int nb, float delta) |
134 | /* float *a lpc coefficients */ | 132 | /* float *a lpc coefficients */ |
135 | /* int order order of LPC coefficients (10) */ | 133 | /* int order order of LPC coefficients (10) */ |
136 | /* float *freq LSP frequencies in radians */ | 134 | /* float *freq LSP frequencies in radians */ |
137 | /* int nb number of sub-intervals (4) */ | 135 | /* int nb number of sub-intervals (4) */ |
138 | /* float delta grid spacing interval (0.02) */ | 136 | /* float delta grid spacing interval (0.02) */ |
139 | { | 137 | { |
140 | float psuml,psumr,psumm,temp_xr,xl,xr,xm = 0; | 138 | float psuml, psumr, psumm, temp_xr, xl, xr, xm = 0; |
141 | float temp_psumr; | 139 | float temp_psumr; |
142 | int i,j,m,flag,k; | 140 | int i, j, m, flag, k; |
143 | float *px; /* ptrs of respective P'(z) & Q'(z) */ | 141 | float *px; /* ptrs of respective P'(z) & Q'(z) */ |
144 | float *qx; | 142 | float *qx; |
145 | float *p; | 143 | float *p; |
146 | float *q; | 144 | float *q; |
147 | float *pt; /* ptr used for cheb_poly_eval() | 145 | float *pt; /* ptr used for cheb_poly_eval() |
148 | whether P' or Q' */ | 146 | whether P' or Q' */ |
149 | int roots=0; /* number of roots found */ | 147 | int roots = 0; /* number of roots found */ |
150 | float Q[order + 1]; | 148 | float Q[order + 1]; |
151 | float P[order + 1]; | 149 | float P[order + 1]; |
152 | 150 | ||
151 | flag = 1; | ||
152 | m = order / 2; /* order of P'(z) & Q'(z) polynimials */ | ||
153 | |||
154 | /* Allocate memory space for polynomials */ | ||
155 | |||
156 | /* determine P'(z)'s and Q'(z)'s coefficients where | ||
157 | P'(z) = P(z)/(1 + z^(-1)) and Q'(z) = Q(z)/(1-z^(-1)) */ | ||
158 | |||
159 | px = P; /* initilaise ptrs */ | ||
160 | qx = Q; | ||
161 | p = px; | ||
162 | q = qx; | ||
163 | *px++ = 1.0; | ||
164 | *qx++ = 1.0; | ||
165 | for (i = 1; i <= m; i++) { | ||
166 | *px++ = a[i] + a[order + 1 - i] - *p++; | ||
167 | *qx++ = a[i] - a[order + 1 - i] + *q++; | ||
168 | } | ||
169 | px = P; | ||
170 | qx = Q; | ||
171 | for (i = 0; i < m; i++) { | ||
172 | *px = 2 * *px; | ||
173 | *qx = 2 * *qx; | ||
174 | px++; | ||
175 | qx++; | ||
176 | } | ||
177 | px = P; /* re-initialise ptrs */ | ||
178 | qx = Q; | ||
179 | |||
180 | /* Search for a zero in P'(z) polynomial first and then alternate to Q'(z). | ||
181 | Keep alternating between the two polynomials as each zero is found */ | ||
182 | |||
183 | xr = 0; /* initialise xr to zero */ | ||
184 | xl = 1.0; /* start at point xl = 1 */ | ||
185 | |||
186 | for (j = 0; j < order; j++) { | ||
187 | if (j % 2) /* determines whether P' or Q' is eval. */ | ||
188 | pt = qx; | ||
189 | else | ||
190 | pt = px; | ||
191 | |||
192 | psuml = cheb_poly_eva(pt, xl, order); /* evals poly. at xl */ | ||
153 | flag = 1; | 193 | flag = 1; |
154 | m = order/2; /* order of P'(z) & Q'(z) polynimials */ | 194 | while (flag && (xr >= -1.0)) { |
155 | 195 | xr = xl - delta; /* interval spacing */ | |
156 | /* Allocate memory space for polynomials */ | 196 | psumr = cheb_poly_eva(pt, xr, order); /* poly(xl-delta_x) */ |
157 | 197 | temp_psumr = psumr; | |
158 | /* determine P'(z)'s and Q'(z)'s coefficients where | 198 | temp_xr = xr; |
159 | P'(z) = P(z)/(1 + z^(-1)) and Q'(z) = Q(z)/(1-z^(-1)) */ | 199 | |
160 | 200 | /* if no sign change increment xr and re-evaluate | |
161 | px = P; /* initilaise ptrs */ | 201 | poly(xr). Repeat til sign change. if a sign change has |
162 | qx = Q; | 202 | occurred the interval is bisected and then checked again |
163 | p = px; | 203 | for a sign change which determines in which interval the |
164 | q = qx; | 204 | zero lies in. If there is no sign change between poly(xm) |
165 | *px++ = 1.0; | 205 | and poly(xl) set interval between xm and xr else set |
166 | *qx++ = 1.0; | 206 | interval between xl and xr and repeat till root is located |
167 | for(i=1;i<=m;i++){ | 207 | within the specified limits */ |
168 | *px++ = a[i]+a[order+1-i]-*p++; | 208 | |
169 | *qx++ = a[i]-a[order+1-i]+*q++; | 209 | if (((psumr * psuml) < 0.0) || (psumr == 0.0)) { |
170 | } | 210 | roots++; |
171 | px = P; | 211 | |
172 | qx = Q; | 212 | psumm = psuml; |
173 | for(i=0;i<m;i++){ | 213 | for (k = 0; k <= nb; k++) { |
174 | *px = 2**px; | 214 | xm = (xl + xr) / 2; /* bisect the interval */ |
175 | *qx = 2**qx; | 215 | psumm = cheb_poly_eva(pt, xm, order); |
176 | px++; | 216 | if (psumm * psuml > 0.) { |
177 | qx++; | 217 | psuml = psumm; |
178 | } | 218 | xl = xm; |
179 | px = P; /* re-initialise ptrs */ | 219 | } else { |
180 | qx = Q; | 220 | psumr = psumm; |
181 | 221 | xr = xm; | |
182 | /* Search for a zero in P'(z) polynomial first and then alternate to Q'(z). | 222 | } |
183 | Keep alternating between the two polynomials as each zero is found */ | 223 | } |
184 | 224 | ||
185 | xr = 0; /* initialise xr to zero */ | 225 | /* once zero is found, reset initial interval to xr */ |
186 | xl = 1.0; /* start at point xl = 1 */ | 226 | freq[j] = (xm); |
187 | 227 | xl = xm; | |
188 | 228 | flag = 0; /* reset flag for next search */ | |
189 | for(j=0;j<order;j++){ | 229 | } else { |
190 | if(j%2) /* determines whether P' or Q' is eval. */ | 230 | psuml = temp_psumr; |
191 | pt = qx; | 231 | xl = temp_xr; |
192 | else | 232 | } |
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 | } | 233 | } |
234 | } | ||
240 | 235 | ||
241 | /* convert from x domain to radians */ | 236 | /* convert from x domain to radians */ |
242 | 237 | ||
243 | for(i=0; i<order; i++) { | 238 | for (i = 0; i < order; i++) { |
244 | freq[i] = acosf(freq[i]); | 239 | freq[i] = acosf(freq[i]); |
245 | } | 240 | } |
246 | 241 | ||
247 | return(roots); | 242 | return (roots); |
248 | } | 243 | } |
249 | 244 | ||
250 | /*---------------------------------------------------------------------------*\ | 245 | /*---------------------------------------------------------------------------*\ |
@@ -263,59 +258,56 @@ void lsp_to_lpc(float *lsp, float *ak, int order) | |||
263 | /* float *ak array of LPC coefficients */ | 258 | /* float *ak array of LPC coefficients */ |
264 | /* int order order of LPC coefficients */ | 259 | /* int order order of LPC coefficients */ |
265 | 260 | ||
266 | |||
267 | { | 261 | { |
268 | int i,j; | 262 | int i, j; |
269 | float xout1,xout2,xin1,xin2; | 263 | float xout1, xout2, xin1, xin2; |
270 | float *pw,*n1,*n2,*n3,*n4 = 0; | 264 | float *pw, *n1, *n2, *n3, *n4 = 0; |
271 | float freq[order]; | 265 | float freq[order]; |
272 | float Wp[(order * 4) + 2]; | 266 | float Wp[(order * 4) + 2]; |
273 | 267 | ||
274 | /* convert from radians to the x=cos(w) domain */ | 268 | /* convert from radians to the x=cos(w) domain */ |
275 | 269 | ||
276 | for(i=0; i<order; i++) | 270 | for (i = 0; i < order; i++) freq[i] = cosf(lsp[i]); |
277 | freq[i] = cosf(lsp[i]); | 271 | |
278 | 272 | pw = Wp; | |
279 | pw = Wp; | 273 | |
280 | 274 | /* initialise contents of array */ | |
281 | /* initialise contents of array */ | 275 | |
282 | 276 | for (i = 0; i <= 4 * (order / 2) + 1; i++) { /* set contents of buffer to 0 */ | |
283 | for(i=0;i<=4*(order/2)+1;i++){ /* set contents of buffer to 0 */ | 277 | *pw++ = 0.0; |
284 | *pw++ = 0.0; | 278 | } |
285 | } | 279 | |
286 | 280 | /* Set pointers up */ | |
287 | /* Set pointers up */ | 281 | |
288 | 282 | pw = Wp; | |
289 | pw = Wp; | 283 | xin1 = 1.0; |
290 | xin1 = 1.0; | 284 | xin2 = 1.0; |
291 | xin2 = 1.0; | 285 | |
292 | 286 | /* reconstruct P(z) and Q(z) by cascading second order polynomials | |
293 | /* reconstruct P(z) and Q(z) by cascading second order polynomials | 287 | in form 1 - 2xz(-1) +z(-2), where x is the LSP coefficient */ |
294 | in form 1 - 2xz(-1) +z(-2), where x is the LSP coefficient */ | 288 | |
295 | 289 | for (j = 0; j <= order; j++) { | |
296 | for(j=0;j<=order;j++){ | 290 | for (i = 0; i < (order / 2); i++) { |
297 | for(i=0;i<(order/2);i++){ | 291 | n1 = pw + (i * 4); |
298 | n1 = pw+(i*4); | 292 | n2 = n1 + 1; |
299 | n2 = n1 + 1; | 293 | n3 = n2 + 1; |
300 | n3 = n2 + 1; | 294 | n4 = n3 + 1; |
301 | n4 = n3 + 1; | 295 | xout1 = xin1 - 2 * (freq[2 * i]) * *n1 + *n2; |
302 | xout1 = xin1 - 2*(freq[2*i]) * *n1 + *n2; | 296 | xout2 = xin2 - 2 * (freq[2 * i + 1]) * *n3 + *n4; |
303 | xout2 = xin2 - 2*(freq[2*i+1]) * *n3 + *n4; | 297 | *n2 = *n1; |
304 | *n2 = *n1; | 298 | *n4 = *n3; |
305 | *n4 = *n3; | 299 | *n1 = xin1; |
306 | *n1 = xin1; | 300 | *n3 = xin2; |
307 | *n3 = xin2; | 301 | xin1 = xout1; |
308 | xin1 = xout1; | 302 | xin2 = xout2; |
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 | } | 303 | } |
304 | xout1 = xin1 + *(n4 + 1); | ||
305 | xout2 = xin2 - *(n4 + 2); | ||
306 | ak[j] = (xout1 + xout2) * 0.5; | ||
307 | *(n4 + 1) = xin1; | ||
308 | *(n4 + 2) = xin2; | ||
309 | |||
310 | xin1 = 0.0; | ||
311 | xin2 = 0.0; | ||
312 | } | ||
320 | } | 313 | } |
321 | |||