Chris@9: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@9: 
Chris@9: /*
Chris@9:     pYIN - A fundamental frequency estimator for monophonic audio
Chris@9:     Centre for Digital Music, Queen Mary, University of London.
Chris@9:     
Chris@9:     This program is free software; you can redistribute it and/or
Chris@9:     modify it under the terms of the GNU General Public License as
Chris@9:     published by the Free Software Foundation; either version 2 of the
Chris@9:     License, or (at your option) any later version.  See the file
Chris@9:     COPYING included with this distribution for more information.
Chris@9: */
Chris@9: 
matthiasm@0: #include "SparseHMM.h"
matthiasm@0: #include <vector>
matthiasm@0: #include <cstdio>
matthiasm@0: #include <iostream>
matthiasm@0: 
matthiasm@0: using std::vector;
matthiasm@0: using std::pair;
matthiasm@0: 
mail@132: SparseHMM::SparseHMM(int fixedLag) :
mail@132:     m_fixedLag(fixedLag),
mail@130:     m_nState(0),
mail@130:     m_nTrans(0),
mail@130:     m_init(0),
mail@130:     m_from(0),
mail@130:     m_to(0),
mail@130:     m_transProb(0),
mail@130:     m_scale(0),
mail@130:     m_psi(0),
mail@130:     m_delta(0),
mail@130:     m_oldDelta(0)
mail@130: {
mail@130: 
mail@130: }
mail@130: 
mail@130: void
mail@130: SparseHMM::build()
mail@130: { }
mail@130: 
Chris@145: std::vector<int> 
mail@130: SparseHMM::decodeViterbi(std::vector<vector<double> > obsProb) 
matthiasm@0: {
Chris@142:     int nFrame = obsProb.size();
mail@130:     if (nFrame < 1) {
Chris@38:         return vector<int>();
Chris@38:     }
Chris@38: 
mail@130:     initialise(obsProb[0]);
Chris@38: 
Chris@38:     // rest of forward step
Chris@142:     for (int iFrame = 1; iFrame < nFrame; ++iFrame)
Chris@38:     {
mail@130:         process(obsProb[iFrame]);
mail@130:     }
Chris@38: 
mail@132:     vector<int> path = track();
mail@130:     return(path);
mail@130: }
mail@130: 
mail@131: void
mail@131: SparseHMM::reset()
mail@131: {
mail@131:     m_scale.clear();
mail@131:     m_psi.clear();
Chris@142:     for (int i = 0; i < int(m_delta.size()); ++i) m_delta[i] = 0;
Chris@142:     for (int i = 0; i < int(m_oldDelta.size()); ++i) m_oldDelta[i] = 0;
mail@131: }
mail@130: 
mail@130: void
mail@130: SparseHMM::initialise(vector<double> firstObs)
mail@130: {
mail@131:     reset();
matthiasm@0: 
matthiasm@0:     double deltasum = 0;
matthiasm@0: 
matthiasm@0:     // initialise first frame
Chris@142:     for (int iState = 0; iState < m_nState; ++iState)
matthiasm@0:     {
mail@130:         m_oldDelta[iState] = m_init[iState] * firstObs[iState];
mail@130:         deltasum += m_oldDelta[iState];
matthiasm@0:     }
matthiasm@0: 
Chris@142:     for (int iState = 0; iState < m_nState; ++iState)
matthiasm@0:     {
mail@130:         m_oldDelta[iState] /= deltasum; // normalise (scale)
matthiasm@0:     }
matthiasm@0: 
mail@130:     m_scale.push_back(1.0/deltasum);
mail@130:     m_psi.push_back(vector<int>(m_nState,0));
mail@130: }
matthiasm@0: 
mail@130: int
mail@130: SparseHMM::process(vector<double> newObs)
mail@130: {
mail@130:     vector<int> tempPsi = vector<int>(m_nState,0);
mail@130: 
mail@130:     // calculate best previous state for every current state
Chris@142:     int fromState;
Chris@142:     int toState;
mail@130:     double currentTransProb;
mail@130:     double currentValue;
mail@130:     
mail@130:     // this is the "sparse" loop
Chris@142:     for (int iTrans = 0; iTrans < m_nTrans; ++iTrans)
matthiasm@0:     {
mail@130:         fromState = m_from[iTrans];
mail@130:         toState = m_to[iTrans];
mail@130:         currentTransProb = m_transProb[iTrans];
matthiasm@0:         
mail@130:         currentValue = m_oldDelta[fromState] * currentTransProb;
mail@130:         if (currentValue > m_delta[toState])
matthiasm@0:         {
mail@130:             // will be multiplied by the right obs later!
mail@130:             m_delta[toState] = currentValue;
mail@130:             tempPsi[toState] = fromState;
matthiasm@0:         }
matthiasm@0:     }
mail@130:     m_psi.push_back(tempPsi);
matthiasm@0: 
mail@130:     double deltasum = 0;
Chris@142:     for (int jState = 0; jState < m_nState; ++jState)
mail@130:     {
mail@130:         m_delta[jState] *= newObs[jState];
mail@130:         deltasum += m_delta[jState];
mail@130:     }
mail@130: 
mail@130:     if (deltasum > 0)
mail@130:     {
Chris@142:         for (int iState = 0; iState < m_nState; ++iState)
mail@130:         {
mail@130:             m_oldDelta[iState] = m_delta[iState] / deltasum;// normalise (scale)
mail@130:             m_delta[iState] = 0;
mail@130:         }
mail@130:         m_scale.push_back(1.0/deltasum);
mail@130:     } else
mail@130:     {
Chris@149: //        std::cerr << "WARNING: Viterbi has been fed some zero "
Chris@149: //            "probabilities, at least they become zero "
Chris@149: //            "in combination with the model." << std::endl;
Chris@142:         for (int iState = 0; iState < m_nState; ++iState)
mail@130:         {
mail@130:             m_oldDelta[iState] = 1.0/m_nState;
mail@130:             m_delta[iState] = 0;
mail@130:         }
mail@130:         m_scale.push_back(1.0);
mail@130:     }
mail@132: 
Chris@142:     if (m_fixedLag > 0 && int(m_psi.size()) > m_fixedLag)
mail@132:     {
mail@132:         m_psi.pop_front();
mail@132:         m_scale.pop_front();
mail@132:     }
mail@132: 
mail@132:     // std::cerr << m_fixedLag << " " << m_psi.size() << std::endl;
mail@132: 
mail@130:     return 0;
mail@130: }
mail@130: 
Chris@145: vector<int>
mail@132: SparseHMM::track()
mail@130: {
matthiasm@0:     // initialise backward step
Chris@142:     int nFrame = m_psi.size();
mail@130: 
mail@130:     // The final output path (current assignment arbitrary, makes sense only for 
mail@130:     // Chordino, where nChord-1 is the "no chord" label)
mail@130:     vector<int> path = vector<int>(nFrame, m_nState-1);
mail@130: 
matthiasm@0:     double bestValue = 0;
Chris@142:     for (int iState = 0; iState < m_nState; ++iState)
matthiasm@0:     {
mail@130:         double currentValue = m_oldDelta[iState];
matthiasm@0:         if (currentValue > bestValue)
matthiasm@0:         {
matthiasm@0:             bestValue = currentValue;            
matthiasm@0:             path[nFrame-1] = iState;
matthiasm@0:         }
matthiasm@0:     }
matthiasm@0: 
mail@130:     // Rest of backward step
matthiasm@0:     for (int iFrame = nFrame-2; iFrame != -1; --iFrame)
matthiasm@0:     {
mail@130:         path[iFrame] = m_psi[iFrame+1][path[iFrame+1]];
matthiasm@0:     }
mail@130:         
matthiasm@0:     return path;
Chris@142: }