libs/vamp-pyin/MonoPitchHMM.cpp

   1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
   2
   3 /*
   4     pYIN - A fundamental frequency estimator for monophonic audio
   5     Centre for Digital Music, Queen Mary, University of London.
   6
   7     This program is free software; you can redistribute it and/or
   8     modify it under the terms of the GNU General Public License as
   9     published by the Free Software Foundation; either version 2 of the
  10     License, or (at your option) any later version.  See the file
  11     COPYING included with this distribution for more information.
  12 */
  13
  14 #include "MonoPitchHMM.h"
  15
  16 #include <boost/math/distributions.hpp>
  17
  18 #include <cstdio>
  19 #include <cmath>
  20
  21 using std::vector;
  22 using std::pair;
  23
  24 MonoPitchHMM::MonoPitchHMM() :
  25 m_minFreq(61.735),
  26 m_nBPS(5),
  27 m_nPitch(0),
  28 m_transitionWidth(0),
  29 m_selfTrans(0.99),
  30 m_yinTrust(.5),
  31 m_freqs(0)
  32 {
  33     m_transitionWidth = 5*(m_nBPS/2) + 1;
  34     m_nPitch = 69 * m_nBPS;
  35     m_freqs = vector<double>(2*m_nPitch);
  36     for (size_t iPitch = 0; iPitch < m_nPitch; ++iPitch)
  37     {
  38         m_freqs[iPitch] = m_minFreq * std::pow(2, iPitch * 1.0 / (12 * m_nBPS));
  39         m_freqs[iPitch+m_nPitch] = -m_freqs[iPitch];
  40     }
  41     build();
  42 }
  43
  44 const vector<double>
  45 MonoPitchHMM::calculateObsProb(const vector<pair<double, double> > pitchProb)
  46 {
  47     vector<double> out = vector<double>(2*m_nPitch+1);
  48     double probYinPitched = 0;
  49     // BIN THE PITCHES
  50     for (size_t iPair = 0; iPair < pitchProb.size(); ++iPair)
  51     {
  52         double freq = 440. * std::pow(2, (pitchProb[iPair].first - 69)/12);
  53         if (freq <= m_minFreq) continue;
  54         double d = 0;
  55         double oldd = 1000;
  56         for (size_t iPitch = 0; iPitch < m_nPitch; ++iPitch)
  57         {
  58             d = std::abs(freq-m_freqs[iPitch]);
  59             if (oldd < d && iPitch > 0)
  60             {
  61                 // previous bin must have been the closest
  62                 out[iPitch-1] = pitchProb[iPair].second;
  63                 probYinPitched += out[iPitch-1];
  64                 break;
  65             }
  66             oldd = d;
  67         }
  68     }
  69
  70     double probReallyPitched = m_yinTrust * probYinPitched;
  71     // std::cerr << probReallyPitched << " " << probYinPitched << std::endl;
  72     // damn, I forget what this is all about...
  73     for (size_t iPitch = 0; iPitch < m_nPitch; ++iPitch)
  74     {
  75         if (probYinPitched > 0) out[iPitch] *= (probReallyPitched/probYinPitched) ;
  76         out[iPitch+m_nPitch] = (1 - probReallyPitched) / m_nPitch;
  77     }
  78     // out[2*m_nPitch] = m_yinTrust * (1 - probYinPitched);
  79     return(out);
  80 }
  81
  82 void
  83 MonoPitchHMM::build()
  84 {
  85     // INITIAL VECTOR
  86     init = vector<double>(2*m_nPitch, 1.0 / 2*m_nPitch);
  87
  88     // TRANSITIONS
  89     for (size_t iPitch = 0; iPitch < m_nPitch; ++iPitch)
  90     {
  91         int theoreticalMinNextPitch = static_cast<int>(iPitch)-static_cast<int>(m_transitionWidth/2);
  92         size_t minNextPitch = iPitch>m_transitionWidth/2 ? iPitch-m_transitionWidth/2 : 0;
  93         size_t maxNextPitch = iPitch<m_nPitch-m_transitionWidth/2 ? iPitch+m_transitionWidth/2 : m_nPitch-1;
  94
  95         // WEIGHT VECTOR
  96         double weightSum = 0;
  97         vector<double> weights;
  98         for (size_t i = minNextPitch; i <= maxNextPitch; ++i)
  99         {
 100             if (i <= iPitch)
 101             {
 102                 weights.push_back(i-theoreticalMinNextPitch+1);
 103                 // weights.push_back(i-theoreticalMinNextPitch+1+m_transitionWidth/2);
 104             } else {
 105                 weights.push_back(iPitch-theoreticalMinNextPitch+1-(i-iPitch));
 106                 // weights.push_back(iPitch-theoreticalMinNextPitch+1-(i-iPitch)+m_transitionWidth/2);
 107             }
 108             weightSum += weights[weights.size()-1];
 109         }
 110
 111         // std::cerr << minNextPitch << "  " << maxNextPitch << std::endl;
 112         // TRANSITIONS TO CLOSE PITCH
 113         for (size_t i = minNextPitch; i <= maxNextPitch; ++i)
 114         {
 115             from.push_back(iPitch);
 116             to.push_back(i);
 117             transProb.push_back(weights[i-minNextPitch] / weightSum * m_selfTrans);
 118
 119             from.push_back(iPitch);
 120             to.push_back(i+m_nPitch);
 121             transProb.push_back(weights[i-minNextPitch] / weightSum * (1-m_selfTrans));
 122
 123             from.push_back(iPitch+m_nPitch);
 124             to.push_back(i+m_nPitch);
 125             transProb.push_back(weights[i-minNextPitch] / weightSum * m_selfTrans);
 126             // transProb.push_back(weights[i-minNextPitch] / weightSum * 0.5);
 127
 128             from.push_back(iPitch+m_nPitch);
 129             to.push_back(i);
 130             transProb.push_back(weights[i-minNextPitch] / weightSum * (1-m_selfTrans));
 131             // transProb.push_back(weights[i-minNextPitch] / weightSum * 0.5);
 132         }
 133
 134         // TRANSITION TO UNVOICED
 135         // from.push_back(iPitch+m_nPitch);
 136         // to.push_back(2*m_nPitch);
 137         // transProb.push_back(1-m_selfTrans);
 138
 139         // TRANSITION FROM UNVOICED TO PITCH
 140         // from.push_back(2*m_nPitch);
 141         // to.push_back(iPitch+m_nPitch);
 142         // transProb.push_back(1.0/m_nPitch);
 143     }
 144     // UNVOICED SELFTRANSITION
 145     // from.push_back(2*m_nPitch);
 146     // to.push_back(2*m_nPitch);
 147     // transProb.push_back(m_selfTrans);
 148
 149     // for (size_t i = 0; i < from.size(); ++i) {
 150     //     std::cerr << "P(["<< from[i] << " --> " << to[i] << "]) = " << transProb[i] << std::endl;
 151     // }
 152
 153 }