Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/TriangularFilterBank.cpp
Go to the documentation of this file.
00001 /*
00002 ** Copyright (C) 1998-2012 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 "TriangularFilterBank.h"
00020 
00021 
00022 using std::ostringstream;
00023 
00024 
00025 using namespace Marsyas;
00026 
00027 TriangularFilterBank::TriangularFilterBank(mrs_string name):MarSystem("TriangularFilterBank",name)
00028 {
00029   addControls();
00030   pfftSize_ = 0;
00031   psamplingRate_ = 0;
00032   mfcc_offsets_ = NULL;
00033   pcepstralCoefs_ = 0;
00034   cepstralCoefs_ = TriangularFilterBank::cepstralCoefs_default;
00035 }
00036 
00037 TriangularFilterBank::TriangularFilterBank(const TriangularFilterBank& 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_ = TriangularFilterBank::cepstralCoefs_default;
00045 }
00046 
00047 TriangularFilterBank::~TriangularFilterBank()
00048 {
00049   if (mfcc_offsets_!=NULL)
00050     delete [] mfcc_offsets_;
00051 }
00052 
00053 MarSystem*
00054 TriangularFilterBank::clone() const
00055 {
00056   return new TriangularFilterBank(*this);
00057 }
00058 
00059 
00060 void
00061 TriangularFilterBank::addControls() {
00063   addControl("mrs_natural/coefficients", TriangularFilterBank::cepstralCoefs_default, ctrl_cepstralCoefs_);
00064   setControlState("mrs_natural/coefficients", true);
00065 }
00066 
00067 void
00068 TriangularFilterBank::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 "TriangularFilterBankxx"
00088   mrs_string inObsName = stringSplit(ctrl_inObsNames_->to<mrs_string>(), ",")[0];
00089   ostringstream oss;
00090   for (i=0; i < cepstralCoefs_; ++i)
00091   {
00092     oss << "TriangularFilterBank" << 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 
00113     ctrl_onObservations_->setValue((mrs_natural)totalFilters_, NOUPDATE);
00114 
00115     // Linear filter boundaries
00116     for (i=0; i< linearFilters_; ++i)
00117     {
00118       freqs_(i) = lowestFrequency_ + i * linearSpacing_;
00119     }
00120 
00121     // Logarithmic filter boundaries
00122     mrs_real first_log = freqs_(linearFilters_-1);
00123     for (i=1; i<=logFilters_+2; ++i)
00124     {
00125       freqs_(linearFilters_-1+i) = first_log * pow(logSpacing_, (mrs_real)i);
00126     }
00127 
00128     // Triangles information
00129     for (i=0; i<totalFilters_; ++i)
00130     {
00131       lower_(i) = freqs_(i);
00132     }
00133 
00134     for (i=1; i<= totalFilters_; ++i)
00135     {
00136       center_(i-1) = freqs_(i);
00137     }
00138 
00139     for (i=2; i<= totalFilters_+1; ++i)
00140     {
00141       upper_(i-2) = freqs_(i);
00142     }
00143 
00144     for (i=0; i<totalFilters_; ++i)
00145     {
00146       triangle_heights_(i) = (mrs_real)(2.0 / (upper_(i) - lower_(i)));
00147     }
00148 
00149     fftFreqs_.stretch(fftSize_);
00150 
00151     for (i=0; i< fftSize_; ++i)
00152     {
00153       fftFreqs_(i) = (float)i / (float)fftSize_ * (float)samplingRate_;
00154     }
00155 
00156     mfccFilterWeights_.create(totalFilters_, fftSize_);
00157     mfccDCT_.create(cepstralCoefs_, totalFilters_);
00158 
00159     mrs_natural chan;
00160 
00161     // NEIL's filter weight speedup
00162     if (pfftSize_!=fftSize_)
00163     {
00164       if (mfcc_offsets_!=NULL)
00165       {
00166         delete [] mfcc_offsets_;
00167       }
00168       mfcc_offsets_ = new int[totalFilters_*fftSize_*2];
00169     }
00170     // Initialize mfccFilterWeights
00171     for (chan = 0; chan < totalFilters_; chan++)
00172     {
00173       // NEIL's filter weight speedup
00174       int len=-1;
00175       int pos=0;
00176       for (i=0; i< fftSize_; ++i)
00177       {
00178         if ((fftFreqs_(i) > lower_(chan))&& (fftFreqs_(i) <= center_(chan)))
00179         {
00180           mfccFilterWeights_(chan, i) = triangle_heights_(chan) *
00181                                         ((fftFreqs_(i) - lower_(chan))/(center_(chan) - lower_(chan)));
00182           // NEIL's filter weight speedup
00183           if (len==-1)
00184           {
00185             pos=i;
00186           }
00187           len=i;
00188         }
00189         if ((fftFreqs_(i) > center_(chan)) && (fftFreqs_(i) <= upper_(chan)))
00190         {
00191           mfccFilterWeights_(chan, i) = triangle_heights_(chan) *
00192                                         ((upper_(chan) - fftFreqs_(i))/(upper_(chan) - center_(chan)));
00193           // NEIL's filter weight speedup
00194           if (len==-1)
00195           {
00196             pos=i;
00197           }
00198           len=i;
00199         }
00200       }
00201       // NEIL's filter weight speedup
00202       mfcc_offsets_[chan] = pos;
00203       mfcc_offsets_[chan+totalFilters_] = len;
00204     }
00205 
00206     // Initialize TriangularFilterBank_DCT
00207     mrs_real scale_fac = (mrs_real)(1.0/ sqrt((mrs_real)(totalFilters_/2)));
00208     for (j = 0; j<cepstralCoefs_; j++)
00209     {
00210       for (i=0; i< totalFilters_; ++i)
00211       {
00212         mfccDCT_(j, i) = scale_fac * cos(j * (2*i +1) * PI/2/totalFilters_);
00213         if (j == 0)
00214         {
00215           mfccDCT_(j,i) *= (mrs_real)(sqrt(2.0)/2.0);
00216         }
00217       }
00218     }
00219   }
00220 
00221 
00222   pfftSize_ = fftSize_;
00223   psamplingRate_ = samplingRate_;
00224 
00225   fmagnitude_.stretch(ctrl_inObservations_->to<mrs_natural>() * 2);
00226   earMagnitude_.stretch(totalFilters_);
00227 }
00228 
00229 void
00230 TriangularFilterBank::myProcess(realvec& in, realvec& out)
00231 {
00232   mrs_natural i,k,o;
00233 
00234   // mirror the spectrum
00235   for (o=0; o < inObservations_; o++)
00236   {
00237     fmagnitude_(o) = in(o,0);
00238   }
00239 
00240   for (o=0; o< inObservations_; o++)
00241   {
00242     fmagnitude_(o + inObservations_) = fmagnitude_(inObservations_ - o -1);
00243   }
00244 
00245   mrs_real sum =0.0;
00246   // Calculate the filterbank responce
00247   for (i=0; i<totalFilters_; ++i)
00248   {
00249     sum = 0.0;
00250     // NEIL's filter weight speedup
00251     for (k=mfcc_offsets_[i]; k<=mfcc_offsets_[i+totalFilters_]; k++)
00252     {
00253       sum += (mfccFilterWeights_(i, k) * fmagnitude_(k));
00254     }
00255     if (sum != 0.0)
00256     {
00257       earMagnitude_(i) = log10(20.0 * sum + 1.0);
00258     }
00259     else
00260     {
00261       earMagnitude_(i) = 0.0;
00262     }
00263   }
00264 
00265   for (o=0; o < onObservations_; o++)
00266     out(o,0) = earMagnitude_(o);
00267 
00268 //  MRSMSG(out);
00269 
00270 }
00271