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@276: #include <complex>
tomwalters@276: 
tomwalters@276: #include "Modules/BMM/ModuleGammatone.h"
tomwalters@276: 
tomwalters@276: namespace aimc {
tomwalters@276: using std::complex;
tomwalters@276: ModuleGammatone::ModuleGammatone(Parameters *params) : Module(params) {
tomwalters@276:   module_identifier_ = "gtfb";
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@276: bool ModuleGammatone::InitializeInternal(const SignalBank& input) {
tomwalters@276:   // Calculate number of channels, and centre frequencies
tomwalters@276: 
tomwalters@276:   forward_.resize(num_channels_);
tomwalters@276:   feedback_.resize(num_channels_);
tomwalters@276: 
tomwalters@276:   for (int ch = 0; ch < num_channels_; ++ch) {
tomwalters@276:     float erb = Freq2ERBw(cf)
tomwalters@276: 
tomwalters@276:     // Sample interval
tomwalters@276:     float dt = 1.0f / fs;
tomwalters@276: 
tomwalters@276:     // Bandwidth parameter
tomwalters@276:     float b = 1.019f * 2.0f * M_PI * erb;
tomwalters@276: 
tomwalters@276:     // All of the following expressions are derived in Apple TR #35, "An 
tomwalters@276:     // Efficient Implementation of the Patterson-Holdsworth Cochlear 
tomwalters@276:     // Filter Bank".
tomwalters@276: 
tomwalters@276:     // Calculate the gain:
tomwalters@276:     complex<float> exponent(0.0f, 2.0f * cf * M_PI * T);
tomwalters@276:     complex<float> ec = exp(2.0f * complex_exponent);
tomwalters@276:     complex<float> two_cf_pi_t(2.0f * cf * M_PI * T, 0.0f);
tomwalters@276:     complex<float> two_pow(pow(2.0f, (3.0f/2.0f)), 0.0f);
tomwalters@276: 
tomwalters@276:     complex<float> p = -2.0f * ec * T + 2.0f * exp(-(B * T) + exponent) * dt;
tomwalters@276: 
tomwalters@276:     float gain = abs(
tomwalters@276:       (part1 * (cos(two_cf_pi_t) - sqrt(3.0f - twopow) * sin(two_cf_pi_t)))
tomwalters@276:       * (part1 * (cos(two_cf_pi_t) + sqrt(3.0f - twopow) * sin(two_cf_pi_t)))
tomwalters@276:       * (part1 * (cos(two_cf_pi_t) - sqrt(3.0f + twopow) * sin(two_cf_pi_t)))
tomwalters@276:       * (part1 * (cos(two_cf_pi_t) + sqrt(3.0f + twopow) * sin(two_cf_pi_t)))
tomwalters@276:       / pow(-2.0f / exp(2.0f * b * dt) - 2.0f * ec + 2.0f * (1.0f + ec)
tomwalters@276:             / exp(b * dt), 4.0f));
tomwalters@276: 
tomwalters@276:     // The filter coefficients themselves:
tomwalters@276:     forward[ch].resize(5, 0.0f);
tomwalters@276:     feedback[ch].resize(9, 0.0f);
tomwalters@276: 
tomwalters@276:     forward[ch][0] = pow(T, 4.0f) / gain;
tomwalters@276:     forward[ch][1] = (-4.0f * pow(T, 4.0f) * cos(2.0f * cf * M_PI * dt)
tomwalters@276:                       / exp(b * dt) / gain);
tomwalters@276:     forward[ch][2] = (6.0f * pow(T, 4.0f) * cos(4.0f * cf * M_PI * dt)
tomwalters@276:                       / exp(2.0f * b * dt) / gain);
tomwalters@276:     forward[ch][3] = (-4.0f * pow(T, 4.0f) * cos(6.0f * cf * M_PI * dt)
tomwalters@276:                       / exp(3.0f * b * dt) / gain);
tomwalters@276:     forward[ch][4] = (pow(T,4.0f) * cos(8.0f * cf * M_PI * dt)
tomwalters@276:                       / exp(4.0f * b * dt) / gain);
tomwalters@276: 
tomwalters@276:     feedback[ch][0] = 1.0f;
tomwalters@276:     feedback[ch][1] = -8.0f * cos(2.0f * cf * M_PI * T) / exp(b * dt);
tomwalters@276:     feedback[ch][2] = (4.0f * (4.0f + 3.0f * cos(4.0f * cf * M_PI * dt))
tomwalters@276:                        / exp(2.0f * b * dt));
tomwalters@276:     feedback[ch][3] = (-8.0f * (6.0f * cos(2.0f * cf * M_PI * dt)
tomwalters@276:                                 + cos(6.0f * cf * M_PI * dt))
tomwalters@276:                        / exp(3.0f * b * dt));
tomwalters@276:     feedback[ch][4] = (2.0f * (18.0f + 16.0f * cos(4.0f * cf * M_PI * dt)
tomwalters@276:                                + cos(8.0f * cf * M_PI * dt))
tomwalters@276:                        / exp(4.0f * b * dt));
tomwalters@276:     feedback[ch][5] = (-8.0f * (6.0f * cos(2.0f * cf * M_PI * T)
tomwalters@276:                                 + cos(6.0f * cf * M_PI * T))
tomwalters@276:                        / exp(5.0f * B * T));
tomwalters@276:     feedback[ch][6] = (4.0f * (4.0f + 3.0f * cos(4.0f * cf * M_PI * T))
tomwalters@276:                        / exp(6.0f * B * T));
tomwalters@276:     feedback[ch][7] = -8.0f * cos(2.0f * cf * M_PI * dt) / exp(7.0f * b * dt);
tomwalters@276:     feedback[ch][8] = exp(-8.0f * b * dt);
tomwalters@276:   }
tomwalters@276: }
tomwalters@276: } // namespace aimc