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 "../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