tomwalters@276: // Copyright 2009-2010, Thomas Walters
tomwalters@276: //
tomwalters@276: // AIM-C: A C++ implementation of the Auditory Image Model
tomwalters@276: // http://www.acousticscale.org/AIMC
tomwalters@276: //
tomwalters@276: // This program is free software: you can redistribute it and/or modify
tomwalters@276: // it under the terms of the GNU General Public License as published by
tomwalters@276: // the Free Software Foundation, either version 3 of the License, or
tomwalters@276: // (at your option) any later version.
tomwalters@276: //
tomwalters@276: // This program is distributed in the hope that it will be useful,
tomwalters@276: // but WITHOUT ANY WARRANTY; without even the implied warranty of
tomwalters@276: // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
tomwalters@276: // GNU General Public License for more details.
tomwalters@276: //
tomwalters@276: // You should have received a copy of the GNU General Public License
tomwalters@276: // along with this program.  If not, see <http://www.gnu.org/licenses/>.
tomwalters@276: 
tomwalters@276: /*! \file
tomwalters@276:  *  \brief Slaney's gammatone filterbank
tomwalters@276:  *
tomwalters@276:  *  \author Thomas Walters <tom@acousticscale.org>
tomwalters@276:  *  \date created 2009/11/13
tomwalters@276:  *  \version \$Id$
tomwalters@276:  */
tomwalters@276: 
tomwalters@287: #include <cmath>
tomwalters@276: #include <complex>
tomwalters@277: #include "Support/ERBTools.h"
tomwalters@276: 
tomwalters@276: #include "Modules/BMM/ModuleGammatone.h"
tomwalters@276: 
tomwalters@276: namespace aimc {
tomwalters@277: using std::vector;
tomwalters@276: using std::complex;
tomwalters@276: ModuleGammatone::ModuleGammatone(Parameters *params) : Module(params) {
tomwalters@277:   module_identifier_ = "gt";
tomwalters@276:   module_type_ = "bmm";
tomwalters@276:   module_description_ = "Gammatone filterbank (Slaney's IIR gammatone)";
tomwalters@276:   module_version_ = "$Id$";
tomwalters@276: 
tomwalters@276:   num_channels_ = parameters_->DefaultInt("gtfb.channel_count", 200);
tomwalters@276:   min_frequency_ = parameters_->DefaultFloat("gtfb.min_frequency", 86.0f);
tomwalters@276:   max_frequency_ = parameters_->DefaultFloat("gtfb.max_frequency", 16000.0f);
tomwalters@276: }
tomwalters@276: 
tomwalters@277: ModuleGammatone::~ModuleGammatone() {
tomwalters@277: }
tomwalters@277: 
tomwalters@277: void ModuleGammatone::ResetInternal() {
tomwalters@277:   state_.resize(num_channels_);
tomwalters@277:   for (int i = 0; i < num_channels_; ++i) {
tomwalters@277:     state_[i].resize(9, 0.0f);
tomwalters@277:   }
tomwalters@277: }
tomwalters@277: 
tomwalters@276: bool ModuleGammatone::InitializeInternal(const SignalBank& input) {
tomwalters@276:   // Calculate number of channels, and centre frequencies
tomwalters@277:   float erb_max = ERBTools::Freq2ERB(max_frequency_);
tomwalters@277:   float erb_min = ERBTools::Freq2ERB(min_frequency_);
tomwalters@277:   float delta_erb = (erb_max - erb_min) / (num_channels_ - 1);
tomwalters@277: 
tomwalters@277:   centre_frequencies_.resize(num_channels_);
tomwalters@277:   float erb_current = erb_min;
tomwalters@277: 
tomwalters@277:   for (int i = 0; i < num_channels_; ++i) {
tomwalters@280:     centre_frequencies_[i] = ERBTools::ERB2Freq(erb_current);
tomwalters@280:     erb_current += delta_erb;
tomwalters@277:   }
tomwalters@276: 
tomwalters@276:   forward_.resize(num_channels_);
tomwalters@277:   back_.resize(num_channels_);
tomwalters@277:   state_.resize(num_channels_);
tomwalters@276: 
tomwalters@276:   for (int ch = 0; ch < num_channels_; ++ch) {
tomwalters@277:     float cf = centre_frequencies_[ch];
tomwalters@277:     float erb = ERBTools::Freq2ERBw(cf);
tomwalters@276: 
tomwalters@276:     // Sample interval
tomwalters@277:     float dt = 1.0f / input.sample_rate();
tomwalters@276: 
tomwalters@276:     // Bandwidth parameter
tomwalters@276:     float b = 1.019f * 2.0f * M_PI * erb;
tomwalters@276: 
tomwalters@280:     // All of the following expressions are derived in Apple TR #35, "An
tomwalters@280:     // Efficient Implementation of the Patterson-Holdsworth Cochlear
tomwalters@276:     // Filter Bank".
tomwalters@276: 
tomwalters@276:     // Calculate the gain:
tomwalters@277:     float cpt = cf * M_PI * dt;
tomwalters@277:     complex<float> exponent(0.0f, 2.0f * cpt);
tomwalters@277:     complex<float> ec = exp(2.0f * exponent);
tomwalters@277:     complex<float> two_cf_pi_t(2.0f * cpt, 0.0f);
tomwalters@277:     complex<float> two_pow(pow(2.0f, (3.0f / 2.0f)), 0.0f);
tomwalters@277:     complex<float> p = -2.0f * ec * dt
tomwalters@277:                        + 2.0f * exp(-(b * dt) + exponent) * dt;
tomwalters@277:     complex<float> b_dt(b * dt, 0.0f);
tomwalters@276: 
tomwalters@276:     float gain = abs(
tomwalters@277:       (p * (cos(two_cf_pi_t) - sqrt(3.0f - two_pow) * sin(two_cf_pi_t)))
tomwalters@277:       * (p * (cos(two_cf_pi_t) + sqrt(3.0f - two_pow) * sin(two_cf_pi_t)))
tomwalters@277:       * (p * (cos(two_cf_pi_t) - sqrt(3.0f + two_pow) * sin(two_cf_pi_t)))
tomwalters@277:       * (p * (cos(two_cf_pi_t) + sqrt(3.0f + two_pow) * sin(two_cf_pi_t)))
tomwalters@277:       / pow(-2.0f / exp(2.0f * b_dt) - 2.0f * ec + 2.0f * (1.0f + ec)
tomwalters@277:             / exp(b_dt), 4.0f));
tomwalters@276: 
tomwalters@276:     // The filter coefficients themselves:
tomwalters@277:     const int coeff_count = 9;
tomwalters@277:     forward_[ch].resize(coeff_count, 0.0f);
tomwalters@277:     back_[ch].resize(coeff_count, 0.0f);
tomwalters@277:     state_[ch].resize(coeff_count, 0.0f);
tomwalters@276: 
tomwalters@277:     forward_[ch][0] = pow(dt, 4.0f) / gain;
tomwalters@277:     forward_[ch][1] = (-4.0f * pow(dt, 4.0f) * cos(2.0f * cpt)
tomwalters@276:                       / exp(b * dt) / gain);
tomwalters@277:     forward_[ch][2] = (6.0f * pow(dt, 4.0f) * cos(4.0f * cpt)
tomwalters@276:                       / exp(2.0f * b * dt) / gain);
tomwalters@277:     forward_[ch][3] = (-4.0f * pow(dt, 4.0f) * cos(6.0f * cpt)
tomwalters@276:                       / exp(3.0f * b * dt) / gain);
tomwalters@277:     forward_[ch][4] = (pow(dt, 4.0f) * cos(8.0f * cpt)
tomwalters@276:                       / exp(4.0f * b * dt) / gain);
tomwalters@277:     // Note: the remainder of the forward vector is zero-padded
tomwalters@276: 
tomwalters@277:     back_[ch][0] = 1.0f;
tomwalters@277:     back_[ch][1] = -8.0f * cos(2.0f * cpt) / exp(b * dt);
tomwalters@277:     back_[ch][2] = (4.0f * (4.0f + 3.0f * cos(4.0f * cpt))
tomwalters@277:                     / exp(2.0f * b * dt));
tomwalters@277:     back_[ch][3] = (-8.0f * (6.0f * cos(2.0f * cpt) + cos(6.0f * cpt))
tomwalters@277:                     / exp(3.0f * b * dt));
tomwalters@277:     back_[ch][4] = (2.0f * (18.0f + 16.0f * cos(4.0f * cpt) + cos(8.0f * cpt))
tomwalters@277:                     / exp(4.0f * b * dt));
tomwalters@277:     back_[ch][5] = (-8.0f * (6.0f * cos(2.0f * cpt) + cos(6.0f * cpt))
tomwalters@277:                     / exp(5.0f * b * dt));
tomwalters@277:     back_[ch][6] = (4.0f * (4.0f + 3.0f * cos(4.0f * cpt))
tomwalters@277:                     / exp(6.0f * b * dt));
tomwalters@277:     back_[ch][7] = -8.0f * cos(2.0f * cpt) / exp(7.0f * b * dt);
tomwalters@277:     back_[ch][8] = exp(-8.0f * b * dt);
tomwalters@276:   }
tomwalters@277:   output_.Initialize(num_channels_,
tomwalters@277:                      input.buffer_length(),
tomwalters@277:                      input.sample_rate());
tomwalters@277:   return true;
tomwalters@276: }
tomwalters@277: 
tomwalters@277: void ModuleGammatone::Process(const SignalBank &input) {
tomwalters@277:   output_.set_start_time(input.start_time());
tomwalters@277:   int audio_channel = 0;
tomwalters@277: 
tomwalters@277:   vector<vector<float> >::iterator b = forward_.begin();
tomwalters@277:   vector<vector<float> >::iterator a = back_.begin();
tomwalters@277:   vector<vector<float> >::iterator s = state_.begin();
tomwalters@277: 
tomwalters@277:   for (int ch = 0; ch < num_channels_; ++ch, ++a, ++b, ++s) {
tomwalters@277:     for (int i = 0; i < input.buffer_length(); ++i) {
tomwalters@280:       // Direct-form-II IIR filter
tomwalters@277:       float in = input.sample(audio_channel, i);
tomwalters@277:       float out = (*b)[0] * in + (*s)[0];
tomwalters@277:       for (unsigned int stage = 1; stage < s->size(); ++stage)
tomwalters@277:         (*s)[stage - 1] = (*b)[stage] * in - (*a)[stage] * out + (*s)[stage];
tomwalters@277:       output_.set_sample(ch, i, out);
tomwalters@277:     }
tomwalters@277:   }
tomwalters@277:   PushOutput();
tomwalters@277: }
tomwalters@277: 
tomwalters@280: }  // namespace aimc