Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/HWPS.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 "HWPS.h"
00020 #include <marsyas/NumericLib.h>
00021 #include <algorithm>
00022 
00023 using std::ostringstream;
00024 using std::min;
00025 using std::max;
00026 using std::abs;
00027 
00028 using namespace Marsyas;
00029 
00030 HWPS::HWPS(mrs_string name):MarSystem("HWPS", name)
00031 {
00032   addControls();
00033 }
00034 
00035 HWPS::HWPS(const HWPS& a) : MarSystem(a)
00036 {
00037   ctrl_histSize_ = getctrl("mrs_natural/histSize");
00038   ctrl_calcDistance_ = getctrl("mrs_bool/calcDistance");
00039 }
00040 
00041 HWPS::~HWPS()
00042 {
00043 }
00044 
00045 MarSystem*
00046 HWPS::clone() const
00047 {
00048   return new HWPS(*this);
00049 }
00050 
00051 void
00052 HWPS::addControls()
00053 {
00054   addctrl("mrs_bool/calcDistance", false, ctrl_calcDistance_);
00055   addctrl("mrs_natural/histSize", 20, ctrl_histSize_);
00056 }
00057 
00058 void
00059 HWPS::harmonicWrap(mrs_real peak1Freq, mrs_real peak2Freq, realvec& peak1SetFreqs, realvec& peak2SetFreqs)
00060 {
00061 
00062   // fundamental frequency estimate
00063   mrs_real hF;
00064 
00065   // Use the lowest in frequency highest amplitude
00066   // peak of the frames under consideration
00067   hF = min(peak1SetFreqs(0), peak2SetFreqs(0));
00068 
00069   // Original HWPS using the considered peaks for folding
00070   // hF = min(peak1Freq, peak2Freq);
00071 
00072   // mrs_real mhF = min(hF, abs(peak1Freq-peak2Freq));
00073 
00074   // shift frequencies
00075   peak1SetFreqs -= peak1Freq;
00076   peak2SetFreqs -= peak2Freq;
00077 
00078 
00079   /*
00080     MATLAB_PUT(peak1SetFreqs, "P1");
00081     MATLAB_PUT(peak2SetFreqs, "P2");
00082     MATLAB_EVAL("clf ; subplot(3, 1, 1);  hold ; stem(P1, A1); stem(P2, A2, 'r')");
00083   */
00084 
00085   // wrap frequencies around fundamental freq estimate
00086   peak1SetFreqs /= hF;
00087   peak2SetFreqs /= hF;
00088 
00089   for (mrs_natural k=0 ; k<peak1SetFreqs.getSize() ; k++)
00090   {
00091     peak1SetFreqs(k)=fmod(peak1SetFreqs(k), 1);
00092     //if(peak1SetFreqs(k)<0)
00093     while(peak1SetFreqs(k)<0)//replacing "if" in case of strongly negative (=> multiple wraps)
00094       peak1SetFreqs(k)+=1;
00095   }
00096   for (mrs_natural k=0 ; k<peak2SetFreqs.getSize() ; k++)
00097   {
00098     peak2SetFreqs(k)=fmod(peak2SetFreqs(k), 1);
00099     //if(peak2SetFreqs(k)<0)
00100     while(peak2SetFreqs(k)<0) //replacing "if" in case of strongly negative (=> multiple wraps)
00101       peak2SetFreqs(k)+=1;
00102   }
00103 }
00104 
00105 void
00106 HWPS::discretize(const realvec& peakSetWrapFreqs, const realvec& peakAmps,
00107                  const mrs_natural& histSize, realvec& resultHistogram)
00108 {
00109   mrs_natural index;
00110 
00111   resultHistogram.create(histSize);
00112 
00113   for (mrs_natural i=0 ; i<peakSetWrapFreqs.getSize() ; ++i)
00114   {
00115     index = (mrs_natural) fmod(floor(peakSetWrapFreqs(i)*histSize+.5), histSize);
00116     resultHistogram(index) += peakAmps(i);
00117   }
00118 }
00119 
00120 void
00121 HWPS::myUpdate(MarControlPtr sender)
00122 {
00123   (void) sender;  //suppress warning of unused parameter(s)
00124   if(inSamples_ > 1) {
00125     MRSWARN("HWPS::myUpdate - inSamples > 1 : only first column will be processed!");
00126   }
00127 
00128   ctrl_onObservations_->setValue(1, NOUPDATE);
00129   ctrl_onSamples_->setValue(1, NOUPDATE);
00130   ctrl_osrate_->setValue(ctrl_israte_, NOUPDATE); //[?]
00131   ctrl_onObsNames_->setValue("HWPS", NOUPDATE);
00132 
00133   //the input has the two vectors/matrices to process stacked vertically
00134   if(inObservations_ % 2 != 0) {
00135     MRSWARN("HWPS::myUpdate - input flow controls do not seem to be in a valid format!");
00136   }
00137 
00138   vec_i_.create(ctrl_inObservations_->to<mrs_natural>()/2);
00139   vec_j_.create(ctrl_inObservations_->to<mrs_natural>()/2);
00140 }
00141 
00142 void
00143 HWPS::myProcess(realvec& in, realvec& out)
00144 {
00145   mrs_natural o;
00146   //get the two stacked vectors from the input
00147   for(o=0; o < inObservations_/2; ++o)
00148   {
00149     vec_i_(o) = in(o,0);
00150     vec_j_(o) = in(o+inObservations_/2,0);
00151   }
00152 
00153   //parse data from the input vectors.
00154   //format:
00155   //[peakFreq, frameNumPeaks, pkSetFreq0,...,pkSetFreqN, pkSetAmp0,..., pkSetAmpN, dummyPad, dummyPad, ...]'
00156   //where N = frameNumPeaks-1, and dummyPad are dummy values used for filling the vector when frameNumPeaks < maxFrameNumPeaks.
00157   pk_i_freq_ = vec_i_(HWPS::pkFreqIdx); //peak i frequency
00158   pk_j_freq_ = vec_j_(HWPS::pkFreqIdx); //peak j frequency
00159   i_frameNumPeaks_ = (mrs_natural)vec_i_(HWPS::frameNumPeaksIdx); //peakSet i frameNumPeaks
00160   j_frameNumPeaks_ = (mrs_natural)vec_j_(HWPS::frameNumPeaksIdx); //peakSet j frameNumPeaks
00161 
00162   //get i peaksets freqs and amplitudes
00163   pkSet_i_Freqs_.stretch(i_frameNumPeaks_);
00164   pkSet_i_Amps_.stretch(i_frameNumPeaks_);
00165   for(o=0; o < i_frameNumPeaks_; ++o)
00166   {
00167     pkSet_i_Freqs_(o) = vec_i_(o+pkSetFeatsIdx);
00168     pkSet_i_Amps_(o) = vec_i_(o+pkSetFeatsIdx+i_frameNumPeaks_);
00169   }
00170   //get j peaksets freqs and amplitudes
00171   pkSet_j_Freqs_.stretch(j_frameNumPeaks_);
00172   pkSet_j_Amps_.stretch(j_frameNumPeaks_);
00173   for(o=0; o < j_frameNumPeaks_; ++o)
00174   {
00175     pkSet_j_Freqs_(o) = vec_j_(o+pkSetFeatsIdx);
00176     pkSet_j_Amps_(o) = vec_j_(o+pkSetFeatsIdx+j_frameNumPeaks_);
00177   }
00178 
00179   //perform harmonic wrapping
00180   pkSet_i_WrapFreqs_ = pkSet_i_Freqs_;
00181   pkSet_j_WrapFreqs_ = pkSet_j_Freqs_;
00182   harmonicWrap(pk_i_freq_, pk_j_freq_, pkSet_i_WrapFreqs_, pkSet_j_WrapFreqs_);
00183 
00184   /*
00185     MATLAB_PUT(pkSet_i_WrapFreqs_, "P1");
00186     MATLAB_PUT(pkSet_j_WrapFreqs_, "P2");
00187     MATLAB_PUT(pkSet_i_Amps_, "A1");
00188     MATLAB_PUT(pkSet_j_Amps_, "A2");
00189     MATLAB_EVAL("subplot(3, 1, 2);  hold ; stem(P1, A1); stem(P2, A2, 'r')");
00190   */
00191 
00192   //create histograms for both peaks
00193   histSize_ = ctrl_histSize_->to<mrs_natural>();
00194   discretize(pkSet_i_WrapFreqs_, pkSet_i_Amps_, histSize_, histogram_i_);
00195   discretize(pkSet_j_WrapFreqs_, pkSet_j_Amps_, histSize_, histogram_j_);
00196 
00197   /*
00198     MATLAB_PUT(histogram_i_, "H1");
00199     MATLAB_PUT(histogram_j_, "H2");
00200     MATLAB_EVAL("subplot(3, 1, 3); bar([H1; H2]')");
00201   */
00202 
00203   if(ctrl_calcDistance_->isTrue())
00204   {
00205     //return the cosine DISTANCE between the two histograms and we get the HWPDistance!
00206     out(0) = NumericLib::cosineDistance(histogram_i_, histogram_j_);
00207   }
00208   else
00209   {
00210     //return the cosine SIMILARITY between the two histograms and we get the HWPSimilarity!
00211     out(0) = 1.0 - NumericLib::cosineDistance(histogram_i_, histogram_j_);
00212   }
00213 }