tomwalters@268: // Copyright 2008-2010, Thomas Walters
tomwalters@268: //
tomwalters@268: // AIM-C: A C++ implementation of the Auditory Image Model
tomwalters@268: // http://www.acousticscale.org/AIMC
tomwalters@268: //
tomwalters@318: // Licensed under the Apache License, Version 2.0 (the "License");
tomwalters@318: // you may not use this file except in compliance with the License.
tomwalters@318: // You may obtain a copy of the License at
tomwalters@268: //
tomwalters@318: //     http://www.apache.org/licenses/LICENSE-2.0
tomwalters@268: //
tomwalters@318: // Unless required by applicable law or agreed to in writing, software
tomwalters@318: // distributed under the License is distributed on an "AS IS" BASIS,
tomwalters@318: // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
tomwalters@318: // See the License for the specific language governing permissions and
tomwalters@318: // limitations under the License.
tomwalters@268: 
tomwalters@268: /*! \file
tomwalters@268:  *  \brief Gaussian features - based on MATLAB code by Christian Feldbauer
tomwalters@268:  */
tomwalters@268: 
tomwalters@268: /*!
tomwalters@273:  * \author Thomas Walters <tom@acousticscale.org>
tomwalters@268:  * \date created 2008/06/23
tomwalters@296:  * \version \$Id$
tomwalters@268:  */
tomwalters@268: 
tomwalters@268: #include <math.h>
tomwalters@268: 
tom@440: #ifdef _WINDOWS
tomwalters@427: #include <float.h>
tomwalters@427: #endif
tomwalters@427: 
tomwalters@268: #include "Modules/Features/ModuleGaussians.h"
tomwalters@268: #include "Support/Common.h"
tomwalters@268: 
tomwalters@268: namespace aimc {
tomwalters@278: ModuleGaussians::ModuleGaussians(Parameters *params) : Module(params) {
tomwalters@268:   // Set module metadata
tomwalters@268:   module_description_ = "Gaussian Fitting to SSI profile";
tomwalters@273:   module_identifier_ = "gaussians";
tomwalters@268:   module_type_ = "features";
tomwalters@296:   module_version_ = "$Id$";
tomwalters@268: 
tomwalters@365:   m_iParamNComp = parameters_->DefaultInt("gaussians.ncomp", 4);
tomwalters@365:   m_fParamVar = parameters_->DefaultFloat("gaussians.var", 115.0);
tomwalters@365:   m_fParamPosteriorExp = parameters_->DefaultFloat("gaussians.posterior_exp",
tomwalters@365:                                                    6.0);
tomwalters@365:   m_iParamMaxIt = parameters_->DefaultInt("gaussians.maxit", 250);
tomwalters@268: 
tomwalters@273:   // The parameters system doesn't support tiny numbers well, to define this
tomwalters@273:   // variable as a string, then convert it to a float afterwards
tomwalters@365:   parameters_->DefaultString("gaussians.priors_converged", "1e-7");
tomwalters@365:   priors_converged_ = parameters_->GetFloat("gaussians.priors_converged");
tomwalters@365:   output_positions_ = parameters_->DefaultBool("gaussians.positions", false);
tomwalters@268: }
tomwalters@268: 
tomwalters@268: ModuleGaussians::~ModuleGaussians() {
tomwalters@268: }
tomwalters@268: 
tomwalters@268: bool ModuleGaussians::InitializeInternal(const SignalBank &input) {
tomwalters@268:   m_pA.resize(m_iParamNComp, 0.0f);
tomwalters@268:   m_pMu.resize(m_iParamNComp, 0.0f);
tomwalters@268: 
tomwalters@268:   // Assuming the number of channels is greater than twice the number of
tomwalters@268:   // Gaussian components, this is ok
tomwalters@365:   output_component_count_ = 1; // Energy component
tomwalters@268:   if (input.channel_count() >= 2 * m_iParamNComp) {
tomwalters@365:     output_component_count_ += (m_iParamNComp - 1);
tomwalters@268:   } else {
tomwalters@268:     LOG_ERROR(_T("Too few channels in filterbank to produce sensible "
tomwalters@268:                  "Gaussian features. Either increase the number of filterbank"
tomwalters@268:                  " channels, or decrease the number of Gaussian components"));
tomwalters@268:     return false;
tomwalters@268:   }
tomwalters@268: 
tomwalters@365:   if (output_positions_) {
tomwalters@365:     output_component_count_ += m_iParamNComp;
tomwalters@365:   }
tomwalters@365: 
tomwalters@365:   output_.Initialize(output_component_count_, 1, input.sample_rate());
tomwalters@365: 
tomwalters@268:   m_iNumChannels = input.channel_count();
tomwalters@268:   m_pSpectralProfile.resize(m_iNumChannels, 0.0f);
tomwalters@268: 
tomwalters@268:   return true;
tomwalters@268: }
tomwalters@268: 
tomwalters@275: void ModuleGaussians::ResetInternal() {
tomwalters@268:   m_pSpectralProfile.clear();
tomwalters@268:   m_pSpectralProfile.resize(m_iNumChannels, 0.0f);
tomwalters@292:   m_pA.clear();
tomwalters@292:   m_pA.resize(m_iParamNComp, 0.0f);
tomwalters@292:   m_pMu.clear();
tomwalters@292:   m_pMu.resize(m_iParamNComp, 0.0f);
tomwalters@268: }
tomwalters@268: 
tomwalters@268: void ModuleGaussians::Process(const SignalBank &input) {
tomwalters@273:   if (!initialized_) {
tomwalters@273:     LOG_ERROR(_T("Module ModuleGaussians not initialized."));
tomwalters@273:     return;
tomwalters@273:   }
tom@421:   output_.set_start_time(input.start_time());
tomwalters@268:   // Calculate spectral profile
tomwalters@365:   for (int ch = 0; ch < input.channel_count(); ++ch) {
tomwalters@365:     m_pSpectralProfile[ch] = 0.0f;
tomwalters@365:     for (int i = 0; i < input.buffer_length(); ++i) {
tomwalters@365:       m_pSpectralProfile[ch] += input[ch][i];
tomwalters@268:     }
tomwalters@365:     m_pSpectralProfile[ch] /= static_cast<float>(input.buffer_length());
tomwalters@273:   }
tomwalters@273: 
tomwalters@280:   float spectral_profile_sum = 0.0f;
tomwalters@273:   for (int i = 0; i < input.channel_count(); ++i) {
tomwalters@273:     spectral_profile_sum += m_pSpectralProfile[i];
tomwalters@273:   }
tomwalters@273: 
tomwalters@365:   // Set the last component of the feature vector to be the log energy
tomwalters@280:   float logsum = log(spectral_profile_sum);
tomwalters@273:   if (!isinf(logsum)) {
tomwalters@365:     output_.set_sample(output_component_count_ - 1, 0, logsum);
tomwalters@273:   } else {
tomwalters@365:     output_.set_sample(output_component_count_ - 1, 0, -1000.0);
tomwalters@268:   }
tomwalters@268: 
tomwalters@365:   for (int ch = 0; ch < input.channel_count(); ++ch) {
tomwalters@427:     m_pSpectralProfile[ch] = pow(m_pSpectralProfile[ch], 0.8f);
tomwalters@268:   }
tomwalters@268: 
tomwalters@268:   RubberGMMCore(2, true);
tomwalters@268: 
tomwalters@365:   float mean1 = m_pMu[0];
tomwalters@365:   float mean2 = m_pMu[1];
tomwalters@280:   // LOG_INFO(_T("Orig. mean 0 = %f"), m_pMu[0]);
tomwalters@280:   // LOG_INFO(_T("Orig. mean 1 = %f"), m_pMu[1]);
tomwalters@280:   // LOG_INFO(_T("Orig. prob 0 = %f"), m_pA[0]);
tomwalters@280:   // LOG_INFO(_T("Orig. prob 1 = %f"), m_pA[1]);
tomwalters@268: 
tomwalters@365:   float a1 = 0.05 * m_pA[0];
tomwalters@365:   float a2 = 1.0 - 0.25 * m_pA[1];
tomwalters@268: 
tomwalters@280:   // LOG_INFO(_T("fA1 = %f"), fA1);
tomwalters@280:   // LOG_INFO(_T("fA2 = %f"), fA2);
tomwalters@274: 
tomwalters@365:   float gradient = (mean2 - mean1) / (a2 - a1);
tomwalters@365:   float intercept = mean2 - gradient * a2;
tomwalters@274: 
tomwalters@280:   // LOG_INFO(_T("fGradient = %f"), fGradient);
tomwalters@280:   // LOG_INFO(_T("fIntercept = %f"), fIntercept);
tomwalters@268: 
tomwalters@268:   for (int i = 0; i < m_iParamNComp; ++i) {
tomwalters@280:     m_pMu[i] = (static_cast<float>(i)
tomwalters@280:                 / (static_cast<float>(m_iParamNComp) - 1.0f))
tomwalters@365:                 * gradient + intercept;
tomwalters@280:                 // LOG_INFO(_T("mean %d = %f"), i, m_pMu[i]);
tomwalters@268:   }
tomwalters@268: 
tomwalters@268:   for (int i = 0; i < m_iParamNComp; ++i) {
tomwalters@280:     m_pA[i] = 1.0f / static_cast<float>(m_iParamNComp);
tomwalters@268:   }
tomwalters@268: 
tomwalters@268:   RubberGMMCore(m_iParamNComp, false);
tomwalters@268: 
tomwalters@365:   // Amplitudes first
tomwalters@268:   for (int i = 0; i < m_iParamNComp - 1; ++i) {
tomwalters@268:     if (!isnan(m_pA[i])) {
tomwalters@268:       output_.set_sample(i, 0, m_pA[i]);
tomwalters@268:     } else {
tomwalters@268:       output_.set_sample(i, 0, 0.0f);
tomwalters@268:     }
tomwalters@268:   }
tomwalters@273: 
tomwalters@365:   // Then means if required
tomwalters@365:   if (output_positions_) {
tomwalters@365:     int idx = 0;
tomwalters@365:     for (int i = m_iParamNComp - 1; i < 2 * m_iParamNComp - 1; ++i) {
tomwalters@365:       if (!isnan(m_pMu[i])) {
tomwalters@365:         output_.set_sample(i, 0, m_pMu[idx]);
tomwalters@365:       } else {
tomwalters@365:         output_.set_sample(i, 0, 0.0f);
tomwalters@365:       }
tomwalters@365:       ++idx;
tomwalters@365:     }
tomwalters@365:   }
tomwalters@365: 
tomwalters@268:   PushOutput();
tomwalters@268: }
tomwalters@268: 
tomwalters@268: bool ModuleGaussians::RubberGMMCore(int iNComponents, bool bDoInit) {
tomwalters@268:   int iSizeX = m_iNumChannels;
tomwalters@268: 
tomwalters@268:   // Normalise the spectral profile
tomwalters@365:   float SpectralProfileTotal = 0.0f;
tomwalters@268:   for (int iCount = 0; iCount < iSizeX; iCount++) {
tomwalters@365:     SpectralProfileTotal += m_pSpectralProfile[iCount];
tomwalters@268:   }
tomwalters@268:   for (int iCount = 0; iCount < iSizeX; iCount++) {
tomwalters@365:     m_pSpectralProfile[iCount] /= SpectralProfileTotal;
tomwalters@268:   }
tomwalters@268: 
tomwalters@268:   if (bDoInit) {
tomwalters@268:     // Uniformly spaced components
tomwalters@280:     float dd = (iSizeX - 1.0f) / iNComponents;
tomwalters@268:     for (int i = 0; i < iNComponents; i++) {
tomwalters@268:       m_pMu[i] = dd / 2.0f + (i * dd);
tomwalters@268:       m_pA[i] = 1.0f / iNComponents;
tomwalters@268:     }
tomwalters@268:   }
tomwalters@268: 
tomwalters@280:   vector<float> pA_old;
tomwalters@268:   pA_old.resize(iNComponents);
tomwalters@280:   vector<float> pP_mod_X;
tomwalters@268:   pP_mod_X.resize(iSizeX);
tomwalters@280:   vector<float> pP_comp;
tomwalters@268:   pP_comp.resize(iSizeX * iNComponents);
tomwalters@268: 
tomwalters@268:   for (int iIteration = 0; iIteration < m_iParamMaxIt; iIteration++) {
tomwalters@268:     // (re)calculate posteriors (component probability given observation)
tomwalters@268:     // denominator: the model density at all observation points X
tomwalters@268:     for (int i = 0; i < iSizeX; ++i) {
tomwalters@268:       pP_mod_X[i] = 0.0f;
tomwalters@268:     }
tomwalters@268: 
tomwalters@365:     for (int c = 0; c < iNComponents; c++) {
tomwalters@268:       for (int iCount = 0; iCount < iSizeX; iCount++) {
tomwalters@268:         pP_mod_X[iCount] += 1.0f / sqrt(2.0f * M_PI * m_fParamVar)
tomwalters@280:                             * exp((-0.5f)
tomwalters@365:                             * pow(static_cast<float>(iCount+1) - m_pMu[c], 2)
tomwalters@365:                             / m_fParamVar) * m_pA[c];
tomwalters@268:       }
tomwalters@268:     }
tomwalters@268: 
tomwalters@268:     for (int i = 0; i < iSizeX * iNComponents; ++i) {
tomwalters@268:       pP_comp[i] = 0.0f;
tomwalters@268:     }
tomwalters@268: 
tomwalters@268:     for (int i = 0; i < iNComponents; i++) {
tomwalters@268:       for (int iCount = 0; iCount < iSizeX; iCount++) {
tomwalters@268:         pP_comp[iCount + i * iSizeX] =
tomwalters@268:           1.0f / sqrt(2.0f * M_PI * m_fParamVar)
tomwalters@280:           * exp((-0.5f) * pow((static_cast<float>(iCount + 1) - m_pMu[i]), 2)
tomwalters@280:           / m_fParamVar);
tomwalters@268:         pP_comp[iCount + i * iSizeX] =
tomwalters@268:           pP_comp[iCount + i * iSizeX] * m_pA[i] / pP_mod_X[iCount];
tomwalters@268:       }
tomwalters@268:     }
tomwalters@268: 
tomwalters@268:     for (int iCount = 0; iCount < iSizeX; ++iCount) {
tomwalters@280:       float fSum = 0.0f;
tomwalters@268:       for (int i = 0; i < iNComponents; ++i) {
tomwalters@268:         pP_comp[iCount+i*iSizeX] = pow(pP_comp[iCount + i * iSizeX],
tomwalters@280:                                        m_fParamPosteriorExp);  // expansion
tomwalters@268:         fSum += pP_comp[iCount+i*iSizeX];
tomwalters@268:       }
tomwalters@268:       for (int i = 0; i < iNComponents; ++i)
tomwalters@268:         pP_comp[iCount+i*iSizeX] = pP_comp[iCount + i * iSizeX] / fSum;
tomwalters@268:         // renormalisation
tomwalters@268:     }
tomwalters@268: 
tomwalters@268:     for (int i = 0; i < iNComponents; ++i) {
tomwalters@268:       pA_old[i] = m_pA[i];
tomwalters@268:       m_pA[i] = 0.0f;
tomwalters@268:       for (int iCount = 0; iCount < iSizeX; ++iCount) {
tomwalters@268:         m_pA[i] += pP_comp[iCount + i * iSizeX] * m_pSpectralProfile[iCount];
tomwalters@268:       }
tomwalters@268:     }
tomwalters@268: 
tomwalters@268:     // finish when already converged
tomwalters@280:     float fPrdist = 0.0f;
tomwalters@268:     for (int i = 0; i < iNComponents; ++i) {
tomwalters@268:       fPrdist += pow((m_pA[i] - pA_old[i]), 2);
tomwalters@268:     }
tomwalters@268:     fPrdist /= iNComponents;
tomwalters@268: 
tomwalters@365:     if (fPrdist < priors_converged_) {
tomwalters@280:       // LOG_INFO("Converged!");
tomwalters@268:       break;
tomwalters@268:     }
tomwalters@280:     // LOG_INFO("Didn't converge!");
tomwalters@274: 
tomwalters@268: 
tomwalters@268:     // update means (positions)
tomwalters@268:     for (int i = 0 ; i < iNComponents; ++i) {
tomwalters@280:       float mu_old = m_pMu[i];
tomwalters@268:       if (m_pA[i] > 0.0f) {
tomwalters@268:         m_pMu[i] = 0.0f;
tomwalters@268:         for (int iCount = 0; iCount < iSizeX; ++iCount) {
tomwalters@268:           m_pMu[i] += m_pSpectralProfile[iCount]
tomwalters@280:                       * pP_comp[iCount + i * iSizeX]
tomwalters@280:                       * static_cast<float>(iCount + 1);
tomwalters@268:         }
tomwalters@268:         m_pMu[i] /= m_pA[i];
tomwalters@268:         if (isnan(m_pMu[i])) {
tomwalters@268:           m_pMu[i] = mu_old;
tomwalters@268:         }
tomwalters@268:       }
tomwalters@268:     }
tomwalters@280:   }  // loop over iterations
tomwalters@268: 
tomwalters@268:   // Ensure they are sorted, using a really simple bubblesort
tomwalters@268:   bool bSorted = false;
tomwalters@268:   while (!bSorted) {
tomwalters@268:     bSorted = true;
tomwalters@268:     for (int i = 0; i < iNComponents - 1; ++i) {
tomwalters@268:       if (m_pMu[i] > m_pMu[i + 1]) {
tomwalters@280:         float fTemp = m_pMu[i];
tomwalters@268:         m_pMu[i] = m_pMu[i + 1];
tomwalters@268:         m_pMu[i + 1] = fTemp;
tomwalters@268:         fTemp = m_pA[i];
tomwalters@268:         m_pA[i] = m_pA[i + 1];
tomwalters@268:         m_pA[i + 1] = fTemp;
tomwalters@268:         bSorted = false;
tomwalters@268:       }
tomwalters@268:     }
tomwalters@268:   }
tomwalters@268:   return true;
tomwalters@268: }
tomwalters@280: }  // namespace aimc
tomwalters@268: