IT++ Logo Newcom Logo

lpcfunc.cpp

Go to the documentation of this file.
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
SourceForge Logo

Generated on Thu Apr 19 14:43:45 2007 for IT++ by Doxygen 1.5.1