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 "Yin.h"
matthiasm@0: 
matthiasm@0: #include "vamp-sdk/FFT.h"
matthiasm@0: #include "MeanFilter.h"
matthiasm@0: #include "YinUtil.h"
matthiasm@0: 
matthiasm@0: #include <vector>
matthiasm@0: #include <cstdlib>
matthiasm@0: #include <cstdio>
matthiasm@0: #include <cmath>
matthiasm@0: #include <complex>
matthiasm@0: 
matthiasm@0: using std::vector;
matthiasm@0: 
matthiasm@0: Yin::Yin(size_t frameSize, size_t inputSampleRate, double thresh) : 
matthiasm@0:     m_frameSize(frameSize),
matthiasm@0:     m_inputSampleRate(inputSampleRate),
matthiasm@0:     m_thresh(thresh),
matthiasm@0:     m_threshDistr(2),
matthiasm@0:     m_yinBufferSize(frameSize/2)
matthiasm@0: {
matthiasm@0:     if (frameSize & (frameSize-1)) {
matthiasm@0:         throw "N must be a power of two";
matthiasm@0:     }
matthiasm@0: }
matthiasm@0: 
matthiasm@0: Yin::~Yin() 
matthiasm@0: {
matthiasm@0: }
matthiasm@0: 
matthiasm@0: Yin::YinOutput
matthiasm@0: Yin::process(const double *in) const {
matthiasm@0:     
matthiasm@0:     double* yinBuffer = new double[m_yinBufferSize];
matthiasm@0: 
matthiasm@0:     // calculate aperiodicity function for all periods
matthiasm@0:     YinUtil::fastDifference(in, yinBuffer, m_yinBufferSize);    
matthiasm@0:     YinUtil::cumulativeDifference(yinBuffer, m_yinBufferSize);
matthiasm@0: 
matthiasm@0:     int tau = 0;
matthiasm@0:     tau = YinUtil::absoluteThreshold(yinBuffer, m_yinBufferSize, m_thresh);
matthiasm@0:         
matthiasm@0:     double interpolatedTau;
matthiasm@0:     double aperiodicity;
matthiasm@0:     double f0;
matthiasm@0:     
matthiasm@0:     if (tau!=0)
matthiasm@0:     {
matthiasm@0:         interpolatedTau = YinUtil::parabolicInterpolation(yinBuffer, abs(tau), m_yinBufferSize);
matthiasm@0:         f0 = m_inputSampleRate * (1.0 / interpolatedTau);
matthiasm@0:     } else {
matthiasm@0:         interpolatedTau = 0;
matthiasm@0:         f0 = 0;
matthiasm@0:     }
matthiasm@0:     double rms = std::sqrt(YinUtil::sumSquare(in, 0, m_yinBufferSize)/m_yinBufferSize);
matthiasm@0:     aperiodicity = yinBuffer[abs(tau)];
matthiasm@0:     // std::cerr << aperiodicity << std::endl;
matthiasm@0:     if (tau < 0) f0 = -f0;
matthiasm@0: 
matthiasm@0:     Yin::YinOutput yo(f0, 1-aperiodicity, rms);
matthiasm@0:     for (size_t iBuf = 0; iBuf < m_yinBufferSize; ++iBuf)
matthiasm@0:     {
matthiasm@0:         yo.salience.push_back(yinBuffer[iBuf] < 1 ? 1-yinBuffer[iBuf] : 0); // why are the values sometimes < 0 if I don't check?
matthiasm@0:     }
matthiasm@0:     
matthiasm@0:     delete [] yinBuffer;
matthiasm@0:     return yo;
matthiasm@0: }
matthiasm@0: 
matthiasm@0: Yin::YinOutput
matthiasm@0: Yin::processProbabilisticYin(const double *in) const {
matthiasm@0:     
matthiasm@0:     double* yinBuffer = new double[m_yinBufferSize];
matthiasm@0: 
matthiasm@0:     // calculate aperiodicity function for all periods
matthiasm@0:     YinUtil::fastDifference(in, yinBuffer, m_yinBufferSize);    
matthiasm@0:     YinUtil::cumulativeDifference(yinBuffer, m_yinBufferSize);
matthiasm@0: 
matthiasm@0:     vector<double> peakProbability = YinUtil::yinProb(yinBuffer, m_threshDistr, m_yinBufferSize);
matthiasm@0:     
matthiasm@0:     // calculate overall "probability" from peak probability
matthiasm@0:     double probSum = 0;
matthiasm@0:     for (size_t iBin = 0; iBin < m_yinBufferSize; ++iBin)
matthiasm@0:     {
matthiasm@0:         probSum += peakProbability[iBin];
matthiasm@0:     }
matthiasm@0:         
matthiasm@0:     Yin::YinOutput yo(0,0,0);
matthiasm@0:     for (size_t iBuf = 0; iBuf < m_yinBufferSize; ++iBuf)
matthiasm@0:     {
matthiasm@0:         yo.salience.push_back(peakProbability[iBuf]);
matthiasm@0:         if (peakProbability[iBuf] > 0)
matthiasm@0:         {
matthiasm@0:             double currentF0 = 
matthiasm@0:                 m_inputSampleRate * (1.0 /
matthiasm@0:                 YinUtil::parabolicInterpolation(yinBuffer, iBuf, m_yinBufferSize));
matthiasm@0:             yo.freqProb.push_back(pair<double, double>(currentF0, peakProbability[iBuf]));
matthiasm@0:         }
matthiasm@0:     }
matthiasm@0:     
matthiasm@0:     // std::cerr << yo.freqProb.size() << std::endl;
matthiasm@0:     
matthiasm@0:     delete [] yinBuffer;
matthiasm@0:     return yo;
matthiasm@0: }
matthiasm@0: 
matthiasm@0: 
matthiasm@0: int
matthiasm@0: Yin::setThreshold(double parameter)
matthiasm@0: {
matthiasm@0:     m_thresh = static_cast<float>(parameter);
matthiasm@0:     return 0;
matthiasm@0: }
matthiasm@0: 
matthiasm@0: int
matthiasm@0: Yin::setThresholdDistr(float parameter)
matthiasm@0: {
matthiasm@0:     m_threshDistr = static_cast<size_t>(parameter);
matthiasm@0:     return 0;
matthiasm@0: }
matthiasm@0: 
matthiasm@0: int
matthiasm@0: Yin::setFrameSize(size_t parameter)
matthiasm@0: {
matthiasm@0:     m_frameSize = parameter;
matthiasm@0:     m_yinBufferSize = m_frameSize/2;
matthiasm@0:     return 0;
matthiasm@0: }
matthiasm@0: 
matthiasm@0: // int
matthiasm@0: // Yin::setRemoveUnvoiced(bool parameter)
matthiasm@0: // {
matthiasm@0: //     m_removeUnvoiced = parameter;
matthiasm@0: //     return 0;
matthiasm@0: // }