Marsyas
0.6.0-alpha
|
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 }