Cepstral features: improve documentation

Author: Laurent El Shafey

include/bob/ap/Ceps.h

@@ -29,9 +29,6 @@
 #include <vector>
 #include "bob/sp/FFT1D.h"
 
-const double ENERGY_FLOOR = 1.0;
-const double FBANK_OUT_FLOOR = 1.0;
-
 namespace bob {
 /**
  * \ingroup libap_api
@@ -55,21 +52,21 @@ class Ceps
 {
   public:
     /**
-     * @brief Constructor: Initialize working arrays
+     * @brief Constructor. Initializes working arrays
      */
     Ceps(double sampling_frequency, double win_length_ms=20., double win_shift_ms=10.,
       size_t n_filters=24, size_t n_ceps=19, double f_min=0., 
       double f_max=4000., size_t delta_win=2, double pre_emphasis_coef=0.95,
       bool mel_scale=true, bool dct_norm=false);
 
     /**
-     * @brief Get the Cepstral Shape
+     * @brief Gets the Cepstral features shape for a given input/input length
      */
     blitz::TinyVector<int,2> getCepsShape(const size_t input_length) const;
     blitz::TinyVector<int,2> getCepsShape(const blitz::Array<double,1>& input) const;
 
     /**
-     * @brief Compute Cepstral features
+     * @brief Computes Cepstral features
      */
     void operator()(const blitz::Array<double,1>& input, blitz::Array<double,2>& output);
 
@@ -245,7 +242,7 @@ class Ceps
      * @brief Computes the first order derivative from the given input. 
      * This methods is used to compute both the delta's and double delta's.
      */
-    void addDerivative(const blitz::Array<double,2>& input, blitz::Array<double,2>& output);
+    void addDerivative(const blitz::Array<double,2>& input, blitz::Array<double,2>& output) const;
     /**
      * @brief Converts a frequency in Herz to the corresponding one in Mel
      */
@@ -258,27 +255,32 @@ class Ceps
      * @brief Pre-emphasises the signal by applying the first order equation
      * \f$data_{n} := data_{n} − a*data_{n−1}\f$
      */
-    void pre_emphasis(blitz::Array<double,1> &data);
+    void pre_emphasis(blitz::Array<double,1> &data) const;
     /**
      * @brief Applies the Hamming window to the signal
      */
-    void hammingWindow(blitz::Array<double,1> &data);
+    void hammingWindow(blitz::Array<double,1> &data) const;
 
+    /**
+     * @brief Computes the power-spectrum of the FFT of the input frame and
+     * applies the triangular filter bank
+     */
     void logFilterBank(blitz::Array<double,1>& x);
     /**
-     * @brief Apply triangular filter bank to the input array and return the log
-     * of the energy in each band.
+     * @brief Applies the triangular filter bank to the input array and 
+     * returns the logarithm of the energy in each band.
      */
-    void logTriangularFBank(blitz::Array<double,1>& data);
-    double logEnergy(blitz::Array<double,1> &data);
+    void logTriangularFilterBank(blitz::Array<double,1>& data) const;
     /**
-     * @brief Apply a p order DCT to vector v1.
-     * Results are returned through v2.
-     * If {m[1],...,m[N]} are the output of the filters, then
-     *    c[i]=sqrt(2/N)*sum for j=1 to N of(m[j]cos(M_PI*i*(j-0.5)/N) i=1,...,p
-     * This is what is implemented here with arrays indexed from 0 to N-1.
+     * @brief Computes the logarithm of the energy
      */
-    void transformDCT(blitz::Array<double,1>& ceps_row);
+    double logEnergy(blitz::Array<double,1> &data) const;
+    /**
+     * @brief Applies the DCT to the cepstral features:
+     * \f$out[i]=sqrt(2/N)*sum_{j=1}^{N} (in[j]cos(M_PI*i*(j-0.5)/N)\f$
+     */
+    void applyDct(blitz::Array<double,1>& ceps_row) const;
+
     void initWinSize();
     void initWinLength();
     void initWinShift();
@@ -287,8 +289,8 @@ class Ceps
     void initCacheDctKernel();
     void initCacheFilterBank();
     /**
-     * @brief Initialize the table m_p_index, which contains the indices of the
-     * cut-off frequencies. It looks like something like this:
+     * @brief Initialize the table m_p_index, which contains the indices of
+     * the cut-off frequencies of the triangular filters.. It looks like:
      *
      *                      filter 2
      *                   <------------->
@@ -298,8 +300,8 @@ class Ceps
      *         0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9  ..........
      *             ^     ^     ^       ^             ^
      *             |     |     |       |             |
-     *          p_in[0]  |  p_in[2]    |          p_in[4]
-     *                p_in[1]       p_in[3]
+     *            t[0]   |    t[2]     |           t[4]
+     *                  t[1]          t[3]
      *
      */
     void initCachePIndex();
@@ -322,16 +324,21 @@ class Ceps
     bool m_with_energy;
     bool m_with_delta;
     bool m_with_delta_delta;
+    double m_energy_floor;
+    double m_fb_out_floor;
+    double m_log_energy_floor;
+    double m_log_fb_out_floor;
+
     blitz::Array<double,2> m_dct_kernel;
     blitz::Array<double,1> m_hamming_kernel;
-    blitz::Array<double,1> m_filters;
     blitz::Array<int,1>  m_p_index;
     std::vector<blitz::Array<double,1> > m_filter_bank;
     bob::sp::FFT1D m_fft;
 
-    mutable blitz::Array<double,1> m_cache_frame;
-    mutable blitz::Array<std::complex<double>,1>  m_cache_complex1;
-    mutable blitz::Array<std::complex<double>,1>  m_cache_complex2;
+    mutable blitz::Array<double,1> m_cache_frame_d;
+    mutable blitz::Array<std::complex<double>,1>  m_cache_frame_c1;
+    mutable blitz::Array<std::complex<double>,1>  m_cache_frame_c2;
+    mutable blitz::Array<double,1> m_cache_filters;
 
     friend class TestCeps;
 };
@@ -349,16 +356,16 @@ class TestCeps
     { return m_ceps.getCepsShape(input_length); }
     blitz::TinyVector<int,2> getCepsShape(const blitz::Array<double,1>& input) const
     { return m_ceps.getCepsShape(input); }
-    blitz::Array<double,1> getFilter(void) { return m_ceps.m_filters; }
+    blitz::Array<double,1> getFilterOutput() { return m_ceps.m_cache_filters; }
 
     void operator()(const blitz::Array<double,1>& input, blitz::Array<double,2>& ceps_2D)
     { m_ceps(input, ceps_2D);}
     void hammingWindow(blitz::Array<double,1>& data){ m_ceps.hammingWindow(data); }
     void pre_emphasis(blitz::Array<double,1>& data){ m_ceps.pre_emphasis(data); }
     void logFilterBank(blitz::Array<double,1>& x){ m_ceps.logFilterBank(x); }
-    void logTriangularFBank(blitz::Array<double,1>& data){ m_ceps.logTriangularFBank(data); }
+    void logTriangularFilterBank(blitz::Array<double,1>& data){ m_ceps.logTriangularFilterBank(data); }
     double logEnergy(blitz::Array<double,1> &data){ return m_ceps.logEnergy(data); }
-    void transformDCT(blitz::Array<double,1>& ceps_row) { m_ceps.transformDCT(ceps_row); }
+    void applyDct(blitz::Array<double,1>& ceps_row) { m_ceps.applyDct(ceps_row); }
 };
 
 }

python/bob/ap/test/test_ceps.py

@@ -248,7 +248,7 @@ def cepstral_features_extraction(obj, rate_wavsample, win_length_ms, win_shift_m
       obj.assertAlmostEqual(filters[kk], filt2[kk], 7, "Error in log Filtering")
 
     ceps = dct_transform(filters, n_filters, dct_kernel, n_ceps, dct_norm)
-    ceps2 = ct.dct_transform(n_ceps)
+    ceps2 = ct.apply_dct(n_ceps)
 
 
     if(with_energy):
@@ -453,4 +453,4 @@ def test_cepstral(self):
     with_delta_delta = False
     cepstral_comparison_run(self,rate_wavsample, win_length_ms, win_shift_ms, n_filters, n_ceps, dct_norm, f_min, f_max, delta_win,
                                pre_emphasis_coef, mel_scale, with_energy, with_delta, with_delta_delta)
-    
+    

src/ap/cxx/Ceps.cc

@@ -34,13 +34,16 @@ bob::ap::Ceps::Ceps( double sampling_frequency, double win_length_ms, double win
   m_delta_win(delta_win), m_pre_emphasis_coeff(pre_emphasis_coeff),
   m_mel_scale(mel_scale), m_dct_norm(dct_norm),
   m_with_energy(true), m_with_delta(true), m_with_delta_delta(true),
-  m_filter_bank(), m_fft(1)
+  m_energy_floor(1.), m_fb_out_floor(1.), m_fft(1)
 {
   initWinLength();
   initWinShift();
 
-  m_filters.resize(m_n_filters);
+  // Initializes logarithm of flooring values
+  m_log_energy_floor = log(m_energy_floor);
+  m_log_fb_out_floor = log(m_fb_out_floor);
 
+  m_cache_filters.resize(m_n_filters);
   initCacheDctKernel();
 }
 
@@ -70,7 +73,7 @@ void bob::ap::Ceps::setWinShiftMs(double win_shift_ms)
 void bob::ap::Ceps::setNFilters(size_t n_filters)
 { 
   m_n_filters = n_filters; 
-  m_filters.resize(m_n_filters); 
+  m_cache_filters.resize(m_n_filters); 
   initCacheFilterBank(); 
   initCacheDctKernel(); 
 }
@@ -132,10 +135,10 @@ void bob::ap::Ceps::initWinShift()
 void bob::ap::Ceps::initWinSize()
 {
   m_win_size = (size_t)pow(2.0,ceil(log((double)m_win_length)/log(2)));
-  m_cache_frame.resize(m_win_size);
+  m_cache_frame_d.resize(m_win_size);
   m_fft.reset(m_win_size);
-  m_cache_complex1.resize(m_win_size);
-  m_cache_complex2.resize(m_win_size);
+  m_cache_frame_c1.resize(m_win_size);
+  m_cache_frame_c2.resize(m_win_size);
 }
 
 void bob::ap::Ceps::initCacheHammingKernel()
@@ -262,26 +265,26 @@ void bob::ap::Ceps::operator()(const blitz::Array<double,1>& input,
   for(int i=0; i<n_frames; ++i) 
   {
     // Set padded frame to zero
-    m_cache_frame = 0.;
+    m_cache_frame_d = 0.;
     // Extract frame input vector
     blitz::Range ri(i*(int)m_win_shift,i*(int)m_win_shift+(int)m_win_length-1);
-    m_cache_frame(rf) = input(ri);
+    m_cache_frame_d(rf) = input(ri);
     // Substract mean value
-    m_cache_frame -= blitz::mean(m_cache_frame);
+    m_cache_frame_d -= blitz::mean(m_cache_frame_d);
 
     // Update output with energy if required
     if(m_with_energy)
-      ceps_matrix(i,(int)m_n_ceps) = logEnergy(m_cache_frame);
+      ceps_matrix(i,(int)m_n_ceps) = logEnergy(m_cache_frame_d);
 
     // Apply pre-emphasis
-    pre_emphasis(m_cache_frame);
+    pre_emphasis(m_cache_frame_d);
     // Apply the Hamming window
-    hammingWindow(m_cache_frame);
+    hammingWindow(m_cache_frame_d);
     // Filter with the triangular filter bank (either in linear or Mel domain)
-    logFilterBank(m_cache_frame);
+    logFilterBank(m_cache_frame_d);
     // Apply DCT kernel and update the output 
     blitz::Array<double,1> ceps_matrix_row(ceps_matrix(i,r1));
-    transformDCT(ceps_matrix_row);
+    applyDct(ceps_matrix_row);
   }
 
   blitz::Range rall = blitz::Range::all();
@@ -302,7 +305,7 @@ void bob::ap::Ceps::operator()(const blitz::Array<double,1>& input,
   }
 }
 
-void bob::ap::Ceps::pre_emphasis(blitz::Array<double,1> &data)
+void bob::ap::Ceps::pre_emphasis(blitz::Array<double,1> &data) const
 {
   if(m_pre_emphasis_coeff!=0.)
   { 
@@ -315,7 +318,7 @@ void bob::ap::Ceps::pre_emphasis(blitz::Array<double,1> &data)
   }
 }
 
-void bob::ap::Ceps::hammingWindow(blitz::Array<double,1> &data)
+void bob::ap::Ceps::hammingWindow(blitz::Array<double,1> &data) const
 {
   blitz::Range r(0,(int)m_win_length-1);
   data(r) *= m_hamming_kernel;
@@ -324,46 +327,44 @@ void bob::ap::Ceps::hammingWindow(blitz::Array<double,1> &data)
 void bob::ap::Ceps::logFilterBank(blitz::Array<double,1>& x)
 {
   // Apply the FFT
-  m_cache_complex1 = bob::core::cast<std::complex<double> >(x);
-  m_fft(m_cache_complex1, m_cache_complex2);
+  m_cache_frame_c1 = bob::core::cast<std::complex<double> >(x);
+  m_fft(m_cache_frame_c1, m_cache_frame_c2);
 
   // Take the the power spectrum of the first part of the output of the FFT
   blitz::Range r(0,(int)m_win_size/2);
   blitz::Array<double,1> x_half(x(r));
-  blitz::Array<std::complex<double>,1> complex_half(m_cache_complex2(r));
+  blitz::Array<std::complex<double>,1> complex_half(m_cache_frame_c2(r));
   x_half = blitz::abs(complex_half);
 
   // Apply the Triangular filter bank to this power spectrum
-  logTriangularFBank(x);
+  logTriangularFilterBank(x);
 }
 
-void bob::ap::Ceps::logTriangularFBank(blitz::Array<double,1>& data)
+void bob::ap::Ceps::logTriangularFilterBank(blitz::Array<double,1>& data) const
 {
   for(int i=0; i<(int)m_n_filters; ++i)
   {
     blitz::Array<double,1> data_slice(data(blitz::Range(m_p_index(i),m_p_index(i+2))));
     double res = blitz::sum(data_slice * m_filter_bank[i]);
-    m_filters(i)=(res < FBANK_OUT_FLOOR)?(double)log(FBANK_OUT_FLOOR):(double)log(res);
+    m_cache_filters(i)= (res < m_fb_out_floor ? m_log_fb_out_floor : log(res));
   }
 }
 
-double bob::ap::Ceps::logEnergy(blitz::Array<double,1> &data)
+double bob::ap::Ceps::logEnergy(blitz::Array<double,1> &data) const
 {
   blitz::Array<double,1> data_p(data(blitz::Range(0,(int)m_win_length-1)));
   double gain = blitz::sum(blitz::pow2(data_p));
-  gain = gain < ENERGY_FLOOR ?
-      (double)(log(ENERGY_FLOOR)) : (double)(log(gain));
-  return (gain);
+  return (gain < m_energy_floor ? m_log_energy_floor : log(gain)); 
 }
 
-void bob::ap::Ceps::transformDCT(blitz::Array<double,1>& ceps_row)
+void bob::ap::Ceps::applyDct(blitz::Array<double,1>& ceps_row) const
 {
   blitz::firstIndex i;
   blitz::secondIndex j;
-  ceps_row = blitz::sum(m_filters(j) * m_dct_kernel(i,j), j);
+  ceps_row = blitz::sum(m_cache_filters(j) * m_dct_kernel(i,j), j);
 }
 
-void bob::ap::Ceps::addDerivative(const blitz::Array<double,2>& input, blitz::Array<double,2>& output)
+void bob::ap::Ceps::addDerivative(const blitz::Array<double,2>& input, blitz::Array<double,2>& output) const
 {
   // Initialize output to zero
   output = 0.;

src/ap/python/ceps.cc

@@ -31,7 +31,7 @@ using namespace boost::python;
 static const char* CEPS_DOC = "Objects of this class, after configuration, can extract Cepstral Features from a 1D array/signal.";
 static const char* TESTCEPS_DOC = "Objects of this class, after configuration, can be used to test the private methods of bob.ap.Ceps.";
 
-static object py_ceps_analysis(bob::ap::Ceps& ceps, bob::python::const_ndarray input)
+static object py_forward(bob::ap::Ceps& ceps, bob::python::const_ndarray input)
 {
   // Gets the shape of the feature
   const blitz::Array<double,1> input_ = input.bz<double,1>();
@@ -64,8 +64,6 @@ static boost::python::tuple py_get_ceps_shape(bob::ap::Ceps& ceps, object input_
 static double py_logEnergy(bob::ap::TestCeps& ceps, bob::python::ndarray data)
 {
   blitz::Array<double,1> data_ = data.bz<double,1>();
-
-  // Get the logEnergy
   return ceps.logEnergy(data_);
 }
 
@@ -88,17 +86,17 @@ static object py_logFilterBank(bob::ap::TestCeps& ceps, bob::python::ndarray dat
   ceps.logFilterBank(data_);
   bob::python::ndarray filter(bob::core::array::t_float64, n_filters);
   blitz::Array<double,1> filter_ = filter.bz<double,1>();
-  // Gets the filter Bank
-  filter_ = ceps.getFilter();
+  // Gets the filter bank output
+  filter_ = ceps.getFilterOutput();
   return filter.self();
 }
 
-static object py_transformDCT(bob::ap::TestCeps& ceps, int n_ceps)
+static object py_applyDct(bob::ap::TestCeps& ceps, int n_ceps)
 {
   bob::python::ndarray ceps_row(bob::core::array::t_float64, n_ceps);
   blitz::Array<double,1> ceps_row_ = ceps_row.bz<double,1>();
-  // Get the Cepstral features
-  ceps.transformDCT(ceps_row_);
+  // Gets the Cepstral features
+  ceps.applyDct(ceps_row_);
   return ceps_row.self();
 }
 
@@ -107,7 +105,7 @@ void bind_ap_ceps()
 {
   class_<bob::ap::Ceps, boost::shared_ptr<bob::ap::Ceps> >("Ceps", CEPS_DOC, init<double, optional<double, double, size_t, size_t, double, double, int, double, bool, bool> >
   ((arg("sampling_frequency"), arg("win_length_ms")=20., arg("win_shift_ms")=10., arg("n_filters")=24, arg("n_ceps")=19, arg("f_min")=0., arg("f_max")=4000., arg("delta_win"), arg("pre_emphasis_coeff")=0.95, arg("mel_scale")=true, arg("dct_norm")=false)))
-        .add_property("sampling_frequency", &bob::ap::Ceps::getSamplingFrequency, &bob::ap::Ceps::setSamplingFrequency, "The sample frequency of the input data")
+        .add_property("sampling_frequency", &bob::ap::Ceps::getSamplingFrequency, &bob::ap::Ceps::setSamplingFrequency, "The sampling frequency of the input data")
         .add_property("win_length_ms", &bob::ap::Ceps::getWinLengthMs, &bob::ap::Ceps::setWinLengthMs, "The window length of the cepstral analysis in milliseconds")
         .add_property("win_length", &bob::ap::Ceps::getWinLength, "The normalized window length wrt. to the sample frequency")
         .add_property("win_shift_ms", &bob::ap::Ceps::getWinShiftMs, &bob::ap::Ceps::setWinShiftMs, "The window shift of the cepstral analysis in milliseconds")
@@ -123,8 +121,8 @@ void bind_ap_ceps()
         .add_property("with_energy", &bob::ap::Ceps::getWithEnergy, &bob::ap::Ceps::setWithEnergy, "Tells if we add the energy to the output feature")
         .add_property("with_delta", &bob::ap::Ceps::getWithDelta, &bob::ap::Ceps::setWithDelta, "Tells if we add the first derivatives to the output feature")
         .add_property("with_delta_delta", &bob::ap::Ceps::getWithDeltaDelta, &bob::ap::Ceps::setWithDeltaDelta, "Tells if we add the second derivatives to the output feature")
-        .def("__call__", &py_ceps_analysis, (arg("input")), "Compute the features")
-        .def("get_ceps_shape", &py_get_ceps_shape, (arg("n_size"), arg("input_data")), "Compute the shape of the output features")
+        .def("__call__", &py_forward, (arg("input")), "Computes the cepstral features")
+        .def("get_ceps_shape", &py_get_ceps_shape, (arg("n_size"), arg("input_data")), "Computes the shape of the output features")
         ;
 
   class_<bob::ap::TestCeps, boost::shared_ptr<bob::ap::TestCeps> >("TestCeps", TESTCEPS_DOC, init<bob::ap::Ceps&>((arg("ceps"))))
@@ -134,7 +132,7 @@ void bind_ap_ceps()
         .def("pre_emphasis", &py_emphasis, (arg("data")), "compute pre-emphasis")
         .def("hamming_window", &py_hammingWindow, (arg("data")), "compute the wraped signal on a hamming Window")
         .def("log_filter_bank", &py_logFilterBank, (arg("data"), arg("m_win_size"), arg("n_filters")), "compute log Filter Bank")
-        .def("dct_transform", &py_transformDCT, (arg("n_ceps")), "DCT Transform")
+        .def("apply_dct", &py_applyDct, (arg("n_ceps")), "DCT Transform")
       ;
 }