Marsyas  0.5.0-beta1
/Users/jleben/code/marsyas/src/marsyas/marsystems/PowerSpectrum.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 "PowerSpectrum.h"
00020 
00021 using std::ostringstream;
00022 using namespace Marsyas;
00023 
00024 #define PSD_POWER 1
00025 #define PSD_MAG 2
00026 #define PSD_DB  3
00027 #define PSD_WDB  4
00028 #define PSD_PD  5
00029 #define PSD_LOGMAG 6
00030 #define PSD_LOGMAG2 7
00031 
00032 PowerSpectrum::PowerSpectrum(mrs_string name):MarSystem("PowerSpectrum",name)
00033 {
00034   ntype_ = PSD_POWER;
00035   N2_ = 0;
00036   re_ = 0.0;
00037   im_ = 0.0;
00038   dB_ = 0.0;
00039   pwr_ = 0.0;
00040 
00041   addControls();
00042 }
00043 
00044 PowerSpectrum::PowerSpectrum(const PowerSpectrum& a):MarSystem(a)
00045 {
00046   ctrl_spectrumType_ = getctrl("mrs_string/spectrumType");
00047 }
00048 
00049 
00050 PowerSpectrum::~PowerSpectrum()
00051 {
00052 }
00053 
00054 void
00055 PowerSpectrum::addControls()
00056 {
00057   addctrl("mrs_string/spectrumType", "power", ctrl_spectrumType_);
00058   setctrlState("mrs_string/spectrumType", true);
00059 }
00060 
00061 MarSystem*
00062 PowerSpectrum::clone() const
00063 {
00064   return new PowerSpectrum(*this);
00065 }
00066 
00067 void
00068 PowerSpectrum::myUpdate(MarControlPtr sender)
00069 {
00070   (void) sender;  //suppress warning of unused parameter(s)
00071 
00072   //Spectrum outputs N values, corresponding to N/2+1
00073   //complex and unique spectrum points - see Spectrum.h documentation
00074   N2_ = ctrl_inObservations_->to<mrs_natural>()/2 + 1;
00075 
00076   ctrl_onSamples_->setValue(ctrl_inSamples_, NOUPDATE);
00077   ctrl_onObservations_->setValue(N2_, NOUPDATE); //outputs N/2+1 real values
00078   ctrl_osrate_->setValue(ctrl_israte_->to<mrs_real>(), NOUPDATE);
00079 
00080   stype_ = ctrl_spectrumType_->to<mrs_string>();
00081   if (stype_ == "power")
00082     ntype_ = PSD_POWER;
00083   else if (stype_ == "magnitude")
00084     ntype_ = PSD_MAG;
00085   else if (stype_ == "decibels")
00086     ntype_ = PSD_DB; //PUT DB!! (DB->NEW)
00087   else if (stype_ == "wrongdBonsets")
00088     ntype_ = PSD_WDB;
00089   else if (stype_ == "powerdensity")
00090     ntype_ = PSD_PD;
00091   else if (stype_ == "logmagnitude")
00092     ntype_ = PSD_LOGMAG;
00093   else if (stype_ == "logmagnitude2")
00094     ntype_ = PSD_LOGMAG2;
00095 
00096   // Add prefix to the observation names.
00097   mrs_string inObsNames = ctrl_inObsNames_->to<mrs_string>();
00098   ctrl_onObsNames_->setValue("Power_" + stype_ + ctrl_inObsNames_->to<mrs_string>(), NOUPDATE);
00099 
00100   // gtzan: removed cluttered output
00101   // ctrl_onObsNames_->setValue(obsNamesAddPrefix(inObsNames, "Power" + stype_ + "_"), NOUPDATE);
00102 }
00103 
00104 void
00105 PowerSpectrum::myProcess(realvec& in, realvec& out)
00106 {
00107   mrs_natural t,o;
00108   for (t=0; t < inSamples_; ++t)
00109   {
00110     for (o=0; o < N2_; o++)
00111     {
00112       if (o==0) //DC bin (i.e. 0)
00113       {
00114         re_ = in(0,t);
00115         im_ = 0.0;
00116       }
00117       else if (o == N2_-1) //Nyquist bin (i.e. N/2)
00118       {
00119         re_ = in(1,t);
00120         im_ = 0.0;
00121       }
00122       else //all other bins
00123       {
00124         re_ = in(2*o, t);
00125         im_ = in(2*o+1, t);
00126       }
00127 
00128       switch (ntype_)
00129       {
00130       case PSD_POWER:
00131         out(o, t) = re_*re_ + im_*im_;
00132         break;
00133       case PSD_MAG:
00134         out(o,t) = sqrt(re_ * re_ + im_ * im_);
00135         break;
00136       case PSD_DB:
00137         // TODO(sness) - Check validity of this magic number.  Should probably be FLT_MIN
00138         dB_ = (mrs_real)(10*log10(re_ * re_ + im_ * im_ + 0.000000001));
00139         out(o,t) = dB_;
00140         break;
00141       case PSD_WDB:
00142         // NOTE : FIXME!!!
00143         // 20*log10() seems to work better for toy_with_onsets
00144         // This is not good, and should be fixed inside the code that
00145         // looks for onsets.
00146         dB_ = (mrs_real)(20*log10(re_ * re_ + im_ * im_ + 0.000000001));
00147         if (dB_ < -100) dB_ = -100;
00148         out(o,t) = dB_;
00149         break;
00150       case PSD_PD:
00151         pwr_ = re_ * re_ + im_ * im_;
00152         out(o,t) = (mrs_real)(2.0 * pwr_) / N2_;
00153         break;
00154       case PSD_LOGMAG2:
00155         out(o,t) = (mrs_real)log10(1+sqrt(re_ * re_ + im_ * im_));
00156         break;
00157       case PSD_LOGMAG:
00158         out(o,t) = log(1+1000.0 * sqrt(re_ * re_ + im_ * im_));
00159         //  out(o,t) = pow(sqrt(re_ * re_ + im_ * im_), 0.5);
00160         // out(o,t) = asin(sqrt(re_ * re_ + im_ * im_));
00161 
00162       }
00163     }
00164   }
00165 
00166   //MATLAB_PUT(out, "PowerSpectrum");
00167   //MATLAB_EVAL("plot(PowerSpectrum)");
00168 }