IT++ 4.3.1
lpcfunc.cpp
Go to the documentation of this file.
1
29
31#include <itpp/base/matfunc.h>
32#include <itpp/signal/sigfun.h>
33#include <itpp/stat/misc_stat.h>
34#include <iostream>
35
37
38using std::cout;
39using std::endl;
40
41namespace itpp
42{
43
44// Autocorrelation sequence to reflection coefficients conversion.
45vec ac2rc(const vec &ac);
46// Autocorrelation sequence to prediction polynomial conversion.
47vec ac2poly(const vec &ac);
48// Inverse sine parameters to reflection coefficients conversion.
49vec is2rc(const vec &is);
50// Reflection coefficients to autocorrelation sequence conversion.
51vec rc2ac(const vec &rc);
52// Reflection coefficients to inverse sine parameters conversion.
53vec rc2is(const vec &rc);
54
55vec autocorr(const vec &x, int order)
56{
57 if (order < 0) order = x.size();
58
59 vec R(order + 1);
60 double sum;
61 int i, j;
62
63 for (i = 0;i < order + 1;i++) {
64 sum = 0;
65 for (j = 0;j < x.size() - i;j++) {
66 sum += x[j] * x[j+i];
67 }
68 R[i] = sum;
69 }
70 return R;
71}
72
73vec levinson(const vec &R2, int order)
74{
75 vec R = R2;
76 R[0] = R[0] * (1. + 1.e-9);
77
78 if (order < 0) order = R.length() - 1;
79 double k, alfa, s;
80 double *any = new double[order+1];
81 double *a = new double[order+1];
82 int j, m;
83 vec out(order + 1);
84
85 a[0] = 1;
86 alfa = R[0];
87 if (alfa <= 0) {
88 out.clear();
89 out[0] = 1;
90 return out;
91 }
92 for (m = 1;m <= order;m++) {
93 s = 0;
94 for (j = 1;j < m;j++) {
95 s = s + a[j] * R[m-j];
96 }
97
98 k = -(R[m] + s) / alfa;
99 if (fabs(k) >= 1.0) {
100 cout << "levinson : panic! abs(k)>=1, order " << m << ". Aborting..." << endl ;
101 for (j = m;j <= order;j++) {
102 a[j] = 0;
103 }
104 break;
105 }
106 for (j = 1;j < m;j++) {
107 any[j] = a[j] + k * a[m-j];
108 }
109 for (j = 1;j < m;j++) {
110 a[j] = any[j];
111 }
112 a[m] = k;
113 alfa = alfa * (1 - k * k);
114 }
115 for (j = 0;j < out.length();j++) {
116 out[j] = a[j];
117 }
118 delete any;
119 delete a;
120 return out;
121}
122
123vec lpc(const vec &x, int order)
124{
125 return levinson(autocorr(x, order), order);
126}
127
128vec poly2ac(const vec &poly)
129{
130 vec a = poly;
131 int order = a.length() - 1;
132 double alfa, s, *any = new double[order+1];
133 int j, m;
134 vec r(order + 1);
135 vec k = poly2rc(a);
136
137 it_error_if(a[0] != 1, "poly2ac : not an lpc filter");
138 r[0] = 1;
139 alfa = 1;
140 for (m = 1;m <= order;m++) {
141 s = 0;
142 for (j = 1;j < m;j++) {
143 s = s + a[j] * r[m-j];
144 }
145 r[m] = -s - alfa * k[m-1];
146 for (j = 1;j < m;j++) {
147 any[j] = a[j] + k[m-1] * a[m-j];
148 }
149 for (j = 1;j < m;j++) {
150 a[j] = any[j];
151 }
152 a[m] = k[m-1];
153 alfa = alfa * (1 - sqr(k[m-1]));
154 }
155 delete any;
156 return r;
157}
158
159vec poly2rc(const vec &a)
160{
161 // a is [1 xx xx xx], a.size()=order+1
162 int m, i;
163 int order = a.size() - 1;
164 vec k(order);
165 vec any(order + 1), aold(a);
166
167 for (m = order - 1;m > 0;m--) {
168 k[m] = aold[m+1] ;
169 if (fabs(k[m]) > 1) k[m] = 1.0 / k[m];
170 for (i = 0;i < m;i++) {
171 any[i+1] = (aold[i+1] - aold[m-i] * k[m]) / (1 - k[m] * k[m]);
172 }
173 aold = any;
174 }
175 k[0] = any[1];
176 if (fabs(k[0]) > 1) k[0] = 1.0 / k[0];
177 return k;
178}
179
180vec rc2poly(const vec &k)
181{
182 int m, i;
183 vec a(k.length() + 1), any(k.length() + 1);
184
185 a[0] = 1;
186 any[0] = 1;
187 a[1] = k[0];
188 for (m = 1;m < k.size();m++) {
189 any[m+1] = k[m];
190 for (i = 0;i < m;i++) {
191 any[i+1] = a[i+1] + a[m-i] * k[m];
192 }
193 a = any;
194 }
195 return a;
196}
197
198vec rc2lar(const vec &k)
199{
200 short m;
201 vec LAR(k.size());
202
203 for (m = 0;m < k.size();m++) {
204 LAR[m] = std::log((1 + k[m]) / (1 - k[m]));
205 }
206 return LAR;
207}
208
209vec lar2rc(const vec &LAR)
210{
211 short m;
212 vec k(LAR.size());
213
214 for (m = 0;m < LAR.size();m++) {
215 k[m] = (std::exp(LAR[m]) - 1) / (std::exp(LAR[m]) + 1);
216 }
217 return k;
218}
219
220double FNevChebP_double(double x, const double c[], int n)
221{
222 int i;
223 double b0 = 0.0, b1 = 0.0, b2 = 0.0;
224
225 for (i = n - 1; i >= 0; --i) {
226 b2 = b1;
227 b1 = b0;
228 b0 = 2.0 * x * b1 - b2 + c[i];
229 }
230 return (0.5 * (b0 - b2 + c[0]));
231}
232
233double FNevChebP(double x, const double c[], int n)
234{
235 int i;
236 double b0 = 0.0, b1 = 0.0, b2 = 0.0;
237
238 for (i = n - 1; i >= 0; --i) {
239 b2 = b1;
240 b1 = b0;
241 b0 = 2.0 * x * b1 - b2 + c[i];
242 }
243 return (0.5 * (b0 - b2 + c[0]));
244}
245
246vec poly2lsf(const vec &pc)
247{
248 int np = pc.length() - 1;
249 vec lsf(np);
250
251 vec fa((np + 1) / 2 + 1), fb((np + 1) / 2 + 1);
252 vec ta((np + 1) / 2 + 1), tb((np + 1) / 2 + 1);
253 double *t;
254 double xlow, xmid, xhigh;
255 double ylow, ymid, yhigh;
256 double xroot;
257 double dx;
258 int i, j, nf;
259 int odd;
260 int na, nb, n;
261 double ss, aa;
262 double DW = (0.02 * pi);
263 int NBIS = 4;
264
265 odd = (np % 2 != 0);
266 if (odd) {
267 nb = (np + 1) / 2;
268 na = nb + 1;
269 }
270 else {
271 nb = np / 2 + 1;
272 na = nb;
273 }
274
275 fa[0] = 1.0;
276 for (i = 1, j = np; i < na; ++i, --j)
277 fa[i] = pc[i] + pc[j];
278
279 fb[0] = 1.0;
280 for (i = 1, j = np; i < nb; ++i, --j)
281 fb[i] = pc[i] - pc[j];
282
283 if (odd) {
284 for (i = 2; i < nb; ++i)
285 fb[i] = fb[i] + fb[i-2];
286 }
287 else {
288 for (i = 1; i < na; ++i) {
289 fa[i] = fa[i] - fa[i-1];
290 fb[i] = fb[i] + fb[i-1];
291 }
292 }
293
294 ta[0] = fa[na-1];
295 for (i = 1, j = na - 2; i < na; ++i, --j)
296 ta[i] = 2.0 * fa[j];
297
298 tb[0] = fb[nb-1];
299 for (i = 1, j = nb - 2; i < nb; ++i, --j)
300 tb[i] = 2.0 * fb[j];
301
302 nf = 0;
303 t = ta._data();
304 n = na;
305 xroot = 2.0;
306 xlow = 1.0;
307 ylow = FNevChebP_double(xlow, t, n);
308
309
310 ss = std::sin(DW);
311 aa = 4.0 - 4.0 * std::cos(DW) - ss;
312 while (xlow > -1.0 && nf < np) {
313 xhigh = xlow;
314 yhigh = ylow;
315 dx = aa * xhigh * xhigh + ss;
316 xlow = xhigh - dx;
317 if (xlow < -1.0)
318 xlow = -1.0;
319 ylow = FNevChebP_double(xlow, t, n);
320 if (ylow * yhigh <= 0.0) {
321 dx = xhigh - xlow;
322 for (i = 1; i <= NBIS; ++i) {
323 dx = 0.5 * dx;
324 xmid = xlow + dx;
325 ymid = FNevChebP_double(xmid, t, n);
326 if (ylow * ymid <= 0.0) {
327 yhigh = ymid;
328 xhigh = xmid;
329 }
330 else {
331 ylow = ymid;
332 xlow = xmid;
333 }
334 }
335 if (yhigh != ylow)
336 xmid = xlow + dx * ylow / (ylow - yhigh);
337 else
338 xmid = xlow + dx;
339 lsf[nf] = std::acos((double) xmid);
340 ++nf;
341 if (xmid >= xroot) {
342 xmid = xlow - dx;
343 }
344 xroot = xmid;
345 if (t == ta._data()) {
346 t = tb._data();
347 n = nb;
348 }
349 else {
350 t = ta._data();
351 n = na;
352 }
353 xlow = xmid;
354 ylow = FNevChebP_double(xlow, t, n);
355 }
356 }
357 if (nf != np) {
358 cout << "poly2lsf: WARNING: failed to find all lsfs" << endl ;
359 }
360 return lsf;
361}
362
363vec lsf2poly(const vec &f)
364{
365 int m = f.length();
366 vec pc(m + 1);
367 double c1, c2, *a;
368 vec p(m + 1), q(m + 1);
369 int mq, n, i, nor;
370
371 it_error_if(m % 2 != 0, "lsf2poly: THIS ROUTINE WORKS ONLY FOR EVEN m");
372 pc[0] = 1.0;
373 a = pc._data() + 1;
374 mq = m >> 1;
375 for (i = 0 ; i <= m ; i++) {
376 q[i] = 0.;
377 p[i] = 0.;
378 }
379 p[0] = q[0] = 1.;
380 for (n = 1; n <= mq; n++) {
381 nor = 2 * n;
382 c1 = 2 * std::cos(f[nor-1]);
383 c2 = 2 * std::cos(f[nor-2]);
384 for (i = nor; i >= 2; i--) {
385 q[i] += q[i-2] - c1 * q[i-1];
386 p[i] += p[i-2] - c2 * p[i-1];
387 }
388 q[1] -= c1;
389 p[1] -= c2;
390 }
391 a[0] = 0.5 * (p[1] + q[1]);
392 for (i = 1, n = 2; i < m ; i++, n++)
393 a[i] = 0.5 * (p[i] + p[n] + q[n] - q[i]);
394
395 return pc;
396}
397
398vec poly2cepstrum(const vec &a)
399{
400 vec c(a.length() - 1);
401
402 for (int n = 1;n <= c.length();n++) {
403 c[n-1] = a[n];
404 for (int k = 1;k < n;k++) {
405 c[n-1] -= double(k) / n * a[n-k] * c[k-1];
406 }
407 }
408 return c;
409}
410
411vec poly2cepstrum(const vec &a, int num)
412{
413 it_error_if(num < a.length(), "a2cepstrum : not allowed cepstrum length");
414 vec c(num);
415 int n;
416
417 for (n = 1;n < a.length();n++) {
418 c[n-1] = a[n];
419 for (int k = 1;k < n;k++) {
420 c[n-1] -= double(k) / n * a[n-k] * c[k-1];
421 }
422 }
423 for (n = a.length();n <= c.length();n++) {
424 c[n-1] = 0;
425 for (int k = n - a.length() + 1;k < n;k++) {
426 c[n-1] -= double(k) / n * a[n-k] * c[k-1];
427 }
428 }
429 return c;
430}
431
432vec cepstrum2poly(const vec &c)
433{
434 vec a(c.length() + 1);
435
436 a[0] = 1;
437 for (int n = 1;n <= c.length();n++) {
438 a[n] = c[n-1];
439 for (int k = 1;k < n;k++) {
440 a[n] += double(k) / n * a[n-k] * c[k-1];
441 }
442 }
443 return a;
444}
445
446vec chirp(const vec &a, double factor)
447{
448 vec temp(a.length());
449 int i;
450 double f = factor;
451
452 it_error_if(a[0] != 1, "chirp : a[0] should be 1");
453 temp[0] = a[0];
454 for (i = 1;i < a.length();i++) {
455 temp[i] = a[i] * f;
456 f *= factor;
457 }
458 return temp;
459}
460
461vec schurrc(const vec &R, int order)
462{
463 if (order == -1) order = R.length() - 1;
464
465 vec k(order), scratch(2*order + 2);
466
467 int m;
468 int h;
469 double ex;
470 double *ep;
471 double *en;
472
473 ep = scratch._data();
474 en = scratch._data() + order + 1;
475
476 m = 0;
477 while (m < order) {
478 m++;
479 ep[m] = R[m];
480 en[m] = R[m-1];
481 }
482 if (en[1] < 1.0) en[1] = 1.0;
483 h = -1;
484 while (h < order) {
485 h++;
486 k[h] = -ep[h+1] / en[1];
487 en[1] = en[1] + k[h] * ep[h+1];
488 if (h == (order - 1)) {
489 // cout << "k: " << k << endl ;
490 return k;
491 }
492 ep[order] = ep[order] + k[h] * en[order-h];
493 m = h + 1;
494 while (m < (order - 1)) {
495 m++;
496 ex = ep[m] + k[h] * en[m-h];
497 en[m-h] = en[m-h] + k[h] * ep[m];
498 ep[m] = ex;
499 }
500 }
501 return k; // can never come here
502}
503
504vec lerouxguegenrc(const vec &R, int order)
505{
506 vec k(order);
507
508 double *r, *rny;
509 int j, m;
510 int M = order;
511
512 r = new double[2*M+1];
513 rny = new double[2*M+1];
514
515 for (j = 0;j <= M;j++) {
516 r[M-j] = r[M+j] = R[j];
517 }
518 for (m = 1;m <= M;m++) {
519 k[m-1] = -r[M+m] / r[M];
520 for (j = -M;j <= M;j++) {
521 rny[M+j] = r[M+j] + k[m-1] * r[M+m-j];
522 }
523 for (j = -M;j <= M;j++) {
524 r[M+j] = rny[M+j];
525 }
526 }
527 delete r;
528 delete rny;
529 return k;
530}
531
532double sd(const vec &In1, const vec &In2)
533{
534 return std::sqrt(37.722339402*energy(poly2cepstrum(In1, 32) - poly2cepstrum(In2, 32)));
535}
536
537// highestfreq=1 gives entire band
538double sd(const vec &In1, const vec &In2, double highestfreq)
539{
540 vec Diff = sqr(abs(log10(filter_spectrum(In1, In2))));
541 double S = 0;
542 for (int i = 0;i < round(highestfreq*129);i++) {
543 S = S + Diff(i);
544 }
545 S = S * 100 / round(highestfreq * 129);
546 return std::sqrt(S);
547}
548
549} // namespace itpp
550
#define it_error_if(t, s)
Abort if t is true.
Definition itassert.h:117
vec log10(const vec &x)
log-10 of the elements
Definition log_exp.h:271
ITPP_EXPORT vec autocorr(const vec &x, int order)
Computes the autocorrelation function.
ITPP_EXPORT vec chirp(const vec &a, double factor)
Returns a chirped version of the input vector.
ITPP_EXPORT vec ac2rc(const vec &ac)
ac2rc - Autocorrelation sequence to reflection coefficients conversion.
ITPP_EXPORT vec rc2lar(const vec &rc)
rc2lar - Reflection coefficients to log area ratios conversion.
ITPP_EXPORT vec lar2rc(const vec &lar)
lar2rc - Log area ratios to reflection coefficients conversion.
ITPP_EXPORT vec schurrc(const vec &R, int order)
schurrc - Schur algorithm.
ITPP_EXPORT vec poly2lsf(const vec &poly)
poly2lsf - Prediction polynomial to line spectral frequencies conversion.
ITPP_EXPORT vec poly2cepstrum(const vec &a)
poly2cepstrum - Prediction polynomial to cepstrum conversion.
ITPP_EXPORT vec rc2poly(const vec &rc)
rc2poly - Reflection coefficients to prediction polynomial conversion.
ITPP_EXPORT vec ac2poly(const vec &ac)
ac2poly - Autocorrelation sequence to prediction polynomial conversion.
ITPP_EXPORT vec rc2is(const vec &rc)
rc2is - Reflection coefficients to inverse sine parameters conversion.
ITPP_EXPORT vec lpc(const vec &x, int order)
lpc - Linear Predictive Coefficients using autocorrelation method.
ITPP_EXPORT vec levinson(const vec &R2, int order)
Levinson - Levinson-Durbin recursion.
ITPP_EXPORT vec lerouxguegenrc(const vec &R, int order)
Computes reflection coefficients from autocorrelation, using the Le-Roux-Guegen algorithm.
ITPP_EXPORT vec poly2rc(const vec &poly)
poly2rc - Prediction polynomial to reflection coefficients conversion.
ITPP_EXPORT double sd(const vec &In1, const vec &In2)
Spectral distortion between two vectors, in dB.
ITPP_EXPORT vec poly2ac(const vec &poly)
poly2ac - Prediction polynomial to autocorrelation sequence conversion.
ITPP_EXPORT vec is2rc(const vec &is)
is2rc - Inverse sine parameters to reflection coefficients conversion.
ITPP_EXPORT vec cepstrum2poly(const vec &c)
cepstrum2poly - Cepstrum to prediction polynomial conversion.
ITPP_EXPORT vec rc2ac(const vec &rc)
rc2ac - Reflection coefficients to autocorrelation sequence conversion.
ITPP_EXPORT vec lsf2poly(const vec &lsf)
lsf2poly - Line spectral frequencies to prediction polynomial conversion.
T sum(const Vec< T > &v)
Sum of all elements in the vector.
Definition matfunc.h:59
vec sqr(const cvec &data)
Absolute square of elements.
Definition elem_math.cpp:36
void poly(const vec &r, vec &p)
Create a polynomial of the given roots.
Definition poly.cpp:40
vec filter_spectrum(const vec &a, int nfft)
Power spectrum calculation of a filter.
Definition sigfun.cpp:325
double energy(const Vec< T > &v)
Calculate the energy: squared 2-norm. energy(v)=sum(abs(v).^2)
Definition misc_stat.h:253
Implementation of linear prediction functions, and conversion between common representations of linea...
Various functions on vectors and matrices - header file.
Miscellaneous statistics functions and classes - header file.
itpp namespace
Definition itmex.h:37
const double pi
Constant Pi.
Definition misc.h:103
ITPP_EXPORT bool any(const bvec &testvec)
Returns true if any element is one and false otherwise.
ITPP_EXPORT double round(double x)
Round to nearest integer, return result in double.
bin abs(const bin &inbin)
absolute value of bin
Definition binary.h:174
Definitions of signal processing functions.