Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/LPC.cpp
Go to the documentation of this file.
00001 /*
00002 ** Copyright (C) 1998-2006 George Tzanetakis <gtzan@cs.uvic.ca>
00003 **
00004 ** This program is free software; you can redistribute it and/or modify
00005 ** it under the terms of the GNU General Public License as published by
00006 ** the Free Software Foundation; either version 2 of the License, or
00007 ** (at your option) any later version.
00008 **
00009 ** This program is distributed in the hope that it will be useful,
00010 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012 ** GNU General Public License for more details.
00013 **
00014 ** You should have received a copy of the GNU General Public License
00015 ** along with this program; if not, write to the Free Software
00016 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00017 */
00018 
00019 #include "LPC.h"
00020 #include "../common_source.h"
00021 
00022 // #define _MATLAB_LPC_
00023 
00024 using std::ostringstream;
00025 using namespace Marsyas;
00026 
00027 LPC::LPC(mrs_string name):MarSystem("LPC",name)
00028 {
00029   addControls();
00030 }
00031 
00032 LPC::LPC(const LPC& a):MarSystem(a)
00033 {
00034   ctrl_coeffs_ = getctrl("mrs_realvec/coeffs");
00035   ctrl_pitch_ = getctrl("mrs_real/pitch");
00036   ctrl_power_ = getctrl("mrs_real/power");
00037 }
00038 
00039 LPC::~LPC()
00040 {
00041 }
00042 
00043 MarSystem*
00044 LPC::clone() const
00045 {
00046   return new LPC(*this);
00047 }
00048 
00049 void
00050 LPC::addControls()
00051 {
00052   addctrl("mrs_natural/order", (mrs_natural)10);
00053   addctrl("mrs_realvec/coeffs", realvec(), ctrl_coeffs_);
00054   addctrl("mrs_real/pitch", 0.0, ctrl_pitch_);
00055   addctrl("mrs_real/power", 0.0, ctrl_power_);
00056   setctrlState("mrs_natural/order", true);
00057   addctrl("mrs_real/lambda", (mrs_real)0.0);
00058   addctrl("mrs_real/gamma", (mrs_real)1.0);
00059 }
00060 
00061 void
00062 LPC::myUpdate(MarControlPtr sender)
00063 {
00064   (void) sender;  //suppress warning of unused parameter(s)
00065   MRSDIAG("LPC.cpp - LPC:myUpdate");
00066 
00067   order_ = getctrl("mrs_natural/order")->to<mrs_natural>();
00068 
00069   setctrl("mrs_natural/onObservations", (mrs_natural)(order_+2)); // <order_> coefs + pitch value + power value
00070   setctrl("mrs_natural/onSamples", (mrs_natural)1);
00071   setctrl("mrs_real/osrate", getctrl("mrs_real/israte"));
00072 
00073   //LPC features names
00074   ostringstream oss;
00075   for (mrs_natural i = 0; i < order_; ++i)
00076     oss << "LPC_" << i+1 << ",";
00077   oss << "LPC_Pitch," << "LPC_Gain,";
00078   setctrl("mrs_string/onObsNames", oss.str());
00079 
00080   temp_.create(order_ ,order_);
00081   Zs_.create(order_);
00082 
00083   {
00084     MarControlAccessor acc(ctrl_coeffs_, NOUPDATE);
00085     realvec& coeffs = acc.to<mrs_realvec>();
00086     coeffs.stretch(order_+1);
00087   }
00088 
00089   //MATLAB engine stuff - for testing and validation purposes only!
00090 #ifdef _MATLAB_LPC_
00091   MATLAB_EVAL("LPC_pitch = [];");
00092 #endif
00093 
00094 }
00095 
00096 //perhaps this method could become an independent MarSystem in the future...
00097 void
00098 LPC::autocorrelationWarped(const realvec& in, realvec& r, mrs_real& pitch, mrs_real lambda)
00099 {
00100   //*Based on the code from: http://www.musicdsp.org/showone.php?id=137
00101 
00102   //find the order-P autocorrelation array, R, for the sequence x
00103   //of length L and using a warping factor of lambda
00104 
00105   mrs_real* x = in.getData();
00106   mrs_natural L = in.getSize();
00107   mrs_real* R = r.getData();
00108   mrs_natural P = in.getSize()/2;//order_;
00109 
00110   mrs_real* dl = new mrs_real[L];
00111   mrs_real* Rt = new mrs_real[L];
00112   mrs_real r1,r2,r1t;
00113   R[0]=0;
00114   Rt[0]=0;
00115   r1=0;
00116   r2=0;
00117   r1t=0;
00118   for(mrs_natural k=0; k<L; k++)
00119   {
00120     Rt[0]+=x[k]*x[k];
00121 
00122     dl[k]=r1-lambda*(x[k]-r2);
00123     r1 = x[k];
00124     r2 = dl[k];
00125   }
00126   for(mrs_natural i=1; i<=P; ++i)
00127   {
00128     Rt[i]=0;
00129     r1=0;
00130     r2=0;
00131     for(mrs_natural k=0; k<L; k++)
00132     {
00133       Rt[i]+=dl[k]*x[k];
00134 
00135       r1t = dl[k];
00136       dl[k]=r1-lambda*(r1t-r2);
00137       r1 = r1t;
00138       r2 = dl[k];
00139     }
00140   }
00141 
00142   for(long i=0; i<=P; ++i)
00143     R[i]=Rt[i]/in.getSize(); // [ML] change /
00144 
00145   delete[] dl;
00146   delete[] Rt;
00147 
00148   //----------------------------------------------------
00149   // estimate pitch
00150   //----------------------------------------------------
00151   mrs_real temp = r(0);
00152   //set peak searching start point to 2% of total window size [?]
00153   mrs_real j = (mrs_real)in.getSize() * 0.02;
00154   //detect first local minimum...
00155   while (r((mrs_natural)j) < temp && j < in.getSize()/2)
00156   {
00157     temp = r((mrs_natural)j);
00158     j++;
00159   }
00160   //... and then from there, detect higher peak => period estimate!
00161   temp = 0.0;
00162   for (mrs_natural i=(mrs_natural)j; i < in.getSize() * 0.5; ++i)
00163   {
00164     if (r(i) > temp)
00165     {
00166       j = i;
00167       temp = r(i);
00168     }
00169   }
00170 
00171   // This code is from the original Marsyas0.1 AutoCorrelation class
00172   //this seems like some sort of Narrowed Autocorrelation (NAC)... [?][!]
00173   //however, it does not seem to fit the NAC definition...
00174   //so, what is this for??
00175   //    mrs_real norm = 1.0 / (mrs_real)in.getSize();
00176   //    mrs_natural k = in.getSize()/2;
00177   //    for (mrs_natural i=0; i < in.getSize()/2; ++i)
00178   //    r(i) *= (k-i) * norm;
00179 
00180   //if autocorr peak not very prominent => not a good period estimate
00181   //so discard it...
00182   if ((r((mrs_natural)j) / r(0)) < 0.4) j=0;
00183   //avoid detection of too high fundamental freqs (e.g. harmonics)?
00184   if (j > in.getSize()/4) j = 0;
00185 
00186   //pitch gets in fact the period (i.e. autocorr lag)
00187   //measured in time samples... maybe this should be converted
00188   //to a more convenient representation (e.g. pitch in Hz).
00189   //Such a conversion would have to depend on the warp factor lambda... [!]
00190   pitch  = (mrs_real)j;
00191 
00192   //MATLAB engine stuff - for testing and validation purposes only!
00193 #ifdef _MATLAB_LPC_
00194   MATLAB_PUT(pitch, "pitch");
00195   MATLAB_EVAL("LPC_pitch = [LPC_pitch pitch];");
00196 #endif
00197 }
00198 
00199 //Perhaps this method could become a member of realvec...
00200 void
00201 LPC::LevinsonDurbin(const realvec& r, realvec& a, realvec& kVec, mrs_real& e)
00202 {
00203   //*Based on the code from: http://www.musicdsp.org/showone.php?id=137
00204 
00205   mrs_natural P = order_;
00206   mrs_real* R = r.getData();
00207   mrs_real* A = a.getData();
00208   mrs_real* K = kVec.getData();
00209   e = 0.0; //prediction error;
00210 
00211   mrs_real Am1[62];
00212 
00213   if(R[0]==0.0)
00214   {
00215     for(mrs_natural i=1; i<=P; ++i)
00216     {
00217       K[i]=0.0;
00218       A[i]=0.0;
00219     }
00220   }
00221   else
00222   {
00223     mrs_real km,Em1,Em;
00224     mrs_natural k,s,m;
00225 
00226     for (k=0; k<=P; k++) {
00227       A[k]=0;
00228       Am1[k]=0;
00229     }
00230 
00231     A[0]=1;
00232     Am1[0]=1;
00233     km=0;
00234     Em1=R[0];
00235 
00236     for (m=1; m<=P; m++)              //m=2:N+1
00237     {
00238       mrs_real err=0.0f;                //err = 0;
00239 
00240       for (k=1; k<=m-1; k++)        //for k=2:m-1
00241         err += Am1[k]*R[m-k];     // err = err + am1(k)*R(m-k+1);
00242 
00243       km = (R[m]-err)/Em1;          //km=(R(m)-err)/Em1;
00244       K[m-1] = -km;
00245       A[m]=km;                      //am(m)=km;
00246 
00247       for (k=1; k<=m-1; k++)        //for k=2:m-1
00248         A[k]=Am1[k]-km*Am1[m-k];  // am(k)=am1(k)-km*am1(m-k+1);
00249 
00250       Em=(1-km*km)*Em1;             //Em=(1-km*km)*Em1;
00251 
00252       for(s=0; s<=P; s++)           //for s=1:N+1
00253         Am1[s] = A[s];            // am1(s) = am(s)
00254 
00255       Em1 = Em;                     //Em1 = Em;
00256       //e = Em1; //prediction error
00257       e = Em1*Em1; //RMS prediction error
00258     }
00259     e = sqrt(e/(mrs_real)P); //RMS prediction error
00260   }
00261 }
00262 
00263 mrs_real
00264 LPC::predictionError(const realvec& data, const realvec& coeffs)
00265 {
00266   mrs_natural i, j;
00267   mrs_real power = 0.0;
00268   mrs_real error, tmp;
00269 
00270   //load delay line with input data...
00271   for (i=0; i< order_; ++i)
00272   {
00273     Zs_(i) = data(order_-i-1);
00274   }
00275   //apply LPC filter and estimate RMS of the error (=~ LPC Gain)
00276   mrs_natural count = 0;
00277   for (i=order_; i<(mrs_natural)data.getSize() ; ++i)
00278   {
00279     tmp = 0.0;
00280     for (j=0; j< order_; j++)
00281       tmp += Zs_(j) * coeffs(j);
00282     for (j=order_-1; j>0; j--)
00283       Zs_(j) = Zs_(j-1);
00284 
00285     Zs_(0) = data(i);
00286     error = data(i) - tmp;
00287     power += error * error;
00288     count ++;
00289   }
00290   return sqrt(power/(mrs_real)count);//=~ RMS LPC Gain
00291 }
00292 
00293 
00294 
00295 
00296 /*
00297 
00298 Computation of the autocorrelation coefficients and the prediction error gain
00299 using the implementation of peter Kabal
00300 http://www-mmsp.ece.mcgill.ca/Documents/Software/index.html
00301 
00302 */
00303 
00304 mrs_real
00305 LPC::VRfDotProd (mrs_real * x1, mrs_real * x2, mrs_natural N)
00306 {
00307   mrs_natural i;
00308   mrs_real sum;
00309 
00310   sum = 0.0;
00311   for (i = 0; i < N; ++i)
00312     sum +=  x1[i] * x2[i];
00313 
00314   return sum;
00315 }
00316 
00317 void
00318 LPC::SPautoc (mrs_real * x, mrs_natural Nx, mrs_real * cor, mrs_natural Nt)
00319 {
00320   mrs_natural i;
00321   mrs_natural N;
00322 
00323   N = Nt;
00324   if (Nt > Nx)
00325     N = Nx;
00326 
00327   for (i = 0; i < N; ++i)
00328     cor[i] =  VRfDotProd (x, &x[i], Nx - i);
00329 
00330   for (i = N; i < Nt; ++i)
00331     cor[i] = 0.0;
00332 
00333   return;
00334 }
00335 
00336 
00337 mrs_real
00338 LPC::SPcorXpc (mrs_real * rxx, mrs_real * pc, mrs_natural Np)
00339 {
00340   mrs_natural i, j, k;
00341   mrs_real rc, tp;
00342   mrs_real perr, t, sum;
00343 
00344   perr = rxx[0];
00345 
00346 
00347   for (k = 0; k < Np; ++k) {
00348 
00349     /* Calculate the next reflection coefficient */
00350     sum = rxx[k+1];
00351     for (i = 0; i < k; ++i)
00352       sum = sum - rxx[k-i] * pc[i];
00353 
00354     /* Check for zero prediction error */
00355     if (perr == 0.0 && sum == 0.0)
00356       rc = 0.0;
00357     else
00358       rc =  (-sum / perr);
00359 
00360     /* Calculate the error (equivalent to perr = perr * (1-rc^2))
00361     A change in sign of perr means that rc (reflection coefficient for stage k)
00362     has a magnitude greater than unity (corresponding to an unstable synthesis
00363     filter)
00364     */
00365     t = perr + rc * sum;
00366 
00367     perr = t;
00368 
00369     /* Convert from reflection coefficients to predictor coefficients */
00370     pc[k] = -rc;
00371     for (i = 0, j = k - 1; i < j; ++i, --j) {
00372       tp = pc[i] + rc * pc[j];
00373       pc[j] = pc[j] + rc * pc[i];
00374       pc[i] = tp;
00375     }
00376     if (i == j)
00377       pc[i] = pc[i] + rc * pc[i];
00378   }
00379 
00380   return perr;
00381 }
00382 
00383 
00384 
00385 void
00386 LPC::myProcess(realvec& in, realvec& out)
00387 {
00388   mrs_natural i;
00389   mrs_real LevinsonError = 0.0;
00390   // FIXME This variable is defined but (possibly) never used.
00391   // mrs_real PredictionError = 0.0;
00392   mrs_real pitch = 0.0, lag = 0.0;
00393 
00394   //MATLAB engine stuff - for testing and validation purposes only!
00395 #ifdef _MATLAB_LPC_
00396   MATLAB_PUT(featureMode_, "featureMode");
00397   MATLAB_PUT(in, "LPC_in");
00398   MATLAB_PUT(order_, "LPC_order");
00399   MATLAB_PUT(getctrl("mrs_real/gamma")->to<mrs_real>(), "LPC_gamma");
00400 #endif
00401 
00402   //-------------------------
00403   // warped autocorrelation
00404   //-------------------------
00405   //calculate warped autocorrelation coeffs
00406   realvec r(in.getSize());
00407 
00408   //--------------------------
00409   // Levison-Durbin recursion
00410   //--------------------------
00411   //Calculate LPC alpha and reflections coeffs
00412   //using Levison-Durbin recursion
00413   //output format for LPC coeffs:
00414   // y(n) = x(n) - a(1)x(n-1) - ... - a(order_-1)x(n-order_-1)
00415   // a = [1 a(1) a(2) ... a(order_-1)]
00416   realvec a(order_+1);
00417   realvec k(order_+1); //reflection coeffs
00418 
00419   // previous implementation from musicDSP without correct gain estimation
00420   // LevinsonDurbin(r, a, k, LevinsonError);
00421 
00422   // implementation adapted from peter Kabal
00423   //SPautoc (in.getData(), in.getSize(), k.getData(), k.getSize());
00424   //LevinsonError = SPcorXpc (k.getData(), a.getData(), a.getSize()-1);
00425   //cout << LevinsonError << endl;
00426 
00427   //this also estimates the pitch - does it work if lambda != 0 [?] [ML] normalization issue
00428   autocorrelationWarped(in, r, lag, getctrl("mrs_real/lambda")->to<mrs_real>());
00429   LevinsonError = SPcorXpc (r.getData(), a.getData(), a.getSize()-1);
00430 
00431   // LevinsonError /= in.getSize(); [ML] add this if SPautoc used
00432   LevinsonError = sqrt(LevinsonError);
00433   // pitch in Hz
00434   pitch = getctrl("mrs_real/israte")->to<mrs_real>()/lag ;
00435 
00436   //--------------------------
00437   // LPC coeffs
00438   //--------------------------
00439   //output LPC coeffs
00440   // y(n) = x(n) - a(1)x(n-1) - ... - a(order_-1)x(n-order_-1)
00441   // a = [1 a(1) a(2) ... a(order_-1)]
00442   // out = [a(1) a(2) ... a(order_-1)]
00443 
00444 
00445   for(i=0; i < order_; ++i)
00446     // out(i) = a(i+1); // musicDsp implementation
00447     out(i) = -a(i); // musicDsp implementation
00448 
00449   //------------------------------
00450   //output pitch and power values
00451   //------------------------------
00452   out(order_) = pitch; // lag in samples <= from ::autocorrelationWarped(...) - does it work if lambda != 0 [?]
00453   out(order_+1) = LevinsonError; //prediction error (= gain? [?])
00454 
00455   //--------------------------
00456   // LPC Pole-shifting
00457   //--------------------------
00458   //verify if Z-Plane pole-shifting should be performed...
00459   mrs_real gamma = getctrl("mrs_real/gamma")->to<mrs_real>();
00460   if(gamma != 1.0)
00461   {
00462     for(mrs_natural j = 0; j < order_; j++)
00463     {
00464       out(j) = (out(j) * pow((mrs_real)gamma, (mrs_real)j+1));
00465     }
00466   }
00467 
00468   {
00469     MarControlAccessor acc(ctrl_coeffs_);
00470     realvec& coeffs = acc.to<mrs_realvec>();
00471     coeffs(0) = 1.0;
00472     for(i=1; i < order_+1; ++i)
00473     {
00474       // coeffs(i) = -a(i); // musicDsp implementation
00475       coeffs(i) = out(i-1);
00476       ctrl_pitch_->setValue(pitch);
00477       ctrl_power_->setValue(LevinsonError);
00478     }
00479   }
00480 
00481   //MATLAB engine stuff - for testing and validation purposes only!
00482 #ifdef _MATLAB_LPC_
00483   MATLAB_PUT(out, "LPC_out");
00484   MATLAB_PUT(getctrl("mrs_real/pitch")->to<mrs_real>(), "pitch_MARS");
00485   MATLAB_PUT(getctrl("mrs_real/power")->to<mrs_real>(), "g_MARS");
00486   MATLAB_PUT(ctrl_coeffs_->to<mrs_realvec>(), "a_MARS");
00487   MATLAB_EVAL("LPC_test");
00488   mrs_real matlabGain;
00489   MATLAB_GET("LPCgain", matlabGain);
00490 #endif
00491 
00492 }