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: 
chris@52: Yin::Yin(size_t frameSize, size_t inputSampleRate, double thresh) : 
chris@52:     m_frameSize(frameSize),
chris@52:     m_inputSampleRate(inputSampleRate),
chris@52:     m_thresh(thresh),
chris@52:     m_threshDistr(2),
chris@52:     m_yinBufferSize(frameSize/2)
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@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
matthiasm@60:     YinUtil::slowDifference(in, yinBuffer, m_yinBufferSize);    
chris@52:     YinUtil::cumulativeDifference(yinBuffer, m_yinBufferSize);
chris@52: 
chris@52:     int tau = 0;
chris@52:     tau = YinUtil::absoluteThreshold(yinBuffer, m_yinBufferSize, m_thresh);
chris@52:         
chris@52:     double interpolatedTau;
chris@52:     double aperiodicity;
chris@52:     double f0;
chris@52:     
chris@54:     if (tau!=0 && tau!=m_yinBufferSize-1)
chris@52:     {
chris@52:         interpolatedTau = YinUtil::parabolicInterpolation(yinBuffer, abs(tau), m_yinBufferSize);
chris@52:         f0 = m_inputSampleRate * (1.0 / interpolatedTau);
chris@52:     } else {
chris@52:         interpolatedTau = 0;
chris@52:         f0 = 0;
chris@52:     }
chris@52:     double rms = std::sqrt(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@54: 
chris@52:     double* yinBuffer = new double[m_yinBufferSize];
chris@52: 
chris@52:     // calculate aperiodicity function for all periods
matthiasm@60:     YinUtil::slowDifference(in, yinBuffer, m_yinBufferSize);    
chris@52:     YinUtil::cumulativeDifference(yinBuffer, m_yinBufferSize);
chris@52: 
chris@52:     vector<double> peakProbability = YinUtil::yinProb(yinBuffer, m_threshDistr, m_yinBufferSize);
chris@54:         
chris@54:     // basic yin output
chris@54:     Yin::YinOutput yo(0,0,0);
chris@54:     for (size_t iBuf = 1; iBuf < m_yinBufferSize-1; ++iBuf)
chris@52:     {
chris@52:         if (peakProbability[iBuf] > 0)
chris@52:         {
chris@52:             double currentF0 = 
chris@52:                 m_inputSampleRate * (1.0 /
chris@52:                 YinUtil::parabolicInterpolation(yinBuffer, iBuf, m_yinBufferSize));
chris@52:             yo.freqProb.push_back(pair<double, double>(currentF0, peakProbability[iBuf]));
chris@52:         }
chris@52:     }
chris@52:     
chris@54:     // add salience
chris@54:     for (size_t iBuf = 0; iBuf < m_yinBufferSize; ++iBuf) {
chris@54:         yo.salience.push_back(peakProbability[iBuf]);
chris@54:     }
chris@54: 
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: 
chris@52: // int
chris@52: // Yin::setRemoveUnvoiced(bool parameter)
chris@52: // {
chris@52: //     m_removeUnvoiced = parameter;
chris@52: //     return 0;
chris@52: // }
chris@54: 
chris@54: float
chris@54: Yin::constrainedMinPick(const double *in, const float minFreq, const int maxFreq) const {
chris@54:     
chris@54:     double* yinBuffer = new double[m_yinBufferSize];
chris@54: 
chris@54:     // calculate aperiodicity function for all periods
matthiasm@60:     YinUtil::slowDifference(in, yinBuffer, m_yinBufferSize);    
chris@54:     YinUtil::cumulativeDifference(yinBuffer, m_yinBufferSize);
chris@54:     
chris@54:     int minPeriod = m_inputSampleRate / maxFreq;
chris@54:     int maxPeriod = m_inputSampleRate / minFreq;
chris@54:     
chris@54:     if (minPeriod < 0 || maxPeriod > m_yinBufferSize || minPeriod > maxPeriod) {
chris@54:         delete [] yinBuffer;
chris@54:         return 0.f;
chris@54:     }
chris@54:     
chris@54:     float bestVal = 1000;
chris@54:     int   bestTau = 0;
chris@54:     for (int tau = minPeriod; tau <= maxPeriod; ++tau)
chris@54:     {
chris@54:         if (yinBuffer[tau] < bestVal) 
chris@54:         {
chris@54:             bestVal = yinBuffer[tau];
chris@54:             bestTau = tau;
chris@54:         }
chris@54:     }
chris@54:     
chris@54:     float interpolatedTau =
chris@54:         YinUtil::parabolicInterpolation(yinBuffer, bestTau, m_yinBufferSize);
chris@54:     
chris@54:     delete [] yinBuffer;
chris@54:     return m_inputSampleRate * (1.0 / interpolatedTau);
matthiasm@49: }