Marsyas
0.6.0-alpha
|
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 "../common_source.h" 00020 #include "AutoCorrelation.h" 00021 #include "Windowing.h" 00022 00023 using std::cout; 00024 using std::endl; 00025 00026 using std::ostringstream; 00027 using namespace Marsyas; 00028 00029 AutoCorrelation::AutoCorrelation(mrs_string name):MarSystem("AutoCorrelation",name) 00030 { 00031 myfft_ = NULL; 00032 addControls(); 00033 } 00034 00035 AutoCorrelation::~AutoCorrelation() 00036 { 00037 delete myfft_; 00038 } 00039 00040 // copy constructor 00041 AutoCorrelation::AutoCorrelation(const AutoCorrelation& a):MarSystem(a) 00042 { 00043 myfft_ = NULL; 00044 00045 ctrl_magcompress_ = getctrl("mrs_real/magcompress"); 00046 ctrl_normalize_ = getctrl("mrs_natural/normalize"); 00047 ctrl_octaveCost_ = getctrl("mrs_real/octaveCost"); 00048 ctrl_voicingThreshold_ = getctrl("mrs_real/voicingThreshold"); 00049 ctrl_aliasedOutput_ = getctrl("mrs_bool/aliasedOutput"); 00050 ctrl_makePositive_ = getctrl("mrs_bool/makePositive"); 00051 ctrl_setr0to1_ = getctrl("mrs_bool/setr0to1"); 00052 ctrl_setr0to0_ = getctrl("mrs_bool/setr0to0"); 00053 ctrl_lowCutoff_ = getctrl("mrs_real/lowCutoff"); 00054 ctrl_highCutoff_ = getctrl("mrs_real/highCutoff"); 00055 } 00056 00057 void 00058 AutoCorrelation::addControls() 00059 { 00060 addctrl("mrs_real/magcompress", 2.0, ctrl_magcompress_); 00061 addctrl("mrs_natural/normalize", 0, ctrl_normalize_); 00062 addctrl("mrs_real/octaveCost", 0.0, ctrl_octaveCost_); 00063 addctrl("mrs_real/voicingThreshold", 0.1, ctrl_voicingThreshold_); 00064 addctrl("mrs_bool/aliasedOutput", false, ctrl_aliasedOutput_); 00065 addctrl("mrs_bool/makePositive", false, ctrl_makePositive_); 00066 addctrl("mrs_bool/setr0to1", false, ctrl_setr0to1_); 00067 addctrl("mrs_bool/setr0to0", true, ctrl_setr0to0_); 00068 addctrl("mrs_real/lowCutoff", 0.0, ctrl_lowCutoff_); 00069 addctrl("mrs_real/highCutoff", 1.0, ctrl_highCutoff_); 00070 00071 00072 ctrl_normalize_->setState(true); 00073 ctrl_octaveCost_->setState(true); 00074 ctrl_voicingThreshold_->setState(true); 00075 ctrl_aliasedOutput_->setState(true); 00076 ctrl_lowCutoff_->setState(true); 00077 ctrl_highCutoff_->setState(true); 00078 } 00079 00080 MarSystem* 00081 AutoCorrelation::clone() const 00082 { 00083 return new AutoCorrelation(*this); 00084 } 00085 00086 void 00087 AutoCorrelation::myUpdate(MarControlPtr sender) 00088 { 00089 (void) sender; //suppress warning of unused parameter(s) 00090 00091 if(!myfft_) 00092 myfft_ = new fft(); 00093 00094 00095 setctrl("mrs_natural/onSamples", getctrl("mrs_natural/inSamples")); 00096 setctrl("mrs_natural/onObservations", getctrl("mrs_natural/inObservations")); 00097 setctrl("mrs_real/osrate", getctrl("mrs_real/israte")); 00098 00099 // round up with ceil() 00100 lowSamples_ = (mrs_natural) ceil(inSamples_ 00101 * getctrl("mrs_real/lowCutoff")->to<mrs_real>()); 00102 numSamples_ = (mrs_natural) ceil( inSamples_ 00103 * getctrl("mrs_real/highCutoff")->to<mrs_real>() 00104 ) - lowSamples_; 00105 00106 if(ctrl_aliasedOutput_->to<mrs_bool>()) 00107 fftSize_ = inSamples_; //will create aliasing! 00108 else 00109 { 00110 //compute fft with a size of the next power of 2 of 2*inSamples-1 samples 00111 fftSize_ = (mrs_natural)pow(2.0, ceil(log(2.0*numSamples_-1.0)/log(2.0))); 00112 } 00113 00114 scratch_.create(fftSize_); 00115 00116 00117 // only working for hanning window 00118 normalize_ = 0; 00119 if(getctrl("mrs_natural/normalize")->to<mrs_natural>()) 00120 { 00121 cout << "NORM INIT" << endl; 00122 realvec tmp(getctrl("mrs_natural/onSamples")->to<mrs_natural>()); 00123 normalize_ = 1; 00124 norm_.create(getctrl("mrs_natural/onSamples")->to<mrs_natural>()); 00125 norm_.setval(1); 00126 Windowing win("Windowing"); 00127 win.updControl("mrs_string/type", "Hanning"); 00128 win.updControl("mrs_natural/inSamples", norm_.getCols()); 00129 win.updControl("mrs_natural/inObservations", norm_.getRows()); 00130 win.process(norm_, tmp); 00131 00132 AutoCorrelation autocorr("Autocorrelation"); 00133 autocorr.updControl("mrs_natural/inSamples", norm_.getCols()); 00134 autocorr.updControl("mrs_natural/inObservations", norm_.getRows()); 00135 autocorr.update(); 00136 autocorr.process(tmp, norm_); 00137 00138 for (mrs_natural i = 0 ; i < norm_.getSize() ; ++i) 00139 norm_(i) = 1/norm_(i); 00140 } 00141 00142 octaveCost_ = getctrl("mrs_real/octaveCost")->to<mrs_real>(); 00143 voicing_ = getctrl("mrs_real/voicingThreshold")->to<mrs_real>(); 00144 if(octaveCost_) 00145 { 00146 octaveCost_ *= octaveCost_; 00147 octaveMax_ = octaveCost_*log(36.0*inSamples_); 00148 } 00149 } 00150 00151 void 00152 AutoCorrelation::myProcess(realvec& in, realvec& out) 00153 { 00154 00155 00156 00157 00158 00159 mrs_natural t,o; 00160 k_ = ctrl_magcompress_->to<mrs_real>(); 00161 00162 // Copy to output to perform inplace fft and zeropad to double size 00163 00164 scratch_.create(fftSize_); //scratch_ needs to be reset every time 00165 for (o=0; o < inObservations_; o++) 00166 { 00167 for (t=lowSamples_; t < (lowSamples_+numSamples_); t++) 00168 { 00169 scratch_(t-(lowSamples_)) = in(o,t); 00170 } 00171 00172 00173 00174 00175 00176 //zeropad 00177 for(t=(lowSamples_+numSamples_); t < fftSize_; t++) 00178 scratch_(t) = 0.0; 00179 00180 00181 //get pointer to data (THIS BREAKS ENCAPSULATION! FIXME [!]) 00182 mrs_real *tmp = scratch_.getData(); 00183 00184 //compute forward FFT (of size fftSize_) 00185 myfft_->rfft(tmp, fftSize_/2, FFT_FORWARD); //rfft() takes as second argument half of the desired FFT size (see fft.cpp) 00186 00187 // Special case for zero and Nyquist/2, 00188 // which only have real part 00189 if (k_ == 2.0) 00190 { 00191 re_ = tmp[0]; 00192 tmp[0] = re_ * re_; 00193 re_ = tmp[1]; 00194 tmp[1] = re_ * re_; 00195 } 00196 else 00197 { 00198 re_ = tmp[0]; 00199 re_ = sqrt(re_ * re_); 00200 tmp[0] = pow(re_, k_); 00201 re_ = tmp[1]; 00202 re_ = sqrt(re_ * re_); 00203 tmp[1] = pow(re_, k_); 00204 } 00205 00206 // Compress the magnitude spectrum and zero 00207 // the imaginary part. 00208 for (t=1; t < fftSize_/2; t++) 00209 { 00210 re_ = tmp[2*t]; 00211 im_ = tmp[2*t+1]; 00212 if (k_ == 2.0) 00213 am_ = re_ * re_ + im_ * im_; 00214 else 00215 { 00216 am_ = sqrt(re_ * re_ + im_ * im_); 00217 am_ = pow(am_, k_); 00218 } 00219 tmp[2*t] = am_; 00220 tmp[2*t+1] = 0; 00221 } 00222 00223 // Take the inverse Fourier Transform (of size fftSize_) 00224 myfft_->rfft(tmp, fftSize_/2, FFT_INVERSE); 00225 00226 // Copy result to output 00227 if(normalize_) 00228 { 00229 cout << "NORM Normalization happening" << endl; 00230 for (t=0; t < onSamples_; t++) 00231 { 00232 out(o,t) = scratch_(t)*norm_(t); 00233 } 00234 } 00235 else 00236 for (t=0; t < onSamples_; t++) 00237 { 00238 // out(o,t) = 0.1 * scratch_(t) + 0.99 * out(o,t); 00239 out(o,t) = 1.0 * scratch_(t) + 0.0 * out(o,t); 00240 // out(o,t) = 0.5 * scratch_(t) + 0.5 * out(o,t); 00241 // out(o,t) += scratch_(t); 00242 00243 } 00244 00245 } 00246 00247 00248 if (ctrl_makePositive_->to<mrs_bool>()) 00249 { 00250 out -= out.minval(); 00251 } 00252 00253 if(octaveCost_) //is there a reference for this octaveCost computation [?] 00254 { 00255 for (o=0; o < inObservations_; o++) 00256 { 00257 mrs_real maxOut = 0; 00258 for (t=1 ; t<onSamples_/2 ; t++) 00259 if (out(o, t)> out(o, t+1) && out(o, t) > out(o, t-1) && out(o, t)>maxOut) 00260 maxOut = out(o, t) ; 00261 //cout << maxOut/out(o, 0)<< " " << 1+voicing_ << << endl; 00262 00263 if(maxOut && maxOut/out(o, 0) > 1-voicing_) 00264 for (t=1; t < onSamples_; t++) 00265 out(o, t) += octaveMax_-octaveCost_*log(36.0*t); 00266 else 00267 out.setval(0); 00268 } 00269 } 00270 00271 if (ctrl_setr0to1_->to<mrs_bool>()) 00272 { 00273 // out -= out.minval(); 00274 00275 /* for (o=0; o < onObservations_; o++) 00276 for (t=0; t< onSamples_-1; t++) 00277 { 00278 out(o,t) = out(o,t) / (onSamples_ - 1 - t); 00279 if (t > onSamples_-1-100) 00280 out(o,t) = 0.0; 00281 } 00282 */ 00283 00284 00285 00286 // mrs_real myNorm = out(0,0); 00287 // if (myNorm > 0) 00288 // out /= myNorm; 00289 } 00290 00291 00292 00293 00294 00295 if (ctrl_setr0to0_->to<mrs_bool>()) 00296 { 00297 00298 // for (o=0; o < onObservations_; o++) 00299 // out(o,0) = 0.0; 00300 00301 00302 00303 for (o=0; o < onObservations_; o++) 00304 for (t=0; t < onSamples_; t++) 00305 { 00306 out(o,t) = out(o,t); 00307 } 00308 } 00309 00310 /* 00311 MATLAB_PUT(in, "corr_in"); 00312 MATLAB_PUT(out, "corr"); 00313 00314 MATLAB_EVAL("subplot(211)"); 00315 MATLAB_EVAL("plot(corr_in)"); 00316 MATLAB_EVAL("subplot(212)"); 00317 MATLAB_EVAL("plot(corr)"); 00318 */ 00319 }