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