tomwalters@0: // Copyright 2008-2010, Thomas Walters
tomwalters@0: //
tomwalters@0: // AIM-C: A C++ implementation of the Auditory Image Model
tomwalters@0: // http://www.acousticscale.org/AIMC
tomwalters@0: //
tomwalters@0: // This program is free software: you can redistribute it and/or modify
tomwalters@0: // it under the terms of the GNU General Public License as published by
tomwalters@0: // the Free Software Foundation, either version 3 of the License, or
tomwalters@0: // (at your option) any later version.
tomwalters@0: //
tomwalters@0: // This program is distributed in the hope that it will be useful,
tomwalters@0: // but WITHOUT ANY WARRANTY; without even the implied warranty of
tomwalters@0: // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
tomwalters@0: // GNU General Public License for more details.
tomwalters@0: //
tomwalters@0: // You should have received a copy of the GNU General Public License
tomwalters@0: // along with this program.  If not, see <http://www.gnu.org/licenses/>.
tomwalters@0: 
tomwalters@0: /*! \file
tomwalters@0:  *  \brief Gaussian features - based on MATLAB code by Christian Feldbauer
tomwalters@0:  */
tomwalters@0: 
tomwalters@0: /*!
tomwalters@0:  * \author Tom Walters <tcw24@cam.ac.uk>
tomwalters@0:  * \date created 2008/06/23
tomwalters@0:  * \version \$Id: ModuleGaussians.cc 2 2010-02-02 12:59:50Z tcw $
tomwalters@0:  */
tomwalters@0: 
tomwalters@0: #include <math.h>
tomwalters@0: 
tomwalters@0: #include "Modules/Features/ModuleGaussians.h"
tomwalters@0: #include "Support/Common.h"
tomwalters@0: 
tomwalters@0: namespace aimc {
tomwalters@0: ModuleGaussians::ModuleGaussians(Parameters *pParam)
tomwalters@0: : Module(pParam) {
tomwalters@0:   // Set module metadata
tomwalters@0:   module_description_ = "Gaussian Fitting to SSI profile";
tomwalters@0:   module_identifier_ = "gaussians"; // unique identifier for the module
tomwalters@0:   module_type_ = "features";
tomwalters@0:   module_version_ = "$Id: ModuleGaussians.cc 2 2010-02-02 12:59:50Z tcw $";
tomwalters@0: 
tomwalters@0:   parameters_->SetDefault("features.gaussians.ncomp", "4");
tomwalters@0:   m_iParamNComp = parameters_->GetInt("features.gaussians.ncomp");
tomwalters@0: 
tomwalters@0:   parameters_->SetDefault("features.gaussians.var", "115.0");
tomwalters@0:   m_fParamVar = parameters_->GetFloat("features.gaussians.var");
tomwalters@0: 
tomwalters@0:   parameters_->SetDefault("features.gaussians.posterior_exp", "6.0");
tomwalters@0:   m_fParamPosteriorExp =
tomwalters@0:     parameters_->GetFloat("features.gaussians.posterior_exp");
tomwalters@0: 
tomwalters@0:   parameters_->SetDefault("features.gaussians.maxit", "250");
tomwalters@0:   m_iParamMaxIt = parameters_->GetInt("features.gaussians.maxit");
tomwalters@0: 
tomwalters@0:   parameters_->SetDefault("features.gaussians.priors_converged", "1e-7");
tomwalters@0:   m_fParamPriorsConverged =
tomwalters@0:     parameters_->GetInt("features.gaussians.priors_converged");
tomwalters@0: }
tomwalters@0: 
tomwalters@0: ModuleGaussians::~ModuleGaussians() {
tomwalters@0: }
tomwalters@0: 
tomwalters@0: bool ModuleGaussians::InitializeInternal(const SignalBank &input) {
tomwalters@0:   m_pA.resize(m_iParamNComp, 0.0f);
tomwalters@0:   m_pMu.resize(m_iParamNComp, 0.0f);
tomwalters@0: 
tomwalters@0:   // Assuming the number of channels is greater than twice the number of
tomwalters@0:   // Gaussian components, this is ok
tomwalters@0:   if (input.channel_count() >= 2 * m_iParamNComp) {
tomwalters@0:     output_.Initialize(1, m_iParamNComp, input.sample_rate());
tomwalters@0:   } else {
tomwalters@0:     LOG_ERROR(_T("Too few channels in filterbank to produce sensible "
tomwalters@0:                  "Gaussian features. Either increase the number of filterbank"
tomwalters@0:                  " channels, or decrease the number of Gaussian components"));
tomwalters@0:     return false;
tomwalters@0:   }
tomwalters@0: 
tomwalters@0:   m_iNumChannels = input.channel_count();
tomwalters@0:   m_pSpectralProfile.resize(m_iNumChannels, 0.0f);
tomwalters@0: 
tomwalters@0:   return true;
tomwalters@0: }
tomwalters@0: 
tomwalters@0: void ModuleGaussians::Reset() {
tomwalters@0:   m_pSpectralProfile.clear();
tomwalters@0:   m_pSpectralProfile.resize(m_iNumChannels, 0.0f);
tomwalters@0: }
tomwalters@0: 
tomwalters@0: void ModuleGaussians::Process(const SignalBank &input) {
tomwalters@0:   int iAudCh = 0;
tomwalters@0: 
tomwalters@0:   // Calculate spectral profile
tomwalters@0:   for (int iChannel = 0;
tomwalters@0:        iChannel < input.channel_count();
tomwalters@0:        ++iChannel) {
tomwalters@0:     m_pSpectralProfile[iChannel] = 0.0f;
tomwalters@0:     for (int iSample = 0;
tomwalters@0:          iSample < input.buffer_length();
tomwalters@0:          ++iSample) {
tomwalters@0:       m_pSpectralProfile[iChannel] += input[iChannel][iSample];
tomwalters@0:     }
tomwalters@0:   }
tomwalters@0: 
tomwalters@0:   for (int iChannel = 0;
tomwalters@0:        iChannel < input.channel_count();
tomwalters@0:        ++iChannel) {
tomwalters@0:     m_pSpectralProfile[iChannel] = pow(m_pSpectralProfile[iChannel], 0.8);
tomwalters@0:   }
tomwalters@0: 
tomwalters@0:   float spectral_profile_sum = 0.0f;
tomwalters@0:   for (int i = 0; i < input.channel_count(); ++i) {
tomwalters@0:     spectral_profile_sum += m_pSpectralProfile[i];
tomwalters@0:   }
tomwalters@0: 
tomwalters@0:   RubberGMMCore(2, true);
tomwalters@0: 
tomwalters@0:   float fMean1 = m_pMu[0];
tomwalters@0:   float fMean2 = m_pMu[1];
tomwalters@0: 
tomwalters@0:   float fA1 = 0.05 * m_pA[0];
tomwalters@0:   float fA2 = 1.0 - 0.25 * m_pA[1];
tomwalters@0: 
tomwalters@0:   float fGradient = (fMean2 - fMean1) / (fA2 - fA1);
tomwalters@0:   float fIntercept = fMean2 - fGradient * fA2;
tomwalters@0: 
tomwalters@0:   for (int i = 0; i < m_iParamNComp; ++i) {
tomwalters@0:     m_pMu[i] = ((float)i / (float)m_iParamNComp - 1.0f)
tomwalters@0:                   * -fGradient + fIntercept;
tomwalters@0:   }
tomwalters@0: 
tomwalters@0:   for (int i = 0; i < m_iParamNComp; ++i) {
tomwalters@0:     m_pA[i] = 1.0f / (float)m_iParamNComp;
tomwalters@0:   }
tomwalters@0: 
tomwalters@0:   RubberGMMCore(m_iParamNComp, false);
tomwalters@0: 
tomwalters@0:   for (int i = 0; i < m_iParamNComp - 1; ++i) {
tomwalters@0:     if (!isnan(m_pA[i])) {
tomwalters@0:       output_.set_sample(i, 0, m_pA[i]);
tomwalters@0:     } else {
tomwalters@0:       output_.set_sample(i, 0, 0.0f);
tomwalters@0:     }
tomwalters@0:   }
tomwalters@0:   /*for (int i = m_iParamNComp; i < m_iParamNComp * 2; ++i) {
tomwalters@0:     m_pOutputData->getSignal(i)->setSample(iAudCh, 0, m_pMu[i-m_iParamNComp]);
tomwalters@0:   }*/
tomwalters@0:   double logsum = log(spectral_profile_sum);
tomwalters@0:   if (!isinf(logsum)) {
tomwalters@0:     output_.set_sample(m_iParamNComp - 1, 0, logsum);
tomwalters@0:   } else {
tomwalters@0:     output_.set_sample(m_iParamNComp - 1, 0, -1000.0);
tomwalters@0:   }
tomwalters@0:   PushOutput();
tomwalters@0: }
tomwalters@0: 
tomwalters@0: bool ModuleGaussians::RubberGMMCore(int iNComponents, bool bDoInit) {
tomwalters@0:   int iSizeX = m_iNumChannels;
tomwalters@0: 
tomwalters@0:   // Normalise the spectral profile
tomwalters@0:   float fSpectralProfileTotal = 0.0f;
tomwalters@0:   for (int iCount = 0; iCount < iSizeX; iCount++) {
tomwalters@0:     fSpectralProfileTotal += m_pSpectralProfile[iCount];
tomwalters@0:   }
tomwalters@0:   for (int iCount = 0; iCount < iSizeX; iCount++) {
tomwalters@0:     m_pSpectralProfile[iCount] /= fSpectralProfileTotal;
tomwalters@0:   }
tomwalters@0: 
tomwalters@0:   if (bDoInit) {
tomwalters@0:     // Uniformly spaced components
tomwalters@0:     float dd = (iSizeX - 1.0f) / iNComponents;
tomwalters@0:     for (int i = 0; i < iNComponents; i++) {
tomwalters@0:       m_pMu[i] = dd / 2.0f + (i * dd);
tomwalters@0:       m_pA[i] = 1.0f / iNComponents;
tomwalters@0:     }
tomwalters@0:   }
tomwalters@0: 
tomwalters@0:   vector<float> pA_old;
tomwalters@0:   pA_old.resize(iNComponents);
tomwalters@0:   vector<float> pP_mod_X;
tomwalters@0:   pP_mod_X.resize(iSizeX);
tomwalters@0:   vector<float> pP_comp;
tomwalters@0:   pP_comp.resize(iSizeX * iNComponents);
tomwalters@0: 
tomwalters@0:   for (int iIteration = 0; iIteration < m_iParamMaxIt; iIteration++) {
tomwalters@0:     // (re)calculate posteriors (component probability given observation)
tomwalters@0:     // denominator: the model density at all observation points X
tomwalters@0:     for (int i = 0; i < iSizeX; ++i) {
tomwalters@0:       pP_mod_X[i] = 0.0f;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     for (int i = 0; i < iNComponents; i++) {
tomwalters@0:       for (int iCount = 0; iCount < iSizeX; iCount++) {
tomwalters@0:         pP_mod_X[iCount] += 1.0f / sqrt(2.0f * M_PI * m_fParamVar)
tomwalters@0:                             * exp((-0.5f) * pow(((float)iCount-m_pMu[i]), 2)
tomwalters@0:                                   / m_fParamVar) * m_pA[i];
tomwalters@0:       }
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     for (int i = 0; i < iSizeX * iNComponents; ++i) {
tomwalters@0:       pP_comp[i] = 0.0f;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     for (int i = 0; i < iNComponents; i++) {
tomwalters@0:       for (int iCount = 0; iCount < iSizeX; iCount++) {
tomwalters@0:         pP_comp[iCount + i * iSizeX] =
tomwalters@0:           1.0f / sqrt(2.0f * M_PI * m_fParamVar)
tomwalters@0:           * exp((-0.5f) * pow(((float)iCount - m_pMu[i]), 2) / m_fParamVar);
tomwalters@0:         pP_comp[iCount + i * iSizeX] =
tomwalters@0:           pP_comp[iCount + i * iSizeX] * m_pA[i] / pP_mod_X[iCount];
tomwalters@0:       }
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     for (int iCount = 0; iCount < iSizeX; ++iCount) {
tomwalters@0:       float fSum = 0.0f;
tomwalters@0:       for (int i = 0; i < iNComponents; ++i) {
tomwalters@0:         pP_comp[iCount+i*iSizeX] = pow(pP_comp[iCount + i * iSizeX],
tomwalters@0:                                        m_fParamPosteriorExp); // expansion
tomwalters@0:         fSum += pP_comp[iCount+i*iSizeX];
tomwalters@0:       }
tomwalters@0:       for (int i = 0; i < iNComponents; ++i)
tomwalters@0:         pP_comp[iCount+i*iSizeX] = pP_comp[iCount + i * iSizeX] / fSum;
tomwalters@0:         // renormalisation
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     for (int i = 0; i < iNComponents; ++i) {
tomwalters@0:       pA_old[i] = m_pA[i];
tomwalters@0:       m_pA[i] = 0.0f;
tomwalters@0:       for (int iCount = 0; iCount < iSizeX; ++iCount) {
tomwalters@0:         m_pA[i] += pP_comp[iCount + i * iSizeX] * m_pSpectralProfile[iCount];
tomwalters@0:       }
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     // finish when already converged
tomwalters@0:     float fPrdist = 0.0f;
tomwalters@0:     for (int i = 0; i < iNComponents; ++i) {
tomwalters@0:       fPrdist += pow((m_pA[i] - pA_old[i]), 2);
tomwalters@0:     }
tomwalters@0:     fPrdist /= iNComponents;
tomwalters@0: 
tomwalters@0:     if (fPrdist < m_fParamPriorsConverged) {
tomwalters@0:       LOG_INFO("Converged!");
tomwalters@0:       break;
tomwalters@0:     }
tomwalters@0: 
tomwalters@0:     // update means (positions)
tomwalters@0:     for (int i = 0 ; i < iNComponents; ++i) {
tomwalters@0:       float mu_old = m_pMu[i];
tomwalters@0:       if (m_pA[i] > 0.0f) {
tomwalters@0:         m_pMu[i] = 0.0f;
tomwalters@0:         for (int iCount = 0; iCount < iSizeX; ++iCount) {
tomwalters@0:           m_pMu[i] += m_pSpectralProfile[iCount]
tomwalters@0:                       * pP_comp[iCount + i * iSizeX] * (float)iCount;
tomwalters@0:         }
tomwalters@0:         m_pMu[i] /= m_pA[i];
tomwalters@0:         if (isnan(m_pMu[i])) {
tomwalters@0:           m_pMu[i] = mu_old;
tomwalters@0:         }
tomwalters@0:       }
tomwalters@0:     }
tomwalters@0:   } // loop over iterations
tomwalters@0: 
tomwalters@0:   // Ensure they are sorted, using a really simple bubblesort
tomwalters@0:   bool bSorted = false;
tomwalters@0:   while (!bSorted) {
tomwalters@0:     bSorted = true;
tomwalters@0:     for (int i = 0; i < iNComponents - 1; ++i) {
tomwalters@0:       if (m_pMu[i] > m_pMu[i + 1]) {
tomwalters@0:         float fTemp = m_pMu[i];
tomwalters@0:         m_pMu[i] = m_pMu[i + 1];
tomwalters@0:         m_pMu[i + 1] = fTemp;
tomwalters@0:         fTemp = m_pA[i];
tomwalters@0:         m_pA[i] = m_pA[i + 1];
tomwalters@0:         m_pA[i + 1] = fTemp;
tomwalters@0:         bSorted = false;
tomwalters@0:       }
tomwalters@0:     }
tomwalters@0:   }
tomwalters@0:   return true;
tomwalters@0: }
tomwalters@0: }  //namespace aimc
tomwalters@0: