Marsyas
0.6.0-alpha
|
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 }