Marsyas
0.6.0-alpha
|
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 }