Marsyas  0.5.0-beta1
/Users/jleben/code/marsyas/src/marsyas/marsystems/MFCC.cpp
Go to the documentation of this file.
00001 /*
00002 ** Copyright (C) 1998-2011 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 "MFCC.h"
00020 
00021 
00022 using std::ostringstream;
00023 
00024 
00025 using namespace Marsyas;
00026 
00027 MFCC::MFCC(mrs_string name):MarSystem("MFCC",name)
00028 {
00029   addControls();
00030   pfftSize_ = 0;
00031   psamplingRate_ = 0;
00032   mfcc_offsets_ = NULL;
00033   pcepstralCoefs_ = 0;
00034   cepstralCoefs_ = MFCC::cepstralCoefs_default;
00035 }
00036 
00037 MFCC::MFCC(const MFCC& a) : MarSystem(a)
00038 {
00039   ctrl_cepstralCoefs_ = getctrl("mrs_natural/coefficients");
00040   pfftSize_ = 0;
00041   psamplingRate_ = 0;
00042   mfcc_offsets_ = NULL;
00043   pcepstralCoefs_ = 0;
00044   cepstralCoefs_ = MFCC::cepstralCoefs_default;
00045 }
00046 
00047 MFCC::~MFCC()
00048 {
00049   if (mfcc_offsets_!=NULL)
00050     delete [] mfcc_offsets_;
00051 }
00052 
00053 MarSystem*
00054 MFCC::clone() const
00055 {
00056   return new MFCC(*this);
00057 }
00058 
00059 
00060 void
00061 MFCC::addControls() {
00063   addControl("mrs_natural/coefficients", MFCC::cepstralCoefs_default, ctrl_cepstralCoefs_);
00064   setControlState("mrs_natural/coefficients", true);
00065 }
00066 
00067 void
00068 MFCC::myUpdate(MarControlPtr sender)
00069 {
00070   (void) sender;  //suppress warning of unused parameter(s)
00071 
00072   // Get the number of cepstral coefficients from the control
00073   cepstralCoefs_ = ctrl_cepstralCoefs_->to<mrs_natural>();
00074 
00075   ctrl_onSamples_->setValue((mrs_natural)1, NOUPDATE);
00076   ctrl_onObservations_->setValue((mrs_natural)cepstralCoefs_, NOUPDATE);
00077   ctrl_osrate_->setValue(ctrl_israte_, NOUPDATE);
00078 
00079   // Initialize frequency boundaries for filters
00080   mrs_natural i,j;
00081   fftSize_ = 2 * (ctrl_inObservations_->to<mrs_natural>()-1); //PowerSpectrum outputs N/2+1 "magnitude" spectral points!
00082   if (fftSize_ == 0) return; // skip unnecessary updates
00083 
00084   samplingRate_ = (mrs_natural) (ctrl_israte_->to<mrs_real>() * fftSize_);
00085 
00086   // Set the observation names: use the first item of the
00087   // input observation names and prefix it with "MFCCxx"
00088   mrs_string inObsName = stringSplit(ctrl_inObsNames_->to<mrs_string>(), ",")[0];
00089   ostringstream oss;
00090   for (i=0; i < cepstralCoefs_; ++i)
00091   {
00092     oss << "MFCC" << i << "_" << inObsName << ",";
00093   }
00094   ctrl_onObsNames_->setValue(oss.str(), NOUPDATE);
00095 
00096   if ((pfftSize_ != fftSize_) || (psamplingRate_ != samplingRate_) || (pcepstralCoefs_ != cepstralCoefs_))
00097   {
00098 
00099     freqs_.create(42);
00100     lowestFrequency_ = 133.3333f;
00101     linearFilters_ = 13;
00102     linearSpacing_ = 66.66666f;
00103     logFilters_ = 27;
00104     logSpacing_ = 1.0711703f;
00105 
00106     totalFilters_ = linearFilters_ + logFilters_;
00107     lower_.create(totalFilters_);
00108     center_.create(totalFilters_);
00109     upper_.create(totalFilters_);
00110     triangle_heights_.create(totalFilters_);
00111 
00112     // Linear filter boundaries
00113     for (i=0; i< linearFilters_; ++i)
00114     {
00115       freqs_(i) = lowestFrequency_ + i * linearSpacing_;
00116     }
00117 
00118     // Logarithmic filter boundaries
00119     mrs_real first_log = freqs_(linearFilters_-1);
00120     for (i=1; i<=logFilters_+2; ++i)
00121     {
00122       freqs_(linearFilters_-1+i) = first_log * pow(logSpacing_, (mrs_real)i);
00123     }
00124 
00125     // Triangles information
00126     for (i=0; i<totalFilters_; ++i)
00127     {
00128       lower_(i) = freqs_(i);
00129       center_(i) = freqs_(i+1);
00130       upper_(i) = freqs_(i+2);
00131       triangle_heights_(i) = (mrs_real)(2.0 / (upper_(i) - lower_(i)));
00132     }
00133 
00134     fftFreqs_.stretch(fftSize_);
00135 
00136     for (i=0; i< fftSize_; ++i)
00137     {
00138       fftFreqs_(i) = (float)i / (float)fftSize_ * (float)samplingRate_;
00139     }
00140 
00141     mfccFilterWeights_.create(totalFilters_, fftSize_);
00142     mfccDCT_.create(cepstralCoefs_, totalFilters_);
00143 
00144     mrs_natural chan;
00145 
00146     // NEIL's filter weight speedup
00147     if (pfftSize_!=fftSize_)
00148     {
00149       if (mfcc_offsets_!=NULL)
00150       {
00151         delete [] mfcc_offsets_;
00152       }
00153       mfcc_offsets_ = new int[totalFilters_*fftSize_*2];
00154     }
00155     // Initialize mfccFilterWeights
00156     for (chan = 0; chan < totalFilters_; chan++)
00157     {
00158       // NEIL's filter weight speedup
00159       int len=-1;
00160       int pos=0;
00161       for (i=0; i< fftSize_; ++i)
00162       {
00163         if ((fftFreqs_(i) > lower_(chan))&& (fftFreqs_(i) <= center_(chan)))
00164         {
00165           mfccFilterWeights_(chan, i) = triangle_heights_(chan) *
00166                                         ((fftFreqs_(i) - lower_(chan))/(center_(chan) - lower_(chan)));
00167           // NEIL's filter weight speedup
00168           if (len==-1)
00169           {
00170             pos=i;
00171           }
00172           len=i;
00173         }
00174         if ((fftFreqs_(i) > center_(chan)) && (fftFreqs_(i) <= upper_(chan)))
00175         {
00176           mfccFilterWeights_(chan, i) = triangle_heights_(chan) *
00177                                         ((upper_(chan) - fftFreqs_(i))/(upper_(chan) - center_(chan)));
00178           // NEIL's filter weight speedup
00179           if (len==-1)
00180           {
00181             pos=i;
00182           }
00183           len=i;
00184         }
00185       }
00186       // NEIL's filter weight speedup
00187       mfcc_offsets_[chan] = pos;
00188       mfcc_offsets_[chan+totalFilters_] = len;
00189     }
00190 
00191     // Initialize MFCC_DCT
00192     mrs_real scale_fac = (mrs_real)(1.0/ sqrt((mrs_real)(totalFilters_/2)));
00193     for (j = 0; j<cepstralCoefs_; j++)
00194     {
00195       for (i=0; i< totalFilters_; ++i)
00196       {
00197         mfccDCT_(j, i) = scale_fac * cos(j * (2*i +1) * PI/2/totalFilters_);
00198         if (j == 0)
00199         {
00200           mfccDCT_(j,i) *= (mrs_real)(sqrt(2.0)/2.0);
00201         }
00202       }
00203     }
00204   }
00205 
00206 
00207   pfftSize_ = fftSize_;
00208   psamplingRate_ = samplingRate_;
00209 
00210   fmagnitude_.stretch(ctrl_inObservations_->to<mrs_natural>() * 2);
00211   earMagnitude_.stretch(totalFilters_);
00212 }
00213 
00214 void
00215 MFCC::myProcess(realvec& in, realvec& out)
00216 {
00217   mrs_natural i,k,o;
00218 
00219   // mirror the spectrum
00220   for (o=0; o < inObservations_; o++)
00221   {
00222     fmagnitude_(o) = in(o,0);
00223   }
00224 
00225   for (o=0; o< inObservations_; o++)
00226   {
00227     fmagnitude_(o + inObservations_) = fmagnitude_(inObservations_ - o -1);
00228   }
00229 
00230   mrs_real sum =0.0;
00231   // Calculate the filterbank responce
00232   for (i=0; i<totalFilters_; ++i)
00233   {
00234     sum = 0.0;
00235     // NEIL's filter weight speedup
00236     for (k=mfcc_offsets_[i]; k<=mfcc_offsets_[i+totalFilters_]; k++)
00237     {
00238       sum += (mfccFilterWeights_(i, k) * fmagnitude_(k));
00239     }
00240     if (sum != 0.0)
00241     {
00242       earMagnitude_(i) = log10(sum);
00243     }
00244     else
00245     {
00246       earMagnitude_(i) = 0.0;
00247     }
00248   }
00249   /* The way it used to be : NEIL
00250   for (i=0; i<totalFilters_; ++i) {
00251   sum = 0.0;
00252   for (k=0; k<fftSize_; k++) {
00253   sum += (mfccFilterWeights_(i, k) * fmagnitude_(k));
00254   }
00255   if (sum != 0.0)
00256   earMagnitude_(i) = log10(sum);
00257   else
00258   earMagnitude_(i) = 0.0;
00259   }
00260   */
00261 
00262   // Take the DCT
00263   for (o=0; o < cepstralCoefs_; o++)
00264   {
00265     sum =0.0;
00266     for (k=0; k < totalFilters_; k++)
00267     {
00268       sum += (mfccDCT_(o,k) * earMagnitude_(k));
00269     }
00270     out(o,0) = sum;
00271   }
00272 }