Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/AutoCorrelation.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 "../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 }