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
Generated on Wed Feb 9 2011 13:47:17 for IT++ by Doxygen 1.7.3