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 "MonoPitchHMM.h"
matthiasm@0: 
matthiasm@0: #include <boost/math/distributions.hpp>
matthiasm@0: 
matthiasm@0: #include <cstdio>
matthiasm@0: #include <cmath>
mail@130: #include <iostream>
matthiasm@0: 
matthiasm@0: using std::vector;
matthiasm@0: using std::pair;
matthiasm@0: 
mail@132: MonoPitchHMM::MonoPitchHMM(int fixedLag) :
Chris@141:     SparseHMM(fixedLag),
Chris@138:     m_minFreq(61.735),
Chris@138:     m_nBPS(5),
Chris@138:     m_nPitch(0),
Chris@138:     m_transitionWidth(0),
Chris@138:     m_selfTrans(0.99),
Chris@138:     m_yinTrust(.5),
Chris@138:     m_freqs(0)
matthiasm@0: {
matthiasm@0:     m_transitionWidth = 5*(m_nBPS/2) + 1;
matthiasm@102:     m_nPitch = 69 * m_nBPS;
mail@130:     m_nState = 2 * m_nPitch; // voiced and unvoiced
matthiasm@25:     m_freqs = vector<double>(2*m_nPitch);
Chris@138:     for (int iPitch = 0; iPitch < m_nPitch; ++iPitch)
matthiasm@0:     {
matthiasm@0:         m_freqs[iPitch] = m_minFreq * std::pow(2, iPitch * 1.0 / (12 * m_nBPS));
matthiasm@0:         m_freqs[iPitch+m_nPitch] = -m_freqs[iPitch];
matthiasm@0:     }
matthiasm@0:     build();
matthiasm@0: }
matthiasm@0: 
Chris@145: vector<double>
Chris@156: MonoPitchHMM::calculateObsProb(const vector<pair<double, double> > &pitchProb)
matthiasm@0: {
matthiasm@0:     vector<double> out = vector<double>(2*m_nPitch+1);
matthiasm@0:     double probYinPitched = 0;
Chris@156:     int nPair = int(pitchProb.size());
Chris@156:     
matthiasm@0:     // BIN THE PITCHES
Chris@156:     for (int iPair = 0; iPair < nPair; ++iPair)
matthiasm@0:     {
matthiasm@0:         double freq = 440. * std::pow(2, (pitchProb[iPair].first - 69)/12);
matthiasm@0:         if (freq <= m_minFreq) continue;
matthiasm@0:         double d = 0;
matthiasm@0:         double oldd = 1000;
Chris@138:         for (int iPitch = 0; iPitch < m_nPitch; ++iPitch)
matthiasm@0:         {
matthiasm@0:             d = std::abs(freq-m_freqs[iPitch]);
matthiasm@0:             if (oldd < d && iPitch > 0)
matthiasm@0:             {
matthiasm@0:                 // previous bin must have been the closest
matthiasm@0:                 out[iPitch-1] = pitchProb[iPair].second;
matthiasm@0:                 probYinPitched += out[iPitch-1];
matthiasm@0:                 break;
matthiasm@0:             }
matthiasm@0:             oldd = d;
matthiasm@0:         }
matthiasm@0:     }
matthiasm@0:     
matthiasm@0:     double probReallyPitched = m_yinTrust * probYinPitched;
matthiasm@58:     // std::cerr << probReallyPitched << " " << probYinPitched << std::endl;
matthiasm@58:     // damn, I forget what this is all about...
Chris@138:     for (int iPitch = 0; iPitch < m_nPitch; ++iPitch)
matthiasm@0:     {
matthiasm@0:         if (probYinPitched > 0) out[iPitch] *= (probReallyPitched/probYinPitched) ;
matthiasm@0:         out[iPitch+m_nPitch] = (1 - probReallyPitched) / m_nPitch;
matthiasm@0:     }
matthiasm@0:     // out[2*m_nPitch] = m_yinTrust * (1 - probYinPitched);
matthiasm@0:     return(out);
matthiasm@0: }
matthiasm@0: 
matthiasm@0: void
matthiasm@0: MonoPitchHMM::build()
matthiasm@0: {
matthiasm@0:     // INITIAL VECTOR
mail@130:     m_init = vector<double>(2*m_nPitch, 1.0 / 2*m_nPitch);
matthiasm@0:     
matthiasm@0:     // TRANSITIONS
Chris@138:     for (int iPitch = 0; iPitch < int(m_nPitch); ++iPitch)
matthiasm@0:     {
Chris@138:         int theoreticalMinNextPitch = iPitch-m_transitionWidth/2;
matthiasm@0:         int minNextPitch = iPitch>m_transitionWidth/2 ? iPitch-m_transitionWidth/2 : 0;
matthiasm@0:         int maxNextPitch = iPitch<m_nPitch-m_transitionWidth/2 ? iPitch+m_transitionWidth/2 : m_nPitch-1;
matthiasm@0:         
matthiasm@0:         // WEIGHT VECTOR
matthiasm@0:         double weightSum = 0;
matthiasm@0:         vector<double> weights;
Chris@138:         for (int i = minNextPitch; i <= maxNextPitch; ++i)
matthiasm@0:         {
matthiasm@0:             if (i <= iPitch)
matthiasm@0:             {
matthiasm@0:                 weights.push_back(i-theoreticalMinNextPitch+1);
matthiasm@0:                 // weights.push_back(i-theoreticalMinNextPitch+1+m_transitionWidth/2);
matthiasm@0:             } else {
matthiasm@0:                 weights.push_back(iPitch-theoreticalMinNextPitch+1-(i-iPitch));
matthiasm@0:                 // weights.push_back(iPitch-theoreticalMinNextPitch+1-(i-iPitch)+m_transitionWidth/2);
matthiasm@0:             }
matthiasm@0:             weightSum += weights[weights.size()-1];
matthiasm@0:         }
matthiasm@0:         
matthiasm@0:         // std::cerr << minNextPitch << "  " << maxNextPitch << std::endl;
matthiasm@0:         // TRANSITIONS TO CLOSE PITCH
Chris@138:         for (int i = minNextPitch; i <= maxNextPitch; ++i)
matthiasm@0:         {
mail@130:             m_from.push_back(iPitch);
mail@130:             m_to.push_back(i);
mail@130:             m_transProb.push_back(weights[i-minNextPitch] / weightSum * m_selfTrans);
matthiasm@0: 
mail@130:             m_from.push_back(iPitch);
mail@130:             m_to.push_back(i+m_nPitch);
mail@130:             m_transProb.push_back(weights[i-minNextPitch] / weightSum * (1-m_selfTrans));
matthiasm@0: 
mail@130:             m_from.push_back(iPitch+m_nPitch);
mail@130:             m_to.push_back(i+m_nPitch);
mail@130:             m_transProb.push_back(weights[i-minNextPitch] / weightSum * m_selfTrans);
matthiasm@0:             // transProb.push_back(weights[i-minNextPitch] / weightSum * 0.5);
matthiasm@0:             
mail@130:             m_from.push_back(iPitch+m_nPitch);
mail@130:             m_to.push_back(i);
mail@130:             m_transProb.push_back(weights[i-minNextPitch] / weightSum * (1-m_selfTrans));
matthiasm@0:             // transProb.push_back(weights[i-minNextPitch] / weightSum * 0.5);
matthiasm@0:         }
matthiasm@0: 
matthiasm@0:         // TRANSITION TO UNVOICED
matthiasm@0:         // from.push_back(iPitch+m_nPitch);
matthiasm@0:         // to.push_back(2*m_nPitch);
matthiasm@0:         // transProb.push_back(1-m_selfTrans);
matthiasm@0:         
matthiasm@0:         // TRANSITION FROM UNVOICED TO PITCH
matthiasm@25:         // from.push_back(2*m_nPitch);
matthiasm@25:         // to.push_back(iPitch+m_nPitch);
matthiasm@25:         // transProb.push_back(1.0/m_nPitch);
matthiasm@0:     }
matthiasm@0:     // UNVOICED SELFTRANSITION
matthiasm@0:     // from.push_back(2*m_nPitch);
matthiasm@0:     // to.push_back(2*m_nPitch);
matthiasm@0:     // transProb.push_back(m_selfTrans);
matthiasm@0:     
Chris@138:     // for (int i = 0; i < from.size(); ++i) {
matthiasm@0:     //     std::cerr << "P(["<< from[i] << " --> " << to[i] << "]) = " << transProb[i] << std::endl;
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);
Chris@9: }
mail@132: 
mail@132: /*
mail@132: Takes a state number and a pitch-prob vector, then finds the pitch that would
mail@132: have been closest to the pitch of the state. Easy to understand? ;)
mail@132: */
Chris@145: float
Chris@156: MonoPitchHMM::nearestFreq(int state, const vector<pair<double, double> > &pitchProb)
mail@132: {
mail@132:     float hmmFreq = m_freqs[state];
mail@132:     // std::cerr << "hmmFreq " << hmmFreq << std::endl;
mail@132:     float bestFreq = 0;
mail@132:     float leastDist = 10000;
mail@132:     if (hmmFreq > 0)
mail@132:     {
mail@132:         // This was a Yin estimate, so try to get original pitch estimate back
mail@132:         // ... a bit hacky, since we could have direclty saved the frequency
mail@132:         // that was assigned to the HMM bin in hmm.calculateObsProb -- but would
mail@132:         // have had to rethink the interface of that method.
mail@132: 
mail@132:         // std::cerr << "pitch prob size " << pitchProb.size() << std::endl;
mail@132: 
mail@132:         for (size_t iPt = 0; iPt < pitchProb.size(); ++iPt)
mail@132:         {
mail@132:             float freq = 440. * 
mail@132:                          std::pow(2, 
mail@132:                                   (pitchProb[iPt].first - 69)/12);
mail@132:             float dist = std::abs(hmmFreq-freq);
mail@132:             if (dist < leastDist)
mail@132:             {
mail@132:                 leastDist = dist;
mail@132:                 bestFreq = freq;
mail@132:             }
mail@132:         }
mail@132:     } else {
mail@132:         bestFreq = hmmFreq;
mail@132:     }
mail@132:     return bestFreq;
Chris@141: }