Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/HarmonicStrength.cpp
Go to the documentation of this file.
00001 /*
00002 ** Copyright (C) 1998-2010 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 "../common_source.h"
00020 #include "HarmonicStrength.h"
00021 
00022 using std::ostringstream;
00023 using namespace Marsyas;
00024 
00025 using namespace std;
00026 
00027 HarmonicStrength::HarmonicStrength(mrs_string name) : MarSystem("HarmonicStrength", name)
00028 {
00029   addControls();
00030 }
00031 
00032 HarmonicStrength::HarmonicStrength(const HarmonicStrength& a) : MarSystem(a)
00033 {
00034   ctrl_base_frequency_ = getctrl("mrs_real/base_frequency");
00035   ctrl_harmonics_ = getctrl("mrs_realvec/harmonics");
00036   ctrl_harmonicsSize_ = getctrl("mrs_natural/harmonicsSize");
00037   ctrl_harmonicsWidth_ = getctrl("mrs_real/harmonicsWidth");
00038   ctrl_inharmonicity_B_ = getctrl("mrs_real/inharmonicity_B");
00039 }
00040 
00041 
00042 HarmonicStrength::~HarmonicStrength()
00043 {
00044 }
00045 
00046 MarSystem*
00047 HarmonicStrength::clone() const
00048 {
00049   return new HarmonicStrength(*this);
00050 }
00051 
00052 void
00053 HarmonicStrength::addControls()
00054 {
00055   addctrl("mrs_real/base_frequency", 440.0, ctrl_base_frequency_);
00056   addctrl("mrs_realvec/harmonics", realvec(0), ctrl_harmonics_);
00057   addctrl("mrs_natural/harmonicsSize", 0, ctrl_harmonicsSize_);
00058   setctrlState("mrs_natural/harmonicsSize", true);
00059   addctrl("mrs_real/harmonicsWidth", 0.05, ctrl_harmonicsWidth_);
00060   addctrl("mrs_natural/type", 0);
00061   addctrl("mrs_real/inharmonicity_B", 0.0, ctrl_inharmonicity_B_);
00062 }
00063 
00064 void
00065 HarmonicStrength::myUpdate(MarControlPtr sender)
00066 {
00067   MRSDIAG("HarmonicStrength.cpp - HarmonicStrength:myUpdate");
00068 
00070   MarSystem::myUpdate(sender);
00071 
00072   mrs_natural num_harmonics = ctrl_harmonicsSize_->to<mrs_natural>();
00073   // setup default harmonics
00074   {
00075     // leave this here for the scope of acc
00076     MarControlAccessor acc(ctrl_harmonics_);
00077     mrs_realvec& harmonics = acc.to<mrs_realvec>();
00078     if ((num_harmonics > 0) && (harmonics.getSize() == 0))
00079     {
00080       harmonics.stretch(num_harmonics);
00081       for (mrs_natural i=0; i < num_harmonics; ++i)
00082       {
00083         harmonics(i) = i+1;
00084       }
00085     }
00086   }
00087   ctrl_onObservations_->setValue(ctrl_harmonicsSize_->to<mrs_natural>(), NOUPDATE);
00088 
00089   //features names
00090   mrs_string orig = ctrl_inObsNames_->to<mrs_string>();
00091   // remove final comma in name
00092   orig = orig.substr(0, orig.size()-1);
00093   ostringstream oss;
00094   for (mrs_natural i = 0; i < num_harmonics; ++i)
00095   {
00096     oss << "HarmonicStrength_" + orig << i+1 << ",";
00097   }
00098   setctrl("mrs_string/onObsNames", oss.str());
00099 }
00100 
00101 
00102 // used inside myProcess
00103 mrs_real
00104 HarmonicStrength::quadratic_interpolation(mrs_real best_bin,
00105     mrs_realvec& in, mrs_natural t)
00106 {
00107   if ((best_bin == 0) || (best_bin == in.getRows()-1))
00108   {
00109     // don't try to interpolate using data that's
00110     // outside of the spectogram
00111     // TODO: find some convincing DSP thing to do in this case
00112     return in( (mrs_natural)best_bin, t);
00113   }
00114   // https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks.html
00115   // a = alpha, b = beta, g = gamma
00116   mrs_real a = in( (mrs_natural)best_bin - 1, t);
00117   mrs_real b = in( (mrs_natural)best_bin + 0, t);
00118   mrs_real g = in( (mrs_natural)best_bin + 1, t);
00119 
00120   mrs_real p = 0.5 * (a-g)/(a-2*b+g);
00121   // avoid some NaNs
00122   if ((p < -0.5) || (p > 0.5))
00123   {
00124     return b;
00125   }
00126   mrs_real yp = b - 0.25*(a-g)*p;
00127   // avoid all NaNs
00128   if (yp < b)
00129   {
00130     // I think this happens because the search window doesn't
00131     // encompass the entire spectrum, so the "highest" bin
00132     // might not actually be the highest one, if it was on the
00133     // edge of search window.
00134     // TODO: find some convincing DSP thing to do in this case
00135     return b;
00136   }
00137   return yp;
00138 }
00139 
00140 mrs_real
00141 HarmonicStrength::find_peak_magnitude(mrs_real central_bin, mrs_realvec& in,
00142                                       mrs_natural t, mrs_real low,
00143                                       mrs_real high)
00144 {
00145   // find peak in 2*radius
00146   // in real-world signals, the harmonic isn't always a
00147   // literal multiple of the base frequency.  We allow a bit
00148   // of margin (the "radius") to find the best bin
00149   mrs_natural best_bin = -1;
00150   mrs_real best_magnitude = 0;
00151   if (low < 0)
00152   {
00153     low = 0;
00154   }
00155   if (high > inObservations_ - 1)
00156   {
00157     high = inObservations_ - 1;
00158   }
00159   for (mrs_natural i = (mrs_natural)low; i < high; i++)
00160   {
00161     if (in(i,t) > best_magnitude)
00162     {
00163       best_bin = i;
00164       best_magnitude = in(i,t);
00165     }
00166   }
00167   if (best_bin >= 0)
00168   {
00169     best_magnitude = quadratic_interpolation(best_bin, in, t);
00170   }
00171   else
00172   {
00173     best_magnitude = in( (mrs_natural)central_bin, t);
00174   }
00175 
00176   return best_magnitude;
00177 }
00178 
00179 void
00180 HarmonicStrength::myProcess(realvec& in, realvec& out)
00181 {
00182   mrs_natural num_harmonics = ctrl_harmonicsSize_->to<mrs_natural>();
00183   mrs_real base_freq = ctrl_base_frequency_->to<mrs_real>();
00184   MarControlAccessor acc(ctrl_harmonics_);
00185   mrs_realvec& harmonics = acc.to<mrs_realvec>();
00186   mrs_real width = ctrl_harmonicsWidth_->to<mrs_real>();
00187 
00188   mrs_real freq2bin = 1.0 / ctrl_israte_->to<mrs_real>();
00189   //mrs_real bin2freq = ctrl_israte_->to<mrs_real>();
00190 
00191   // Iterate over the samples (remember, FFT is vertical)
00192   for (mrs_natural t = 0; t < inSamples_; t++)
00193   {
00194     mrs_real energy_rms = 0.0;
00195     for (mrs_natural o = 0; o < inObservations_; o++)
00196     {
00197       energy_rms += in(o, t) * in(o,t);
00198     }
00199     energy_rms = sqrt(energy_rms);
00200 
00201     for (mrs_natural h = 0; h < num_harmonics; h++)
00202     {
00203       mrs_real n = harmonics(h);
00204       mrs_real B = ctrl_inharmonicity_B_->to<mrs_real>();
00205 
00206       mrs_real freq = n * base_freq * sqrt(1.0 + B*n*n);
00207       mrs_real bin = freq * freq2bin;
00208       //mrs_real freqold= n*base_freq;
00209       //mrs_real binold = freqold*freq2bin;
00210       //cout<<B<<"\t"<<freq<<"\t"<<bin<<endl;
00211       mrs_real low = bin - width * inObservations_;
00212       mrs_real high = bin + width * inObservations_;
00213       //cout<<low<<"\t"<<high<<endl;
00214       //cout<<low*bin2freq<<"\t"<<high*bin2freq<<endl;
00215       mrs_real magnitude = find_peak_magnitude(bin, in, t, low, high);
00216       if (magnitude == 0)
00217       {
00218         magnitude = 1e-60;
00219       }
00220       else
00221       {
00222         switch (getctrl("mrs_natural/type")->to<mrs_natural>())
00223         {
00224         case 0:
00225           //out(h, t) = log(magnitude / energy_rms);
00226           out(h, t) = magnitude / energy_rms;
00227           break;
00228         case 1:
00229           out(h, t) = magnitude;
00230           break;
00231         case 2:
00232           out(h, t) = log(magnitude);
00233           break;
00234         default:
00235           out(h, t) = magnitude;
00236           break;
00237         }
00238       }
00239       /*
00240       if (out(h,t) != out(h,t)) {
00241           cout<<"zzzzzzzzzz NaN detected! zzzzzzz"<<endl;
00242           cout<<magnitude<<"\t"<<energy_rms<<endl;
00243       }
00244       */
00245     }
00246   }
00247 }
00248