IT++ Logo

lpcfunc.cpp

Go to the documentation of this file.
00001 
00031 #include <itpp/srccode/lpcfunc.h>
00032 #include <itpp/base/matfunc.h>
00033 #include <itpp/signal/sigfun.h>
00034 #include <itpp/stat/misc_stat.h>
00035 #include <iostream>
00036 
00038 
00039 using std::cout;
00040 using std::endl;
00041 
00042 namespace itpp
00043 {
00044 
00045 // Autocorrelation sequence to reflection coefficients conversion.
00046 vec ac2rc(const vec &ac);
00047 // Autocorrelation sequence to prediction polynomial conversion.
00048 vec ac2poly(const vec &ac);
00049 // Inverse sine parameters to reflection coefficients conversion.
00050 vec is2rc(const vec &is);
00051 // Reflection coefficients to autocorrelation sequence conversion.
00052 vec rc2ac(const vec &rc);
00053 // Reflection coefficients to inverse sine parameters conversion.
00054 vec rc2is(const vec &rc);
00055 
00056 vec autocorr(const vec &x, int order)
00057 {
00058   if (order < 0) order = x.size();
00059 
00060   vec R(order + 1);
00061   double sum;
00062   int i, j;
00063 
00064   for (i = 0;i < order + 1;i++) {
00065     sum = 0;
00066     for (j = 0;j < x.size() - i;j++) {
00067       sum += x[j] * x[j+i];
00068     }
00069     R[i] = sum;
00070   }
00071   return R;
00072 }
00073 
00074 vec levinson(const vec &R2, int order)
00075 {
00076   vec R = R2;
00077   R[0] = R[0] * (1. + 1.e-9);
00078 
00079   if (order < 0) order = R.length() - 1;
00080   double k, alfa, s;
00081   double *any = new double[order+1];
00082   double *a = new double[order+1];
00083   int j, m;
00084   vec out(order + 1);
00085 
00086   a[0] = 1;
00087   alfa = R[0];
00088   if (alfa <= 0) {
00089     out.clear();
00090     out[0] = 1;
00091     return out;
00092   }
00093   for (m = 1;m <= order;m++) {
00094     s = 0;
00095     for (j = 1;j < m;j++) {
00096       s = s + a[j] * R[m-j];
00097     }
00098 
00099     k = -(R[m] + s) / alfa;
00100     if (fabs(k) >= 1.0) {
00101       cout << "levinson : panic! abs(k)>=1, order " << m << ". Aborting..." << endl ;
00102       for (j = m;j <= order;j++) {
00103         a[j] = 0;
00104       }
00105       break;
00106     }
00107     for (j = 1;j < m;j++) {
00108       any[j] = a[j] + k * a[m-j];
00109     }
00110     for (j = 1;j < m;j++) {
00111       a[j] = any[j];
00112     }
00113     a[m] = k;
00114     alfa = alfa * (1 - k * k);
00115   }
00116   for (j = 0;j < out.length();j++) {
00117     out[j] = a[j];
00118   }
00119   delete any;
00120   delete a;
00121   return out;
00122 }
00123 
00124 vec lpc(const vec &x, int order)
00125 {
00126   return levinson(autocorr(x, order), order);
00127 }
00128 
00129 vec poly2ac(const vec &poly)
00130 {
00131   vec  a = poly;
00132   int order = a.length() - 1;
00133   double alfa, s, *any = new double[order+1];
00134   int j, m;
00135   vec  r(order + 1);
00136   vec  k = poly2rc(a);
00137 
00138   it_error_if(a[0] != 1, "poly2ac : not an lpc filter");
00139   r[0] = 1;
00140   alfa = 1;
00141   for (m = 1;m <= order;m++) {
00142     s = 0;
00143     for (j = 1;j < m;j++) {
00144       s = s + a[j] * r[m-j];
00145     }
00146     r[m] = -s - alfa * k[m-1];
00147     for (j = 1;j < m;j++) {
00148       any[j] = a[j] + k[m-1] * a[m-j];
00149     }
00150     for (j = 1;j < m;j++) {
00151       a[j] = any[j];
00152     }
00153     a[m] = k[m-1];
00154     alfa = alfa * (1 - sqr(k[m-1]));
00155   }
00156   delete any;
00157   return r;
00158 }
00159 
00160 vec poly2rc(const vec &a)
00161 {
00162   // a is [1 xx xx xx], a.size()=order+1
00163   int   m, i;
00164   int    order = a.size() - 1;
00165   vec k(order);
00166   vec any(order + 1), aold(a);
00167 
00168   for (m = order - 1;m > 0;m--) {
00169     k[m] = aold[m+1] ;
00170     if (fabs(k[m]) > 1) k[m] = 1.0 / k[m];
00171     for (i = 0;i < m;i++) {
00172       any[i+1] = (aold[i+1] - aold[m-i] * k[m]) / (1 - k[m] * k[m]);
00173     }
00174     aold = any;
00175   }
00176   k[0] = any[1];
00177   if (fabs(k[0]) > 1) k[0] = 1.0 / k[0];
00178   return k;
00179 }
00180 
00181 vec rc2poly(const vec &k)
00182 {
00183   int  m, i;
00184   vec a(k.length() + 1), any(k.length() + 1);
00185 
00186   a[0] = 1;
00187   any[0] = 1;
00188   a[1] = k[0];
00189   for (m = 1;m < k.size();m++) {
00190     any[m+1] = k[m];
00191     for (i = 0;i < m;i++) {
00192       any[i+1] = a[i+1] + a[m-i] * k[m];
00193     }
00194     a = any;
00195   }
00196   return a;
00197 }
00198 
00199 vec rc2lar(const vec &k)
00200 {
00201   short m;
00202   vec LAR(k.size());
00203 
00204   for (m = 0;m < k.size();m++) {
00205     LAR[m] = std::log((1 + k[m]) / (1 - k[m]));
00206   }
00207   return LAR;
00208 }
00209 
00210 vec lar2rc(const vec &LAR)
00211 {
00212   short m;
00213   vec k(LAR.size());
00214 
00215   for (m = 0;m < LAR.size();m++) {
00216     k[m] = (std::exp(LAR[m]) - 1) / (std::exp(LAR[m]) + 1);
00217   }
00218   return k;
00219 }
00220 
00221 double FNevChebP_double(double  x, const double c[], int n)
00222 {
00223   int i;
00224   double b0 = 0.0, b1 = 0.0, b2 = 0.0;
00225 
00226   for (i = n - 1; i >= 0; --i) {
00227     b2 = b1;
00228     b1 = b0;
00229     b0 = 2.0 * x * b1 - b2 + c[i];
00230   }
00231   return (0.5 * (b0 - b2 + c[0]));
00232 }
00233 
00234 double FNevChebP(double  x, const double c[], int n)
00235 {
00236   int i;
00237   double b0 = 0.0, b1 = 0.0, b2 = 0.0;
00238 
00239   for (i = n - 1; i >= 0; --i) {
00240     b2 = b1;
00241     b1 = b0;
00242     b0 = 2.0 * x * b1 - b2 + c[i];
00243   }
00244   return (0.5 * (b0 - b2 + c[0]));
00245 }
00246 
00247 vec poly2lsf(const vec &pc)
00248 {
00249   int np = pc.length() - 1;
00250   vec lsf(np);
00251 
00252   vec fa((np + 1) / 2 + 1), fb((np + 1) / 2 + 1);
00253   vec ta((np + 1) / 2 + 1), tb((np + 1) / 2 + 1);
00254   double *t;
00255   double xlow, xmid, xhigh;
00256   double ylow, ymid, yhigh;
00257   double xroot;
00258   double dx;
00259   int i, j, nf;
00260   int odd;
00261   int na, nb, n;
00262   double ss, aa;
00263   double DW = (0.02 * pi);
00264   int  NBIS = 4;
00265 
00266   odd = (np % 2 != 0);
00267   if (odd) {
00268     nb = (np + 1) / 2;
00269     na = nb + 1;
00270   }
00271   else {
00272     nb = np / 2 + 1;
00273     na = nb;
00274   }
00275 
00276   fa[0] = 1.0;
00277   for (i = 1, j = np; i < na; ++i, --j)
00278     fa[i] = pc[i] + pc[j];
00279 
00280   fb[0] = 1.0;
00281   for (i = 1, j = np; i < nb; ++i, --j)
00282     fb[i] = pc[i] - pc[j];
00283 
00284   if (odd) {
00285     for (i = 2; i < nb; ++i)
00286       fb[i] = fb[i] + fb[i-2];
00287   }
00288   else {
00289     for (i = 1; i < na; ++i) {
00290       fa[i] = fa[i] - fa[i-1];
00291       fb[i] = fb[i] + fb[i-1];
00292     }
00293   }
00294 
00295   ta[0] = fa[na-1];
00296   for (i = 1, j = na - 2; i < na; ++i, --j)
00297     ta[i] = 2.0 * fa[j];
00298 
00299   tb[0] = fb[nb-1];
00300   for (i = 1, j = nb - 2; i < nb; ++i, --j)
00301     tb[i] = 2.0 * fb[j];
00302 
00303   nf = 0;
00304   t = ta._data();
00305   n = na;
00306   xroot = 2.0;
00307   xlow = 1.0;
00308   ylow = FNevChebP_double(xlow, t, n);
00309 
00310 
00311   ss = std::sin(DW);
00312   aa = 4.0 - 4.0 * std::cos(DW)  - ss;
00313   while (xlow > -1.0 && nf < np) {
00314     xhigh = xlow;
00315     yhigh = ylow;
00316     dx = aa * xhigh * xhigh + ss;
00317     xlow = xhigh - dx;
00318     if (xlow < -1.0)
00319       xlow = -1.0;
00320     ylow = FNevChebP_double(xlow, t, n);
00321     if (ylow * yhigh <= 0.0) {
00322       dx = xhigh - xlow;
00323       for (i = 1; i <= NBIS; ++i) {
00324         dx = 0.5 * dx;
00325         xmid = xlow + dx;
00326         ymid = FNevChebP_double(xmid, t, n);
00327         if (ylow * ymid <= 0.0) {
00328           yhigh = ymid;
00329           xhigh = xmid;
00330         }
00331         else {
00332           ylow = ymid;
00333           xlow = xmid;
00334         }
00335       }
00336       if (yhigh != ylow)
00337         xmid = xlow + dx * ylow / (ylow - yhigh);
00338       else
00339         xmid = xlow + dx;
00340       lsf[nf] = std::acos((double) xmid);
00341       ++nf;
00342       if (xmid >= xroot) {
00343         xmid = xlow - dx;
00344       }
00345       xroot = xmid;
00346       if (t == ta._data()) {
00347         t = tb._data();
00348         n = nb;
00349       }
00350       else {
00351         t = ta._data();
00352         n = na;
00353       }
00354       xlow = xmid;
00355       ylow = FNevChebP_double(xlow, t, n);
00356     }
00357   }
00358   if (nf != np) {
00359     cout << "poly2lsf: WARNING: failed to find all lsfs" << endl ;
00360   }
00361   return lsf;
00362 }
00363 
00364 vec lsf2poly(const vec &f)
00365 {
00366   int m = f.length();
00367   vec  pc(m + 1);
00368   double c1, c2, *a;
00369   vec  p(m + 1), q(m + 1);
00370   int mq, n, i, nor;
00371 
00372   it_error_if(m % 2 != 0, "lsf2poly: THIS ROUTINE WORKS ONLY FOR EVEN m");
00373   pc[0] = 1.0;
00374   a = pc._data() + 1;
00375   mq = m >> 1;
00376   for (i = 0 ; i <= m ; i++) {
00377     q[i] = 0.;
00378     p[i] = 0.;
00379   }
00380   p[0] = q[0] = 1.;
00381   for (n = 1; n <= mq; n++) {
00382     nor = 2 * n;
00383     c1 = 2 * std::cos(f[nor-1]);
00384     c2 = 2 * std::cos(f[nor-2]);
00385     for (i = nor; i >= 2; i--) {
00386       q[i] += q[i-2] - c1 * q[i-1];
00387       p[i] += p[i-2] - c2 * p[i-1];
00388     }
00389     q[1] -= c1;
00390     p[1] -= c2;
00391   }
00392   a[0] = 0.5 * (p[1] + q[1]);
00393   for (i = 1, n = 2; i < m ; i++, n++)
00394     a[i] = 0.5 * (p[i] + p[n] + q[n] - q[i]);
00395 
00396   return pc;
00397 }
00398 
00399 vec poly2cepstrum(const vec &a)
00400 {
00401   vec c(a.length() - 1);
00402 
00403   for (int n = 1;n <= c.length();n++) {
00404     c[n-1] = a[n];
00405     for (int k = 1;k < n;k++) {
00406       c[n-1] -= double(k) / n * a[n-k] * c[k-1];
00407     }
00408   }
00409   return c;
00410 }
00411 
00412 vec poly2cepstrum(const vec &a, int num)
00413 {
00414   it_error_if(num < a.length(), "a2cepstrum : not allowed cepstrum length");
00415   vec c(num);
00416   int n;
00417 
00418   for (n = 1;n < a.length();n++) {
00419     c[n-1] = a[n];
00420     for (int k = 1;k < n;k++) {
00421       c[n-1] -= double(k) / n * a[n-k] * c[k-1];
00422     }
00423   }
00424   for (n = a.length();n <= c.length();n++) {
00425     c[n-1] = 0;
00426     for (int k = n - a.length() + 1;k < n;k++) {
00427       c[n-1] -= double(k) / n * a[n-k] * c[k-1];
00428     }
00429   }
00430   return c;
00431 }
00432 
00433 vec cepstrum2poly(const vec &c)
00434 {
00435   vec a(c.length() + 1);
00436 
00437   a[0] = 1;
00438   for (int n = 1;n <= c.length();n++) {
00439     a[n] = c[n-1];
00440     for (int k = 1;k < n;k++) {
00441       a[n] += double(k) / n * a[n-k] * c[k-1];
00442     }
00443   }
00444   return a;
00445 }
00446 
00447 vec chirp(const vec &a, double factor)
00448 {
00449   vec    temp(a.length());
00450   int    i;
00451   double   f = factor;
00452 
00453   it_error_if(a[0] != 1, "chirp : a[0] should be 1");
00454   temp[0] = a[0];
00455   for (i = 1;i < a.length();i++) {
00456     temp[i] = a[i] * f;
00457     f *= factor;
00458   }
00459   return temp;
00460 }
00461 
00462 vec schurrc(const vec &R, int order)
00463 {
00464   if (order == -1) order = R.length() - 1;
00465 
00466   vec    k(order), scratch(2*order + 2);
00467 
00468   int m;
00469   int h;
00470   double ex;
00471   double *ep;
00472   double *en;
00473 
00474   ep = scratch._data();
00475   en = scratch._data() + order + 1;
00476 
00477   m = 0;
00478   while (m < order) {
00479     m++;
00480     ep[m] = R[m];
00481     en[m] = R[m-1];
00482   }
00483   if (en[1] < 1.0) en[1] = 1.0;
00484   h = -1;
00485   while (h < order) {
00486     h++;
00487     k[h] = -ep[h+1] / en[1];
00488     en[1] = en[1] + k[h] * ep[h+1];
00489     if (h == (order - 1)) {
00490       // cout << "k: " << k << endl ;
00491       return k;
00492     }
00493     ep[order] = ep[order] + k[h] * en[order-h];
00494     m = h + 1;
00495     while (m < (order - 1)) {
00496       m++;
00497       ex = ep[m] + k[h] * en[m-h];
00498       en[m-h] = en[m-h] + k[h] * ep[m];
00499       ep[m] = ex;
00500     }
00501   }
00502   return k;  // can never come here
00503 }
00504 
00505 vec lerouxguegenrc(const vec &R, int order)
00506 {
00507   vec    k(order);
00508 
00509   double  *r, *rny;
00510   int j, m;
00511   int M = order;
00512 
00513   r = new double[2*M+1];
00514   rny = new double[2*M+1];
00515 
00516   for (j = 0;j <= M;j++) {
00517     r[M-j] = r[M+j] = R[j];
00518   }
00519   for (m = 1;m <= M;m++) {
00520     k[m-1] = -r[M+m] / r[M];
00521     for (j = -M;j <= M;j++) {
00522       rny[M+j] = r[M+j] + k[m-1] * r[M+m-j];
00523     }
00524     for (j = -M;j <= M;j++) {
00525       r[M+j] = rny[M+j];
00526     }
00527   }
00528   delete r;
00529   delete rny;
00530   return k;
00531 }
00532 
00533 double sd(const vec &In1, const vec &In2)
00534 {
00535   return std::sqrt(37.722339402*energy(poly2cepstrum(In1, 32) - poly2cepstrum(In2, 32)));
00536 }
00537 
00538 // highestfreq=1 gives entire band
00539 double sd(const vec &In1, const vec &In2, double highestfreq)
00540 {
00541   vec Diff = sqr(abs(log10(filter_spectrum(In1, In2))));
00542   double S = 0;
00543   for (int i = 0;i < round(highestfreq*129);i++) {
00544     S = S + Diff(i);
00545   }
00546   S = S * 100 / round(highestfreq * 129);
00547   return std::sqrt(S);
00548 }
00549 
00550 } // namespace itpp
00551 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Wed Feb 9 2011 13:47:25 for IT++ by Doxygen 1.7.3