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 "MonoNoteHMM.h"
matthiasm@0: 
matthiasm@0: #include <boost/math/distributions.hpp>
matthiasm@0: 
matthiasm@0: #include <cstdio>
matthiasm@0: #include <cmath>
matthiasm@0: 
matthiasm@0: using std::vector;
matthiasm@0: using std::pair;
matthiasm@0: 
mail@132: MonoNoteHMM::MonoNoteHMM(int fixedLag) :
mail@132:     SparseHMM(fixedLag),
matthiasm@0:     par()
matthiasm@0: {
matthiasm@0:     build();
matthiasm@0: }
matthiasm@0: 
Chris@145: vector<double>
Chris@156: MonoNoteHMM::calculateObsProb(const vector<pair<double, double> > &pitchProb)
matthiasm@0: {
matthiasm@0:     // pitchProb is a list of pairs (pitches and their probabilities)
matthiasm@0:     
matthiasm@0:     size_t nCandidate = pitchProb.size();
matthiasm@0:     
matthiasm@0:     // what is the probability of pitched
matthiasm@0:     double pIsPitched = 0;
mail@133:     for (size_t iCand = 0; iCand < nCandidate; ++iCand)
matthiasm@0:     {
mail@133:         pIsPitched += pitchProb[iCand].second;
matthiasm@0:     }
matthiasm@0: 
mail@133:     pIsPitched = pIsPitched * (1-par.priorWeight) + 
mail@133:                      par.priorPitchedProb * par.priorWeight;
matthiasm@0: 
matthiasm@0:     vector<double> out = vector<double>(par.n);
matthiasm@0:     double tempProbSum = 0;
matthiasm@0:     for (size_t i = 0; i < par.n; ++i)
matthiasm@0:     {
matthiasm@103:         if (i % par.nSPP != 2)
matthiasm@0:         {
matthiasm@0:             // std::cerr << getMidiPitch(i) << std::endl;
matthiasm@0:             double tempProb = 0;
matthiasm@0:             if (nCandidate > 0)
matthiasm@0:             {
matthiasm@0:                 double minDist = 10000.0;
matthiasm@0:                 double minDistProb = 0;
matthiasm@0:                 size_t minDistCandidate = 0;
mail@133:                 for (size_t iCand = 0; iCand < nCandidate; ++iCand)
matthiasm@0:                 {
mail@133:                     double currDist = std::abs(getMidiPitch(i)-
mail@133:                                                pitchProb[iCand].first);
matthiasm@0:                     if (currDist < minDist)
matthiasm@0:                     {
matthiasm@0:                         minDist = currDist;
mail@133:                         minDistProb = pitchProb[iCand].second;
mail@133:                         minDistCandidate = iCand;
matthiasm@0:                     }
matthiasm@0:                 }
matthiasm@101:                 tempProb = std::pow(minDistProb, par.yinTrust) * 
matthiasm@101:                            boost::math::pdf(pitchDistr[i], 
matthiasm@101:                                             pitchProb[minDistCandidate].first);
matthiasm@0:             } else {
matthiasm@0:                 tempProb = 1;
matthiasm@0:             }
matthiasm@0:             tempProbSum += tempProb;
matthiasm@0:             out[i] = tempProb;
matthiasm@0:         }
matthiasm@0:     }
matthiasm@0:     
matthiasm@0:     for (size_t i = 0; i < par.n; ++i)
matthiasm@0:     {
matthiasm@103:         if (i % par.nSPP != 2)
matthiasm@0:         {
matthiasm@0:             if (tempProbSum > 0) 
matthiasm@0:             {
matthiasm@0:                 out[i] = out[i] / tempProbSum * pIsPitched;
matthiasm@0:             }
matthiasm@0:         } else {
matthiasm@0:             out[i] = (1-pIsPitched) / (par.nPPS * par.nS);
matthiasm@0:         }
matthiasm@0:     }
Chris@146:     
matthiasm@0:     return(out);
matthiasm@0: }
matthiasm@0: 
matthiasm@0: void
matthiasm@0: MonoNoteHMM::build()
matthiasm@0: {
matthiasm@0:     // the states are organised as follows:
matthiasm@0:     // 0-2. lowest pitch
matthiasm@0:     //    0. attack state
matthiasm@0:     //    1. stable state
matthiasm@0:     //    2. silent state
matthiasm@103:     // 3-5. second-lowest pitch
matthiasm@103:     //    3. attack state
matthiasm@0:     //    ...
matthiasm@0:     
mail@130:     m_nState = par.n;
mail@130: 
matthiasm@0:     // observation distributions
matthiasm@0:     for (size_t iState = 0; iState < par.n; ++iState)
matthiasm@0:     {
matthiasm@0:         pitchDistr.push_back(boost::math::normal(0,1));
matthiasm@103:         if (iState % par.nSPP == 2)
matthiasm@0:         {
matthiasm@0:             // silent state starts tracking
mail@130:             m_init.push_back(1.0/(par.nS * par.nPPS));
matthiasm@0:         } else {
mail@130:             m_init.push_back(0.0);            
matthiasm@0:         }
matthiasm@0:     }
matthiasm@0: 
matthiasm@0:     for (size_t iPitch = 0; iPitch < (par.nS * par.nPPS); ++iPitch)
matthiasm@0:     {
matthiasm@0:         size_t index = iPitch * par.nSPP;
matthiasm@0:         double mu = par.minPitch + iPitch * 1.0/par.nPPS;
matthiasm@0:         pitchDistr[index] = boost::math::normal(mu, par.sigmaYinPitchAttack);
matthiasm@0:         pitchDistr[index+1] = boost::math::normal(mu, par.sigmaYinPitchStable);
matthiasm@0:         pitchDistr[index+2] = boost::math::normal(mu, 1.0); // dummy
matthiasm@0:     }
matthiasm@0:     
matthiasm@0:     boost::math::normal noteDistanceDistr(0, par.sigma2Note);
matthiasm@0: 
matthiasm@0:     for (size_t iPitch = 0; iPitch < (par.nS * par.nPPS); ++iPitch)
matthiasm@0:     {
matthiasm@0:         // loop through all notes and set sparse transition probabilities
matthiasm@0:         size_t index = iPitch * par.nSPP;
matthiasm@0: 
matthiasm@0:         // transitions from attack state
mail@130:         m_from.push_back(index);
mail@130:         m_to.push_back(index);
mail@130:         m_transProb.push_back(par.pAttackSelftrans);
matthiasm@0: 
mail@130:         m_from.push_back(index);
mail@130:         m_to.push_back(index+1);
mail@130:         m_transProb.push_back(1-par.pAttackSelftrans);
matthiasm@0: 
matthiasm@0:         // transitions from stable state
mail@130:         m_from.push_back(index+1);
mail@130:         m_to.push_back(index+1); // to itself
mail@130:         m_transProb.push_back(par.pStableSelftrans);
matthiasm@0:         
mail@130:         m_from.push_back(index+1);
mail@130:         m_to.push_back(index+2); // to silent
mail@130:         m_transProb.push_back(par.pStable2Silent);
matthiasm@0: 
matthiasm@0:         // the "easy" transitions from silent state
mail@130:         m_from.push_back(index+2);
mail@130:         m_to.push_back(index+2);
mail@130:         m_transProb.push_back(par.pSilentSelftrans);
matthiasm@0:         
matthiasm@0:         
matthiasm@106:         // the more complicated transitions from the silent
matthiasm@0:         double probSumSilent = 0;
matthiasm@106: 
matthiasm@0:         vector<double> tempTransProbSilent;
matthiasm@0:         for (size_t jPitch = 0; jPitch < (par.nS * par.nPPS); ++jPitch)
matthiasm@0:         {
matthiasm@0:             int fromPitch = iPitch;
matthiasm@0:             int toPitch = jPitch;
matthiasm@101:             double semitoneDistance = 
matthiasm@101:                 std::abs(fromPitch - toPitch) * 1.0 / par.nPPS;
matthiasm@0:             
mail@133: 
matthiasm@101:             if (semitoneDistance == 0 || 
matthiasm@101:                 (semitoneDistance > par.minSemitoneDistance 
matthiasm@101:                  && semitoneDistance < par.maxJump))
matthiasm@0:             {
matthiasm@0:                 size_t toIndex = jPitch * par.nSPP; // note attack index
matthiasm@0: 
matthiasm@101:                 double tempWeightSilent = boost::math::pdf(noteDistanceDistr, 
matthiasm@101:                                                            semitoneDistance);
matthiasm@0:                 probSumSilent += tempWeightSilent;
matthiasm@0: 
matthiasm@0:                 tempTransProbSilent.push_back(tempWeightSilent);
matthiasm@0: 
mail@130:                 m_from.push_back(index+2);
mail@130:                 m_to.push_back(toIndex);
matthiasm@0:             }
matthiasm@0:         }
matthiasm@0:         for (size_t i = 0; i < tempTransProbSilent.size(); ++i)
matthiasm@0:         {
mail@133:             m_transProb.push_back((1-par.pSilentSelftrans) * 
mail@133:                                   tempTransProbSilent[i]/probSumSilent);
matthiasm@0:         }
matthiasm@0:     }
mail@130:     m_nTrans = m_transProb.size();
mail@130:     m_delta = vector<double>(m_nState);
mail@130:     m_oldDelta = vector<double>(m_nState);
matthiasm@0: }
matthiasm@0: 
matthiasm@0: double
matthiasm@0: MonoNoteHMM::getMidiPitch(size_t index)
matthiasm@0: {
matthiasm@0:     return pitchDistr[index].mean();
matthiasm@0: }
matthiasm@0: 
matthiasm@0: double
matthiasm@0: MonoNoteHMM::getFrequency(size_t index)
matthiasm@0: {
matthiasm@0:     return 440 * pow(2, (pitchDistr[index].mean()-69)/12);
matthiasm@0: }