Marsyas
0.6.0-alpha
|
00001 /* 00002 ** Copyright (C) 1998-2008 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 "PvUnconvert.h" 00020 00021 #include <algorithm> 00022 00023 using std::ostringstream; 00024 00025 using namespace Marsyas; 00026 00027 PvUnconvert::PvUnconvert(mrs_string name):MarSystem("PvUnconvert",name) 00028 { 00029 //type_ = "PvUnconvert"; 00030 //name_ = name; 00031 00032 addControls(); 00033 transient_counter_ = 0; 00034 } 00035 00036 00037 PvUnconvert::PvUnconvert(const PvUnconvert& a):MarSystem(a) 00038 { 00039 ctrl_mode_ = getctrl("mrs_string/mode"); 00040 ctrl_peakPicking_ = getctrl("mrs_string/peakPicking"); 00041 00042 ctrl_lastphases_ = getctrl("mrs_realvec/lastphases"); 00043 ctrl_analysisphases_ = getctrl("mrs_realvec/analysisphases"); 00044 ctrl_regions_ = getctrl("mrs_realvec/regions"); 00045 ctrl_magnitudes_ = getctrl("mrs_realvec/magnitudes"); 00046 ctrl_peaks_ = getctrl("mrs_realvec/peaks"); 00047 00048 ctrl_phaselock_ = getctrl("mrs_bool/phaselock"); 00049 00050 transient_counter_ = 0; 00051 00052 } 00053 00054 00055 PvUnconvert::~PvUnconvert() 00056 { 00057 00058 } 00059 00060 00061 MarSystem* 00062 PvUnconvert::clone() const 00063 { 00064 return new PvUnconvert(*this); 00065 } 00066 00067 00068 void 00069 PvUnconvert::addControls() 00070 { 00071 addctrl("mrs_natural/Interpolation", MRS_DEFAULT_SLICE_NSAMPLES/4); 00072 addctrl("mrs_natural/Decimation", MRS_DEFAULT_SLICE_NSAMPLES/4); 00073 addctrl("mrs_string/mode", "loose_phaselock", ctrl_mode_); 00074 addctrl("mrs_string/peakPicking", "multires", ctrl_peakPicking_); 00075 addctrl("mrs_realvec/lastphases", realvec(), ctrl_lastphases_); 00076 addctrl("mrs_realvec/analysisphases", realvec(), ctrl_analysisphases_); 00077 addctrl("mrs_realvec/regions", realvec(), ctrl_regions_); 00078 addctrl("mrs_realvec/magnitudes", realvec(), ctrl_magnitudes_); 00079 addctrl("mrs_realvec/peaks", realvec(), ctrl_peaks_); 00080 addctrl("mrs_bool/phaselock", false, ctrl_phaselock_); 00081 } 00082 00083 void 00084 PvUnconvert::myUpdate(MarControlPtr sender) 00085 { 00086 (void) sender; //suppress warning of unused parameter(s) 00087 setctrl("mrs_natural/onSamples", getctrl("mrs_natural/inSamples")); 00088 setctrl("mrs_natural/onObservations", getctrl("mrs_natural/inObservations")->to<mrs_natural>() - 2); 00089 setctrl("mrs_real/osrate", getctrl("mrs_real/israte")->to<mrs_real>() / getctrl("mrs_natural/onObservations")->to<mrs_natural>()); 00090 00091 mrs_natural onObservations = getctrl("mrs_natural/onObservations")->to<mrs_natural>(); 00092 mrs_real israte = getctrl("mrs_real/israte")->to<mrs_real>(); 00093 00094 N2_ = onObservations/2; 00095 00096 if (N2_+1 != (mrs_natural)phase_.getSize()) 00097 { 00098 00099 { 00100 MarControlAccessor acc(ctrl_lastphases_); 00101 mrs_realvec& lastphases = acc.to<mrs_realvec>(); 00102 lastphases.create(N2_+1); 00103 } 00104 00105 { 00106 MarControlAccessor acc(ctrl_analysisphases_); 00107 mrs_realvec& analysisphases = acc.to<mrs_realvec>(); 00108 analysisphases.create(N2_+1); 00109 } 00110 00111 00112 { 00113 MarControlAccessor acc(ctrl_regions_); 00114 mrs_realvec& regions = acc.to<mrs_realvec>(); 00115 regions.create(N2_+1); 00116 for (int i=0; i < N2_+1; ++i) 00117 regions(i) = i; 00118 } 00119 00120 { 00121 MarControlAccessor acc(ctrl_magnitudes_); 00122 mrs_realvec& magnitudes = acc.to<mrs_realvec>(); 00123 magnitudes.create(N2_+1); 00124 } 00125 00126 { 00127 MarControlAccessor acc(ctrl_peaks_); 00128 mrs_realvec& peaks = acc.to<mrs_realvec>(); 00129 peaks.create(N2_+1); 00130 } 00131 00132 00133 phase_.create(N2_+1); 00134 lphase_.create(N2_+1); 00135 iphase_.create(N2_+1); 00136 lmag_.create(N2_+1); 00137 } 00138 00139 00140 fundamental_ = (mrs_real) (israte / onObservations); 00141 factor_ = (((getctrl("mrs_natural/Interpolation")->to<mrs_natural>()* TWOPI)/(israte))); 00142 00143 00144 } 00145 00146 00147 int 00148 PvUnconvert::subband(int bin) 00149 { 00150 int si; 00151 si = 0; 00152 00153 if (bin < 16) 00154 si = 0; 00155 else if ((bin >= 16) && (bin < 32)) 00156 si = 1; 00157 else if (bin < 512) 00158 si = (int)(log(bin*1.0) / log(2.0))-3; 00159 else if (bin > 512) 00160 si = 6; 00161 return si; 00162 } 00163 00164 00165 bool 00166 PvUnconvert::isPeak(int bin, mrs_realvec& magnitudes, mrs_real maxAmp) 00167 { 00168 bool res = true; 00169 00170 int h = subband(bin); 00171 h = 2; 00172 00173 if ((bin > 2) && (bin <= N2_-2)) 00174 for (int i = bin-h; i < bin+h; ++i) 00175 { 00176 if (magnitudes(bin) < magnitudes(i)) 00177 res = false; 00178 } 00179 00180 if (magnitudes(bin) < 0.005 * maxAmp) 00181 res = false; 00182 return res; 00183 } 00184 00185 00186 00187 void 00188 PvUnconvert::myProcess(realvec& in, realvec& out) 00189 { 00190 00191 mrs_natural t; 00192 MarControlAccessor acc(ctrl_lastphases_); 00193 mrs_realvec& lastphases = acc.to<mrs_realvec>(); 00194 MarControlAccessor acc1(ctrl_analysisphases_); 00195 mrs_realvec& analysisphases = acc1.to<mrs_realvec>(); 00196 00197 MarControlAccessor acc2(ctrl_regions_); 00198 mrs_realvec& regions = acc2.to<mrs_realvec>(); 00199 00200 MarControlAccessor acc3(ctrl_magnitudes_); 00201 mrs_realvec& magnitudes = acc3.to<mrs_realvec>(); 00202 00203 MarControlAccessor acc4(ctrl_peaks_); 00204 mrs_realvec& peaks = acc4.to<mrs_realvec>(); 00205 00206 peaks.setval(0.0); 00207 00208 00209 static int count = 0; 00210 count++; 00211 00212 mrs_natural re, amp, im, freq; 00213 mrs_real avg_re; 00214 mrs_real avg_im; 00215 00216 const mrs_string& mode = ctrl_mode_->to<mrs_string>(); 00217 00218 mrs_real interpolation = getctrl("mrs_natural/Interpolation")->to<mrs_natural>() * 1.0; 00219 mrs_real decimation = getctrl("mrs_natural/Decimation")->to<mrs_natural>() * 1.0; 00220 mrs_real tratio = interpolation / decimation; 00221 00222 00223 00224 mrs_real beta = 0.66 + tratio/3.0; 00225 if (mode == "identity_phaselock") 00226 beta = 1.0; 00227 00228 mrs_real maxAmp =0.0; 00229 00230 // calculate magnitude 00231 for (t=0; t <= N2_; t++) 00232 { 00233 re = amp = 2*t; 00234 im = freq = 2*t+1; 00235 if (t== N2_) 00236 { 00237 re = 1; 00238 } 00239 magnitudes(t) = in(re,0); 00240 if (t==N2_) 00241 magnitudes(t) = 0.0; 00242 if (magnitudes(t) > maxAmp) 00243 maxAmp = magnitudes(t); 00244 } 00245 00246 00247 if (mode == "loose_phaselock") 00248 { 00249 for (t=0; t <= N2_; t++) 00250 { 00251 re = amp = 2*t; 00252 im = freq = 2*t+1; 00253 if (t == N2_) 00254 re = 1; 00255 00256 phase_(t) = lastphases(t) + interpolation * in(freq,0); 00257 } 00258 } 00259 00260 // propagate phases for peaks and calculate regions of influence 00261 if ((mode == "identity_phaselock")||(mode == "scaled_phaselock")) 00262 { 00263 int previous_peak=0; 00264 int peak = 0; 00265 int peakCount = 0; 00266 00267 for (t=0; t <= N2_; t++) 00268 { 00269 re = amp = 2*t; 00270 im = freq = 2*t+1; 00271 if (t == N2_) 00272 re = 1; 00273 if (isPeak(t, magnitudes, maxAmp)) 00274 { 00275 peakCount++; 00276 if (mode == "identity_phaselock") 00277 iphase_(t) = lastphases(t) + interpolation * in(freq,0); 00278 else if (mode == "scaled_phaselock") 00279 iphase_(t) = lastphases((mrs_natural)regions(t)) + interpolation * in(freq,0); 00280 } 00281 } 00282 00283 for (t=0; t <= N2_; t++) 00284 { 00285 if (isPeak(t, magnitudes, maxAmp)) 00286 { 00287 // calculate significant peaks and corresponding 00288 // non-overlapping intervals 00289 peak = t; 00290 00291 if (peak-previous_peak == 1) 00292 regions(peak) = peak; 00293 else 00294 { 00295 for (int j=previous_peak; j< previous_peak + (int)((peak-previous_peak)/2.0); j++) 00296 { 00297 peaks(j) = magnitudes(previous_peak); 00298 regions(j) = previous_peak; 00299 } 00300 00301 for (int j= previous_peak + (int)((peak-previous_peak)/2.0); j < peak; j++) 00302 { 00303 peaks(j) = magnitudes(peak); 00304 regions(j) = peak; 00305 } 00306 } 00307 previous_peak = peak; 00308 } 00309 00310 } 00311 } 00312 00313 00314 00315 // resynthesis for all bins 00316 for (t=0; t <= N2_; t++) 00317 { 00318 re = amp = 2*t; 00319 im = freq = 2*t+1; 00320 00321 if ((mode == "identity_phaselock")||(mode == "scaled_phaselock")) 00322 { 00323 if (t == N2_) 00324 re = 1; 00325 00326 00327 // unwrap analysis phases 00328 while (analysisphases(t) > PI) 00329 analysisphases(t) -= TWOPI; 00330 while (analysisphases(t) < -PI) 00331 analysisphases(t) += TWOPI; 00332 00333 iphase_(t) = iphase_((mrs_natural)regions(t)) + 00334 beta * (analysisphases(t) - analysisphases((mrs_natural)regions(t))); 00335 00336 // sinusoidal trajectory continuation heuristic 00337 /* if (t - regions(t) > subband(t)) 00338 iphase_(t) = phase_(regions(t)) + beta * (analysisphases(t) - analysisphases(regions(t))); 00339 else 00340 { 00341 iphase_(t) = lastphases(t) + interpolation * in(freq,0); 00342 } 00343 */ 00344 00345 // } 00346 00347 if (ctrl_phaselock_->to<mrs_bool>()) 00348 { 00349 iphase_(t) = analysisphases(t); 00350 } 00351 00352 out(re,0) = magnitudes(t) * cos(iphase_(t)); 00353 if (t != N2_) 00354 out(im,0) = -magnitudes(t) * sin(iphase_(t)); 00355 lastphases(t) = iphase_(t); 00356 } 00357 else if (mode == "loose_phaselock") 00358 { 00359 if (t == N2_) 00360 re = 1; 00361 00362 if ((t >= 2) && (t < N2_-2)) 00363 { 00364 avg_re = magnitudes(t) * cos(phase_(t)) + 00365 0.25 * magnitudes(t-1) * cos(phase_(t-1)) + 00366 0.25 * magnitudes(t+1) * cos(phase_(t)); 00367 avg_im = -magnitudes(t) * sin(phase_(t)) - 00368 -0.25 * magnitudes(t-1) * sin(phase_(t-1)) - 00369 -0.25 * magnitudes(t+1) * sin(phase_(t)); 00370 lphase_(t) = -atan2(avg_im,avg_re); 00371 } 00372 else 00373 lphase_(t) = phase_(t); 00374 00375 lmag_(t) = magnitudes(t); 00376 00377 if (ctrl_phaselock_->to<mrs_bool>()) 00378 { 00379 00380 lphase_ = analysisphases; 00381 ctrl_phaselock_->setValue(false); 00382 } 00383 00384 /* while (lphase_(t) > PI) 00385 lphase_(t) -= TWOPI; 00386 while (lphase_(t) < -PI) 00387 lphase_(t) += TWOPI; 00388 */ 00389 00390 00391 out(re,0) = lmag_(t) * cos(lphase_(t)); 00392 if (t != N2_) 00393 out(im,0) = -lmag_(t) * sin(lphase_(t)); 00394 lastphases(t) = lphase_(t); 00395 } 00396 else // classic 00397 { 00398 if (t == N2_) 00399 { 00400 re = 1; 00401 } 00402 phase_(t) = lastphases(t) + interpolation * in(freq,0); 00403 00404 if (ctrl_phaselock_->to<mrs_bool>()) 00405 { 00406 phase_(t) = analysisphases(t); 00407 } 00408 00409 out(re,0) = magnitudes(t) * cos(phase_(t)); 00410 if (t != N2_) 00411 out(im,0) = -magnitudes(t) * sin(phase_(t)); 00412 lastphases(t) = phase_(t); 00413 } 00414 } 00415 00416 // Reset phaselock control 00417 if (ctrl_phaselock_->to<mrs_bool>()) 00418 { 00419 ctrl_phaselock_->setValue(false); 00420 00421 } 00422 00423 } 00424 00425 00426 00427 00428 00429 00430 00431 00432 00433 00434 00435 00436 00437 00438 00439 00440 00441 00442 00443 00444 00445 00446