changeset 5:3c782dec2fc0

- Ported over HTK file output
- Added some more meat to the Slaney IIR gammatone implementation
- Ported over the AIM-MAT sf2003 parabola strobe algorithm
- Finished making the SAI implementation compile
- Ported over the strobe list class (now uses STL deques internally)

--- a/SConstruct	Tue Feb 16 18:40:55 2010 +0000
+++ b/SConstruct	Thu Feb 18 16:55:40 2010 +0000
@@ -97,10 +97,13 @@
                   'Support/Parameters.cc',
                   'Support/Module.cc',
                   'Modules/Input/ModuleFileInput.cc',
+                  'Modules/BMM/ModuleGammatone.cc',
                   'Modules/BMM/ModulePZFC.cc',
                   'Modules/NAP/ModuleHCL.cc',
-                  #'Modules/SAI/ModuleSAI.cc',
-                  'Modules/Features/ModuleGaussians.cc']
+                  'Modules/Strobes/ModuleParabola.cc',
+                  'Modules/SAI/ModuleSAI.cc',
+                  'Modules/Features/ModuleGaussians.cc',
+                  'Modules/Output/FileOutputHTK.cc']
     
 if not target_platform == 'win32':
   # On windows, utf support is builtin for SimpleIni

--- a/src/Modules/BMM/ModuleGammatone.cc	Tue Feb 16 18:40:55 2010 +0000
+++ b/src/Modules/BMM/ModuleGammatone.cc	Thu Feb 18 16:55:40 2010 +0000
@@ -25,13 +25,16 @@
  */
 
 #include <complex>
+#include <math.h>
+#include "Support/ERBTools.h"
 
 #include "Modules/BMM/ModuleGammatone.h"
 
 namespace aimc {
+using std::vector;
 using std::complex;
 ModuleGammatone::ModuleGammatone(Parameters *params) : Module(params) {
-  module_identifier_ = "gtfb";
+  module_identifier_ = "gt";
   module_type_ = "bmm";
   module_description_ = "Gammatone filterbank (Slaney's IIR gammatone)";
   module_version_ = "$Id$";
@@ -41,17 +44,40 @@
   max_frequency_ = parameters_->DefaultFloat("gtfb.max_frequency", 16000.0f);
 }
 
+ModuleGammatone::~ModuleGammatone() {
+}
+
+void ModuleGammatone::ResetInternal() {
+  state_.resize(num_channels_);
+  for (int i = 0; i < num_channels_; ++i) {
+    state_[i].resize(9, 0.0f);
+  }
+}
+
 bool ModuleGammatone::InitializeInternal(const SignalBank& input) {
   // Calculate number of channels, and centre frequencies
+  float erb_max = ERBTools::Freq2ERB(max_frequency_);
+  float erb_min = ERBTools::Freq2ERB(min_frequency_);
+  float delta_erb = (erb_max - erb_min) / (num_channels_ - 1);
+
+  centre_frequencies_.resize(num_channels_);
+  float erb_current = erb_min;
+
+  for (int i = 0; i < num_channels_; ++i) {
+   centre_frequencies_[i] = ERBTools::ERB2Freq(erb_current);
+   erb_current += delta_erb;
+  }
 
   forward_.resize(num_channels_);
-  feedback_.resize(num_channels_);
+  back_.resize(num_channels_);
+  state_.resize(num_channels_);
 
   for (int ch = 0; ch < num_channels_; ++ch) {
-    float erb = Freq2ERBw(cf)
+    float cf = centre_frequencies_[ch];
+    float erb = ERBTools::Freq2ERBw(cf);
 
     // Sample interval
-    float dt = 1.0f / fs;
+    float dt = 1.0f / input.sample_rate();
 
     // Bandwidth parameter
     float b = 1.019f * 2.0f * M_PI * erb;
@@ -61,52 +87,80 @@
     // Filter Bank".
 
     // Calculate the gain:
-    complex<float> exponent(0.0f, 2.0f * cf * M_PI * T);
-    complex<float> ec = exp(2.0f * complex_exponent);
-    complex<float> two_cf_pi_t(2.0f * cf * M_PI * T, 0.0f);
-    complex<float> two_pow(pow(2.0f, (3.0f/2.0f)), 0.0f);
-
-    complex<float> p = -2.0f * ec * T + 2.0f * exp(-(B * T) + exponent) * dt;
+    float cpt = cf * M_PI * dt;
+    complex<float> exponent(0.0f, 2.0f * cpt);
+    complex<float> ec = exp(2.0f * exponent);
+    complex<float> two_cf_pi_t(2.0f * cpt, 0.0f);
+    complex<float> two_pow(pow(2.0f, (3.0f / 2.0f)), 0.0f);
+    complex<float> p = -2.0f * ec * dt
+                       + 2.0f * exp(-(b * dt) + exponent) * dt;
+    complex<float> b_dt(b * dt, 0.0f);
 
     float gain = abs(
-      (part1 * (cos(two_cf_pi_t) - sqrt(3.0f - twopow) * sin(two_cf_pi_t)))
-      * (part1 * (cos(two_cf_pi_t) + sqrt(3.0f - twopow) * sin(two_cf_pi_t)))
-      * (part1 * (cos(two_cf_pi_t) - sqrt(3.0f + twopow) * sin(two_cf_pi_t)))
-      * (part1 * (cos(two_cf_pi_t) + sqrt(3.0f + twopow) * sin(two_cf_pi_t)))
-      / pow(-2.0f / exp(2.0f * b * dt) - 2.0f * ec + 2.0f * (1.0f + ec)
-            / exp(b * dt), 4.0f));
+      (p * (cos(two_cf_pi_t) - sqrt(3.0f - two_pow) * sin(two_cf_pi_t)))
+      * (p * (cos(two_cf_pi_t) + sqrt(3.0f - two_pow) * sin(two_cf_pi_t)))
+      * (p * (cos(two_cf_pi_t) - sqrt(3.0f + two_pow) * sin(two_cf_pi_t)))
+      * (p * (cos(two_cf_pi_t) + sqrt(3.0f + two_pow) * sin(two_cf_pi_t)))
+      / pow(-2.0f / exp(2.0f * b_dt) - 2.0f * ec + 2.0f * (1.0f + ec)
+            / exp(b_dt), 4.0f));
 
     // The filter coefficients themselves:
-    forward[ch].resize(5, 0.0f);
-    feedback[ch].resize(9, 0.0f);
+    const int coeff_count = 9;
+    forward_[ch].resize(coeff_count, 0.0f);
+    back_[ch].resize(coeff_count, 0.0f);
+    state_[ch].resize(coeff_count, 0.0f);
 
-    forward[ch][0] = pow(T, 4.0f) / gain;
-    forward[ch][1] = (-4.0f * pow(T, 4.0f) * cos(2.0f * cf * M_PI * dt)
+    forward_[ch][0] = pow(dt, 4.0f) / gain;
+    forward_[ch][1] = (-4.0f * pow(dt, 4.0f) * cos(2.0f * cpt)
                       / exp(b * dt) / gain);
-    forward[ch][2] = (6.0f * pow(T, 4.0f) * cos(4.0f * cf * M_PI * dt)
+    forward_[ch][2] = (6.0f * pow(dt, 4.0f) * cos(4.0f * cpt)
                       / exp(2.0f * b * dt) / gain);
-    forward[ch][3] = (-4.0f * pow(T, 4.0f) * cos(6.0f * cf * M_PI * dt)
+    forward_[ch][3] = (-4.0f * pow(dt, 4.0f) * cos(6.0f * cpt)
                       / exp(3.0f * b * dt) / gain);
-    forward[ch][4] = (pow(T,4.0f) * cos(8.0f * cf * M_PI * dt)
+    forward_[ch][4] = (pow(dt, 4.0f) * cos(8.0f * cpt)
                       / exp(4.0f * b * dt) / gain);
+    // Note: the remainder of the forward vector is zero-padded
 
-    feedback[ch][0] = 1.0f;
-    feedback[ch][1] = -8.0f * cos(2.0f * cf * M_PI * T) / exp(b * dt);
-    feedback[ch][2] = (4.0f * (4.0f + 3.0f * cos(4.0f * cf * M_PI * dt))
-                       / exp(2.0f * b * dt));
-    feedback[ch][3] = (-8.0f * (6.0f * cos(2.0f * cf * M_PI * dt)
-                                + cos(6.0f * cf * M_PI * dt))
-                       / exp(3.0f * b * dt));
-    feedback[ch][4] = (2.0f * (18.0f + 16.0f * cos(4.0f * cf * M_PI * dt)
-                               + cos(8.0f * cf * M_PI * dt))
-                       / exp(4.0f * b * dt));
-    feedback[ch][5] = (-8.0f * (6.0f * cos(2.0f * cf * M_PI * T)
-                                + cos(6.0f * cf * M_PI * T))
-                       / exp(5.0f * B * T));
-    feedback[ch][6] = (4.0f * (4.0f + 3.0f * cos(4.0f * cf * M_PI * T))
-                       / exp(6.0f * B * T));
-    feedback[ch][7] = -8.0f * cos(2.0f * cf * M_PI * dt) / exp(7.0f * b * dt);
-    feedback[ch][8] = exp(-8.0f * b * dt);
+    back_[ch][0] = 1.0f;
+    back_[ch][1] = -8.0f * cos(2.0f * cpt) / exp(b * dt);
+    back_[ch][2] = (4.0f * (4.0f + 3.0f * cos(4.0f * cpt))
+                    / exp(2.0f * b * dt));
+    back_[ch][3] = (-8.0f * (6.0f * cos(2.0f * cpt) + cos(6.0f * cpt))
+                    / exp(3.0f * b * dt));
+    back_[ch][4] = (2.0f * (18.0f + 16.0f * cos(4.0f * cpt) + cos(8.0f * cpt))
+                    / exp(4.0f * b * dt));
+    back_[ch][5] = (-8.0f * (6.0f * cos(2.0f * cpt) + cos(6.0f * cpt))
+                    / exp(5.0f * b * dt));
+    back_[ch][6] = (4.0f * (4.0f + 3.0f * cos(4.0f * cpt))
+                    / exp(6.0f * b * dt));
+    back_[ch][7] = -8.0f * cos(2.0f * cpt) / exp(7.0f * b * dt);
+    back_[ch][8] = exp(-8.0f * b * dt);
   }
+  output_.Initialize(num_channels_,
+                     input.buffer_length(),
+                     input.sample_rate());
+  return true;
 }
+
+void ModuleGammatone::Process(const SignalBank &input) {
+  output_.set_start_time(input.start_time());
+  int audio_channel = 0;
+
+  vector<vector<float> >::iterator b = forward_.begin();
+  vector<vector<float> >::iterator a = back_.begin();
+  vector<vector<float> >::iterator s = state_.begin();
+
+  for (int ch = 0; ch < num_channels_; ++ch, ++a, ++b, ++s) {
+    for (int i = 0; i < input.buffer_length(); ++i) {
+      // Direct-form-II IIR filter  
+      float in = input.sample(audio_channel, i);
+      float out = (*b)[0] * in + (*s)[0];
+      for (unsigned int stage = 1; stage < s->size(); ++stage)
+        (*s)[stage - 1] = (*b)[stage] * in - (*a)[stage] * out + (*s)[stage];
+      output_.set_sample(ch, i, out);
+    }
+  }
+  PushOutput();
+}
+
 } // namespace aimc
\ No newline at end of file

--- a/src/Modules/BMM/ModuleGammatone.h	Tue Feb 16 18:40:55 2010 +0000
+++ b/src/Modules/BMM/ModuleGammatone.h	Thu Feb 18 16:55:40 2010 +0000
@@ -26,6 +26,10 @@
 
 #include <vector>
 
+#include "Support/Module.h"
+#include "Support/Parameters.h"
+#include "Support/SignalBank.h"
+
 namespace aimc {
 using std::vector;
 class ModuleGammatone : public Module {
@@ -39,7 +43,8 @@
   virtual bool InitializeInternal(const SignalBank& input);
   virtual void ResetInternal();
   vector<vector<float> > forward_;
-  vector<vector<float> > feedback_;
+  vector<vector<float> > back_;
+  vector<vector<float> > state_;
   vector<float> centre_frequencies_;
   int num_channels_;
   float max_frequency_;

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Modules/BMM/ModuleGammatone_test.py	Thu Feb 18 16:55:40 2010 +0000
@@ -0,0 +1,83 @@
+#!/usr/bin/env python
+# encoding: utf-8
+#
+# AIM-C: A C++ implementation of the Auditory Image Model
+# http://www.acousticscale.org/AIMC
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+ModuleGammatone_test.py
+
+Created by Thomas Walters on 2010-02-15.
+Copyright 2010 Thomas Walters <[e-mail]>
+Test for the Slaney IIR gammatone.
+"""
+
+import aimc
+from scipy import io
+
+def main():
+  data_file = "src/Modules/BMM/testdata/gammatone.mat"
+  data = io.loadmat(data_file)
+  
+  # The margin of error allowed between the returned values from AIM-C and
+  # the stored MATLAB values.
+  epsilon = 0.000001;
+  
+  input_wave = data["input_wave"]
+  sample_rate = data["sample_rate"]
+  centre_frequencies = data["centre_frequencies"]
+  expected_output = data["expected_output"]
+  
+  (channel_count, buffer_length, frame_count) = expected_output.shape
+  
+  input_sig = aimc.SignalBank()
+  input_sig.Initialize(1, buffer_length, 44100)
+  parameters = aimc.Parameters()
+  mod_gt = aimc.ModuleGammatone(parameters)
+  mod_gt.Initialize(input_sig)
+  
+  correct_count = 0;
+  incorrect_count  = 0;
+  for p in range(0, profile_count):
+    profile = given_profiles[p]
+    features = matlab_features[p]
+    for i in range(0, channel_count):
+      profile_sig.set_sample(i, 0, profile[i])
+    mod_gauss.Process(profile_sig)
+    out_sig = mod_gauss.GetOutputBank()
+    error = False;
+    for j in range(0, out_sig.channel_count()):
+      if (abs(out_sig.sample(j, 0) - features[j]) > epsilon):
+        error = True;
+        incorrect_count += 1;
+      else:
+        correct_count += 1;
+    if error:
+      print("Mismatch at profile %d" % (p))
+      print("AIM-C values: %f %f %f %f" % (out_sig.sample(0, 0), out_sig.sample(1, 0), out_sig.sample(2, 0), out_sig.sample(3, 0)))
+      print("MATLAB values: %f %f %f %f" % (features[0], features[1], features[2], features[3]))
+      print("")
+  percent_correct = 100 * correct_count / (correct_count + incorrect_count)
+  print("Total correct: %f percent" % (percent_correct))
+  if percent_correct == 100:
+    print("=== TEST PASSED ===")
+  else:
+    print("=== TEST FAILED! ===")
+
+  pass
+
+
+if __name__ == '__main__':
+  main()

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Modules/Output/FileOutputHTK.cc	Thu Feb 18 16:55:40 2010 +0000
@@ -0,0 +1,194 @@
+// Copyright 2006-2010, Thomas Walters
+//
+// AIM-C: A C++ implementation of the Auditory Image Model
+// http://www.acousticscale.org/AIMC
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+/*!
+ * \file
+ * \brief File output in the HTK format. 
+ *
+ * \author Tom Walters <[e-mail]> 
+ * \author Willem van Engen <[e-mail]>
+ * \date created 2006/10/30
+ * \version \$Id$
+ */
+
+#ifdef _WINDOWS
+#	include <direct.h> // for _mkdir&_rmdir
+#else
+#	include <sys/types.h>
+#	include <dirent.h> // for opendir&friends
+#endif
+#include <stdio.h>
+#include <string.h>
+#include <cmath>
+
+#include "Modules/Output/FileOutputHTK.h"
+
+namespace aimc {
+FileOutputHTK::FileOutputHTK(Parameters *params) : Module(params) {
+	file_handle_ = NULL;
+	header_written_ = false;
+  filename_[0] = '\0';
+  frame_period_ms_ = 0.0f;
+}
+
+FileOutputHTK::~FileOutputHTK() {
+	if (file_handle_ != NULL)
+		CloseFile();
+}
+
+bool FileOutputHTK::OpenFile(const char* filename, float frame_period_ms) {
+	if (file_handle_ != NULL) {
+	  LOG_ERROR(_T("Couldn't open output file. A file is already open."));
+  	return false;
+	}
+
+	// Check that the output file exists and is writeable
+	if ((file_handle_ = fopen(filename, "wb"))==NULL ) {
+		LOG_ERROR(_T("Couldn't open output file '%s' for writing."), filename);
+		return false;
+	}
+	strcpy(filename_, filename);
+	sample_count_ = 0;
+  frame_period_ms_ = frame_period_ms;
+  header_written_ = false;
+	return true;
+}
+
+bool FileOutputHTK::InitializeInternal(const SignalBank &input) {
+  if (file_handle_ == NULL) {
+    LOG_ERROR(_T("Couldn't initialize file output. "
+                 "Please call FileOutputHTK::OpenFile first"));
+		return false;
+  }
+  if (header_written_) {
+    LOG_ERROR(_T("A header has already been written on the output file."
+                 "Please call FileOutputHTK::CloseFile to close that file, "
+                 "and FileOutputHTK::OpenFile to open an new one before "
+                 "calling FileOutputHTK::Initialize again."));
+		return false;
+  }
+  channel_count_ = input.channel_count();
+  buffer_length_ = input.buffer_length();
+  WriteHeader(channel_count_ * buffer_length_, frame_period_ms_);
+  return true;
+}
+
+void FileOutputHTK::ResetInternal() {
+  if (file_handle_ != NULL && !header_written_) {
+    WriteHeader(channel_count_ * buffer_length_, frame_period_ms_);
+  }
+}
+
+void FileOutputHTK::WriteHeader(int num_elements, float period_ms) {
+	if (header_written_)
+		return;
+
+	/* HTK format file: (taken from the HTK book - section 5.10.1)
+	 * Header: 12 bytes in total, contains:
+	 * sample_count - number of samples in file (4-byte integer)(long)
+	 * sample_period - sample period in 100ns units (4-byte integer)(long)
+	 * sample_size - number of bytes per sample (2-byte integer)(short)
+	 * parameter_kind - a code indicating the sample kind (2-byte integer)(short)
+	 */
+
+	 // To be filled in when the file is done
+	int32_t sample_count = 0;
+
+  int32_t sample_period = floor(1e4 * period_ms);
+  int16_t sample_size = num_elements * sizeof(float);
+
+  // User-defined coefficients with energy term
+	int16_t parameter_kind = H_USER + H_E;
+
+	// Fix endianness
+	sample_count = ByteSwap32(sample_count);
+	sample_period = ByteSwap32(sample_period);
+	sample_size = ByteSwap16(sample_size);
+	parameter_kind = ByteSwap16(parameter_kind);
+
+  // Enter header values. sample_count is a dummy value which is filled in on
+  // file close
+	fwrite(&sample_count, sizeof(sample_count), 1, file_handle_);
+	fwrite(&sample_period, sizeof(sample_period), 1, file_handle_);
+	fwrite(&sample_size, sizeof(sample_size), 1, file_handle_);
+	fwrite(&parameter_kind, sizeof(parameter_kind), 1, file_handle_);
+	fflush(file_handle_);
+
+	header_written_ = true;
+}
+
+
+void FileOutputHTK::Process(const SignalBank &input) {
+  if (file_handle_ == NULL) {
+    LOG_ERROR(_T("Couldn't process file output. No file is open."
+                 "Please call FileOutputHTK::OpenFile first"));
+		return;
+  }
+
+  if (!header_written_) {
+    LOG_ERROR(_T("No header has been written on the output file yet. Please"
+                 "call FileOutputHTK::Initialize() before calling "
+                 "FileOutputHTK::Process()"));
+		return;
+  }
+	float s;
+
+	for (int ch = 0; ch < input.channel_count(); ch++) {
+    for (int i = 0; i < input.buffer_length(); i++) {
+      s = input.sample(ch, i);
+      s = ByteSwapFloat(s);
+      fwrite(&s, sizeof(float), 1, file_handle_);
+    }
+  }
+	sample_count_++;
+}
+
+bool FileOutputHTK::CloseFile() {
+	if (file_handle_ == NULL)
+    return false;
+
+	// Write the first 4 bytes of the file
+	// with how many samples there are in the file
+	fflush(file_handle_);
+	rewind(file_handle_);
+	fflush(file_handle_);
+  int32_t samples = sample_count_;
+	samples = ByteSwap32(samples);
+	fwrite(&samples, sizeof(samples), 1, file_handle_);
+
+	// And close the file
+	fclose(file_handle_);
+	file_handle_ = NULL;
+	return true;
+}
+
+float FileOutputHTK::ByteSwapFloat(float d) {
+  // Endianness fix
+  float a;
+  unsigned char *dst = (unsigned char *)&a;
+  unsigned char *src = (unsigned char *)&d;
+
+  dst[0] = src[3];
+  dst[1] = src[2];
+  dst[2] = src[1];
+  dst[3] = src[0];
+
+  return a;
+}
+}  //namespace aimc
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Modules/Output/FileOutputHTK.h	Thu Feb 18 16:55:40 2010 +0000
@@ -0,0 +1,106 @@
+/*!
+ * \file
+ * \brief File output to HTK format class definition
+ *
+ * \author Tom Walters <[e-mail]> and Willem van Engen <[e-mail]>
+ * \date created 2006/10/30
+ * \version \$Header$
+ */
+/* (c) 2006, University of Cambridge, Medical Research Council
+ * http://www.pdn.cam.ac.uk/groups/cnbh/aimmanual
+ */
+#ifndef _AIMC_MODULE_OUTPUT_HTK_H_
+#define _AIMC_MODULE_OUTPUT_HTK_H_
+
+#include "Support/Module.h"
+#include "Support/SignalBank.h"
+
+
+// Defines taken from HTKwrite.c and The HTK Book
+#define H_WAVEFORM		0	//sampled waveform
+#define H_LPC			1	//linear prediction filter coefficients
+#define H_LPREFC		2	//linear prediction reflection coefficients
+#define H_LPCEPSTRA		3	//LPC cepstral coefficients
+#define H_LPDELCEP		4	//LPC cepstra plus delta coefficients
+#define H_IREFC			5	//LPC reflection coef in 16 bit integer format
+#define H_MFCC			6	//mel-frequency cepstral coefficients
+#define H_FBANK			7	//log mel-filter bank channel outputs
+#define H_MELSPEC		8	//linear mel-filter bank channel outputs
+#define H_USER			9	//user defined sample kind
+#define H_DISCRETE		10	//vector quantised data
+#define H_PLP           11  // Perceptual Linear Prediction
+#define H_ANON          12
+
+#define H_E 64 //has energy
+#define H_N 128 //absolute energy suppressed
+#define H_D 256 //has delta coefficients
+#define H_A 512 //has acceleration coefficients
+#define H_C 1024 //is compressed
+#define H_Z 2048 //has zero mean static coef.
+#define H_K 4096 //has CRC checksum
+#define H_O 8192 //has 0th cepstral coef.
+#define H_V 16384 // Attach vq index
+#define H_T 32768 // Attach delta-delta-delta index
+
+// HTK fomat is big-endian...
+#define ByteSwap16(n) \
+( ((((uint16_t) n) << 8) & 0xFF00) | \
+((((uint16_t) n) >> 8) & 0x00FF) )
+
+#define ByteSwap32(n) \
+( ((((uint32_t) n) << 24) & 0xFF000000) | \
+((((uint32_t) n) << 8) & 0x00FF0000) | \
+((((uint32_t) n) >> 8) & 0x0000FF00) | \
+((((uint32_t) n) >> 24) & 0x000000FF) )
+
+/*!
+ * \class FileOutputHTK "Output/FileOutputHTK.h"
+ * \brief File output to HTK class
+ *
+ * This class gives a method for saving either a signal or a profile to HTK format.
+ * \sa Signal, SignalBank
+ */
+namespace aimc {
+class FileOutputHTK : public Module {
+ public:
+	/*! \brief Create a new file output for an HTK format file. Use of this 
+	 *  class only really makes sense for the output of 1-D frames.
+	 */
+	FileOutputHTK(Parameters *pParam);
+	~FileOutputHTK();
+
+	/*! \brief Initialize the output to HTK.
+	 *  \param *filename Filename of the ouptut file to be created.
+	 *  If the file exists it will be overwritten
+	 *  \return Returns true on success of initialization.
+	 */
+	bool OpenFile(const char *filename, float frame_period_ms);
+  bool CloseFile();
+  virtual void Process(const SignalBank &input);
+protected:
+  virtual bool InitializeInternal(const SignalBank &input);
+  virtual void ResetInternal();
+
+  float ByteSwapFloat(float d);
+
+	void WriteHeader(int nelements, float sampPeriod);
+
+	//! \brief Whether initialization is done or not
+	bool header_written_;
+
+	//! \brief Filename
+	char filename_[PATH_MAX];
+	//! \brief Internal pointer to the output file
+	FILE *file_handle_;
+
+	//! \brief Count of the number of samples in the file, written on close
+	int sample_count_;
+
+  int channel_count_;
+  int buffer_length_;
+  float frame_period_ms_;
+};
+}  // namespace aimc
+
+#endif  // _AIMC_MODULE_OUTPUT_HTK_H_
+

--- a/src/Modules/SAI/ModuleSAI.cc	Tue Feb 16 18:40:55 2010 +0000
+++ b/src/Modules/SAI/ModuleSAI.cc	Thu Feb 18 16:55:40 2010 +0000
@@ -29,11 +29,12 @@
 
 #include "Modules/SAI/ModuleSAI.h"
 
+namespace aimc {
 ModuleSAI::ModuleSAI(Parameters *parameters) : Module(parameters) {
+  module_identifier_ = "weighted_sai";
+  module_type_ = "sai";
   module_description_ = "Stabilised auditory image";
-  module_name_ = "sai2007";
-  module_type_ = "sai";
-  module_id_ = "$Id: ModuleSAI.cc 4 2010-02-03 18:44:58Z tcw $";
+  module_version_ = "$Id: ModuleSAI.cc 4 2010-02-03 18:44:58Z tcw $";
 
   min_delay_ms_ = parameters_->DefaultFloat("sai.min_delay_ms", 0.0f);
   max_delay_ms_ = parameters_->DefaultFloat("sai.max_delay_ms", 35.0f);
@@ -43,30 +44,32 @@
                                                    0.03f);
   frame_period_ms_ = parameters_->DefaultFloat("sai.frame_period_ms", 20.0f);
 
-  min_strobe_delay_index_ = 0;
-  max_strobe_delay_index_ = 0;
+  max_concurrent_strobes_
+    = parameters_->DefaultInt("sai.max_concurrent_strobes", 50);
+
+  min_strobe_delay_idx_ = 0;
+  max_strobe_delay_idx_ = 0;
   sai_decay_factor_ = 0.0f;
-  max_concurrent_strobes_ = 0;
-  output_frame_period_ms_ = 0.0f;
-  last_output_frame_time_ms_ = 0.0f;
+  fire_counter_ = 0;
 }
 
 bool ModuleSAI::InitializeInternal(const SignalBank &input) {
   // The SAI output bank must be as long as the SAI's Maximum delay.
   // One sample is added to the SAI buffer length to account for the
   // zero-lag point
-  int sai_buffer_length = 1 + Round(input.sample_rate() * max_delay_ms_);
+  int sai_buffer_length = 1 + floor(input.sample_rate() * max_delay_ms_);
+  channel_count_ = input.channel_count();
 
   // Make an output SignalBank with the same number of channels and centre
   // frequencies as the input, but with a different buffer length
   if (!output_.Initialize(input.channel_count(),
                           sai_buffer_length,
-                          input.sample_rate());) {
+                          input.sample_rate())) {
     LOG_ERROR("Failed to create output buffer in SAI module");
     return false;
   }
   for (int i = 0; i < input.channel_count(); ++i ) {
-    output_.set_centre_frequency(i, input.get_centre_frequency(i));
+    output_.set_centre_frequency(i, input.centre_frequency(i));
   }
 
   // sai_temp_ will be initialized to zero
@@ -75,169 +78,154 @@
     return false;
   }
 
-  last_output_time_ms_ = 0.0f;
-
-  frame_period_samples_ = Round(input.sample_rate()
-                                * frame_period_ms_ / 1000.0f);
-
-  min_strobe_delay_idx_ = Round(input.sample_rate() * min_delay_ms_ / 1000.0f);
-  max_strobe_delay_idx_ = Round(input.sample_rate() * max_delay_ms_ / 1000.0f);
+  frame_period_samples_ = floor(input.sample_rate() * frame_period_ms_
+                                / 1000.0f);
+  min_strobe_delay_idx_ = floor(input.sample_rate() * min_delay_ms_
+                                / 1000.0f);
+  max_strobe_delay_idx_ = floor(input.sample_rate() * max_delay_ms_
+                                / 1000.0f);
 
   // Make sure we don't go past the output buffer's upper bound
-  if (max_strobe_delay_idx_ > output_.buffer_length()))
+  if (max_strobe_delay_idx_ > output_.buffer_length()) {
     max_strobe_delay_idx_ = output_.buffer_length();
+  }
 
   // Define decay factor from time since last sample (see ti2003)
   sai_decay_factor_ = pow(0.5f, 1.0f / (buffer_memory_decay_
                                         * input.sample_rate()));
 
-  // Maximum strobes that can be active at the same time within maxdelay.
-  //! \todo Choose this value in a more principled way
-  max_concurrent_strobes_ = Round(1000.0f * max_delay_ * 5);
-
   // Precompute strobe weights
   strobe_weights_.resize(max_concurrent_strobes_);
   for (int n = 0; n < max_concurrent_strobes_; ++n) {
-    strobe_weights_[n] = pow(1.0f / (n + 1)), strobe_weight_alpha_);
+    strobe_weights_[n] = pow(1.0f / (n + 1), strobe_weight_alpha_);
   }
 
-  // Active Strobes
-  active_strobes_.Resize(input.channel_count());
-  for (int i = 0; i < input.channel_count(); ++i) {
-    active_strobes_[i].Create(max_concurrent_strobes_);
-  }
-  next_strobes_.resize(input.channel_count(), 0);
+  ResetInternal();
 
   return true;
 }
 
 void ModuleSAI::ResetInternal() {
+  // Active Strobes
+  active_strobes_.clear();
+  active_strobes_.resize(channel_count_);
+  fire_counter_ = frame_period_samples_ - 1;
 }
 
 void ModuleSAI::Process(const SignalBank &input) {
-  int s;
-  int c;
-  int output_buffer_length = output_.buffer_length();
 
   // Reset the next strobe times
   next_strobes_.clear();
   next_strobes_.resize(output_.channel_count(), 0);
 
   // Offset the times on the strobes from the previous buffer
-  for (c = 0; c < input.channel_count(), ++c) {
-    active_strobes_[c].shiftStrobes(input.buffer_length());
+  for (int ch = 0; ch < input.channel_count(); ++ch) {
+    active_strobes_[ch].ShiftStrobes(input.buffer_length());
   }
 
-  // Make sure only start time is transferred to the output
-  output_.set_start_time(input.start_time());
-
   // Loop over samples to make the SAI
-  for (s = 0; s < input_buffer_length; ++s) {
+  for (int i = 0; i < input.buffer_length(); ++i) {
     float decay_factor = pow(sai_decay_factor_, fire_counter_);
     // Loop over channels
-    for (c = 0; c < input.channel_count(); ++c) {
+    for (int ch = 0; ch < input.channel_count(); ++ch) {
       // Local convenience variables
-      StrobeList &active_strobes = active_strobes_[c];
-      float centre_frequency = input.get_centre_frequency(c);
-      int next_strobe = next_strobes_[c];
+      StrobeList &active_strobes = active_strobes_[ch];
+      int next_strobe_index = next_strobes_[ch];
 
-      // 1. Update strobes
+      // Update strobes
       // If we are up to or beyond the next strobe...
-      if (next_strobe < input.strobe_count(c)) {
-        if (s == pSigIn->getStrobe(iNextStrobe)) {
-          //A new strobe has arrived
-          // if there aren't already too many strobes active...
-          if ((active_strobes.getStrobeCount() + 1) < max_concurrent_strobes_) {
-            // ...add the active strobe to the list of current strobes
-            // calculate the strobe weight
-            float weight = 1.0f;
-            if (active_strobes.getStrobeCount() > 0) {
-              int last_strobe = active_strobes.getTime(
-                active_strobes.getStrobeCount());
-
-              // If the strobe occured within 10 impulse-response
-              // cycles of the previous strobe, then lower its weight
-              weight = (s - iLastStrobe) / input.sample_rate()
-                       * centre_frequency / 10.0f;
-              if (weight > 1.0f)
-                weight = 1.0f;
-            }
-            pActiveStrobes->addStrobe(iCurrentSample, weight);
-            iNextStrobe++;
-          } else {
-            // We have a problem
-            aimASSERT(0);
+      if (next_strobe_index < input.strobe_count(ch)) {
+        if (i == input.strobe(ch, next_strobe_index)) {
+          // A new strobe has arrived.
+          // If there are too many strobes active, then get rid of the
+          // earliest one
+          if (active_strobes.strobe_count() >= max_concurrent_strobes_) {
+            active_strobes.DeleteFirstStrobe();
           }
 
-          // 2. Having updated the strobes, we now need to update the
-          // strobe weights
+          // Add the active strobe to the list of current strobes and
+          // calculate the strobe weight
+          float weight = 1.0f;
+          if (active_strobes.strobe_count() > 0) {
+            int last_strobe_time = active_strobes.Strobe(
+              active_strobes.strobe_count() - 1).time;
+
+            // If the strobe occured within 10 impulse-response
+            // cycles of the previous strobe, then lower its weight
+            weight = (i - last_strobe_time) / input.sample_rate()
+                     * input.centre_frequency(ch) / 10.0f;
+            if (weight > 1.0f)
+              weight = 1.0f;
+          }
+          active_strobes.AddStrobe(i, weight);
+          next_strobe_index++;
+
+
+          // Update the strobe weights
           float total_strobe_weight = 0.0f;
-          for (int si = 1; si <= pActiveStrobes->getStrobeCount(); ++si) {
-            total_strobe_weight += (pActiveStrobes->getWeight(si) 
-              * m_pStrobeWeights[pActiveStrobes->getStrobeCount() - si]);
+          for (int si = 0; si < active_strobes.strobe_count(); ++si) {
+            total_strobe_weight += (active_strobes.Strobe(si).weight
+              * strobe_weights_[active_strobes.strobe_count() - si - 1]);
           }
-          for (int si = 1; si <= pActiveStrobes->getStrobeCount(); ++si) {
-            active_strobes.setWorkingWeight(si,(active_strobes.getWeight(si)
-              * strobe_weights_[active_strobes.getStrobeCount() - si])
-                / total_strobe_weight);
+          for (int si = 0; si < active_strobes.strobe_count(); ++si) {
+            active_strobes.SetWorkingWeight(si,
+              (active_strobes.Strobe(si).weight
+               * strobe_weights_[active_strobes.strobe_count() - si - 1])
+              / total_strobe_weight);
           }
         }
       }
 
-      // remove inactive strobes...
-      while (pActiveStrobes->getStrobeCount() > 0) {
-        // Get the time of the first strobe (ordering of strobes is
-        // from one, not zero)
-        int iStrobeTime = pActiveStrobes->getTime(1);
-        int iDelay = iCurrentSample - iStrobeTime;
-        // ... do we now need to remove this strobe?
-        if (iDelay > m_maxStrobeDelayIdx)
-          pActiveStrobes->deleteFirstStrobe();
+      // Remove inactive strobes
+      while (active_strobes.strobe_count() > 0) {
+        // Get the relative time of the first strobe, and see if it exceeds
+        // the maximum allowed time.
+        if ((i - active_strobes.Strobe(0).time) > max_strobe_delay_idx_)
+          active_strobes.DeleteFirstStrobe();
         else
           break;
-          // Since the strobes are ordered, we don't need to go
-          // beyond the first still-active strobe
       }
 
-      // 3. Loop over active strobes
-      for (int si = 1; si <= pActiveStrobes->getStrobeCount(); si++) {
-        // 3.1 Add effect of active strobe at correct place in the SAI buffer
-        // Calculate the time from the strobe event to 'now': iDelay
-        int iStrobeTime = pActiveStrobes->getTime(si);
-        int iDelay = iCurrentSample - iStrobeTime;
+      // Update the SAI buffer with the weighted effect of all the active
+      // strobes at the current sample
+      for (int si = 0; si < active_strobes.strobe_count(); ++si) {
+        // Add the effect of active strobe at correct place in the SAI buffer
+        // Calculate 'delay', the time from the strobe event to now
+        int delay = i - active_strobes.Strobe(si).time;
 
         // If the delay is greater than the (user-set)
         // minimum strobe delay, the strobe can be used
-        if (iDelay >= m_minStrobeDelayIdx && iDelay < m_maxStrobeDelayIdx) {
+        if (delay >= min_strobe_delay_idx_ && delay < max_strobe_delay_idx_) {
           // The value at be added to the SAI
-          float sig = pSigIn->getSample(iCurrentSample, audCh);
+          float sig = input.sample(ch, i);
 
           // Weight the sample correctly
-          sig *= pActiveStrobes->getWorkingWeight(si);
+          sig *= active_strobes.Strobe(si).working_weight;
 
           // Adjust the weight acording to the number of samples until the
           // next output frame
-          sig *= fDecayFactor;
+          sig *= decay_factor;
 
           // Update the temporary SAI buffer
-          pSigOut->setSample(iDelay, audCh,
-                             pSigOut->getSample(iDelay, audCh)+sig);
+          output_.set_sample(ch, delay, output_.sample(ch, delay) + sig);
         }
       }
 
-      m_pNextStrobes[bankCh]=iNextStrobe;
+      next_strobes_[ch] = next_strobe_index;
 
     } // End loop over channels
 
+    fire_counter_--;
 
-    //Check to see if we need to output an SAI frame this sample
-    if (m_iFireCounter-- == 0) {
+    // Check to see if we need to output an SAI frame on this sample
+    if (fire_counter_ <= 0) {
       // Decay the SAI by the correct amount and add the current output frame
-      float decay = pow(sai_decay_factor_, fire_period_samples_);
+      float decay = pow(sai_decay_factor_, frame_period_samples_);
 
-      for (c = 0; c < input.channel_count(); ++c) {
+      for (int ch = 0; ch < input.channel_count(); ++ch) {
         for (int i = 0; i < output_.buffer_length(); ++i) {
-          output_.set_sample(c, i, sai_temp_[c][i] + output_[c][i] * decay);
+          output_.set_sample(ch, i,
+                             sai_temp_[ch][i] + output_[ch][i] * decay);
         }
       }
 
@@ -248,10 +236,10 @@
         }
       }
 
-      m_iFireCounter=m_iFirePeriodSamples-1;
+      fire_counter_ = frame_period_samples_ - 1;
 
-      // Make sure the start time is transferred to the output
-      m_pOutputData->setStartTime(m_pInputData->getSignal(0)->getStartTime()+(SIGNAL_SAMPLE)((float)iCurrentSample*1000.0f/(float)m_pInputData->getSamplerate()));
+      // Transfer the current time to the output buffer
+      output_.set_start_time(input.start_time() + i);
       PushOutput();
     }
   } // End loop over samples
@@ -259,3 +247,4 @@
 
 ModuleSAI::~ModuleSAI() {
 }
+}  // namespace aimc

--- a/src/Modules/SAI/ModuleSAI.h	Tue Feb 16 18:40:55 2010 +0000
+++ b/src/Modules/SAI/ModuleSAI.h	Thu Feb 18 16:55:40 2010 +0000
@@ -78,21 +78,18 @@
 
   /*! \brief The maximum number strobes that can be active at the same time.
    *
-   *  A strobe lasts for strobe.maxdelay, there can only be a certain number
-   *  of strobes active at the same time, that's this value. It's used for
-   *  allocating memory buffers, like m_pUnfinishedStrobeCount and
-   *  m_pStrobeWeights.
    */
   int max_concurrent_strobes_;
 
-  int fire_period_samples_;
   int fire_counter_;
 
   //! \brief Period in milliseconds between output frames
-  float output_frame_period_ms_;
+  float frame_period_ms_;
+  int frame_period_samples_;
 
-  //! \brief Time of the last frame output
-  float last_output_frame_time_ms_;
+  float min_delay_ms_;
+  float max_delay_ms_;
+  int channel_count_;
 };
 }  // namespace aimc
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Modules/Strobes/ModuleParabola.cc	Thu Feb 18 16:55:40 2010 +0000
@@ -0,0 +1,151 @@
+// Copyright 2007-2010, Thomas Walters
+//
+// AIM-C: A C++ implementation of the Auditory Image Model
+// http://www.acousticscale.org/AIMC
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+/*!
+ * \file
+ * \brief Parabola strobe detection module - Using the 'parabloa' strobe
+ * criterion from the AIM-MAT sf2003 module
+ *
+ * \author Tom Walters <[e-mail]>
+ * \date created 2007/08/01
+ * \version \$Id:$
+ */
+
+#include <math.h>
+
+#include "Modules/Strobes/ModuleParabola.h"
+
+namespace aimc {
+ModuleParabola::ModuleParabola(Parameters *params) : Module(params) {
+  module_description_ = "sf2003 parabola algorithm";
+  module_identifier_ = "parabola";
+  module_type_ = "strobes";
+  module_version_ = "$Id$";
+
+  // Get data from parameters
+  height_ = parameters_->DefaultFloat("parabola.height", 1.2f);
+  parabw_ = parameters_->DefaultFloat("parabola.width_cycles", 1.5f);
+  strobe_decay_time_ = parameters_->DefaultFloat("parabla.strobe_decay_time",
+                                                 0.02f);
+  max_strobes_ = parameters_->DefaultInt("parabola.max_strobes", 50);
+  channel_count_ = 0;
+}
+
+bool ModuleParabola::InitializeInternal(const SignalBank &input) {
+  output_.Initialize(input);
+  channel_count_ = input.channel_count();
+  sample_rate_ = input.sample_rate();
+
+  // Parameters for the parabola
+  parab_a_.resize(channel_count_);
+  parab_b_.resize(channel_count_);
+  parab_wnull_.resize(channel_count_);
+  parab_var_samples_.resize(channel_count_);
+
+  for (int ch = 0; ch < channel_count_; ++ch) {
+    parab_wnull_[ch] = parabw_ / input.centre_frequency(ch);
+    parab_var_samples_[ch] = floor(parab_wnull_[ch] * sample_rate_);
+    parab_a_[ch] = 4.0f * (1.0f - height_)
+                   / (parab_wnull_[ch] * parab_wnull_[ch]);
+    parab_b_[ch] = -parab_wnull_[ch] / 2.0f;
+  }
+
+  // Number of samples over which the threshold should decay
+  strobe_decay_samples_ = floor(sample_rate_ * strobe_decay_time_);
+
+  // Prepare internal buffers
+  ResetInternal();
+
+  return true;
+}
+
+void ModuleParabola::ResetInternal() {
+  threshold_.resize(channel_count_, 0.0f);
+  last_threshold_.resize(channel_count_, 0.0f);
+  samples_since_last_strobe_.resize(channel_count_, 0);
+
+  prev_sample_.resize(channel_count_, 10000.0f);
+  curr_sample_.resize(channel_count_, 5000.0f);
+  next_sample_.resize(channel_count_, 0.0f);
+}
+
+void ModuleParabola::Process(const SignalBank &input) {
+  float decay_constant;
+
+  for (int ch = 0; ch < output_.channel_count(); ch++) {
+    output_.ResetStrobes(ch);
+  }
+  output_.set_start_time(input.start_time());
+
+  // Loop across samples first, then channels
+  for (int i = 0; i < input.buffer_length(); i++) {
+    // Find strobes in each channel first
+    for (int ch = 0; ch < input.channel_count(); ++ch) {
+      // Shift all the samples by one
+      // curr_sample is the sample at time (i - 1)
+      prev_sample_[ch] = curr_sample_[ch];
+      curr_sample_[ch] = next_sample_[ch];
+      next_sample_[ch] = input.sample(ch, i);
+
+      // Copy input signal to output signal
+      output_.set_sample(ch, i, input.sample(ch, i));
+
+	    if (curr_sample_[ch] >= threshold_[ch]) {
+	      threshold_[ch] = curr_sample_[ch];
+		    if (prev_sample_[ch] < curr_sample_[ch]
+		        && next_sample_[ch] < curr_sample_[ch]) {
+		      // We have a strobe: set threshold and add strobe to the list
+		      output_.AddStrobe(ch, i - 1);
+		      last_threshold_[ch] = threshold_[ch];
+		      parab_var_samples_[ch] =
+		        floor(input.sample_rate()
+		              * (parab_wnull_[ch] - (threshold_[ch]
+		                                     - 2.0f * parab_a_[ch]  *parab_b_[ch])
+		                                    / (2.0f * parab_a_[ch])));
+		    }
+      }
+	    if (output_.strobe_count(ch) > 0) {
+        samples_since_last_strobe_[ch] = (i - 1)
+             - output_.strobe(ch, output_.strobe_count(ch) - 1);
+      } else {
+	      samples_since_last_strobe_[ch] = UINT_MAX;
+      }
+
+	    if (samples_since_last_strobe_[ch] > parab_var_samples_[ch]) {
+	      decay_constant = last_threshold_[ch] / strobe_decay_samples_;
+	      if (threshold_[ch] > decay_constant)
+		      threshold_[ch] -= decay_constant;
+	      else
+		      threshold_[ch] = 0.0f;
+      } else {
+	      threshold_[ch] = last_threshold_[ch]
+	          * (parab_a_[ch] * pow((samples_since_last_strobe_[ch]
+	                                 / input.sample_rate() + parab_b_[ch]),
+	                                2.0f) + height_);
+      }
+    }
+  }
+
+  PushOutput();
+}
+
+
+
+ModuleParabola::~ModuleParabola() {
+}
+}  // namespace aimc

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Modules/Strobes/ModuleParabola.h	Thu Feb 18 16:55:40 2010 +0000
@@ -0,0 +1,106 @@
+// Copyright 2007-2010, Thomas Walters
+//
+// AIM-C: A C++ implementation of the Auditory Image Model
+// http://www.acousticscale.org/AIMC
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+/*!
+ * \file
+ * \brief SAI 2003 module - for any output frame rate
+ *
+ * \author Tom Walters <[e-mail]>
+ * \date created 2007/03/28
+ * \version \$Id$
+ */
+
+#ifndef _AIMC_MODULE_STROBES_PARABOLA_H_
+#define _AIMC_MODULE_STROBES_PARABOLA_H_
+
+#include <vector>
+
+#include "Support/Module.h"
+#include "Support/Parameters.h"
+#include "Support/SignalBank.h"
+
+/*!
+ * \class ModuleParabola "Modules/SAI/ModuleParabola.h"
+ * \brief SF 2003
+ */
+namespace aimc {
+using std::vector;
+class ModuleParabola : public Module {
+ public:
+	ModuleParabola(Parameters *params);
+	virtual ~ModuleParabola();
+	void Process(const SignalBank& input);
+private:
+  /*! \brief Prepare the module
+	 *  \param input Input signal bank
+	 *  \param output true on success false on failure
+	 */
+	virtual bool InitializeInternal(const SignalBank& input);
+
+	virtual void ResetInternal();
+
+
+	//! \brief Number of samples over which the strobe should be decayed to zero
+	int strobe_decay_samples_;
+
+	/*! \brief Current strobe thresholds, one for each bank channel.
+	 *
+	 *  This value is decayed over time.
+	 */
+	vector<float> threshold_;
+
+	/*! \brief Signal value at the last strobe, one for each bank channel.
+	 *
+	 *  This value is not decayed over time.
+	 */
+	vector<float> last_threshold_;
+
+  //! \brief Parabola height parameter
+  float height_;
+
+  //! \brief Parabola width paramters
+  float parabw_;
+
+  //! \brief Parabola a value
+  vector<float> parab_a_;
+
+  //! \brief Parabola b value
+  vector<float> parab_b_;
+
+  //! \brief Parabola calculation variable
+  vector<float> parab_wnull_;
+
+  //! \brief Parabola calculation variable
+  vector<int> parab_var_samples_;
+
+  //! \brief Storage for the number of samples since the last strobe
+  vector<int> samples_since_last_strobe_;
+
+  vector<float> prev_sample_;
+  vector<float> curr_sample_;
+  vector<float> next_sample_;
+
+  float alpha_;
+  int channel_count_;
+  float sample_rate_;
+  float strobe_decay_time_;
+  int max_strobes_;
+};
+}  // namespace aimc
+
+#endif  // _AIM_MODULE_STROBES_PARABOLA_H_

--- a/src/Support/ERBTools.h	Tue Feb 16 18:40:55 2010 +0000
+++ b/src/Support/ERBTools.h	Thu Feb 18 16:55:40 2010 +0000
@@ -17,7 +17,7 @@
 // along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 /*! \file
- *  \brief
+ *  \brief ERB calculations
  */
 
 /*! \author: Thomas Walters <[e-mail]>
@@ -40,6 +40,10 @@
   static float Freq2ERBw(float freq) {
     return 24.7f * (4.37f * freq / 1000.0f + 1.0f);
   }
+
+  static float ERB2Freq(float erb) {
+    return (pow(10, (erb / 21.4f)) - 1.0f) / 4.37f * 1000.0f;
+  }
 };
 }
 

--- a/src/Support/SignalBank.cc	Tue Feb 16 18:40:55 2010 +0000
+++ b/src/Support/SignalBank.cc	Thu Feb 18 16:55:40 2010 +0000
@@ -84,11 +84,12 @@
 
   centre_frequencies_.resize(channel_count_, 0.0f);
   for (int i = 0; i < channel_count_; ++i) {
-    centre_frequencies_[i] = input.get_centre_frequency(i);
+    centre_frequencies_[i] = input.centre_frequency(i);
   }
 
   for (int i = 0; i < channel_count_; ++i) {
     signals_[i].resize(buffer_length_, 0.0f);
+    strobes_[i].resize(0);
   }
   initialized_ = true;
   return true;
@@ -122,8 +123,20 @@
   signals_[channel][index] = value;
 }
 
-const deque<int> &SignalBank::strobes(int channel) const {
-  return strobes_[channel];
+int SignalBank::strobe(int channel, int index) const {
+  return strobes_[channel][index];
+}
+
+int SignalBank::strobe_count(int channel) const {
+  return strobes_[channel].size();
+}
+
+void SignalBank::AddStrobe(int channel, int time) {
+  strobes_[channel].push_back(time);
+}
+
+void SignalBank::ResetStrobes(int channel) {
+  strobes_[channel].resize(0);
 }
 
 float SignalBank::sample_rate() const {
@@ -142,7 +155,7 @@
   start_time_ = start_time;
 }
 
-float SignalBank::get_centre_frequency(int i) const {
+float SignalBank::centre_frequency(int i) const {
   if (i < channel_count_)
     return centre_frequencies_[i];
   else

--- a/src/Support/SignalBank.h	Tue Feb 16 18:40:55 2010 +0000
+++ b/src/Support/SignalBank.h	Thu Feb 18 16:55:40 2010 +0000
@@ -61,12 +61,15 @@
 
   float sample(int channel, int index) const;
   void set_sample(int channel, int index, float value);
-  const deque<int> &strobes(int channel) const;
+  int strobe(int channel, int index) const;
+  int strobe_count(int channel) const;
+  void AddStrobe(int channel, int time);
+  void ResetStrobes(int channel);
   float sample_rate() const;
   int buffer_length() const;
   int start_time() const;
   void set_start_time(int start_time);
-  float get_centre_frequency(int i) const;
+  float centre_frequency(int i) const;
   void set_centre_frequency(int i, float cf);
   bool initialized() const;
   int channel_count() const;
@@ -74,7 +77,7 @@
   int channel_count_;
   int buffer_length_;
   vector<vector<float> > signals_;
-  vector<deque<int> > strobes_;
+  vector<vector<int> > strobes_;
   vector<float> centre_frequencies_;
   float sample_rate_;
   int start_time_;

--- a/src/Support/StrobeList.h	Tue Feb 16 18:40:55 2010 +0000
+++ b/src/Support/StrobeList.h	Thu Feb 18 16:55:40 2010 +0000
@@ -25,25 +25,25 @@
  * \version \$Id: StrobeList.h 1 2010-02-02 11:04:50Z tcw $
  */
 
-#ifndef _AIMC_STROBE_LIST_H_
-#define _AIMC_STROBE_LIST_H_
+#ifndef _AIMC_SUPPORT_STROBE_LIST_H_
+#define _AIMC_SUPPORT_STROBE_LIST_H_
 
+#include <deque>
 #include <math.h>
 
-#include "Support/util.h"
-
-typedef struct StrobePoint {
-  int iTime;
-  float fWeight;
-  float fWorkingWeight;
+namespace aimc {
+using std::deque;
+struct StrobePoint {
+  int time;
+  float weight;
+  float working_weight;
   StrobePoint() {
-   iTime=0;
-   fWeight=0.0f;
-   fWorkingWeight=0.0f;
+   time = 0;
+   weight = 0.0f;
+   working_weight = 0.0f;
   }
 };
 
-class StrobeList;
 /*!
  * \class Signal "Support/StrobeList.h"
  * \brief Modifiable List of Strobe Points, which must be ordered
@@ -54,90 +54,56 @@
   /*! \brief Create a new strobe list
    */
   inline StrobeList() {
-    m_iStrobeCount = 0;
-    m_iMaxStrobes=0;
-    m_iOffset=0;
-    m_pStrobes = NULL;
+    strobes_.resize(0);
   };
 
-  inline void Create(unsigned int iMaxStrobes) {
-    m_iMaxStrobes = iMaxStrobes;
-    m_pStrobes = new StrobePoint[m_iMaxStrobes+1];
-  }
-
   inline ~StrobeList() {
-    if (m_pStrobes != NULL)
-      delete [] m_pStrobes;
   };
 
   //! \brief Return the strobe time (in samples, can be negative)
-  inline unsigned int getTime(unsigned int iStrobeNumber) {
-    return m_pStrobes[GetBufferNumber(iStrobeNumber)].iTime;
-  };
-
-  //! \brief Return the strobe weight
-  inline float getWeight(unsigned int iStrobeNumber) {
-    return m_pStrobes[GetBufferNumber(iStrobeNumber)].fWeight;
-  };
-
-  //! \brief Return the strobe's working weight
-  inline float getWorkingWeight(unsigned int iStrobeNumber) {
-    return m_pStrobes[GetBufferNumber(iStrobeNumber)].fWorkingWeight;
+  inline StrobePoint Strobe(int strobe_number) {
+    return strobes_.at(strobe_number);
   };
 
     //! \brief Set the strobe weight
-  inline void setWeight(unsigned int iStrobeNumber, float fWeight) {
-    m_pStrobes[GetBufferNumber(iStrobeNumber)].fWeight = fWeight;
+  inline void SetWeight(int strobe_number, float weight) {
+    strobes_.at(strobe_number).weight = weight;
   };
 
     //! \brief Set the strobe's working weight
-  inline void setWorkingWeight(unsigned int iStrobeNumber,
-                               float fWorkingWeight) {
-    m_pStrobes[GetBufferNumber(iStrobeNumber)].fWorkingWeight = fWorkingWeight;
+  inline void SetWorkingWeight(int strobe_number, float working_weight) {
+    strobes_.at(strobe_number).working_weight = working_weight;
   };
 
   //! \brief Add a strobe to the list (must be in order)
-  inline void addStrobe(int iTime, float fWeight) {
-    if (m_iStrobeCount + 1 <= m_iMaxStrobes) {
-      m_iStrobeCount++;
-      m_pStrobes[GetBufferNumber(m_iStrobeCount)].iTime=iTime;
-      m_pStrobes[GetBufferNumber(m_iStrobeCount)].fWeight=fWeight;
-    }
+  inline void AddStrobe(int time, float weight) {
+    StrobePoint s;
+    s.time = time;
+    s.weight = weight;
+    strobes_.push_back(s);
   };
 
   //! \brief Delete a strobe from the list
-  inline void deleteFirstStrobe() {
-    if (m_iStrobeCount > 0) {
-      m_iOffset = (m_iOffset+1) % m_iMaxStrobes;
-      m_iStrobeCount--;
-    }
+  inline void DeleteFirstStrobe() {
+    strobes_.pop_front();
   };
 
   //! \brief Get the number of strobes
-  inline unsigned int getStrobeCount() {
-    return m_iStrobeCount;
+  inline int strobe_count() const {
+    return strobes_.size();
   };
 
-  //! \brief Shift the position of all strobes by subtracting iOffset from
+  //! \brief Shift the position of all strobes by subtracting offset from
   //! the time value of each
-  inline void shiftStrobes(unsigned int iOffset) {
-    for (unsigned int i=1; i<m_iStrobeCount+1; i++)
-      m_pStrobes[GetBufferNumber(i)].iTime-=iOffset;
+  inline void ShiftStrobes(int offset) {
+    for (unsigned int i = 0; i < strobes_.size(); ++i)
+      strobes_[i].time -= offset;
   };
 
- protected:
-
  private:
-  unsigned int m_iStrobeCount;
-  unsigned int m_iMaxStrobes;
-  unsigned int m_iOffset;
-  StrobePoint* m_pStrobes;
-
-  inline unsigned int GetBufferNumber(unsigned int iStrobeIndex) {
-    aimASSERT(((iStrobeIndex-1+m_iOffset) % m_iMaxStrobes)<m_iMaxStrobes);
-    return ((iStrobeIndex-1+m_iOffset) % m_iMaxStrobes) ;
-  };
+  deque<StrobePoint> strobes_;
 };
+}  // namespace aimc
 
 #endif /* _AIMC_STROBE_LIST_H_ */
 

--- a/swig/aim_modules.i	Tue Feb 16 18:40:55 2010 +0000
+++ b/swig/aim_modules.i	Thu Feb 18 16:55:40 2010 +0000
@@ -1,4 +1,4 @@
-// Copyright 2008-2010, Thomas Walters
+// Copyright 2010, Thomas Walters
 //
 // AIM-C: A C++ implementation of the Auditory Image Model
 // http://www.acousticscale.org/AIMC
@@ -24,6 +24,7 @@
 #include "Support/Parameters.h"
 #include "Support/SignalBank.h"
 #include "Modules/Features/ModuleGaussians.h"
+#include "Modules/BMM/ModuleGammatone.h"
 %}
 
 namespace aimc {
@@ -100,5 +101,14 @@
 	virtual void Process(const SignalBank &input);
 	void Reset();
 };
+
+class ModuleGammatone : public Module
+{
+ public:
+	ModuleGammatone(Parameters *pParam);
+	virtual ~ModuleGammatone();
+	virtual void Process(const SignalBank &input);
+	void Reset();
+};
 }  // namespace aimc
 

--- a/swig/setup.py	Tue Feb 16 18:40:55 2010 +0000
+++ b/swig/setup.py	Thu Feb 18 16:55:40 2010 +0000
@@ -29,7 +29,8 @@
                                    '../src/Support/Parameters.cc',
                                    '../src/Support/SignalBank.cc', 
                                    '../src/Support/Module.cc',
-                                   '../src/Modules/Features/ModuleGaussians.cc'],
+                                   '../src/Modules/Features/ModuleGaussians.cc',
+                                   '../src/Modules/BMM/ModuleGammatone.cc'],
                         swig_opts = ['-c++','-I../src/'], 
                         include_dirs=['../src/']
                         )