Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/PvUnconvert.cpp
Go to the documentation of this file.
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