view SparseHMM.cpp @ 131:b877df85ad9e fixedlag

mono pitch works now with the refactored HMM implementation

/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */

/*
    pYIN - A fundamental frequency estimator for monophonic audio
    Centre for Digital Music, Queen Mary, University of London.
    
    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 2 of the
    License, or (at your option) any later version.  See the file
    COPYING included with this distribution for more information.
*/

#include "SparseHMM.h"
#include <vector>
#include <cstdio>
#include <iostream>

using std::vector;
using std::pair;

SparseHMM::SparseHMM() :
    m_nState(0),
    m_nTrans(0),
    m_init(0),
    m_from(0),
    m_to(0),
    m_transProb(0),
    m_scale(0),
    m_psi(0),
    m_delta(0),
    m_oldDelta(0)
{

}

const vector<double>
SparseHMM::calculateObsProb(const vector<pair<double, double> > data)
{
    // dummy (virtual?) implementation to be overloaded
    return(vector<double>());
}

void
SparseHMM::build()
{ }

const std::vector<int> 
SparseHMM::decodeViterbi(std::vector<vector<double> > obsProb) 
{
    size_t nFrame = obsProb.size();
    if (nFrame < 1) {
        return vector<int>();
    }

    initialise(obsProb[0]);

    // rest of forward step
    for (size_t iFrame = 1; iFrame < nFrame; ++iFrame)
    {
        process(obsProb[iFrame]);
    }

    vector<int> path = finalise();
    return(path);
}

void
SparseHMM::reset()
{
    m_scale.clear();
    m_psi.clear();
    for (size_t i = 0; i < m_delta.size(); ++i) m_delta[i] = 0;
    for (size_t i = 0; i < m_oldDelta.size(); ++i) m_oldDelta[i] = 0;
}

void
SparseHMM::initialise(vector<double> firstObs)
{
    reset();

    double deltasum = 0;

    // initialise first frame
    for (size_t iState = 0; iState < m_nState; ++iState)
    {
        m_oldDelta[iState] = m_init[iState] * firstObs[iState];
        deltasum += m_oldDelta[iState];
    }

    for (size_t iState = 0; iState < m_nState; ++iState)
    {
        m_oldDelta[iState] /= deltasum; // normalise (scale)
    }

    m_scale.push_back(1.0/deltasum);
    m_psi.push_back(vector<int>(m_nState,0));
}

int
SparseHMM::process(vector<double> newObs)
{
    vector<int> tempPsi = vector<int>(m_nState,0);

    // calculate best previous state for every current state
    size_t fromState;
    size_t toState;
    double currentTransProb;
    double currentValue;
    
    // this is the "sparse" loop
    for (size_t iTrans = 0; iTrans < m_nTrans; ++iTrans)
    {
        fromState = m_from[iTrans];
        toState = m_to[iTrans];
        currentTransProb = m_transProb[iTrans];
        
        currentValue = m_oldDelta[fromState] * currentTransProb;
        if (currentValue > m_delta[toState])
        {
            // will be multiplied by the right obs later!
            m_delta[toState] = currentValue;
            tempPsi[toState] = fromState;
        }
    }
    m_psi.push_back(tempPsi);


    double deltasum = 0;
    for (size_t jState = 0; jState < m_nState; ++jState)
    {
        m_delta[jState] *= newObs[jState];
        deltasum += m_delta[jState];
    }

    if (deltasum > 0)
    {
        for (size_t iState = 0; iState < m_nState; ++iState)
        {
            m_oldDelta[iState] = m_delta[iState] / deltasum;// normalise (scale)
            m_delta[iState] = 0;
        }
        m_scale.push_back(1.0/deltasum);
    } else
    {
        std::cerr << "WARNING: Viterbi has been fed some zero "
            "probabilities, at least they become zero "
            "in combination with the model." << std::endl;
        for (size_t iState = 0; iState < m_nState; ++iState)
        {
            m_oldDelta[iState] = 1.0/m_nState;
            m_delta[iState] = 0;
        }
        m_scale.push_back(1.0);
    }
    return 0;
}

vector<int> 
SparseHMM::finalise()
{
    // initialise backward step
    size_t nFrame = m_psi.size();

    // The final output path (current assignment arbitrary, makes sense only for 
    // Chordino, where nChord-1 is the "no chord" label)
    vector<int> path = vector<int>(nFrame, m_nState-1);

    double bestValue = 0;
    for (size_t iState = 0; iState < m_nState; ++iState)
    {
        double currentValue = m_oldDelta[iState];
        if (currentValue > bestValue)
        {
            bestValue = currentValue;            
            path[nFrame-1] = iState;
        }
    }

    // Rest of backward step
    for (int iFrame = nFrame-2; iFrame != -1; --iFrame)
    {
        path[iFrame] = m_psi[iFrame+1][path[iFrame+1]];
    }
        
    return path;
}