summaryrefslogtreecommitdiff
path: root/lsp.c
diff options
context:
space:
mode:
authorerdgeist <erdgeist@erdgeist.org>2025-08-15 12:42:40 +0200
committererdgeist <erdgeist@erdgeist.org>2025-08-15 12:42:40 +0200
commit30325d24d107dbf133da39f7c96d1510fd1c9449 (patch)
tree932baa5b2a4475821f16dccf9e3e05011daa6d92 /lsp.c
parent9022d768021bbe15c7815cc6f8b64218b46f0e10 (diff)
Bump to codec2 version 1.2.0erdgeist-bump-to-1.2.0
Diffstat (limited to 'lsp.c')
-rw-r--r--lsp.c354
1 files changed, 173 insertions, 181 deletions
diff --git a/lsp.c b/lsp.c
index 05d190e..aa7bac3 100644
--- a/lsp.c
+++ b/lsp.c
@@ -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 89static float cheb_poly_eva(float *coef, float x, int order)
88static float
89cheb_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
133int lpc_to_lsp (float *a, int order, float *freq, int nb, float delta) 131int 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