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
Generated on Sat Jul 9 2011 15:21:33 for IT++ by Doxygen 1.7.4