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: 
matthiasm@0: MonoNoteHMM::MonoNoteHMM() :
matthiasm@0:     par()
matthiasm@0: {
matthiasm@0:     build();
matthiasm@0: }
matthiasm@0: 
matthiasm@0: const vector<double>
matthiasm@0: 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;
matthiasm@0:     for (size_t iCandidate = 0; iCandidate < nCandidate; ++iCandidate)
matthiasm@0:     {
matthiasm@0:         // pIsPitched = pitchProb[iCandidate].second > pIsPitched ? pitchProb[iCandidate].second : pIsPitched;
matthiasm@0:         pIsPitched += pitchProb[iCandidate].second;
matthiasm@0:     }
matthiasm@0: 
matthiasm@0:     // pIsPitched = std::pow(pIsPitched, (1-par.priorWeight)) * std::pow(par.priorPitchedProb, par.priorWeight);
matthiasm@0:     pIsPitched = pIsPitched * (1-par.priorWeight) + 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@0:         if (i % 4 != 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;
matthiasm@0:                 for (size_t iCandidate = 0; iCandidate < nCandidate; ++iCandidate)
matthiasm@0:                 {
matthiasm@0:                     double currDist = std::abs(getMidiPitch(i)-pitchProb[iCandidate].first);
matthiasm@0:                     if (currDist < minDist)
matthiasm@0:                     {
matthiasm@0:                         minDist = currDist;
matthiasm@0:                         minDistProb = pitchProb[iCandidate].second;
matthiasm@0:                         minDistCandidate = iCandidate;
matthiasm@0:                     }
matthiasm@0:                 }
matthiasm@0:                 tempProb = std::pow(minDistProb, par.yinTrust) * boost::math::pdf(pitchDistr[i], 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@0:         if (i % 4 != 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:     }
matthiasm@0: 
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@0:     //    3. inter state
matthiasm@0:     // 4-6. second-lowest pitch
matthiasm@0:     //    4. attack state
matthiasm@0:     //    ...
matthiasm@0:     
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@0:         if (iState % 4 == 2)
matthiasm@0:         {
matthiasm@0:             // silent state starts tracking
matthiasm@0:             init.push_back(1.0/(par.nS * par.nPPS));
matthiasm@0:         } else {
matthiasm@0:             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:         pitchDistr[index+3] = boost::math::normal(mu, par.sigmaYinPitchInter);
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
matthiasm@0:         from.push_back(index);
matthiasm@0:         to.push_back(index);
matthiasm@0:         transProb.push_back(par.pAttackSelftrans);
matthiasm@0: 
matthiasm@0:         from.push_back(index);
matthiasm@0:         to.push_back(index+1);
matthiasm@0:         transProb.push_back(1-par.pAttackSelftrans);
matthiasm@0: 
matthiasm@0:         // transitions from stable state
matthiasm@0:         from.push_back(index+1);
matthiasm@0:         to.push_back(index+1); // to itself
matthiasm@0:         transProb.push_back(par.pStableSelftrans);
matthiasm@0:         
matthiasm@0:         from.push_back(index+1);
matthiasm@0:         to.push_back(index+2); // to silent
matthiasm@0:         transProb.push_back(par.pStable2Silent);
matthiasm@0: 
matthiasm@0:         from.push_back(index+1);
matthiasm@0:         to.push_back(index+3); // to inter-note
matthiasm@0:         transProb.push_back(1-par.pStableSelftrans-par.pStable2Silent);
matthiasm@0: 
matthiasm@0:         // the "easy" transitions from silent state
matthiasm@0:         from.push_back(index+2);
matthiasm@0:         to.push_back(index+2);
matthiasm@0:         transProb.push_back(par.pSilentSelftrans);
matthiasm@0:         
matthiasm@0:         // the "easy" inter state transition
matthiasm@0:         from.push_back(index+3);
matthiasm@0:         to.push_back(index+3);
matthiasm@0:         transProb.push_back(par.pInterSelftrans);
matthiasm@0:         
matthiasm@0:         // the more complicated transitions from the silent and inter state
matthiasm@0:         double probSumSilent = 0;
matthiasm@0:         double probSumInter = 0;
matthiasm@0:         vector<double> tempTransProbInter;
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@0:             double semitoneDistance = std::abs(fromPitch - toPitch) * 1.0 / par.nPPS;
matthiasm@0:             
matthiasm@0:             // if (std::fmod(semitoneDistance, 1) == 0 && semitoneDistance > par.minSemitoneDistance)
matthiasm@0:             if (semitoneDistance == 0 || (semitoneDistance > par.minSemitoneDistance && semitoneDistance < par.maxJump))
matthiasm@0:             {
matthiasm@0:                 size_t toIndex = jPitch * par.nSPP; // note attack index
matthiasm@0: 
matthiasm@0:                 double tempWeightSilent = boost::math::pdf(noteDistanceDistr, semitoneDistance);
matthiasm@0:                 double tempWeightInter = semitoneDistance == 0 ? 0 : tempWeightSilent;
matthiasm@0:                 probSumSilent += tempWeightSilent;
matthiasm@0:                 probSumInter += tempWeightInter;
matthiasm@0: 
matthiasm@0:                 tempTransProbSilent.push_back(tempWeightSilent);
matthiasm@0:                 tempTransProbInter.push_back(tempWeightInter);
matthiasm@0: 
matthiasm@0:                 from.push_back(index+2);
matthiasm@0:                 to.push_back(toIndex);
matthiasm@0:                 from.push_back(index+3);
matthiasm@0:                 to.push_back(toIndex);                
matthiasm@0:             }
matthiasm@0:         }
matthiasm@0:         for (size_t i = 0; i < tempTransProbSilent.size(); ++i)
matthiasm@0:         {
matthiasm@0:             transProb.push_back((1-par.pSilentSelftrans) * tempTransProbSilent[i]/probSumSilent);
matthiasm@0:             transProb.push_back((1-par.pInterSelftrans) * tempTransProbInter[i]/probSumInter);
matthiasm@0:         }
matthiasm@0:     }
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: }