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

Generated on Sat Jul 9 2011 15:21:33 for IT++ by Doxygen 1.7.4