Marsyas  0.6.0-alpha
/usr/src/RPM/BUILD/marsyas-0.6.0/src/marsyas/marsystems/AimGammatone.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 "AimGammatone.h"
00020 #include "../common_source.h"
00021 #include "../ERBTools.h"
00022 
00023 using std::ostringstream;
00024 using std::complex;
00025 using std::abs;
00026 using std::cout;
00027 using std::endl;
00028 
00029 using namespace Marsyas;
00030 
00031 AimGammatone::AimGammatone(mrs_string name):MarSystem("AimGammatone",name)
00032 {
00033   is_initialized = false;
00034   initialized_num_channels = 0;
00035   initialized_min_frequency = 0.0;
00036   initialized_max_frequency = 0.0;
00037   initialized_israte = 0.0;
00038 
00039   is_reset = false;
00040   reset_num_channels = 0;
00041 
00042   addControls();
00043 }
00044 
00045 
00046 AimGammatone::AimGammatone(const AimGammatone& a): MarSystem(a)
00047 {
00048   is_initialized = false;
00049   initialized_num_channels = 0;
00050   initialized_min_frequency = 0.0;
00051   initialized_max_frequency = 0.0;
00052   initialized_israte = 0.0;
00053 
00054   is_reset = false;
00055   reset_num_channels = 0;
00056 
00057   ctrl_num_channels_= getctrl("mrs_natural/num_channels");
00058   ctrl_min_frequency_ = getctrl("mrs_real/min_frequency");
00059   ctrl_max_frequency_ = getctrl("mrs_real/max_frequency");
00060 }
00061 
00062 AimGammatone::~AimGammatone()
00063 {
00064 }
00065 
00066 
00067 MarSystem*
00068 AimGammatone::clone() const
00069 {
00070   return new AimGammatone(*this);
00071 }
00072 
00073 void
00074 AimGammatone::addControls()
00075 {
00076   addControl("mrs_natural/num_channels", 200 , ctrl_num_channels_);
00077   addControl("mrs_real/min_frequency", 86.0 , ctrl_min_frequency_);
00078   addControl("mrs_real/max_frequency", 16000.0 , ctrl_max_frequency_);
00079 }
00080 
00081 void
00082 AimGammatone::myUpdate(MarControlPtr sender)
00083 {
00084   (void) sender;  //suppress warning of unused parameter(s)
00085 
00086   MRSDIAG("AimGammatone.cpp - AimGammatone:myUpdate");
00087   ctrl_onSamples_->setValue(ctrl_inSamples_, NOUPDATE);
00088   ctrl_osrate_->setValue(ctrl_israte_->to<mrs_real>());
00089   ctrl_onObsNames_->setValue("AimGammatone_" + ctrl_inObsNames_->to<mrs_string>() , NOUPDATE);
00090 
00091   // Need to have double the amount of channels, the first set of
00092   // channels are for the signals the second set of channels are for
00093   // the centre frequencies
00094   ctrl_onObservations_->setValue(ctrl_num_channels_->to<mrs_natural>() , NOUPDATE);
00095 
00096   //
00097   // Does the MarSystem need initialization?
00098   //
00099   if (initialized_num_channels != ctrl_num_channels_->to<mrs_natural>() ||
00100       initialized_min_frequency != ctrl_min_frequency_->to<mrs_real>() ||
00101       initialized_max_frequency != ctrl_max_frequency_->to<mrs_real>() ||
00102       initialized_israte != ctrl_israte_->to<mrs_real>()) {
00103     is_initialized = false;
00104   }
00105 
00106   if (!is_initialized) {
00107     InitializeInternal();
00108     is_initialized = true;
00109     initialized_num_channels = ctrl_num_channels_->to<mrs_natural>();
00110     initialized_min_frequency = ctrl_min_frequency_->to<mrs_real>();
00111     initialized_max_frequency = ctrl_max_frequency_->to<mrs_real>();
00112     initialized_israte = ctrl_israte_->to<mrs_real>();
00113   }
00114 
00115   //
00116   // Does the MarSystem need a reset?
00117   //
00118   if (reset_num_channels != ctrl_num_channels_->to<mrs_natural>()) {
00119     is_reset = false;
00120   }
00121 
00122   if (!is_reset) {
00123     ResetInternal();
00124     is_reset = true;
00125     reset_num_channels = ctrl_num_channels_->to<mrs_natural>();
00126   }
00127 
00128 }
00129 
00130 bool
00131 AimGammatone::InitializeInternal() {
00132   mrs_natural num_channels = ctrl_num_channels_->to<mrs_natural>();
00133   double min_frequency = ctrl_min_frequency_->to<mrs_real>();
00134   double max_frequency = ctrl_max_frequency_->to<mrs_real>();
00135 
00136   // Calculate number of channels, and centre frequencies
00137   double erb_max = ERBTools::Freq2ERB(max_frequency);
00138   double erb_min = ERBTools::Freq2ERB(min_frequency);
00139   double delta_erb = (erb_max - erb_min) / (num_channels - 1);
00140 
00141   centre_frequencies_.resize(num_channels);
00142   double erb_current = erb_min;
00143 
00144   for (int i = 0; i < num_channels; ++i) {
00145     centre_frequencies_[i] = ERBTools::ERB2Freq(erb_current);
00146     erb_current += delta_erb;
00147   }
00148 
00149   a_.resize(num_channels);
00150   b1_.resize(num_channels);
00151   b2_.resize(num_channels);
00152   b3_.resize(num_channels);
00153   b4_.resize(num_channels);
00154   state_1_.resize(num_channels);
00155   state_2_.resize(num_channels);
00156   state_3_.resize(num_channels);
00157   state_4_.resize(num_channels);
00158 
00159   for (int ch = 0; ch < num_channels; ++ch) {
00160     double cf = centre_frequencies_[ch];
00161     double erb = ERBTools::Freq2ERBw(cf);
00162 
00163     // Sample interval
00164     double dt = 1.0 / ctrl_israte_->to<mrs_real>();
00165 
00166     // Bandwidth parameter
00167     double b = 1.019 * 2.0 * PI * erb;
00168 
00169     // The following expressions are derived in Apple TR #35, "An
00170     // Efficient Implementation of the Patterson-Holdsworth Cochlear
00171     // Filter Bank" and used in Malcolm Slaney's auditory toolbox, where he
00172     // defines this alternaltive four stage cascade of second-order filters.
00173 
00174     // Calculate the gain:
00175     double cpt = cf * PI * dt;
00176     complex<double> exponent(0.0, 2.0 * cpt);
00177     complex<double> ec = exp(2.0 * exponent);
00178     complex<double> two_cf_pi_t(2.0 * cpt, 0.0);
00179     complex<double> two_pow(pow(2.0, (3.0 / 2.0)), 0.0);
00180     complex<double> p1 = -2.0 * ec * dt;
00181     complex<double> p2 = 2.0 * exp(-(b * dt) + exponent) * dt;
00182     complex<double> b_dt(b * dt, 0.0);
00183 
00184     double gain = abs(
00185                     (p1 + p2 * (cos(two_cf_pi_t) - sqrt(3.0 - two_pow) * sin(two_cf_pi_t)))
00186                     * (p1 + p2 * (cos(two_cf_pi_t) + sqrt(3.0 - two_pow) * sin(two_cf_pi_t)))
00187                     * (p1 + p2 * (cos(two_cf_pi_t) - sqrt(3.0 + two_pow) * sin(two_cf_pi_t)))
00188                     * (p1 + p2 * (cos(two_cf_pi_t) + sqrt(3.0 + two_pow) * sin(two_cf_pi_t)))
00189                     / pow((-2.0 / exp(2.0 * b_dt) - 2.0 * ec + 2.0 * (1.0 + ec)
00190                            / exp(b_dt)), 4));
00191 
00192     // The filter coefficients themselves:
00193     const int coeff_count = 3;
00194     a_[ch].resize(coeff_count, 0.0);
00195     b1_[ch].resize(coeff_count, 0.0);
00196     b2_[ch].resize(coeff_count, 0.0);
00197     b3_[ch].resize(coeff_count, 0.0);
00198     b4_[ch].resize(coeff_count, 0.0);
00199     state_1_[ch].resize(coeff_count, 0.0);
00200     state_2_[ch].resize(coeff_count, 0.0);
00201     state_3_[ch].resize(coeff_count, 0.0);
00202     state_4_[ch].resize(coeff_count, 0.0);
00203 
00204     double B0 = dt;
00205     double B2 = 0.0;
00206 
00207     double B11 = -(2.0 * dt * cos(2.0 * cf * PI * dt) / exp(b * dt)
00208                    + 2.0 * sqrt(3 + pow(2.0, 1.5)) * dt
00209                    * sin(2.0 * cf * PI * dt) / exp(b * dt)) / 2.0;
00210     double B12 = -(2.0 * dt * cos(2.0 * cf * PI * dt) / exp(b * dt)
00211                    - 2.0 * sqrt(3 + pow(2.0, 1.5)) * dt
00212                    * sin(2.0 * cf * PI * dt) / exp(b * dt)) / 2.0;
00213     double B13 = -(2.0 * dt * cos(2.0 * cf * PI * dt) / exp(b * dt)
00214                    + 2.0 * sqrt(3 - pow(2.0, 1.5)) * dt
00215                    * sin(2.0 * cf * PI * dt) / exp(b * dt)) / 2.0;
00216     double B14 = -(2.0 * dt * cos(2.0 * cf * PI * dt) / exp(b * dt)
00217                    - 2.0 * sqrt(3 - pow(2.0, 1.5)) * dt
00218                    * sin(2.0 * cf * PI * dt) / exp(b * dt)) / 2.0;
00219 
00220     a_[ch][0] = 1.0;
00221     a_[ch][1] = -2.0 * cos(2.0 * cf * PI * dt) / exp(b * dt);
00222     a_[ch][2] = exp(-2.0 * b * dt);
00223     b1_[ch][0] = B0 / gain;
00224     b1_[ch][1] = B11 / gain;
00225     b1_[ch][2] = B2 / gain;
00226     b2_[ch][0] = B0;
00227     b2_[ch][1] = B12;
00228     b2_[ch][2] = B2;
00229     b3_[ch][0] = B0;
00230     b3_[ch][1] = B13;
00231     b3_[ch][2] = B2;
00232     b4_[ch][0] = B0;
00233     b4_[ch][1] = B14;
00234     b4_[ch][2] = B2;
00235   }
00236   return true;
00237 }
00238 
00239 void
00240 AimGammatone::ResetInternal() {
00241   mrs_natural num_channels = ctrl_num_channels_->to<mrs_natural>();
00242 
00243   state_1_.resize(num_channels);
00244   state_2_.resize(num_channels);
00245   state_3_.resize(num_channels);
00246   state_4_.resize(num_channels);
00247   for (int i = 0; i < num_channels; ++i) {
00248     state_1_[i].resize(3, 0.0);
00249     state_2_[i].resize(3, 0.0);
00250     state_3_[i].resize(3, 0.0);
00251     state_4_[i].resize(3, 0.0);
00252   }
00253 }
00254 
00255 void
00256 AimGammatone::myProcess(realvec& in, realvec& out)
00257 {
00258   int audio_channel = 0;
00259 
00260   std::vector<std::vector<double> >::iterator b1 = b1_.begin();
00261   std::vector<std::vector<double> >::iterator b2 = b2_.begin();
00262   std::vector<std::vector<double> >::iterator b3 = b3_.begin();
00263   std::vector<std::vector<double> >::iterator b4 = b4_.begin();
00264   std::vector<std::vector<double> >::iterator a = a_.begin();
00265   std::vector<std::vector<double> >::iterator s1 = state_1_.begin();
00266   std::vector<std::vector<double> >::iterator s2 = state_2_.begin();
00267   std::vector<std::vector<double> >::iterator s3 = state_3_.begin();
00268   std::vector<std::vector<double> >::iterator s4 = state_4_.begin();
00269 
00270   // Temporary storage between filter stages
00271   std::vector<double> outbuff(ctrl_inSamples_->to<mrs_natural>());
00272   mrs_natural _channel_count = ctrl_num_channels_->to<mrs_natural>();
00273   mrs_natural _inSamples = ctrl_inSamples_->to<mrs_natural>();
00274 
00275   for (int ch = 0; ch < _channel_count;
00276        ++ch, ++b1, ++b2, ++b3, ++b4, ++a, ++s1, ++s2, ++s3, ++s4) {
00277     for (int i = 0; i < _inSamples; ++i) {
00278       // Direct-form-II IIR filter
00279       double inputsample = in(audio_channel, i);
00280       outbuff[i] = (*b1)[0] * inputsample + (*s1)[0];
00281       for (unsigned int stage = 1; stage < s1->size(); ++stage)
00282         (*s1)[stage - 1] = (*b1)[stage] * inputsample
00283                            - (*a)[stage] * outbuff[i] + (*s1)[stage];
00284     }
00285     for (int i = 0; i < _inSamples; ++i) {
00286       // Direct-form-II IIR filter
00287       double inputsample = outbuff[i];
00288       outbuff[i] = (*b2)[0] * inputsample + (*s2)[0];
00289       for (unsigned int stage = 1; stage < s2->size(); ++stage)
00290         (*s2)[stage - 1] = (*b2)[stage] * inputsample
00291                            - (*a)[stage] * outbuff[i] + (*s2)[stage];
00292     }
00293     for (int i = 0; i < _inSamples; ++i) {
00294       // Direct-form-II IIR filter
00295       double inputsample = outbuff[i];
00296       outbuff[i] = (*b3)[0] * inputsample + (*s3)[0];
00297       for (unsigned int stage = 1; stage < s3->size(); ++stage)
00298         (*s3)[stage - 1] = (*b3)[stage] * inputsample
00299                            - (*a)[stage] * outbuff[i] + (*s3)[stage];
00300     }
00301     for (int i = 0; i < _inSamples; ++i) {
00302       // Direct-form-II IIR filter
00303       double inputsample = outbuff[i];
00304       outbuff[i] = (*b4)[0] * inputsample + (*s4)[0];
00305       for (unsigned int stage = 1; stage < s4->size(); ++stage)
00306         (*s4)[stage - 1] = (*b4)[stage] * inputsample
00307                            - (*a)[stage] * outbuff[i] + (*s4)[stage];
00308       out(ch, i) = outbuff[i];
00309     }
00310   }
00311 
00312   // Copy over the centre frequencies to the second half of the observations
00313   /* for (t = 0; t < ctrl_inSamples_->to<mrs_natural>(); t++) {
00314     for (o = 0; o < _channel_count; o++) {
00315       out(o + _channel_count, t) = centre_frequencies_[o];
00316     }
00317   }
00318   */
00319 
00320 }