00001 00034 #include <itpp/srccode/lpcfunc.h> 00035 #include <itpp/base/stat.h> 00036 #include <itpp/base/sigfun.h> 00037 #include <iostream> 00038 00039 00040 using std::cout; 00041 using std::endl; 00042 00043 namespace itpp { 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; 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;any[0]=1; 00186 a[1]=k[0]; 00187 for (m=1;m<k.size();m++) { 00188 any[m+1]=k[m]; 00189 for (i=0;i<m;i++) { 00190 any[i+1]=a[i+1]+a[m-i]*k[m]; 00191 } 00192 a=any; 00193 } 00194 return a; 00195 } 00196 00197 vec rc2lar(const vec &k) 00198 { 00199 short m; 00200 vec LAR(k.size()); 00201 00202 for (m=0;m<k.size();m++) { 00203 LAR[m]=std::log((1+k[m])/(1-k[m])); 00204 } 00205 return LAR; 00206 } 00207 00208 vec lar2rc(const vec &LAR) 00209 { 00210 short m; 00211 vec k(LAR.size()); 00212 00213 for (m=0;m<LAR.size();m++) { 00214 k[m]=(std::exp(LAR[m])-1)/(std::exp(LAR[m])+1); 00215 } 00216 return k; 00217 } 00218 00219 double FNevChebP_double(double x,const double c[],int n) 00220 { 00221 int i; 00222 double b0=0.0, b1=0.0, b2=0.0; 00223 00224 for (i = n - 1; i >= 0; --i) { 00225 b2 = b1; 00226 b1 = b0; 00227 b0 = 2.0 * x * b1 - b2 + c[i]; 00228 } 00229 return (0.5 * (b0 - b2 + c[0])); 00230 } 00231 00232 double FNevChebP(double x,const double c[],int n) 00233 { 00234 int i; 00235 double b0=0.0, b1=0.0, b2=0.0; 00236 00237 for (i = n - 1; i >= 0; --i) { 00238 b2 = b1; 00239 b1 = b0; 00240 b0 = 2.0 * x * b1 - b2 + c[i]; 00241 } 00242 return (0.5 * (b0 - b2 + c[0])); 00243 } 00244 00245 vec poly2lsf(const vec &pc) 00246 { 00247 int np=pc.length()-1; 00248 vec lsf(np); 00249 00250 vec fa((np+1)/2+1), fb((np+1)/2+1); 00251 vec ta((np+1)/2+1), tb((np+1)/2+1); 00252 double *t; 00253 double xlow, xmid, xhigh; 00254 double ylow, ymid, yhigh; 00255 double xroot; 00256 double dx; 00257 int i, j, nf; 00258 int odd; 00259 int na, nb, n; 00260 double ss, aa; 00261 double DW=(0.02 * pi); 00262 int NBIS=4; 00263 00264 odd = (np % 2 != 0); 00265 if (odd) { 00266 nb = (np + 1) / 2; 00267 na = nb + 1; 00268 } 00269 else { 00270 nb = np / 2 + 1; 00271 na = nb; 00272 } 00273 00274 fa[0] = 1.0; 00275 for (i = 1, j = np; i < na; ++i, --j) 00276 fa[i] = pc[i] + pc[j]; 00277 00278 fb[0] = 1.0; 00279 for (i = 1, j = np; i < nb; ++i, --j) 00280 fb[i] = pc[i] - pc[j]; 00281 00282 if (odd) { 00283 for (i = 2; i < nb; ++i) 00284 fb[i] = fb[i] + fb[i-2]; 00285 } 00286 else { 00287 for (i = 1; i < na; ++i) { 00288 fa[i] = fa[i] - fa[i-1]; 00289 fb[i] = fb[i] + fb[i-1]; 00290 } 00291 } 00292 00293 ta[0] = fa[na-1]; 00294 for (i = 1, j = na - 2; i < na; ++i, --j) 00295 ta[i] = 2.0 * fa[j]; 00296 00297 tb[0] = fb[nb-1]; 00298 for (i = 1, j = nb - 2; i < nb; ++i, --j) 00299 tb[i] = 2.0 * fb[j]; 00300 00301 nf = 0; 00302 t = ta._data(); 00303 n = na; 00304 xroot = 2.0; 00305 xlow = 1.0; 00306 ylow = FNevChebP_double(xlow, t, n); 00307 00308 00309 ss = std::sin (DW); 00310 aa = 4.0 - 4.0 * std::cos (DW) - ss; 00311 while (xlow > -1.0 && nf < np) { 00312 xhigh = xlow; 00313 yhigh = ylow; 00314 dx = aa * xhigh * xhigh + ss; 00315 xlow = xhigh - dx; 00316 if (xlow < -1.0) 00317 xlow = -1.0; 00318 ylow = FNevChebP_double(xlow, t, n); 00319 if (ylow * yhigh <= 0.0) { 00320 dx = xhigh - xlow; 00321 for (i = 1; i <= NBIS; ++i) { 00322 dx = 0.5 * dx; 00323 xmid = xlow + dx; 00324 ymid = FNevChebP_double(xmid, t, n); 00325 if (ylow * ymid <= 0.0) { 00326 yhigh = ymid; 00327 xhigh = xmid; 00328 } 00329 else { 00330 ylow = ymid; 00331 xlow = xmid; 00332 } 00333 } 00334 if (yhigh != ylow) 00335 xmid = xlow + dx * ylow / (ylow - yhigh); 00336 else 00337 xmid = xlow + dx; 00338 lsf[nf] = std::acos((double) xmid); 00339 ++nf; 00340 if (xmid >= xroot) { 00341 xmid = xlow - dx; 00342 } 00343 xroot = xmid; 00344 if (t == ta._data()) { 00345 t = tb._data(); 00346 n = nb; 00347 } 00348 else { 00349 t = ta._data(); 00350 n = na; 00351 } 00352 xlow = xmid; 00353 ylow = FNevChebP_double(xlow, t, n); 00354 } 00355 } 00356 if (nf != np) { 00357 cout << "poly2lsf: WARNING: failed to find all lsfs" << endl ; 00358 } 00359 return lsf; 00360 } 00361 00362 vec lsf2poly(const vec &f) 00363 { 00364 int m=f.length(); 00365 vec pc(m+1); 00366 double c1, c2, *a; 00367 vec p(m+1), q(m+1); 00368 int mq, n, i, nor; 00369 00370 it_error_if(m%2!=0,"lsf2poly: THIS ROUTINE WORKS ONLY FOR EVEN m"); 00371 pc[0] = 1.0; 00372 a = pc._data() + 1; 00373 mq=m>>1; 00374 for(i=0 ; i <= m ; i++) { 00375 q[i]=0.; 00376 p[i]=0.; 00377 } 00378 p[0] = q[0] = 1.; 00379 for(n=1; n <= mq; n++) { 00380 nor=2*n; 00381 c1 = 2*std::cos(f[nor-1]); 00382 c2 = 2*std::cos(f[nor-2]); 00383 for(i=nor; i >= 2; i--) { 00384 q[i] += q[i-2] - c1*q[i-1]; 00385 p[i] += p[i-2] - c2*p[i-1]; 00386 } 00387 q[1] -= c1; 00388 p[1] -= c2; 00389 } 00390 a[0] = 0.5 * (p[1] + q[1]); 00391 for(i=1, n=2; i < m ; i++, n++) 00392 a[i] = 0.5 * (p[i] + p[n] + q[n] - q[i]); 00393 00394 return pc; 00395 } 00396 00397 vec poly2cepstrum(const vec &a) 00398 { 00399 vec c(a.length()-1); 00400 00401 for (int n=1;n<=c.length();n++) { 00402 c[n-1]=a[n]; 00403 for (int k=1;k<n;k++) { 00404 c[n-1]-=double(k)/n*a[n-k]*c[k-1]; 00405 } 00406 } 00407 return c; 00408 } 00409 00410 vec poly2cepstrum(const vec &a, int num) 00411 { 00412 it_error_if(num<a.length(),"a2cepstrum : not allowed cepstrum length"); 00413 vec c(num); 00414 int n; 00415 00416 for (n=1;n<a.length();n++) { 00417 c[n-1]=a[n]; 00418 for (int k=1;k<n;k++) { 00419 c[n-1]-=double(k)/n*a[n-k]*c[k-1]; 00420 } 00421 } 00422 for (n=a.length();n<=c.length();n++) { 00423 c[n-1]=0; 00424 for (int k=n-a.length()+1;k<n;k++) { 00425 c[n-1]-=double(k)/n*a[n-k]*c[k-1]; 00426 } 00427 } 00428 return c; 00429 } 00430 00431 vec cepstrum2poly(const vec &c) 00432 { 00433 vec a(c.length()+1); 00434 00435 a[0]=1; 00436 for (int n=1;n<=c.length();n++) { 00437 a[n]=c[n-1]; 00438 for (int k=1;k<n;k++) { 00439 a[n]+=double(k)/n*a[n-k]*c[k-1]; 00440 } 00441 } 00442 return a; 00443 } 00444 00445 vec chirp(const vec &a, double factor) 00446 { 00447 vec temp(a.length()); 00448 int i; 00449 double f=factor; 00450 00451 it_error_if(a[0]!=1,"chirp : a[0] should be 1"); 00452 temp[0]=a[0]; 00453 for (i=1;i<a.length();i++) { 00454 temp[i]=a[i]*f; 00455 f*=factor; 00456 } 00457 return temp; 00458 } 00459 00460 vec schurrc(const vec &R, int order) 00461 { 00462 if (order==-1) order=R.length()-1; 00463 00464 vec k(order), scratch(2*order+2); 00465 00466 int m; 00467 int h; 00468 double ex; 00469 double *ep; 00470 double *en; 00471 00472 ep = scratch._data(); 00473 en = scratch._data() + order + 1; 00474 00475 m = 0; 00476 while( m < order){ 00477 m++; 00478 ep[m] = R[m]; 00479 en[m] = R[m-1]; 00480 } 00481 if( en[1] < 1.0) en[1] = 1.0; 00482 h = -1; 00483 while( h < order){ 00484 h++; 00485 k[h] = -ep[h+1]/en[1]; 00486 en[1] = en[1] + k[h]*ep[h+1]; 00487 if( h == (order-1)) { 00488 // cout << "k: " << k << endl ; 00489 return k; 00490 } 00491 ep[order] = ep[order] + k[h]*en[order-h]; 00492 m = h+1; 00493 while( m < (order-1)){ 00494 m++; 00495 ex = ep[m] + k[h]*en[m-h]; 00496 en[m-h] = en[m-h] + k[h]*ep[m]; 00497 ep[m] = ex; 00498 } 00499 } 00500 return k; // can never come here 00501 } 00502 00503 vec lerouxguegenrc(const vec &R, int order) 00504 { 00505 vec k(order); 00506 00507 double *r,*rny; 00508 int j,m; 00509 int M=order; 00510 00511 r=new double[2*M+1]; 00512 rny=new double[2*M+1]; 00513 00514 for (j=0;j<=M;j++) { 00515 r[M-j]=r[M+j]=R[j]; 00516 } 00517 for (m=1;m<=M;m++) { 00518 k[m-1]=-r[M+m]/r[M]; 00519 for (j=-M;j<=M;j++) { 00520 rny[M+j]=r[M+j]+k[m-1]*r[M+m-j]; 00521 } 00522 for (j=-M;j<=M;j++) { 00523 r[M+j]=rny[M+j]; 00524 } 00525 } 00526 delete r; 00527 delete rny; 00528 return k; 00529 } 00530 00531 double sd(const vec &In1, const vec &In2) 00532 { 00533 return std::sqrt(37.722339402*energy(poly2cepstrum(In1,32)-poly2cepstrum(In2,32))); 00534 } 00535 00536 // highestfreq=1 gives entire band 00537 double sd(const vec &In1, const vec &In2, double highestfreq) 00538 { 00539 vec Diff=sqr(abs(log10(filter_spectrum(In1,In2)))); 00540 double S=0; 00541 for (int i=0;i<round(highestfreq*129);i++) { 00542 S=S+Diff(i); 00543 } 00544 S=S*100/round(highestfreq*129); 00545 return std::sqrt(S); 00546 } 00547 00548 } // namespace itpp
Generated on Thu Aug 30 02:47:20 2007 for IT++ by Doxygen 1.5.3