IT++ Logo

lpcfunc.cpp

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

Generated on Mon Jan 7 22:28:59 2008 for IT++ by Doxygen 1.5.4