chris@52: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
chris@52: 
chris@52: /*
chris@52:     pYIN - A fundamental frequency estimator for monophonic audio
chris@52:     Centre for Digital Music, Queen Mary, University of London.
chris@52:     
chris@52:     This program is free software; you can redistribute it and/or
chris@52:     modify it under the terms of the GNU General Public License as
chris@52:     published by the Free Software Foundation; either version 2 of the
chris@52:     License, or (at your option) any later version.  See the file
chris@52:     COPYING included with this distribution for more information.
chris@52: */
chris@52: 
chris@52: #include "Yin.h"
chris@52: 
chris@52: #include "vamp-sdk/FFT.h"
chris@52: #include "MeanFilter.h"
chris@52: #include "YinUtil.h"
chris@52: 
chris@52: #include <vector>
chris@52: #include <cstdlib>
chris@52: #include <cstdio>
chris@52: #include <cmath>
chris@52: #include <complex>
chris@52: 
chris@52: using std::vector;
chris@52: 
matthiasm@117: Yin::Yin(size_t frameSize, size_t inputSampleRate, double thresh, bool fast) : 
chris@52:     m_frameSize(frameSize),
chris@52:     m_inputSampleRate(inputSampleRate),
chris@52:     m_thresh(thresh),
chris@52:     m_threshDistr(2),
matthiasm@117:     m_yinBufferSize(frameSize/2),
Chris@136:     m_fast(fast),
Chris@136:     m_yinUtil(new YinUtil(m_yinBufferSize))
chris@52: {
chris@52:     if (frameSize & (frameSize-1)) {
chris@52:       //  throw "N must be a power of two";
chris@52:     }
chris@52: }
chris@52: 
chris@52: Yin::~Yin() 
chris@52: {
Chris@136:     delete m_yinUtil;
chris@52: }
chris@52: 
chris@52: Yin::YinOutput
chris@52: Yin::process(const double *in) const {
chris@52:     
chris@52:     double* yinBuffer = new double[m_yinBufferSize];
chris@52: 
chris@52:     // calculate aperiodicity function for all periods
Chris@136:     if (m_fast) m_yinUtil->fastDifference(in, yinBuffer);
Chris@136:     else m_yinUtil->slowDifference(in, yinBuffer);
matthiasm@117: 
Chris@136:     m_yinUtil->cumulativeDifference(yinBuffer);
chris@52: 
chris@52:     int tau = 0;
Chris@136:     tau = m_yinUtil->absoluteThreshold(yinBuffer, m_thresh);
chris@52:         
chris@52:     double interpolatedTau;
chris@52:     double aperiodicity;
chris@52:     double f0;
chris@52:     
chris@52:     if (tau!=0)
chris@52:     {
Chris@136:         interpolatedTau = m_yinUtil->parabolicInterpolation(yinBuffer, abs(tau));
chris@52:         f0 = m_inputSampleRate * (1.0 / interpolatedTau);
chris@52:     } else {
chris@52:         interpolatedTau = 0;
chris@52:         f0 = 0;
chris@52:     }
Chris@136:     double rms = std::sqrt(m_yinUtil->sumSquare(in, 0, m_yinBufferSize)/m_yinBufferSize);
chris@52:     aperiodicity = yinBuffer[abs(tau)];
chris@52:     // std::cerr << aperiodicity << std::endl;
chris@52:     if (tau < 0) f0 = -f0;
chris@52: 
chris@52:     Yin::YinOutput yo(f0, 1-aperiodicity, rms);
chris@52:     for (size_t iBuf = 0; iBuf < m_yinBufferSize; ++iBuf)
chris@52:     {
chris@52:         yo.salience.push_back(yinBuffer[iBuf] < 1 ? 1-yinBuffer[iBuf] : 0); // why are the values sometimes < 0 if I don't check?
chris@52:     }
chris@52:     
chris@52:     delete [] yinBuffer;
chris@52:     return yo;
chris@52: }
chris@52: 
chris@52: Yin::YinOutput
chris@52: Yin::processProbabilisticYin(const double *in) const {
chris@52:     
chris@52:     double* yinBuffer = new double[m_yinBufferSize];
chris@52: 
chris@52:     // calculate aperiodicity function for all periods
Chris@136:     if (m_fast) m_yinUtil->fastDifference(in, yinBuffer);
Chris@136:     else m_yinUtil->slowDifference(in, yinBuffer);
matthiasm@117: 
Chris@136:     m_yinUtil->cumulativeDifference(yinBuffer);
chris@52: 
Chris@136:     vector<double> peakProbability = m_yinUtil->yinProb(yinBuffer, m_threshDistr);
chris@52:     
chris@52:     // calculate overall "probability" from peak probability
chris@52:     double probSum = 0;
chris@52:     for (size_t iBin = 0; iBin < m_yinBufferSize; ++iBin)
chris@52:     {
chris@52:         probSum += peakProbability[iBin];
chris@52:     }
Chris@136:     double rms = std::sqrt(m_yinUtil->sumSquare(in, 0, m_yinBufferSize)/m_yinBufferSize);
matthiasm@103:     Yin::YinOutput yo(0,0,rms);
chris@52:     for (size_t iBuf = 0; iBuf < m_yinBufferSize; ++iBuf)
chris@52:     {
chris@52:         yo.salience.push_back(peakProbability[iBuf]);
chris@52:         if (peakProbability[iBuf] > 0)
chris@52:         {
chris@52:             double currentF0 = 
chris@52:                 m_inputSampleRate * (1.0 /
Chris@136:                 m_yinUtil->parabolicInterpolation(yinBuffer, iBuf));
chris@52:             yo.freqProb.push_back(pair<double, double>(currentF0, peakProbability[iBuf]));
chris@52:         }
chris@52:     }
chris@52:     
chris@52:     // std::cerr << yo.freqProb.size() << std::endl;
chris@52:     
chris@52:     delete [] yinBuffer;
chris@52:     return yo;
chris@52: }
chris@52: 
chris@52: 
chris@52: int
chris@52: Yin::setThreshold(double parameter)
chris@52: {
chris@52:     m_thresh = static_cast<float>(parameter);
chris@52:     return 0;
chris@52: }
chris@52: 
chris@52: int
chris@52: Yin::setThresholdDistr(float parameter)
chris@52: {
chris@52:     m_threshDistr = static_cast<size_t>(parameter);
chris@52:     return 0;
chris@52: }
chris@52: 
chris@52: int
chris@52: Yin::setFrameSize(size_t parameter)
chris@52: {
chris@52:     m_frameSize = parameter;
chris@52:     m_yinBufferSize = m_frameSize/2;
chris@52:     return 0;
chris@52: }
chris@52: 
matthiasm@117: int
matthiasm@117: Yin::setFast(bool parameter)
matthiasm@117: {
matthiasm@117:     m_fast = parameter;
matthiasm@117:     return 0;
matthiasm@117: }