Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/ERB.cpp
Go to the documentation of this file.
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 
00020 #include "ERB.h"
00021 #include "../common_source.h"
00022 #include "Series.h"
00023 #include "Filter.h"
00024 #include <sstream>
00025 
00026 
00027 using std::ostringstream;
00028 using std::stringstream;
00029 
00030 using namespace Marsyas;
00031 
00032 ERB::ERB(mrs_string name):MarSystem("ERB", name),
00033   filterBank (0)
00034 {
00035   //type_ = "ERB";
00036   //name_ = name;
00037 
00038   fs = 0;
00039   numChannels = 0;
00040   lowFreq = 0;
00041 
00042   addControls();
00043 }
00044 
00045 
00046 ERB::~ERB()
00047 {
00048 }
00049 
00050 MarSystem*
00051 ERB::clone() const
00052 {
00053   return new ERB(*this);
00054 }
00055 
00056 
00057 void
00058 ERB::addControls()
00059 {
00060   addctrl("mrs_natural/numChannels", 1);
00061   addctrl("mrs_real/lowFreq", 100.0f);
00062 
00063   setctrlState("mrs_natural/numChannels", true);
00064   setctrlState("mrs_real/lowFreq", true);
00065 }
00066 
00067 void
00068 ERB::myUpdate(MarControlPtr sender)
00069 {
00070   (void) sender;  //suppress warning of unused parameter(s)
00071   MRSDIAG("ERB.cpp - ERB:myUpdate");
00072 
00073   //FilterBank creation
00074   if (numChannels != getctrl("mrs_natural/numChannels")->to<mrs_natural>()) {
00075     numChannels = getctrl("mrs_natural/numChannels")->to<mrs_natural>();
00076     if (filterBank) delete filterBank;
00077     filterBank = new Fanout("filterBank");
00078     stringstream name;
00079     for(mrs_natural i = 0; i < numChannels; ++i) {
00080       name << "channel_" << i;
00081       Series* channel = new Series(name.str());
00082       name.str("");
00083       name << "filter_" << i << "_" << 0;
00084       Filter* filter_0 = new Filter(name.str());
00085       name.str("");
00086       name << "filter_" << i << "_" << 1;
00087       Filter* filter_1 = new Filter(name.str());
00088       name.str("");
00089       name << "filter_" << i << "_" << 2;
00090       Filter* filter_2 = new Filter(name.str());
00091       name.str("");
00092       name << "filter_" << i << "_" << 3;
00093       Filter* filter_3 = new Filter(name.str());
00094       name.str("");
00095       channel->addMarSystem(filter_0);
00096       channel->addMarSystem(filter_1);
00097       channel->addMarSystem(filter_2);
00098       channel->addMarSystem(filter_3);
00099       filterBank->addMarSystem(channel);
00100     }
00101     centerFreqs.create(numChannels);
00102     fcoefs.create(numChannels, 10);
00103   }
00104 
00105   setctrl("mrs_natural/onSamples", getctrl("mrs_natural/inSamples"));
00106   setctrl("mrs_natural/onObservations", numChannels*getctrl("mrs_natural/inObservations")->to<mrs_natural>());
00107   setctrl("mrs_real/osrate", getctrl("mrs_real/israte"));
00108 
00109   //Coefficients computation
00110   lowFreq = getctrl("mrs_real/lowFreq")->to<mrs_real>();
00111   fs = getctrl("mrs_real/israte")->to<mrs_real>();
00112   highFreq = fs/2;
00113   EarQ = 9.26449f;
00114   minBW = 24.7f;
00115   order = 1;
00116 
00117   for (mrs_natural i = 0; i < numChannels; ++i) {
00118     centerFreqs(i) = - (EarQ*minBW) + exp((i+1)*(-log(highFreq+EarQ*minBW) + log(lowFreq+EarQ*minBW))/numChannels)*(highFreq + EarQ*minBW);
00119   }
00120 
00121   A0 = 1/fs;
00122   A2 = 0.0f;
00123   B0 = 1.0f;
00124 
00125   for (mrs_natural i = 0; i < numChannels; ++i) {
00126     fcoefs(i,0) = A0;
00127     fcoefs(i,1) = A11(centerFreqs(i), B(E(centerFreqs(i))));
00128     fcoefs(i,2) = A12(centerFreqs(i), B(E(centerFreqs(i))));
00129     fcoefs(i,3) = A13(centerFreqs(i), B(E(centerFreqs(i))));
00130     fcoefs(i,4) = A14(centerFreqs(i), B(E(centerFreqs(i))));
00131     fcoefs(i,5) = A2;
00132     fcoefs(i,6) = B0;
00133     fcoefs(i,7) = B1(centerFreqs(i), B(E(centerFreqs(i))));
00134     fcoefs(i,8) = B2( B(E(centerFreqs(i))));
00135     fcoefs(i,9) = gain(centerFreqs(i), B(E(centerFreqs(i))));
00136   }
00137 
00138   //Controls update
00139   stringstream channel,filter,ctrl;
00140   realvec b(1,3),a(1,3);
00141   channel << "Series/channel_0/";
00142   ctrl << channel.str() << "mrs_natural/inSamples";
00143   filterBank->setctrl(ctrl.str(), getctrl("mrs_natural/inSamples"));
00144   ctrl.str("");
00145   ctrl << channel.str() << "mrs_natural/inObservations";
00146   filterBank->setctrl(ctrl.str(), getctrl("mrs_natural/inObservations"));
00147   ctrl.str("");
00148   ctrl << channel.str() << "mrs_natural/onObservations";
00149   filterBank->setctrl(ctrl.str(), getctrl("mrs_natural/inObservations"));
00150   ctrl.str("");
00151   ctrl << channel.str() << "mrs_real/israte";
00152   filterBank->setctrl(ctrl.str(), getctrl("mrs_real/israte"));
00153   for(mrs_natural i = 0; i < numChannels; ++i) {
00154     //filter 0
00155     channel.str("");
00156     channel << "Series/channel_" << i << "/";
00157     filter.str("");
00158     filter << channel.str() << "Filter/filter_" << i << "_" << 0 << "/";
00159     ctrl.str("");
00160     ctrl << filter.str() << "mrs_natural/inSamples";
00161     filterBank->setctrl(ctrl.str(), getctrl("mrs_natural/inSamples"));
00162     ctrl.str("");
00163     ctrl << filter.str() << "mrs_natural/inObservations";
00164     filterBank->setctrl(ctrl.str(), getctrl("mrs_natural/inObservations"));
00165     ctrl.str("");
00166     ctrl << filter.str() << "mrs_natural/onObservations";
00167     filterBank->setctrl(ctrl.str(), getctrl("mrs_natural/inObservations"));
00168     ctrl.str("");
00169     ctrl << filter.str() << "mrs_real/israte";
00170     filterBank->setctrl(ctrl.str(), getctrl("mrs_real/israte"));
00171     a(0) = fcoefs(i,6);
00172     a(1) = fcoefs(i,7);
00173     a(2) = fcoefs(i,8);
00174     b(0) = fcoefs(i,0)/fcoefs(i,9);
00175     b(1) = fcoefs(i,1)/fcoefs(i,9);
00176     b(2) = fcoefs(i,5)/fcoefs(i,9);
00177     ctrl.str("");
00178     ctrl << filter.str() << "mrs_realvec/ncoeffs";
00179     filterBank->setctrl(ctrl.str(), b);
00180     ctrl.str("");
00181     ctrl << filter.str() << "mrs_realvec/dcoeffs";
00182     filterBank->setctrl(ctrl.str(), a);
00183     //filter 1
00184     filter.str("");
00185     filter << channel.str() << "Filter/filter_" << i << "_" << 1 << "/";
00186     a(0) = fcoefs(i,6);
00187     a(1) = fcoefs(i,7);
00188     a(2) = fcoefs(i,8);
00189     b(0) = fcoefs(i,0);
00190     b(1) = fcoefs(i,2);
00191     b(2) = fcoefs(i,5);
00192     ctrl.str("");
00193     ctrl << filter.str() << "mrs_realvec/ncoeffs";
00194     filterBank->setctrl(ctrl.str(), b);
00195     ctrl.str("");
00196     ctrl << filter.str() << "mrs_realvec/dcoeffs";
00197     filterBank->setctrl(ctrl.str(), a);
00198     //filter 2
00199     filter.str("");
00200     filter << channel.str() << "Filter/filter_" << i << "_" << 2 << "/";
00201     a(0) = fcoefs(i,6);
00202     a(1) = fcoefs(i,7);
00203     a(2) = fcoefs(i,8);
00204     b(0) = fcoefs(i,0);
00205     b(1) = fcoefs(i,3);
00206     b(2) = fcoefs(i,5);
00207     ctrl.str("");
00208     ctrl << filter.str() << "mrs_realvec/ncoeffs";
00209     filterBank->setctrl(ctrl.str(), b);
00210     ctrl.str("");
00211     ctrl << filter.str() << "mrs_realvec/dcoeffs";
00212     filterBank->setctrl(ctrl.str(), a);
00213     //filter 3
00214     filter.str("");
00215     filter << channel.str() << "Filter/filter_" << i << "_" << 3 << "/";
00216     a(0) = fcoefs(i,6);
00217     a(1) = fcoefs(i,7);
00218     a(2) = fcoefs(i,8);
00219     b(0) = fcoefs(i,0);
00220     b(1) = fcoefs(i,4);
00221     b(2) = fcoefs(i,5);
00222     ctrl.str("");
00223     ctrl << filter.str() << "mrs_realvec/ncoeffs";
00224     filterBank->setctrl(ctrl.str(), b);
00225     ctrl.str("");
00226     ctrl << filter.str() << "mrs_realvec/dcoeffs";
00227     filterBank->setctrl(ctrl.str(), a);
00228   }
00229   //update the whole filter bank
00230   filterBank->update();
00231 }
00232 
00233 
00234 mrs_real
00235 ERB::E(mrs_real x)
00236 {
00237   return pow(pow(x/EarQ, (mrs_real)order) + pow(minBW, (mrs_real)order), 1/(mrs_real)order);
00238 }
00239 
00240 mrs_real
00241 ERB::B(mrs_real x)
00242 {
00243   return 1.019f*2.0f*PI*x;
00244 }
00245 
00246 mrs_real
00247 ERB::B1(mrs_real x, mrs_real y)
00248 {
00249   return -2.0f*cos(2.0f*x*PI/fs)/exp(y/fs);
00250 }
00251 
00252 mrs_real
00253 ERB::B2(mrs_real x)
00254 {
00255   return exp(-2.0f*x/fs);
00256 }
00257 
00258 mrs_real
00259 ERB::A11(mrs_real x, mrs_real y)
00260 {
00261   return -(2.0f/fs*cos(2.0f*x*PI/fs)/exp(y/fs) + 2.0f*sqrt(3.0f+sqrt(8.0f))/fs*sin(2.0f*x*PI/fs)/exp(y/fs))/2.0f;
00262 }
00263 
00264 
00265 mrs_real
00266 ERB::A12(mrs_real x, mrs_real y)
00267 {
00268   return -(2.0f/fs*cos(2.0f*x*PI/fs)/exp(y/fs) - 2.0f*sqrt(3.0f+sqrt(8.0f))/fs*sin(2.0f*x*PI/fs)/exp(y/fs))/2.0f;
00269 }
00270 
00271 mrs_real
00272 ERB::A13(mrs_real x, mrs_real y)
00273 {
00274   return -(2.0f/fs*cos(2.0f*x*PI/fs)/exp(y/fs) + 2.0f*sqrt(3.0f-sqrt(8.0f))/fs*sin(2.0f*x*PI/fs)/exp(y/fs))/2.0f;
00275 }
00276 
00277 mrs_real
00278 ERB::A14(mrs_real x, mrs_real y)
00279 {
00280   return -(2.0f/fs*cos(2.0f*x*PI/fs)/exp(y/fs) - 2.0f*sqrt(3.0f-sqrt(8.0f))/fs*sin(2.0f*x*PI/fs)/exp(y/fs))/2.0f;
00281 }
00282 
00283 mrs_real
00284 ERB::gain(mrs_real x, mrs_real y)
00285 {
00286   mrs_real r1(2.0f*x*PI/fs);
00287   mrs_real r2(sqrt(3.0f-sqrt(8.0f)));
00288   mrs_real r3(sqrt(3.0f+sqrt(8.0f)));
00289   mrs_real z1r(cos(2.0f*r1)),z1i(sin(2.0f*r1));
00290   mrs_real z2r(cos(r1)),z2i(sin(r1));
00291   mrs_real z3r(-2.0f*z1r/fs),z3i(-2.0f*z1i/fs);
00292   mrs_real z4r(2.0f*exp(-y/fs)*z2r/fs),z4i(2.0f*exp(-y/fs)*z2i/fs);
00293   return abs(z3r + z4r*(cos(r1)-r2*sin(r1)),z3i + z4i*(cos(r1)-r2*sin(r1)))*
00294          abs(z3r + z4r*(cos(r1)+r2*sin(r1)),z3i + z4i*(cos(r1)+r2*sin(r1)))*
00295          abs(z3r + z4r*(cos(r1)-r3*sin(r1)),z3i + z4i*(cos(r1)-r3*sin(r1)))*
00296          abs(z3r + z4r*(cos(r1)+r3*sin(r1)),z3i + z4i*(cos(r1)+r3*sin(r1)))/
00297          pow(abs(-2.0f/exp(2.0f*y/fs) - 2.0f*z1r + 2.0f*(1.0f + z1r)/exp(y/fs),2.0f*z1i*(-1.0f + 1.0f/exp(y/fs))),4);
00298 }
00299 
00300 mrs_real
00301 ERB::abs(mrs_real r1, mrs_real r2)
00302 {
00303   return sqrt(r1*r1 + r2*r2);
00304 }
00305 
00306 void
00307 ERB::myProcess(realvec& in, realvec& out)
00308 {
00309   //checkFlow(in,out);
00310   //lmartins: if (mute_) return;
00311   if(getctrl("mrs_bool/mute")->to<mrs_bool>()) return;
00312 
00313   filterBank->process(in, out);
00314 }