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@288:   state_1_.resize(num_channels_);
tomwalters@288:   state_2_.resize(num_channels_);
tomwalters@288:   state_3_.resize(num_channels_);
tomwalters@288:   state_4_.resize(num_channels_);
tomwalters@277:   for (int i = 0; i < num_channels_; ++i) {
tomwalters@288:     state_1_[i].resize(3, 0.0f);
tomwalters@288:     state_2_[i].resize(3, 0.0f);
tomwalters@288:     state_3_[i].resize(3, 0.0f);
tomwalters@288:     state_4_[i].resize(3, 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@288:   output_.Initialize(num_channels_,
tomwalters@288:                      input.buffer_length(),
tomwalters@288:                      input.sample_rate());
tomwalters@288: 
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@288:     output_.set_centre_frequency(i, centre_frequencies_[i]);
tomwalters@277:   }
tomwalters@276: 
tomwalters@288:   a_.resize(num_channels_);
tomwalters@288:   b1_.resize(num_channels_);
tomwalters@288:   b2_.resize(num_channels_);
tomwalters@288:   b3_.resize(num_channels_);
tomwalters@288:   b4_.resize(num_channels_);
tomwalters@288:   state_1_.resize(num_channels_);
tomwalters@288:   state_2_.resize(num_channels_);
tomwalters@288:   state_3_.resize(num_channels_);
tomwalters@288:   state_4_.resize(num_channels_);
tomwalters@276: 
tomwalters@276:   for (int ch = 0; ch < num_channels_; ++ch) {
tomwalters@288:     double cf = centre_frequencies_[ch];
tomwalters@288:     double erb = ERBTools::Freq2ERBw(cf);
tomwalters@290:     // LOG_INFO("%e", erb);
tomwalters@276: 
tomwalters@276:     // Sample interval
tomwalters@288:     double dt = 1.0f / input.sample_rate();
tomwalters@276: 
tomwalters@276:     // Bandwidth parameter
tomwalters@288:     double b = 1.019f * 2.0f * M_PI * erb;
tomwalters@276: 
tomwalters@288:     // The following expressions are derived in Apple TR #35, "An
tomwalters@280:     // Efficient Implementation of the Patterson-Holdsworth Cochlear
tomwalters@288:     // Filter Bank" and used in Malcolm Slaney's auditory toolbox, where he
tomwalters@288:     // defines this alternaltive four stage cascade of second-order filters.
tomwalters@276: 
tomwalters@276:     // Calculate the gain:
tomwalters@288:     double cpt = cf * M_PI * dt;
tomwalters@288:     complex<double> exponent(0.0, 2.0 * cpt);
tomwalters@288:     complex<double> ec = exp(2.0 * exponent);
tomwalters@288:     complex<double> two_cf_pi_t(2.0 * cpt, 0.0);
tomwalters@288:     complex<double> two_pow(pow(2.0, (3.0 / 2.0)), 0.0);
tomwalters@290:     complex<double> p1 = -2.0 * ec * dt;
tomwalters@290:     complex<double> p2 = 2.0 * exp(-(b * dt) + exponent) * dt;
tomwalters@288:     complex<double> b_dt(b * dt, 0.0);
tomwalters@276: 
tomwalters@288:     double gain = abs(
tomwalters@290:       (p1 + p2 * (cos(two_cf_pi_t) - sqrt(3.0 - two_pow) * sin(two_cf_pi_t)))
tomwalters@290:       * (p1 + p2 * (cos(two_cf_pi_t) + sqrt(3.0 - two_pow) * sin(two_cf_pi_t)))
tomwalters@290:       * (p1 + p2 * (cos(two_cf_pi_t) - sqrt(3.0 + two_pow) * sin(two_cf_pi_t)))
tomwalters@290:       * (p1 + p2 * (cos(two_cf_pi_t) + sqrt(3.0 + two_pow) * sin(two_cf_pi_t)))
tomwalters@290:       / pow((-2.0 / exp(2.0 * b_dt) - 2.0 * ec + 2.0 * (1.0 + ec)
tomwalters@290:             / exp(b_dt)), 4));
tomwalters@290:     LOG_INFO("%e", gain);
tomwalters@276: 
tomwalters@276:     // The filter coefficients themselves:
tomwalters@288:     const int coeff_count = 3;
tomwalters@288:     a_[ch].resize(coeff_count, 0.0f);
tomwalters@288:     b1_[ch].resize(coeff_count, 0.0f);
tomwalters@288:     b2_[ch].resize(coeff_count, 0.0f);
tomwalters@288:     b3_[ch].resize(coeff_count, 0.0f);
tomwalters@288:     b4_[ch].resize(coeff_count, 0.0f);
tomwalters@288:     state_1_[ch].resize(coeff_count, 0.0f);
tomwalters@288:     state_2_[ch].resize(coeff_count, 0.0f);
tomwalters@288:     state_3_[ch].resize(coeff_count, 0.0f);
tomwalters@288:     state_4_[ch].resize(coeff_count, 0.0f);
tomwalters@276: 
tomwalters@288:     double B0 = dt;
tomwalters@288:     double B2 = 0.0f;
tomwalters@276: 
tomwalters@288:     double B11 = -(2.0f * dt * cos(2.0f * cf * M_PI * dt) / exp(b * dt)
tomwalters@288:                    + 2.0f * sqrt(3 + pow(2.0f, 1.5f)) * dt
tomwalters@288:                        * sin(2.0f * cf * M_PI * dt) / exp(b * dt)) / 2.0f;
tomwalters@288:     double B12 = -(2.0f * dt * cos(2.0f * cf * M_PI * dt) / exp(b * dt)
tomwalters@288:                    - 2.0f * sqrt(3 + pow(2.0f, 1.5f)) * dt
tomwalters@288:                        * sin(2.0f * cf * M_PI * dt) / exp(b * dt)) / 2.0f;
tomwalters@288:     double B13 = -(2.0f * dt * cos(2.0f * cf * M_PI * dt) / exp(b * dt)
tomwalters@288:                    + 2.0f * sqrt(3 - pow(2.0f, 1.5f)) * dt
tomwalters@288:                        * sin(2.0f * cf * M_PI * dt) / exp(b * dt)) / 2.0f;
tomwalters@288:     double B14 = -(2.0f * dt * cos(2.0f * cf * M_PI * dt) / exp(b * dt)
tomwalters@288:                    - 2.0f * sqrt(3 - pow(2.0f, 1.5f)) * dt
tomwalters@289:                        * sin(2.0f * cf * M_PI * dt) / exp(b * dt)) / 2.0f;
tomwalters@288: 
tomwalters@288:     a_[ch][0] = 1.0f;
tomwalters@288:     a_[ch][1] = -2.0f * cos(2.0f * cf * M_PI * dt) / exp(b * dt);
tomwalters@288:     a_[ch][2] = exp(-2.0f * b * dt);
tomwalters@288:     b1_[ch][0] = B0 / gain;
tomwalters@288:     b1_[ch][1] = B11 / gain;
tomwalters@288:     b1_[ch][2] = B2 / gain;
tomwalters@288:     b2_[ch][0] = B0;
tomwalters@288:     b2_[ch][1] = B12;
tomwalters@288:     b2_[ch][2] = B2;
tomwalters@288:     b3_[ch][0] = B0;
tomwalters@288:     b3_[ch][1] = B13;
tomwalters@288:     b3_[ch][2] = B2;
tomwalters@288:     b4_[ch][0] = B0;
tomwalters@288:     b4_[ch][1] = B14;
tomwalters@288:     b4_[ch][2] = B2;
tomwalters@276:   }
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@288:   vector<vector<double> >::iterator b1 = b1_.begin();
tomwalters@288:   vector<vector<double> >::iterator b2 = b2_.begin();
tomwalters@288:   vector<vector<double> >::iterator b3 = b3_.begin();
tomwalters@288:   vector<vector<double> >::iterator b4 = b4_.begin();
tomwalters@288:   vector<vector<double> >::iterator a = a_.begin();
tomwalters@288:   vector<vector<double> >::iterator s1 = state_1_.begin();
tomwalters@288:   vector<vector<double> >::iterator s2 = state_2_.begin();
tomwalters@288:   vector<vector<double> >::iterator s3 = state_3_.begin();
tomwalters@288:   vector<vector<double> >::iterator s4 = state_4_.begin();
tomwalters@277: 
tomwalters@288:   // Temporary storage between filter stages
tomwalters@288:   vector<double> out(input.buffer_length());
tomwalters@288:   for (int ch = 0; ch < num_channels_;
tomwalters@288:        ++ch, ++b1, ++b2, ++b3, ++b4, ++a, ++s1, ++s2, ++s3, ++s4) {
tomwalters@277:     for (int i = 0; i < input.buffer_length(); ++i) {
tomwalters@280:       // Direct-form-II IIR filter
tomwalters@288:       double in = input.sample(audio_channel, i);
tomwalters@288:       out[i] = (*b1)[0] * in + (*s1)[0];
tomwalters@288:       for (unsigned int stage = 1; stage < s1->size(); ++stage)
tomwalters@288:         (*s1)[stage - 1] = (*b1)[stage] * in
tomwalters@288:                            - (*a)[stage] * out[i] + (*s1)[stage];
tomwalters@288:     }
tomwalters@288:     for (int i = 0; i < input.buffer_length(); ++i) {
tomwalters@288:       // Direct-form-II IIR filter
tomwalters@288:       double in = out[i];
tomwalters@288:       out[i] = (*b2)[0] * in + (*s2)[0];
tomwalters@288:       for (unsigned int stage = 1; stage < s2->size(); ++stage)
tomwalters@288:         (*s2)[stage - 1] = (*b2)[stage] * in
tomwalters@288:                            - (*a)[stage] * out[i] + (*s2)[stage];
tomwalters@288:     }
tomwalters@288:     for (int i = 0; i < input.buffer_length(); ++i) {
tomwalters@288:       // Direct-form-II IIR filter
tomwalters@288:       double in = out[i];
tomwalters@288:       out[i] = (*b3)[0] * in + (*s3)[0];
tomwalters@288:       for (unsigned int stage = 1; stage < s3->size(); ++stage)
tomwalters@288:         (*s3)[stage - 1] = (*b3)[stage] * in
tomwalters@288:                            - (*a)[stage] * out[i] + (*s3)[stage];
tomwalters@288:     }
tomwalters@288:     for (int i = 0; i < input.buffer_length(); ++i) {
tomwalters@288:       // Direct-form-II IIR filter
tomwalters@288:       double in = out[i];
tomwalters@288:       out[i] = (*b4)[0] * in + (*s4)[0];
tomwalters@288:       for (unsigned int stage = 1; stage < s4->size(); ++stage)
tomwalters@288:         (*s4)[stage - 1] = (*b4)[stage] * in
tomwalters@288:                            - (*a)[stage] * out[i] + (*s4)[stage];
tomwalters@288:       output_.set_sample(ch, i, out[i]);
tomwalters@277:     }
tomwalters@277:   }
tomwalters@277:   PushOutput();
tomwalters@277: }
tomwalters@277: 
tomwalters@280: }  // namespace aimc