IT++ Logo
lpcfunc.cpp
Go to the documentation of this file.
1 
30 #include <itpp/srccode/lpcfunc.h>
31 #include <itpp/base/matfunc.h>
32 #include <itpp/signal/sigfun.h>
33 #include <itpp/stat/misc_stat.h>
34 #include <iostream>
35 
37 
38 using std::cout;
39 using std::endl;
40 
41 namespace itpp
42 {
43 
44 // Autocorrelation sequence to reflection coefficients conversion.
45 vec ac2rc(const vec &ac);
46 // Autocorrelation sequence to prediction polynomial conversion.
47 vec ac2poly(const vec &ac);
48 // Inverse sine parameters to reflection coefficients conversion.
49 vec is2rc(const vec &is);
50 // Reflection coefficients to autocorrelation sequence conversion.
51 vec rc2ac(const vec &rc);
52 // Reflection coefficients to inverse sine parameters conversion.
53 vec rc2is(const vec &rc);
54 
55 vec 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 
73 vec 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 
123 vec lpc(const vec &x, int order)
124 {
125  return levinson(autocorr(x, order), order);
126 }
127 
128 vec 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 
159 vec 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 
180 vec 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 
198 vec 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 
209 vec 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 
220 double 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 
233 double 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 
246 vec 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 
363 vec 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 
398 vec 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 
411 vec 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 
432 vec 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 
446 vec 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 
461 vec 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 
504 vec 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 
532 double 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
538 double 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 
SourceForge Logo

Generated on Sat Jul 6 2013 10:54:25 for IT++ by Doxygen 1.8.2