changeset 0:99bac62ee2da

added PYIN sources, should be compileable

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Makefile.inc	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,82 @@
+
+PLUGIN_EXT	?= .so
+
+CXX	?= g++
+CC	?= gcc
+
+CFLAGS := $(CFLAGS) 
+CXXFLAGS := -I. $(CXXFLAGS)
+
+# LDFLAGS := $(LDFLAGS) -lvamp-sdk
+# PLUGIN_LDFLAGS := $(LDFLAGS) $(PLUGIN_LDFLAGS)
+# TEST_LDFLAGS := $(TEST_LDFLAGS) $(LDFLAGS) -lboost_unit_test_framework
+
+PLUGIN := yintony$(PLUGIN_EXT)
+
+SOURCES := PYIN.cpp \
+           VampYin.cpp \
+           Yin.cpp \
+           YinUtil.cpp \
+           MonoNote.cpp \
+           MonoPitch.cpp \
+           MonoNoteParameters.cpp \
+           SparseHMM.cpp \
+           MonoNoteHMM.cpp \
+           MonoPitchHMM.cpp \
+
+PLUGIN_MAIN := libmain.cpp
+
+TESTS := test/test-meanfilter \
+         test/test-fft \
+         test/test-yin \
+         test/test-mononote
+         
+OBJECTS := $(SOURCES:.cpp=.o)
+OBJECTS := $(OBJECTS:.c=.o)
+
+PLUGIN_OBJECTS := $(OBJECTS) $(PLUGIN_MAIN:.cpp=.o)
+
+all: $(PLUGIN) $(TESTS)
+	for t in $(TESTS); do echo "Running $$t"; ./"$$t" || exit 1; done
+
+plugin: $(PLUGIN)
+
+$(PLUGIN): $(PLUGIN_OBJECTS)
+	$(CXX) -o $@ $^ $(PLUGIN_LDFLAGS)
+
+test/test-meanfilter: test/TestMeanFilter.o $(OBJECTS)
+	$(CXX) -o $@ $^ $(TEST_LDFLAGS)
+
+test/test-fft: test/TestFFT.o $(OBJECTS)
+	$(CXX) -o $@ $^ $(TEST_LDFLAGS)
+	
+test/test-yin: test/TestYin.o $(OBJECTS)
+	$(CXX) -o $@ $^ $(TEST_LDFLAGS)
+	
+test/test-mononote: test/TestMonoNote.o $(OBJECTS)
+	$(CXX) -o $@ $^ $(TEST_LDFLAGS)
+
+clean:		
+		rm -f $(PLUGIN_OBJECTS) test/*.o
+
+distclean:	clean
+		rm -f $(PLUGIN) $(TESTS)
+
+depend:
+		makedepend -Y -fMakefile.inc *.cpp test/*.cpp *.h test/*.h
+
+# DO NOT DELETE
+
+PYIN.o: PYIN.h
+VampYin.o: VampYin.h
+Yin.o: Yin.h
+MonoNoteParameters.o: MonoNoteParameters.h
+MonoNote.o: MonoNote.h
+MonoPitch.o: MonoPitch.h
+MonoPitchHMM.o: MonoPitchHMM.h
+SparseHMM.o: SparseHMM.h
+MonoNoteHMM.o: MonoNoteHMM.h
+libmain.o: PYIN.h VampYin.h
+test/TestMeanFilter.o: MeanFilter.h
+test/TestYin.o: Yin.h
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Makefile.linux64	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,11 @@
+
+CFLAGS := -Wall -O3 -fPIC -I../vamp-plugin-sdk/
+CXXFLAGS := $(CFLAGS)
+
+PLUGIN_LDFLAGS := -shared -Wl,-Bstatic -L../vamp-plugin-sdk -lvamp-sdk -Wl,-Bdynamic 
+TEST_LDFLAGS := -Wl,-Bstatic -L../vamp-plugin-sdk -lvamp-sdk -Wl,-Bdynamic -lboost_unit_test_framework
+
+PLUGIN_EXT := .so
+
+include Makefile.inc
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Makefile.osx	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,11 @@
+ARCHFLAGS :=
+CFLAGS := $(ARCHFLAGS) -O3 -I../inst/include -I/usr/local/boost -Wall -fPIC 
+CXXFLAGS := $(CFLAGS)
+
+LDFLAGS := -lvamp-sdk $(ARCHFLAGS) 
+PLUGIN_LDFLAGS := -dynamiclib -lvamp-sdk
+TEST_LDFLAGS := -lvamp-sdk -lboost_unit_test_framework
+PLUGIN_EXT := .dylib
+
+include Makefile.inc
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/MeanFilter.h	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,74 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    This file is Copyright (c) 2012 Chris Cannam
+  
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
+    ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+
+#ifndef _MEAN_FILTER_H_
+#define _MEAN_FILTER_H_
+
+class MeanFilter
+{
+public:
+    /**
+     * Construct a non-causal mean filter with filter length flen,
+     * that replaces each sample N with the mean of samples
+     * [N-floor(F/2) .. N+floor(F/2)] where F is the filter length.
+     * Only odd F are supported.
+     */
+    MeanFilter(int flen) : m_flen(flen) { }
+    ~MeanFilter() { }
+
+    /**
+     * Filter the n samples in "in" and place the results in "out"
+     */
+    void filter(const double *in, double *out, const int n) {
+	filterSubsequence(in, out, n, n, 0);
+    }
+
+    /**
+     * Filter the n samples starting at the given offset in the
+     * m-element array "in" and place the results in the n-element
+     * array "out"
+     */
+    void filterSubsequence(const double *in, double *out,
+			   const int m, const int n,
+			   const int offset) {
+	int half = m_flen/2;
+	for (int i = 0; i < n; ++i) {
+	    double v = 0;
+	    int n = 0;
+	    for (int j = -half; j <= half; ++j) {
+		int ix = i + j + offset;
+		if (ix >= 0 && ix < m) {
+		    v += in[ix];
+		    ++n;
+		}
+	    }
+	    out[i] = v / n;
+	}
+    }
+
+private:
+    int m_flen;
+};
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/MonoNote.cpp	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,48 @@
+#include "MonoNote.h"
+#include <vector>
+
+#include <cstdio>
+#include <cmath>
+#include <complex>
+
+using std::vector;
+using std::pair;
+
+MonoNote::MonoNote() :
+    hmm()
+{
+}
+
+MonoNote::~MonoNote()
+{
+}
+
+const vector<MonoNote::FrameOutput>
+MonoNote::process(const vector<vector<pair<double, double> > > pitchProb)
+{
+    vector<vector<double> > obsProb;
+    for (size_t iFrame = 0; iFrame < pitchProb.size(); ++iFrame)
+    {
+        obsProb.push_back(hmm.calculateObsProb(pitchProb[iFrame]));
+    }
+    
+    vector<double> *scale = new vector<double>(pitchProb.size());
+    
+    vector<MonoNote::FrameOutput> out; 
+    
+    vector<int> path = hmm.decodeViterbi(obsProb, scale);
+    
+    for (size_t iFrame = 0; iFrame < path.size(); ++iFrame)
+    {
+        double currPitch = -1;
+        int stateKind = 0;
+
+        currPitch = hmm.par.minPitch + (path[iFrame]/hmm.par.nSPP) * 1.0/hmm.par.nPPS;
+        stateKind = (path[iFrame]) % hmm.par.nSPP + 1;
+
+        out.push_back(FrameOutput(iFrame, currPitch, stateKind));
+        // std::cerr << path[iFrame] << " -- "<< pitchProb[iFrame][0].first << " -- "<< currPitch << " -- " << stateKind << std::endl;
+    }
+    delete scale;
+    return(out);
+}
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/MonoNote.h	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,33 @@
+#ifndef _MONONOTE_H_
+#define _MONONOTE_H_
+
+#include "MonoNoteHMM.h"
+#include "MonoNoteParameters.h"
+
+#include <iostream>
+#include <vector>
+#include <exception>
+
+using std::vector;
+using std::pair;
+
+class MonoNote {
+public:
+    MonoNote();
+    virtual ~MonoNote();
+    
+    struct FrameOutput {
+        size_t frameNumber;
+        double pitch;
+        size_t noteState; // unvoiced, attack, stable, release, inter
+        FrameOutput() :  frameNumber(0), pitch(-1.0), noteState(0) { }
+        FrameOutput(size_t _frameNumber, double _pitch, size_t _noteState) :
+            frameNumber(_frameNumber), pitch(_pitch), noteState(_noteState) { }
+    };
+    // pitchProb is a frame-wise vector carrying a vector of pitch-probability pairs
+    const vector<FrameOutput> process(const vector<vector<pair<double, double> > > pitchProb);
+private:
+    MonoNoteHMM hmm;
+};
+
+#endif
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/MonoNoteHMM.cpp	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,206 @@
+#include "MonoNoteHMM.h"
+
+#include <boost/math/distributions.hpp>
+
+#include <cstdio>
+#include <cmath>
+
+using std::vector;
+using std::pair;
+
+MonoNoteHMM::MonoNoteHMM() :
+    par()
+{
+    build();
+}
+
+const vector<double>
+MonoNoteHMM::calculateObsProb(const vector<pair<double, double> > pitchProb)
+{
+    // pitchProb is a list of pairs (pitches and their probabilities)
+    
+    size_t nCandidate = pitchProb.size();
+    
+    // what is the probability of pitched
+    double pIsPitched = 0;
+    for (size_t iCandidate = 0; iCandidate < nCandidate; ++iCandidate)
+    {
+        // pIsPitched = pitchProb[iCandidate].second > pIsPitched ? pitchProb[iCandidate].second : pIsPitched;
+        pIsPitched += pitchProb[iCandidate].second;
+    }
+
+    // pIsPitched = std::pow(pIsPitched, (1-par.priorWeight)) * std::pow(par.priorPitchedProb, par.priorWeight);
+    pIsPitched = pIsPitched * (1-par.priorWeight) + par.priorPitchedProb * par.priorWeight;
+
+    vector<double> out = vector<double>(par.n);
+    double tempProbSum = 0;
+    for (size_t i = 0; i < par.n; ++i)
+    {
+        if (i % 4 != 2)
+        {
+            // std::cerr << getMidiPitch(i) << std::endl;
+            double tempProb = 0;
+            if (nCandidate > 0)
+            {
+                double minDist = 10000.0;
+                double minDistProb = 0;
+                size_t minDistCandidate = 0;
+                for (size_t iCandidate = 0; iCandidate < nCandidate; ++iCandidate)
+                {
+                    double currDist = std::abs(getMidiPitch(i)-pitchProb[iCandidate].first);
+                    if (currDist < minDist)
+                    {
+                        minDist = currDist;
+                        minDistProb = pitchProb[iCandidate].second;
+                        minDistCandidate = iCandidate;
+                    }
+                }
+                tempProb = std::pow(minDistProb, par.yinTrust) * boost::math::pdf(pitchDistr[i], pitchProb[minDistCandidate].first);
+            } else {
+                tempProb = 1;
+            }
+            tempProbSum += tempProb;
+            out[i] = tempProb;
+        }
+    }
+    
+    for (size_t i = 0; i < par.n; ++i)
+    {
+        if (i % 4 != 2)
+        {
+            if (tempProbSum > 0) 
+            {
+                out[i] = out[i] / tempProbSum * pIsPitched;
+            }
+        } else {
+            out[i] = (1-pIsPitched) / (par.nPPS * par.nS);
+        }
+    }
+
+    return(out);
+}
+
+void
+MonoNoteHMM::build()
+{
+    // the states are organised as follows:
+    // 0-2. lowest pitch
+    //    0. attack state
+    //    1. stable state
+    //    2. silent state
+    //    3. inter state
+    // 4-6. second-lowest pitch
+    //    4. attack state
+    //    ...
+    
+    // observation distributions
+    for (size_t iState = 0; iState < par.n; ++iState)
+    {
+        pitchDistr.push_back(boost::math::normal(0,1));
+        if (iState % 4 == 2)
+        {
+            // silent state starts tracking
+            init.push_back(1.0/(par.nS * par.nPPS));
+        } else {
+            init.push_back(0.0);            
+        }
+    }
+
+    for (size_t iPitch = 0; iPitch < (par.nS * par.nPPS); ++iPitch)
+    {
+        size_t index = iPitch * par.nSPP;
+        double mu = par.minPitch + iPitch * 1.0/par.nPPS;
+        pitchDistr[index] = boost::math::normal(mu, par.sigmaYinPitchAttack);
+        pitchDistr[index+1] = boost::math::normal(mu, par.sigmaYinPitchStable);
+        pitchDistr[index+2] = boost::math::normal(mu, 1.0); // dummy
+        pitchDistr[index+3] = boost::math::normal(mu, par.sigmaYinPitchInter);
+    }
+    
+    boost::math::normal noteDistanceDistr(0, par.sigma2Note);
+
+    for (size_t iPitch = 0; iPitch < (par.nS * par.nPPS); ++iPitch)
+    {
+        // loop through all notes and set sparse transition probabilities
+        size_t index = iPitch * par.nSPP;
+
+        // transitions from attack state
+        from.push_back(index);
+        to.push_back(index);
+        transProb.push_back(par.pAttackSelftrans);
+
+        from.push_back(index);
+        to.push_back(index+1);
+        transProb.push_back(1-par.pAttackSelftrans);
+
+        // transitions from stable state
+        from.push_back(index+1);
+        to.push_back(index+1); // to itself
+        transProb.push_back(par.pStableSelftrans);
+        
+        from.push_back(index+1);
+        to.push_back(index+2); // to silent
+        transProb.push_back(par.pStable2Silent);
+
+        from.push_back(index+1);
+        to.push_back(index+3); // to inter-note
+        transProb.push_back(1-par.pStableSelftrans-par.pStable2Silent);
+
+        // the "easy" transitions from silent state
+        from.push_back(index+2);
+        to.push_back(index+2);
+        transProb.push_back(par.pSilentSelftrans);
+        
+        // the "easy" inter state transition
+        from.push_back(index+3);
+        to.push_back(index+3);
+        transProb.push_back(par.pInterSelftrans);
+        
+        // the more complicated transitions from the silent and inter state
+        double probSumSilent = 0;
+        double probSumInter = 0;
+        vector<double> tempTransProbInter;
+        vector<double> tempTransProbSilent;
+        for (size_t jPitch = 0; jPitch < (par.nS * par.nPPS); ++jPitch)
+        {
+            int fromPitch = iPitch;
+            int toPitch = jPitch;
+            double semitoneDistance = std::abs(fromPitch - toPitch) * 1.0 / par.nPPS;
+            
+            // if (std::fmod(semitoneDistance, 1) == 0 && semitoneDistance > par.minSemitoneDistance)
+            if (semitoneDistance == 0 || (semitoneDistance > par.minSemitoneDistance && semitoneDistance < par.maxJump))
+            {
+                size_t toIndex = jPitch * par.nSPP; // note attack index
+
+                double tempWeightSilent = boost::math::pdf(noteDistanceDistr, semitoneDistance);
+                double tempWeightInter = semitoneDistance == 0 ? 0 : tempWeightSilent;
+                probSumSilent += tempWeightSilent;
+                probSumInter += tempWeightInter;
+
+                tempTransProbSilent.push_back(tempWeightSilent);
+                tempTransProbInter.push_back(tempWeightInter);
+
+                from.push_back(index+2);
+                to.push_back(toIndex);
+                from.push_back(index+3);
+                to.push_back(toIndex);                
+            }
+        }
+        for (size_t i = 0; i < tempTransProbSilent.size(); ++i)
+        {
+            transProb.push_back((1-par.pSilentSelftrans) * tempTransProbSilent[i]/probSumSilent);
+            transProb.push_back((1-par.pInterSelftrans) * tempTransProbInter[i]/probSumInter);
+        }
+    }
+}
+
+double
+MonoNoteHMM::getMidiPitch(size_t index)
+{
+    return pitchDistr[index].mean();
+}
+
+double
+MonoNoteHMM::getFrequency(size_t index)
+{
+    return 440 * pow(2, (pitchDistr[index].mean()-69)/12);
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/MonoNoteHMM.h	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,32 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    This file is Copyright (c) 2012 Matthias Mauch
+
+*/
+
+#ifndef _MONONOTEHMM_H_
+#define _MONONOTEHMM_H_
+
+#include "MonoNoteParameters.h"
+#include "SparseHMM.h"
+
+#include <boost/math/distributions.hpp>
+
+#include <vector>
+#include <cstdio>
+
+using std::vector;
+
+class MonoNoteHMM : public SparseHMM
+{
+public:
+    MonoNoteHMM();
+    const std::vector<double> calculateObsProb(const vector<pair<double, double> >);
+    double getMidiPitch(size_t index);
+    double getFrequency(size_t index);
+    void build();
+    MonoNoteParameters par;
+    vector<boost::math::normal> pitchDistr;
+};
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/MonoNoteParameters.cpp	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,31 @@
+#include "MonoNoteParameters.h"
+
+MonoNoteParameters::MonoNoteParameters() :
+    minPitch(36), 
+    nPPS(3), 
+    nS(43), 
+    nSPP(4), // states per pitch
+    n(0),
+    initPi(0), 
+    pAttackSelftrans(0.99),
+    pStableSelftrans(0.99),
+    pStable2Silent(0.005),
+    pSilentSelftrans(0.99), 
+    sigma2Note(1.5),
+    maxJump(13),
+    pInterSelftrans(0.99),
+    priorPitchedProb(.7),
+    priorWeight(0.5),
+    minSemitoneDistance(.5),
+    sigmaYinPitchAttack(4), 
+    sigmaYinPitchStable(0.9),
+    sigmaYinPitchInter(4),
+    yinTrust(0.1)
+{
+    // just in case someone put in a silly value for pRelease2Unvoiced
+    n = nPPS * nS * nSPP;
+}
+
+MonoNoteParameters::~MonoNoteParameters()
+{
+}
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/MonoNoteParameters.h	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,48 @@
+#ifndef _MONONOTEPARAMETERS_H_
+#define _MONONOTEPARAMETERS_H_
+
+#include <iostream>
+#include <vector>
+#include <exception>
+
+using std::vector;
+
+class MonoNoteParameters
+{
+public:
+    MonoNoteParameters();
+    virtual ~MonoNoteParameters();
+    
+    // model architecture parameters
+    size_t minPitch; // lowest pitch in MIDI notes
+    size_t nPPS; // number of pitches per semitone
+    size_t nS; // number of semitones
+    size_t nSPP; // number of states per pitch
+    size_t n; // number of states (will be calcualted from other parameters)
+    
+    // initial state probabilities
+    vector<double> initPi; 
+    
+    // transition parameters
+    double pAttackSelftrans;
+    double pStableSelftrans;
+    double pStable2Silent;
+    double pSilentSelftrans;
+    double sigma2Note; // standard deviation of next note Gaussian distribution
+    double maxJump;
+    double pInterSelftrans;
+    
+    double priorPitchedProb;
+    double priorWeight;
+
+    double minSemitoneDistance; // minimum distance for a transition
+    
+    double sigmaYinPitchAttack;
+    double sigmaYinPitchStable;
+    double sigmaYinPitchInter;
+    
+    double yinTrust;
+    
+};
+
+#endif
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/MonoPitch.cpp	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,74 @@
+#include "MonoPitch.h"
+#include "MonoPitchHMM.h"
+#include <vector>
+
+#include <cstdio>
+#include <cmath>
+#include <complex>
+
+using std::vector;
+using std::pair;
+
+MonoPitch::MonoPitch() :
+    hmm()
+{
+}
+
+MonoPitch::~MonoPitch()
+{
+}
+
+const vector<float>
+MonoPitch::process(const vector<vector<pair<double, double> > > pitchProb)
+{
+    // std::cerr << "before observation prob calculation" << std::endl;
+    vector<vector<double> > obsProb;
+    for (size_t iFrame = 0; iFrame < pitchProb.size(); ++iFrame)
+    {
+        obsProb.push_back(hmm.calculateObsProb(pitchProb[iFrame]));
+    }
+    for (size_t i = 0; i < obsProb[0].size(); ++i) {
+        obsProb[0][i] = 0;
+    }
+    obsProb[0][obsProb[0].size()-1] = 1;
+    
+    // std::cerr << "after observation prob calculation" << std::endl;
+    
+    vector<double> *scale = new vector<double>(pitchProb.size());
+    
+    vector<float> out; 
+    
+    // std::cerr << "before Viterbi decoding" << obsProb.size() << "ng" << obsProb[1].size() << std::endl;
+    vector<int> path = hmm.decodeViterbi(obsProb, scale);
+    // std::cerr << "after Viterbi decoding" << std::endl;
+    
+    for (size_t iFrame = 0; iFrame < path.size(); ++iFrame)
+    {
+        // std::cerr << path[iFrame] << " " << hmm.m_freqs[path[iFrame]] << std::endl;
+        float hmmFreq = hmm.m_freqs[path[iFrame]];
+        float bestFreq = 0;
+        float leastDist = 10000;
+        if (hmmFreq > 0)
+        {
+            // This was a Yin estimate, so try to get original pitch estimate back
+            // ... a bit hacky, since we could have direclty saved the frequency
+            // that was assigned to the HMM bin in hmm.calculateObsProb -- but would
+            // have had to rethink the interface of that method.
+            for (size_t iPitch = 0; iPitch < pitchProb[iFrame].size(); ++iPitch)
+            {
+                float freq = 440. * std::pow(2, (pitchProb[iFrame][iPitch].first - 69)/12);
+                float dist = std::abs(hmmFreq-freq);
+                if (dist < leastDist)
+                {
+                    leastDist = dist;
+                    bestFreq = freq;
+                }
+            }
+        } else {
+            bestFreq = hmmFreq;
+        }
+        out.push_back(bestFreq);
+    }
+    delete scale;
+    return(out);
+}
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/MonoPitch.h	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,24 @@
+#ifndef _MONOPITCH_H_
+#define _MONOPITCH_H_
+
+#include "MonoPitchHMM.h"
+
+#include <iostream>
+#include <vector>
+#include <exception>
+
+using std::vector;
+using std::pair;
+
+class MonoPitch {
+public:
+    MonoPitch();
+    virtual ~MonoPitch();
+    
+    // pitchProb is a frame-wise vector carrying a vector of pitch-probability pairs
+    const vector<float> process(const vector<vector<pair<double, double> > > pitchProb);
+private:
+    MonoPitchHMM hmm;
+};
+
+#endif
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/MonoPitchHMM.cpp	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,142 @@
+#include "MonoPitchHMM.h"
+
+#include <boost/math/distributions.hpp>
+
+#include <cstdio>
+#include <cmath>
+
+using std::vector;
+using std::pair;
+
+MonoPitchHMM::MonoPitchHMM() :
+m_minFreq(55),
+m_nBPS(10),
+m_nPitch(0),
+m_transitionWidth(0),
+m_selfTrans(0.99),
+m_yinTrust(.5),
+m_freqs(0)
+{
+    m_transitionWidth = 5*(m_nBPS/2) + 1;
+    m_nPitch = 48 * m_nBPS;
+    m_freqs = vector<double>(2*m_nPitch+1);
+    for (size_t iPitch = 0; iPitch < m_nPitch; ++iPitch)
+    {
+        m_freqs[iPitch] = m_minFreq * std::pow(2, iPitch * 1.0 / (12 * m_nBPS));
+        // m_freqs[iPitch+m_nPitch] = -2;
+        m_freqs[iPitch+m_nPitch] = -m_freqs[iPitch];
+        // std::cerr << "pitch " << iPitch << " = " << m_freqs[iPitch] << std::endl;
+    }
+    m_freqs[2*m_nPitch] = -1;
+    build();
+}
+
+const vector<double>
+MonoPitchHMM::calculateObsProb(const vector<pair<double, double> > pitchProb)
+{
+    vector<double> out = vector<double>(2*m_nPitch+1);
+    double probYinPitched = 0;
+    // BIN THE PITCHES
+    for (size_t iPair = 0; iPair < pitchProb.size(); ++iPair)
+    {
+        double freq = 440. * std::pow(2, (pitchProb[iPair].first - 69)/12);
+        if (freq <= m_minFreq) continue;
+        double d = 0;
+        double oldd = 1000;
+        for (size_t iPitch = 0; iPitch < m_nPitch; ++iPitch)
+        {
+            d = std::abs(freq-m_freqs[iPitch]);
+            if (oldd < d && iPitch > 0)
+            {
+                // previous bin must have been the closest
+                out[iPitch-1] = pitchProb[iPair].second;
+                probYinPitched += out[iPitch-1];
+                break;
+            }
+            oldd = d;
+        }
+    }
+    
+    double probReallyPitched = m_yinTrust * probYinPitched;
+    for (size_t iPitch = 0; iPitch < m_nPitch; ++iPitch)
+    {
+        if (probYinPitched > 0) out[iPitch] *= (probReallyPitched/probYinPitched) ;
+        out[iPitch+m_nPitch] = (1 - probReallyPitched) / m_nPitch;
+    }
+    // out[2*m_nPitch] = m_yinTrust * (1 - probYinPitched);
+    return(out);
+}
+
+void
+MonoPitchHMM::build()
+{
+    // INITIAL VECTOR
+    init = vector<double>(2*m_nPitch+1);
+    init[2*m_nPitch] = 1; // force first frame to be unvoiced.
+    
+    // TRANSITIONS
+    for (size_t iPitch = 0; iPitch < m_nPitch; ++iPitch)
+    {
+        int theoreticalMinNextPitch = static_cast<int>(iPitch)-static_cast<int>(m_transitionWidth/2);
+        int minNextPitch = iPitch>m_transitionWidth/2 ? iPitch-m_transitionWidth/2 : 0;
+        int maxNextPitch = iPitch<m_nPitch-m_transitionWidth/2 ? iPitch+m_transitionWidth/2 : m_nPitch-1;
+        
+        // WEIGHT VECTOR
+        double weightSum = 0;
+        vector<double> weights;
+        for (size_t i = minNextPitch; i <= maxNextPitch; ++i)
+        {
+            if (i <= iPitch)
+            {
+                weights.push_back(i-theoreticalMinNextPitch+1);
+                // weights.push_back(i-theoreticalMinNextPitch+1+m_transitionWidth/2);
+            } else {
+                weights.push_back(iPitch-theoreticalMinNextPitch+1-(i-iPitch));
+                // weights.push_back(iPitch-theoreticalMinNextPitch+1-(i-iPitch)+m_transitionWidth/2);
+            }
+            weightSum += weights[weights.size()-1];
+        }
+        
+        // std::cerr << minNextPitch << "  " << maxNextPitch << std::endl;
+        // TRANSITIONS TO CLOSE PITCH
+        for (size_t i = minNextPitch; i <= maxNextPitch; ++i)
+        {
+            from.push_back(iPitch);
+            to.push_back(i);
+            transProb.push_back(weights[i-minNextPitch] / weightSum * m_selfTrans);
+
+            from.push_back(iPitch);
+            to.push_back(i+m_nPitch);
+            transProb.push_back(weights[i-minNextPitch] / weightSum * (1-m_selfTrans));
+
+            from.push_back(iPitch+m_nPitch);
+            to.push_back(i+m_nPitch);
+            transProb.push_back(weights[i-minNextPitch] / weightSum * m_selfTrans);
+            // transProb.push_back(weights[i-minNextPitch] / weightSum * 0.5);
+            
+            from.push_back(iPitch+m_nPitch);
+            to.push_back(i);
+            transProb.push_back(weights[i-minNextPitch] / weightSum * (1-m_selfTrans));
+            // transProb.push_back(weights[i-minNextPitch] / weightSum * 0.5);
+        }
+
+        // TRANSITION TO UNVOICED
+        // from.push_back(iPitch+m_nPitch);
+        // to.push_back(2*m_nPitch);
+        // transProb.push_back(1-m_selfTrans);
+        
+        // TRANSITION FROM UNVOICED TO PITCH
+        from.push_back(2*m_nPitch);
+        to.push_back(iPitch+m_nPitch);
+        transProb.push_back(1.0/m_nPitch);
+    }
+    // UNVOICED SELFTRANSITION
+    // from.push_back(2*m_nPitch);
+    // to.push_back(2*m_nPitch);
+    // transProb.push_back(m_selfTrans);
+    
+    // for (size_t i = 0; i < from.size(); ++i) {
+    //     std::cerr << "P(["<< from[i] << " --> " << to[i] << "]) = " << transProb[i] << std::endl;
+    // }
+    
+}
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/MonoPitchHMM.h	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,36 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    This file is Copyright (c) 2012 Matthias Mauch
+
+*/
+
+#ifndef _MONOPITCHHMM_H_
+#define _MONOPITCHHMM_H_
+
+#include "SparseHMM.h"
+
+#include <boost/math/distributions.hpp>
+
+#include <vector>
+#include <cstdio>
+
+using std::vector;
+
+class MonoPitchHMM : public SparseHMM
+{
+public:
+    MonoPitchHMM();
+    const std::vector<double> calculateObsProb(const vector<pair<double, double> >);
+    // double getMidiPitch(size_t index);
+    // double getFrequency(size_t index);
+    void build();
+    double m_minFreq; // 82.40689f/2
+    size_t m_nBPS;
+    size_t m_nPitch;
+    size_t m_transitionWidth;
+    double m_selfTrans;
+    double m_yinTrust;
+    vector<double> m_freqs;
+};
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/PYIN.cpp	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,547 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    This file is Copyright (c) 2012 Chris Cannam
+  
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
+    ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+
+#include "PYIN.h"
+#include "MonoNote.h"
+#include "MonoPitch.h"
+
+#include "vamp-sdk/FFT.h"
+
+#include <vector>
+#include <algorithm>
+
+#include <cstdio>
+#include <cmath>
+#include <complex>
+
+using std::string;
+using std::vector;
+using Vamp::RealTime;
+
+
+PYIN::PYIN(float inputSampleRate) :
+    Plugin(inputSampleRate),
+    m_channels(0),
+    m_stepSize(256),
+    m_blockSize(2048),
+    m_fmin(40),
+    m_fmax(700),
+    m_yin(2048, inputSampleRate, 0.0),
+    m_oF0Candidates(0),
+    m_oF0Probs(0),
+    m_oVoicedProb(0),
+    m_oCandidateSalience(0),
+    m_oSmoothedPitchTrack(0),
+    m_oNotes(0),
+    m_threshDistr(2.0f),
+    m_outputUnvoiced(0),
+    m_pitchProb(0),
+    m_timestamp(0)
+{
+}
+
+PYIN::~PYIN()
+{
+}
+
+string
+PYIN::getIdentifier() const
+{
+    return "yintony";
+}
+
+string
+PYIN::getName() const
+{
+    return "Yintony";
+}
+
+string
+PYIN::getDescription() const
+{
+    return "Monophonic pitch and note tracking based on a probabilistic Yin extension.";
+}
+
+string
+PYIN::getMaker() const
+{
+    return "Matthias Mauch";
+}
+
+int
+PYIN::getPluginVersion() const
+{
+    // Increment this each time you release a version that behaves
+    // differently from the previous one
+    return 1;
+}
+
+string
+PYIN::getCopyright() const
+{
+    return "GPL";
+}
+
+PYIN::InputDomain
+PYIN::getInputDomain() const
+{
+    return TimeDomain;
+}
+
+size_t
+PYIN::getPreferredBlockSize() const
+{
+    return 2048;
+}
+
+size_t 
+PYIN::getPreferredStepSize() const
+{
+    return 256;
+}
+
+size_t
+PYIN::getMinChannelCount() const
+{
+    return 1;
+}
+
+size_t
+PYIN::getMaxChannelCount() const
+{
+    return 1;
+}
+
+PYIN::ParameterList
+PYIN::getParameterDescriptors() const
+{
+    ParameterList list;
+    
+    ParameterDescriptor d;
+
+    d.identifier = "threshdistr";
+    d.name = "Yin threshold distribution";
+    d.description = ".";
+    d.unit = "";
+    d.minValue = 0.0f;
+    d.maxValue = 7.0f;
+    d.defaultValue = 2.0f;
+    d.isQuantized = true;
+    d.quantizeStep = 1.0f;
+    d.valueNames.push_back("Uniform");
+    d.valueNames.push_back("Beta (mean 0.10)");
+    d.valueNames.push_back("Beta (mean 0.15)");
+    d.valueNames.push_back("Beta (mean 0.20)");
+    d.valueNames.push_back("Beta (mean 0.30)");
+    d.valueNames.push_back("Single Value 0.10");
+    d.valueNames.push_back("Single Value 0.15");
+    d.valueNames.push_back("Single Value 0.20");
+    list.push_back(d);
+
+    d.identifier = "outputunvoiced";
+    d.valueNames.clear();
+    d.name = "Output estimates classified as unvoiced?";
+    d.description = ".";
+    d.unit = "";
+    d.minValue = 0.0f;
+    d.maxValue = 2.0f;
+    d.defaultValue = 2.0f;
+    d.isQuantized = true;
+    d.quantizeStep = 1.0f;
+    d.valueNames.push_back("No");
+    d.valueNames.push_back("Yes");
+    d.valueNames.push_back("Yes, as negative frequencies");
+    list.push_back(d);
+
+    return list;
+}
+
+float
+PYIN::getParameter(string identifier) const
+{
+    if (identifier == "threshdistr") {
+            return m_threshDistr;
+    }
+    if (identifier == "outputunvoiced") {
+            return m_outputUnvoiced;
+    }
+    return 0.f;
+}
+
+void
+PYIN::setParameter(string identifier, float value) 
+{
+    if (identifier == "threshdistr")
+    {
+        m_threshDistr = value;
+    }
+    if (identifier == "outputunvoiced")
+    {
+        m_outputUnvoiced = value;
+    }
+    
+}
+
+PYIN::ProgramList
+PYIN::getPrograms() const
+{
+    ProgramList list;
+    return list;
+}
+
+string
+PYIN::getCurrentProgram() const
+{
+    return ""; // no programs
+}
+
+void
+PYIN::selectProgram(string name)
+{
+}
+
+PYIN::OutputList
+PYIN::getOutputDescriptors() const
+{
+    OutputList outputs;
+
+    OutputDescriptor d;
+    
+    int outputNumber = 0;
+
+    d.identifier = "f0candidates";
+    d.name = "F0 Candidates";
+    d.description = "Estimated fundamental frequency candidates.";
+    d.unit = "Hz";
+    d.hasFixedBinCount = false;
+    // d.binCount = 1;
+    d.hasKnownExtents = true;
+    d.minValue = m_fmin;
+    d.maxValue = 500;
+    d.isQuantized = false;
+    d.sampleType = OutputDescriptor::FixedSampleRate;
+    d.sampleRate = (m_inputSampleRate / m_stepSize);
+    d.hasDuration = false;
+    outputs.push_back(d);
+    m_oF0Candidates = outputNumber++;
+
+    d.identifier = "f0probs";
+    d.name = "Candidate Probabilities";
+    d.description = "Probabilities  of estimated fundamental frequency candidates.";
+    d.unit = "";
+    d.hasFixedBinCount = false;
+    // d.binCount = 1;
+    d.hasKnownExtents = true;
+    d.minValue = 0;
+    d.maxValue = 1;
+    d.isQuantized = false;
+    d.sampleType = OutputDescriptor::FixedSampleRate;
+    d.sampleRate = (m_inputSampleRate / m_stepSize);
+    d.hasDuration = false;
+    outputs.push_back(d);
+    m_oF0Probs = outputNumber++;
+    
+    d.identifier = "voicedprob";
+    d.name = "Voiced Probability";
+    d.description = "Probability that the signal is voiced according to Probabilistic Yin.";
+    d.unit = "";
+    d.hasFixedBinCount = true;
+    d.binCount = 1;
+    d.hasKnownExtents = true;
+    d.minValue = 0;
+    d.maxValue = 1;
+    d.isQuantized = false;
+    d.sampleType = OutputDescriptor::FixedSampleRate;
+    d.sampleRate = (m_inputSampleRate / m_stepSize);
+    d.hasDuration = false;
+    outputs.push_back(d);
+    m_oVoicedProb = outputNumber++;
+
+    d.identifier = "candidatesalience";
+    d.name = "Candidate Salience";
+    d.description = "Candidate Salience";
+    d.hasFixedBinCount = true;
+    d.binCount = m_blockSize / 2;
+    d.hasKnownExtents = true;
+    d.minValue = 0;
+    d.maxValue = 1;
+    d.isQuantized = false;
+    d.sampleType = OutputDescriptor::FixedSampleRate;
+    d.sampleRate = (m_inputSampleRate / m_stepSize);
+    d.hasDuration = false;
+    outputs.push_back(d);
+    m_oCandidateSalience = outputNumber++;
+    
+    d.identifier = "smoothedpitchtrack";
+    d.name = "Smoothed Pitch Track";
+    d.description = ".";
+    d.unit = "Hz";
+    d.hasFixedBinCount = true;
+    d.binCount = 1;
+    d.hasKnownExtents = false;
+    // d.minValue = 0;
+    // d.maxValue = 1;
+    d.isQuantized = false;
+    d.sampleType = OutputDescriptor::FixedSampleRate;
+    d.sampleRate = (m_inputSampleRate / m_stepSize);
+    d.hasDuration = false;
+    outputs.push_back(d);
+    m_oSmoothedPitchTrack = outputNumber++;
+
+    d.identifier = "notes";
+    d.name = "Notes";
+    d.description = "Derived fixed-pitch note frequencies";
+    // d.unit = "MIDI unit";
+    d.unit = "Hz";
+    d.hasFixedBinCount = true;
+    d.binCount = 1;
+    d.hasKnownExtents = false;
+    d.isQuantized = false;
+    d.sampleType = OutputDescriptor::VariableSampleRate;
+    d.sampleRate = (m_inputSampleRate / m_stepSize);
+    d.hasDuration = true;
+    outputs.push_back(d);
+    m_oNotes = outputNumber++;
+
+    d.identifier = "notepitchtrack";
+    d.name = "Note pitch track.";
+    d.description = "Estimated fundamental frequencies during notes.";
+    d.unit = "Hz";
+    d.hasFixedBinCount = true;
+    d.binCount = 1;
+    d.hasKnownExtents = true;
+    d.minValue = m_fmin;
+    d.maxValue = 75;
+    d.isQuantized = false;
+    d.sampleType = OutputDescriptor::VariableSampleRate;
+    // d.sampleRate = (m_inputSampleRate / m_stepSize);
+    d.hasDuration = false;
+    outputs.push_back(d);
+    m_oNotePitchTrack = outputNumber++;
+
+    return outputs;
+}
+
+bool
+PYIN::initialise(size_t channels, size_t stepSize, size_t blockSize)
+{
+    if (channels < getMinChannelCount() ||
+	channels > getMaxChannelCount()) return false;
+
+    std::cerr << "PYIN::initialise: channels = " << channels
+          << ", stepSize = " << stepSize << ", blockSize = " << blockSize
+          << std::endl;
+
+    m_channels = channels;
+    m_stepSize = stepSize;
+    m_blockSize = blockSize;
+    
+    reset();
+
+    return true;
+}
+
+void
+PYIN::reset()
+{    
+    m_yin.setThresholdDistr(m_threshDistr);
+    m_yin.setFrameSize(m_blockSize);
+    
+    m_pitchProb.clear();
+    m_timestamp.clear();
+    
+    std::cerr << "PYIN::reset"
+          << ", blockSize = " << m_blockSize
+          << std::endl;
+}
+
+PYIN::FeatureSet
+PYIN::process(const float *const *inputBuffers, RealTime timestamp)
+{
+    timestamp = timestamp + Vamp::RealTime::frame2RealTime(m_blockSize/4, lrintf(m_inputSampleRate));
+    FeatureSet fs;
+    
+    double *dInputBuffers = new double[m_blockSize];
+    for (size_t i = 0; i < m_blockSize; ++i) dInputBuffers[i] = inputBuffers[0][i];
+    
+    Yin::YinOutput yo = m_yin.processProbabilisticYin(dInputBuffers);
+    
+    Feature f;
+    f.hasTimestamp = true;
+    f.timestamp = timestamp;
+    for (size_t i = 0; i < yo.freqProb.size(); ++i)
+    {
+        f.values.push_back(yo.freqProb[i].first);
+    }
+    fs[m_oF0Candidates].push_back(f);
+    
+    f.values.clear();
+    float voicedProb = 0;
+    for (size_t i = 0; i < yo.freqProb.size(); ++i)
+    {
+        f.values.push_back(yo.freqProb[i].second);
+        voicedProb += yo.freqProb[i].second;
+    }
+    fs[m_oF0Probs].push_back(f);
+    
+    f.values.clear();
+    f.values.push_back(voicedProb);
+    fs[m_oVoicedProb].push_back(f);
+
+    f.values.clear();
+    float salienceSum = 0;
+    for (size_t iBin = 0; iBin < yo.salience.size(); ++iBin)
+    {
+        f.values.push_back(yo.salience[iBin]);
+        salienceSum += yo.salience[iBin];
+    }
+    fs[m_oCandidateSalience].push_back(f);
+
+    delete [] dInputBuffers;
+
+    vector<pair<double, double> > tempPitchProb;
+    for (size_t iCandidate = 0; iCandidate < yo.freqProb.size(); ++iCandidate)
+    {
+        double tempPitch = 12 * std::log(yo.freqProb[iCandidate].first/440)/std::log(2.) + 69;
+        tempPitchProb.push_back(pair<double, double>
+            (tempPitch, yo.freqProb[iCandidate].second));
+    }
+    m_pitchProb.push_back(tempPitchProb);
+        
+    m_timestamp.push_back(timestamp);
+
+    return fs;
+}
+
+PYIN::FeatureSet
+PYIN::getRemainingFeatures()
+{
+    FeatureSet fs;
+    Feature f;
+    f.hasTimestamp = true;
+    f.hasDuration = false;
+    
+    // MONO-PITCH STUFF
+    MonoPitch mp;
+    std::cerr << "before viterbi" << std::endl;
+    vector<float> mpOut = mp.process(m_pitchProb);
+    // std::cerr << "after viterbi " << mpOut.size() << " "<< m_timestamp.size() << std::endl;
+    for (size_t iFrame = 0; iFrame < mpOut.size(); ++iFrame)
+    {
+        if (mpOut[iFrame] < 0 && (m_outputUnvoiced==0)) continue;
+        f.timestamp = m_timestamp[iFrame];
+        f.values.clear();
+        if (m_outputUnvoiced == 1)
+        {
+            f.values.push_back(abs(mpOut[iFrame]));
+        } else {
+            f.values.push_back(mpOut[iFrame]);
+        }
+        
+        fs[m_oSmoothedPitchTrack].push_back(f);
+    }
+    
+    // // MONO-NOTE STUFF
+    // MonoNote mn;
+    // 
+    // vector<MonoNote::FrameOutput> mnOut = mn.process(m_pitchProb);
+    // 
+    // f.hasTimestamp = true;
+    // f.hasDuration = true;
+    // f.values.clear();
+    // f.values.push_back(0);
+    // 
+    // Feature fNoteFreqTrack;
+    // fNoteFreqTrack.hasTimestamp = true;
+    // fNoteFreqTrack.hasDuration = false;
+    // 
+    // 
+    // int oldState = -1;
+    // int onsetFrame = 0;
+    // int framesInNote = 0;
+    // double pitchSum = 0;
+    // 
+    // for (size_t iFrame = 0; iFrame < mnOut.size(); ++iFrame)
+    // {
+    //     if (std::pow(2,(mnOut[onsetFrame].pitch - 69) / 12) * 440 >= m_fmin && mnOut[iFrame].noteState != 2 && oldState == 2)
+    //     {
+    //         // greedy pitch track
+    //         vector<double> notePitchTrack;
+    //         for (size_t i = onsetFrame; i <= iFrame; ++i) 
+    //         {
+    //             fNoteFreqTrack.timestamp = m_timestamp[i];
+    //             
+    //             bool hasPitch = 0;
+    //             double bestProb = 0;
+    //             double bestPitch = 0;
+    //             
+    //             for (int iCandidate = (int(m_pitchProb[i].size())) - 1; iCandidate != -1; --iCandidate)
+    //             {
+    //                 // std::cerr << "writing :( " << std::endl;
+    //                 double tempPitch = m_pitchProb[i][iCandidate].first;
+    //                 if (std::abs(tempPitch-mnOut[onsetFrame].pitch) < 5 && m_pitchProb[i][iCandidate].second > bestProb)
+    //                 {
+    //                     bestProb = m_pitchProb[i][iCandidate].second;
+    //                     bestPitch = m_pitchProb[i][iCandidate].first;
+    //                     hasPitch = 1;
+    //                 }
+    //             }
+    //             if (hasPitch) {
+    //                 double tempFreq = std::pow(2,(bestPitch - 69) / 12) * 440;
+    //                 notePitchTrack.push_back(bestPitch);
+    //                 fNoteFreqTrack.values.clear();
+    //                 fNoteFreqTrack.values.push_back(tempFreq); // convert back to Hz (Vamp hosts prefer that)
+    //                 fs[m_oNotePitchTrack].push_back(fNoteFreqTrack);
+    //             }
+    //         }
+    //         
+    //         // closing old note
+    //         f.duration = m_timestamp[iFrame]-m_timestamp[onsetFrame];
+    //         double tempPitch = mnOut[onsetFrame].pitch;
+    //         size_t notePitchTrackSize = notePitchTrack.size();
+    //         if (notePitchTrackSize > 2) {
+    //             std::sort(notePitchTrack.begin(), notePitchTrack.end());
+    //             tempPitch = notePitchTrack[notePitchTrackSize/2]; // median
+    //         }
+    //         f.values[0] = std::pow(2,(tempPitch - 69) / 12) * 440; // convert back to Hz (Vamp hosts prefer that)
+    //         fs[m_oNotes].push_back(f);
+    //         
+    //    
+    //     }
+    // 
+    //     if (mnOut[iFrame].noteState == 1 && oldState != 1)
+    //     {
+    //         // open note
+    //         onsetFrame = iFrame;
+    //         f.timestamp = m_timestamp[iFrame];
+    //         pitchSum = 0;
+    //         framesInNote = 0;
+    //     }
+    //     
+    //     oldState = mnOut[iFrame].noteState;
+    // }
+    
+    
+    return fs;
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/PYIN.h	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,72 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    This file is Copyright (c) 2012 Matthias Mauch
+
+*/
+
+#ifndef _PYIN_H_
+#define _PYIN_H_
+
+#include <vamp-sdk/Plugin.h>
+
+#include "Yin.h"
+
+class PYIN : public Vamp::Plugin
+{
+public:
+    PYIN(float inputSampleRate);
+    virtual ~PYIN();
+
+    std::string getIdentifier() const;
+    std::string getName() const;
+    std::string getDescription() const;
+    std::string getMaker() const;
+    int getPluginVersion() const;
+    std::string getCopyright() const;
+
+    InputDomain getInputDomain() const;
+    size_t getPreferredBlockSize() const;
+    size_t getPreferredStepSize() const;
+    size_t getMinChannelCount() const;
+    size_t getMaxChannelCount() const;
+
+    ParameterList getParameterDescriptors() const;
+    float getParameter(std::string identifier) const;
+    void setParameter(std::string identifier, float value);
+
+    ProgramList getPrograms() const;
+    std::string getCurrentProgram() const;
+    void selectProgram(std::string name);
+
+    OutputList getOutputDescriptors() const;
+
+    bool initialise(size_t channels, size_t stepSize, size_t blockSize);
+    void reset();
+
+    FeatureSet process(const float *const *inputBuffers,
+                       Vamp::RealTime timestamp);
+
+    FeatureSet getRemainingFeatures();
+
+protected:
+    size_t m_channels;
+    size_t m_stepSize;
+    size_t m_blockSize;
+    float m_fmin;
+    float m_fmax;
+    Yin m_yin;
+    
+    mutable int m_oF0Candidates;
+    mutable int m_oF0Probs;
+    mutable int m_oVoicedProb;
+    mutable int m_oCandidateSalience;
+    mutable int m_oSmoothedPitchTrack;
+    mutable int m_oNotes;
+    mutable int m_oNotePitchTrack;
+    mutable float m_threshDistr;
+    mutable float m_outputUnvoiced;
+    mutable vector<vector<pair<double, double> > > m_pitchProb;
+    mutable vector<Vamp::RealTime> m_timestamp;
+};
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SparseHMM.cpp	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,128 @@
+#include "SparseHMM.h"
+#include <vector>
+#include <cstdio>
+#include <iostream>
+
+using std::vector;
+using std::pair;
+
+const vector<double>
+SparseHMM::calculateObsProb(const vector<pair<double, double> > data)
+{
+    // dummy (virtual?) implementation to be overloaded
+    return(vector<double>());
+}
+
+const std::vector<int> 
+SparseHMM::decodeViterbi(std::vector<vector<double> > obsProb,
+                         vector<double> *scale) 
+{
+    size_t nState = init.size();
+    size_t nFrame = obsProb.size();
+    
+    // check for consistency    
+    size_t nTrans = transProb.size();
+    
+    // declaring variables
+    std::vector<double> delta = std::vector<double>(nState);
+    std::vector<double> oldDelta = std::vector<double>(nState);
+    vector<vector<int> > psi; //  "matrix" of remembered indices of the best transitions
+    vector<int> path = vector<int>(nFrame, nState-1); // the final output path (current assignment arbitrary, makes sense only for Chordino, where nChord-1 is the "no chord" label)
+
+    double deltasum = 0;
+
+    // initialise first frame
+    for (size_t iState = 0; iState < nState; ++iState)
+    {
+        oldDelta[iState] = init[iState] * obsProb[0][iState];
+        // std::cerr << iState << " ----- " << init[iState] << std::endl;
+        deltasum += oldDelta[iState];
+    }
+
+    for (size_t iState = 0; iState < nState; ++iState)
+    {
+        oldDelta[iState] /= deltasum; // normalise (scale)
+        // std::cerr << oldDelta[iState] << std::endl;
+    }
+
+    scale->push_back(1.0/deltasum);
+    psi.push_back(vector<int>(nState,0));
+
+    // rest of forward step
+    for (size_t iFrame = 1; iFrame < nFrame; ++iFrame)
+    {
+        deltasum = 0;
+        psi.push_back(vector<int>(nState,0));
+
+        // calculate best previous state for every current state
+        size_t fromState;
+        size_t toState;
+        double currentTransProb;
+        double currentValue;
+        
+        // this is the "sparse" loop
+        for (size_t iTrans = 0; iTrans < nTrans; ++iTrans)
+        {
+            fromState = from[iTrans];
+            toState = to[iTrans];
+            currentTransProb = transProb[iTrans];
+            
+            currentValue = oldDelta[fromState] * currentTransProb;
+            if (currentValue > delta[toState])
+            {
+                delta[toState] = currentValue; // will be multiplied by the right obs later!
+                psi[iFrame][toState] = fromState;
+            }            
+        }
+        
+        for (size_t jState = 0; jState < nState; ++jState)
+        {
+            delta[jState] *= obsProb[iFrame][jState];
+            deltasum += delta[jState];
+        }
+
+        if (deltasum > 0)
+        {
+            for (size_t iState = 0; iState < nState; ++iState)
+            {
+                oldDelta[iState] = delta[iState] / deltasum; // normalise (scale)
+                delta[iState] = 0;
+            }
+            scale->push_back(1.0/deltasum);
+        } else
+        {
+            std::cerr << "WARNING: Viterbi has been fed some zero probabilities, at least they become zero in combination with the model." << std::endl;
+            for (size_t iState = 0; iState < nState; ++iState)
+            {
+                oldDelta[iState] = 1.0/nState;
+                delta[iState] = 0;
+            }
+            scale->push_back(1.0);
+        }
+    }
+
+    // initialise backward step
+    double bestValue = 0;
+    for (size_t iState = 0; iState < nState; ++iState)
+    {
+        double currentValue = oldDelta[iState];
+        if (currentValue > bestValue)
+        {
+            bestValue = currentValue;            
+            path[nFrame-1] = iState;
+        }
+    }
+
+    // rest of backward step
+    for (int iFrame = nFrame-2; iFrame != -1; --iFrame)
+    {
+        path[iFrame] = psi[iFrame+1][path[iFrame+1]];
+    }
+    
+    for (size_t iState = 0; iState < nState; ++iState)
+    {
+        // std::cerr << psi[2][iState] << std::endl;
+    }
+    
+    return path;
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SparseHMM.h	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,28 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    This file is Copyright (c) 2012 Matthias Mauch
+
+*/
+
+#ifndef _SPARSEHMM_H_
+#define _SPARSEHMM_H_
+
+#include <vector>
+#include <cstdio>
+
+using std::vector;
+using std::pair;
+
+class SparseHMM
+{
+public:
+    virtual const std::vector<double> calculateObsProb(const vector<pair<double, double> >);
+    const std::vector<int> decodeViterbi(std::vector<vector<double> > obs, 
+                                   vector<double> *scale);
+    vector<double> init;
+    vector<size_t> from;
+    vector<size_t> to;
+    vector<double> transProb;
+};
+
+#endif
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/VampYin.cpp	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,366 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#include "VampYin.h"
+#include "MonoNote.h"
+
+#include "vamp-sdk/FFT.h"
+
+#include <vector>
+#include <algorithm>
+
+#include <cstdio>
+#include <cmath>
+#include <complex>
+
+using std::string;
+using std::vector;
+using Vamp::RealTime;
+
+
+VampYin::VampYin(float inputSampleRate) :
+    Plugin(inputSampleRate),
+    m_channels(0),
+    m_stepSize(256),
+    m_blockSize(2048),
+    m_fmin(40),
+    m_fmax(1000),
+    m_yin(2048, inputSampleRate, 0.0),
+    m_outNoF0(0),
+    m_outNoPeriodicity(0),
+    m_outNoRms(0),
+    m_outNoSalience(0),
+    m_yinParameter(0.15f),
+    m_outputUnvoiced(0.0f)
+{
+}
+
+VampYin::~VampYin()
+{
+}
+
+string
+VampYin::getIdentifier() const
+{
+    return "yin";
+}
+
+string
+VampYin::getName() const
+{
+    return "Yin";
+}
+
+string
+VampYin::getDescription() const
+{
+    return "A vamp implementation of the Yin algorithm for monophonic frequency estimation.";
+}
+
+string
+VampYin::getMaker() const
+{
+    return "Matthias Mauch";
+}
+
+int
+VampYin::getPluginVersion() const
+{
+    // Increment this each time you release a version that behaves
+    // differently from the previous one
+    return 1;
+}
+
+string
+VampYin::getCopyright() const
+{
+    return "GPL";
+}
+
+VampYin::InputDomain
+VampYin::getInputDomain() const
+{
+    return TimeDomain;
+}
+
+size_t
+VampYin::getPreferredBlockSize() const
+{
+    return 2048;
+}
+
+size_t 
+VampYin::getPreferredStepSize() const
+{
+    return 256;
+}
+
+size_t
+VampYin::getMinChannelCount() const
+{
+    return 1;
+}
+
+size_t
+VampYin::getMaxChannelCount() const
+{
+    return 1;
+}
+
+VampYin::ParameterList
+VampYin::getParameterDescriptors() const
+{
+    ParameterList list;
+    
+    ParameterDescriptor d;
+    d.identifier = "yinThreshold";
+    d.name = "Yin threshold";
+    d.description = "The greedy Yin search for a low value difference function is done once a dip lower than this threshold is reached.";
+    d.unit = "";
+    d.minValue = 0.025f;
+    d.maxValue = 1.0f;
+    d.defaultValue = 0.15f;
+    d.isQuantized = true;
+    d.quantizeStep = 0.025f;
+        
+    list.push_back(d);
+
+    // d.identifier = "removeunvoiced";
+    // d.name = "Remove pitches classified as unvoiced.";
+    // d.description = "If ticked, then the pitch estimator will return the most likely pitch, even if it 'thinks' there isn't any.";
+    // d.unit = "";
+    // d.minValue = 0.0f;
+    // d.maxValue = 1.0f;
+    // d.defaultValue = 0.0f;
+    // d.isQuantized = true;
+    // d.quantizeStep = 1.0f;
+    // d.valueNames.clear();
+    // list.push_back(d);
+
+    d.identifier = "outputunvoiced";
+    d.valueNames.clear();
+    d.name = "Output estimates classified as unvoiced?";
+    d.description = ".";
+    d.unit = "";
+    d.minValue = 0.0f;
+    d.maxValue = 2.0f;
+    d.defaultValue = 2.0f;
+    d.isQuantized = true;
+    d.quantizeStep = 1.0f;
+    d.valueNames.push_back("No");
+    d.valueNames.push_back("Yes");
+    d.valueNames.push_back("Yes, as negative frequencies");
+    list.push_back(d);
+
+    return list;
+}
+
+float
+VampYin::getParameter(string identifier) const
+{
+    if (identifier == "yinThreshold") {
+        return m_yinParameter;
+    }
+    if (identifier == "outputunvoiced") {
+        return m_outputUnvoiced;
+    }
+    return 0.f;
+}
+
+void
+VampYin::setParameter(string identifier, float value) 
+{
+    if (identifier == "yinThreshold")
+    {
+        m_yinParameter = value;
+    }
+    if (identifier == "outputunvoiced")
+    {
+        m_outputUnvoiced = value;
+    }
+}
+
+VampYin::ProgramList
+VampYin::getPrograms() const
+{
+    ProgramList list;
+    return list;
+}
+
+string
+VampYin::getCurrentProgram() const
+{
+    return ""; // no programs
+}
+
+void
+VampYin::selectProgram(string name)
+{
+}
+
+VampYin::OutputList
+VampYin::getOutputDescriptors() const
+{
+    OutputList outputs;
+
+    OutputDescriptor d;
+    
+    int outputNumber = 0;
+
+    d.identifier = "f0";
+    d.name = "Estimated f0";
+    d.description = "Estimated fundamental frequency";
+    d.unit = "Hz";
+    d.hasFixedBinCount = true;
+    d.binCount = 1;
+    d.hasKnownExtents = true;
+    d.minValue = m_fmin;
+    d.maxValue = 500;
+    d.isQuantized = false;
+    d.sampleType = OutputDescriptor::FixedSampleRate;
+    d.sampleRate = (m_inputSampleRate / m_stepSize);
+    d.hasDuration = false;
+    outputs.push_back(d);
+    m_outNoF0 = outputNumber++;
+
+    d.identifier = "periodicity";
+    d.name = "Periodicity";
+    d.description = "by-product of Yin f0 estimation";
+    d.unit = "";
+    d.hasFixedBinCount = true;
+    d.binCount = 1;
+    d.hasKnownExtents = true;
+    d.minValue = 0;
+    d.maxValue = 1;
+    d.isQuantized = false;
+    d.sampleType = OutputDescriptor::FixedSampleRate;
+    d.sampleRate = (m_inputSampleRate / m_stepSize);
+    d.hasDuration = false;
+    outputs.push_back(d);
+    m_outNoPeriodicity = outputNumber++;
+
+    d.identifier = "rms";
+    d.name = "root mean square";
+    d.description = "Root mean square of the waveform.";
+    d.unit = "";
+    d.hasFixedBinCount = true;
+    d.binCount = 1;
+    d.hasKnownExtents = true;
+    d.minValue = 0;
+    d.maxValue = 1;
+    d.isQuantized = false;
+    d.sampleType = OutputDescriptor::FixedSampleRate;
+    d.sampleRate = (m_inputSampleRate / m_stepSize);
+    d.hasDuration = false;
+    outputs.push_back(d);
+    m_outNoRms = outputNumber++;
+
+    d.identifier = "salience";
+    d.name = "Salience";
+    d.description = "Yin Salience";
+    d.hasFixedBinCount = true;
+    d.binCount = m_blockSize / 2;
+    d.hasKnownExtents = true;
+    d.minValue = 0;
+    d.maxValue = 1;
+    d.isQuantized = false;
+    d.sampleType = OutputDescriptor::FixedSampleRate;
+    d.sampleRate = (m_inputSampleRate / m_stepSize);
+    d.hasDuration = false;
+    outputs.push_back(d);
+    m_outNoSalience = outputNumber++;
+
+    return outputs;
+}
+
+bool
+VampYin::initialise(size_t channels, size_t stepSize, size_t blockSize)
+{
+    if (channels < getMinChannelCount() ||
+	channels > getMaxChannelCount()) return false;
+
+    std::cerr << "VampYin::initialise: channels = " << channels
+          << ", stepSize = " << stepSize << ", blockSize = " << blockSize
+          << std::endl;
+
+    m_channels = channels;
+    m_stepSize = stepSize;
+    m_blockSize = blockSize;
+    
+    reset();
+
+    return true;
+}
+
+void
+VampYin::reset()
+{    
+    m_yin.setThreshold(m_yinParameter);
+    m_yin.setFrameSize(m_blockSize);
+        
+    std::cerr << "VampYin::reset: yin threshold set to " << (m_yinParameter)
+          << ", blockSize = " << m_blockSize
+          << std::endl;
+}
+
+VampYin::FeatureSet
+VampYin::process(const float *const *inputBuffers, RealTime timestamp)
+{
+    timestamp = timestamp + Vamp::RealTime::frame2RealTime(m_blockSize/4, lrintf(m_inputSampleRate));
+    FeatureSet fs;
+    
+    double *dInputBuffers = new double[m_blockSize];
+    for (size_t i = 0; i < m_blockSize; ++i) dInputBuffers[i] = inputBuffers[0][i];
+    
+    Yin::YinOutput yo = m_yin.process(dInputBuffers);
+    // std::cerr << "f0 in VampYin: " << yo.f0 << std::endl;
+    Feature f;
+    f.hasTimestamp = true;
+    f.timestamp = timestamp;
+    if (m_outputUnvoiced == 0.0f)
+    {
+        // std::cerr << "f0 in VampYin: " << yo.f0 << std::endl;
+        if (yo.f0 > 0 && yo.f0 < m_fmax && yo.f0 > m_fmin) {
+            f.values.push_back(yo.f0);
+            fs[m_outNoF0].push_back(f);
+        }
+    } else if (m_outputUnvoiced == 1.0f)
+    {
+        if (abs(yo.f0) < m_fmax && abs(yo.f0) > m_fmin) {
+            f.values.push_back(abs(yo.f0));
+            fs[m_outNoF0].push_back(f);
+        }
+    } else
+    {
+        if (abs(yo.f0) < m_fmax && abs(yo.f0) > m_fmin) {
+            f.values.push_back(yo.f0);
+            fs[m_outNoF0].push_back(f);
+        }
+    }
+
+    f.values.clear();
+    f.values.push_back(yo.rms);
+    fs[m_outNoRms].push_back(f);
+    
+    f.values.clear();
+    for (size_t iBin = 0; iBin < yo.salience.size(); ++iBin)
+    {
+        f.values.push_back(yo.salience[iBin]);
+    }
+    fs[m_outNoSalience].push_back(f);
+    
+    f.values.clear();
+    // f.values[0] = yo.periodicity;
+    f.values.push_back(yo.periodicity);
+    fs[m_outNoPeriodicity].push_back(f);
+    
+    delete [] dInputBuffers;
+
+    return fs;
+}
+
+VampYin::FeatureSet
+VampYin::getRemainingFeatures()
+{
+    FeatureSet fs;
+    return fs;
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/VampYin.h	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,67 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    This file is Copyright (c) 2012 Matthias Mauch
+
+*/
+
+#ifndef _VAMPYIN_H_
+#define _VAMPYIN_H_
+
+#include <vamp-sdk/Plugin.h>
+
+#include "Yin.h"
+
+class VampYin : public Vamp::Plugin
+{
+public:
+    VampYin(float inputSampleRate);
+    virtual ~VampYin();
+
+    std::string getIdentifier() const;
+    std::string getName() const;
+    std::string getDescription() const;
+    std::string getMaker() const;
+    int getPluginVersion() const;
+    std::string getCopyright() const;
+
+    InputDomain getInputDomain() const;
+    size_t getPreferredBlockSize() const;
+    size_t getPreferredStepSize() const;
+    size_t getMinChannelCount() const;
+    size_t getMaxChannelCount() const;
+
+    ParameterList getParameterDescriptors() const;
+    float getParameter(std::string identifier) const;
+    void setParameter(std::string identifier, float value);
+
+    ProgramList getPrograms() const;
+    std::string getCurrentProgram() const;
+    void selectProgram(std::string name);
+
+    OutputList getOutputDescriptors() const;
+
+    bool initialise(size_t channels, size_t stepSize, size_t blockSize);
+    void reset();
+
+    FeatureSet process(const float *const *inputBuffers,
+                       Vamp::RealTime timestamp);
+
+    FeatureSet getRemainingFeatures();
+
+protected:
+    size_t m_channels;
+    size_t m_stepSize;
+    size_t m_blockSize;
+    float m_fmin;
+    float m_fmax;
+    Yin m_yin;
+    
+    mutable int m_outNoF0;
+    mutable int m_outNoPeriodicity;
+    mutable int m_outNoRms;
+    mutable int m_outNoSalience;
+    mutable float m_yinParameter;
+    mutable float m_outputUnvoiced;
+};
+
+#endif

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

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Yin.h	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,80 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    This file is Copyright (c) 2012 Chris Cannam
+  
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
+    ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+
+#ifndef _YIN_H_
+#define _YIN_H_
+
+#include "vamp-sdk/FFT.h"
+#include "MeanFilter.h"
+
+#include <cmath>
+
+#include <iostream>
+#include <vector>
+#include <exception>
+
+using std::vector;
+using std::pair;
+
+
+
+class Yin
+{
+public:
+    Yin(size_t frameSize, size_t inputSampleRate, double thresh = 0.2);
+    virtual ~Yin();
+
+    struct YinOutput {
+        double f0;
+        double periodicity;
+        double rms;
+        vector<double> salience;
+        vector<pair<double, double> > freqProb;
+        YinOutput() :  f0(0), periodicity(0), rms(0), 
+            salience(vector<double>(0)), freqProb(vector<pair<double, double> >(0)) { }
+        YinOutput(double _f, double _p, double _r) :
+            f0(_f), periodicity(_p), rms(_r), 
+            salience(vector<double>(0)), freqProb(vector<pair<double, double> >(0)) { }
+        YinOutput(double _f, double _p, double _r, vector<double> _salience) :
+            f0(_f), periodicity(_p), rms(_r), salience(_salience), 
+            freqProb(vector<pair<double, double> >(0)) { }
+    };
+    
+    int setThreshold(double parameter);
+    int setThresholdDistr(float parameter);
+    int setFrameSize(size_t frameSize);
+    // int setRemoveUnvoiced(bool frameSize);
+    YinOutput process(const double *in) const;
+    YinOutput processProbabilisticYin(const double *in) const;
+
+private:
+    mutable size_t m_frameSize;
+    mutable size_t m_inputSampleRate;
+    mutable double m_thresh;
+    mutable size_t m_threshDistr;
+    mutable size_t m_yinBufferSize;
+    // mutable bool m_removeUnvoiced;
+};
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/YinUtil.cpp	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,343 @@
+#include "YinUtil.h"
+
+#include <vector>
+
+#include <cstdio>
+#include <cmath>
+#include <algorithm>
+
+#include <boost/math/distributions.hpp>
+
+void 
+YinUtil::fastDifference(const double *in, double *yinBuffer, const size_t yinBufferSize) 
+{
+    
+    // DECLARE AND INITIALISE
+    // initialisation of most of the arrays here was done in a separate function,
+    // with all the arrays as members of the class... moved them back here.
+    
+    size_t frameSize = 2 * yinBufferSize;
+    
+    for (size_t j = 0; j < yinBufferSize; ++j)
+    {
+        yinBuffer[j] = 0.;
+    }
+
+    double *audioTransformedReal = new double[frameSize];
+    double *audioTransformedImag = new double[frameSize];
+    double *nullImag = new double[frameSize];
+    double *kernel = new double[frameSize];
+    double *kernelTransformedReal = new double[frameSize];
+    double *kernelTransformedImag = new double[frameSize];
+    double *yinStyleACFReal = new double[frameSize];
+    double *yinStyleACFImag = new double[frameSize];
+    double *powerTerms = new double[yinBufferSize];
+    
+    for (size_t j = 0; j < yinBufferSize; ++j)
+    {
+        powerTerms[j] = 0.;
+    }
+    
+    for (size_t j = 0; j < frameSize; ++j)
+    {
+        nullImag[j] = 0.;
+        audioTransformedReal[j] = 0.;
+        audioTransformedImag[j] = 0.;
+        kernel[j] = 0.;
+        kernelTransformedReal[j] = 0.;
+        kernelTransformedImag[j] = 0.;
+        yinStyleACFReal[j] = 0.;
+        yinStyleACFImag[j] = 0.;
+    }
+    
+    // POWER TERM CALCULATION
+    // ... for the power terms in equation (7) in the Yin paper
+    powerTerms[0] = 0.0;
+    for (size_t j = 0; j < yinBufferSize; ++j) {
+        powerTerms[0] += in[j] * in[j];
+    }
+
+    // now iteratively calculate all others (saves a few multiplications)
+    for (size_t tau = 1; tau < yinBufferSize; ++tau) {
+        powerTerms[tau] = powerTerms[tau-1] - in[tau-1] * in[tau-1] + in[tau+yinBufferSize] * in[tau+yinBufferSize];  
+    }
+
+    // YIN-STYLE AUTOCORRELATION via FFT
+    // 1. data
+    Vamp::FFT::forward(frameSize, in, nullImag, audioTransformedReal, audioTransformedImag);
+    
+    // 2. half of the data, disguised as a convolution kernel
+    for (size_t j = 0; j < yinBufferSize; ++j) {
+        kernel[j] = in[yinBufferSize-1-j];
+        kernel[j+yinBufferSize] = 0;
+    }
+    Vamp::FFT::forward(frameSize, kernel, nullImag, kernelTransformedReal, kernelTransformedImag);
+
+    // 3. convolution via complex multiplication -- written into
+    for (size_t j = 0; j < frameSize; ++j) {
+        yinStyleACFReal[j] = audioTransformedReal[j]*kernelTransformedReal[j] - audioTransformedImag[j]*kernelTransformedImag[j]; // real
+        yinStyleACFImag[j] = audioTransformedReal[j]*kernelTransformedImag[j] + audioTransformedImag[j]*kernelTransformedReal[j]; // imaginary
+    }
+    Vamp::FFT::inverse(frameSize, yinStyleACFReal, yinStyleACFImag, audioTransformedReal, audioTransformedImag);
+    
+    // CALCULATION OF difference function
+    // ... according to (7) in the Yin paper.
+    for (size_t j = 0; j < yinBufferSize; ++j) {
+        // taking only the real part
+        yinBuffer[j] = powerTerms[0] + powerTerms[j] - 2 * audioTransformedReal[j+yinBufferSize-1];
+    }
+    delete [] audioTransformedReal;
+    delete [] audioTransformedImag;
+    delete [] nullImag;
+    delete [] kernel;
+    delete [] kernelTransformedReal;
+    delete [] kernelTransformedImag;
+    delete [] yinStyleACFReal;
+    delete [] yinStyleACFImag;
+    delete [] powerTerms;
+}
+
+void 
+YinUtil::cumulativeDifference(double *yinBuffer, const size_t yinBufferSize)
+{    
+    size_t tau;
+    
+    yinBuffer[0] = 1;
+    
+    double runningSum = 0;
+    
+    for (tau = 1; tau < yinBufferSize; ++tau) {
+        runningSum += yinBuffer[tau];
+        if (runningSum == 0)
+        {
+            yinBuffer[tau] = 1;
+        } else {
+            yinBuffer[tau] *= tau / runningSum;
+        }
+    }    
+}
+
+int 
+YinUtil::absoluteThreshold(const double *yinBuffer, const size_t yinBufferSize, const double thresh)
+{
+    size_t tau;
+    size_t minTau = 0;
+    double minVal = 1000.;
+    
+    // using Joren Six's "loop construct" from TarsosDSP
+    tau = 2;
+    while (tau < yinBufferSize)
+    {
+        if (yinBuffer[tau] < thresh)
+        {
+            while (tau+1 < yinBufferSize && yinBuffer[tau+1] < yinBuffer[tau])
+            {
+                ++tau;
+            }
+            return tau;
+        } else {
+            if (yinBuffer[tau] < minVal)
+            {
+                minVal = yinBuffer[tau];
+                minTau = tau;
+            }
+        }
+        ++tau;
+    }
+    if (minTau > 0)
+    {
+        return -minTau;
+    }
+    return 0;
+}
+
+
+std::vector<double>
+YinUtil::yinProb(const double *yinBuffer, const size_t prior, const size_t yinBufferSize) 
+{
+    double minWeight = 0.01;
+    size_t tau;
+    std::vector<float> thresholds;
+    std::vector<float> distribution;
+    std::vector<double> peakProb = std::vector<double>(yinBufferSize);
+    float uniformDist[100] = {0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000};
+    float betaDist1[100] = {0.028911,0.048656,0.061306,0.068539,0.071703,0.071877,0.069915,0.066489,0.062117,0.057199,0.052034,0.046844,0.041786,0.036971,0.032470,0.028323,0.024549,0.021153,0.018124,0.015446,0.013096,0.011048,0.009275,0.007750,0.006445,0.005336,0.004397,0.003606,0.002945,0.002394,0.001937,0.001560,0.001250,0.000998,0.000792,0.000626,0.000492,0.000385,0.000300,0.000232,0.000179,0.000137,0.000104,0.000079,0.000060,0.000045,0.000033,0.000024,0.000018,0.000013,0.000009,0.000007,0.000005,0.000003,0.000002,0.000002,0.000001,0.000001,0.000001,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000};
+    float betaDist2[100] = {0.012614,0.022715,0.030646,0.036712,0.041184,0.044301,0.046277,0.047298,0.047528,0.047110,0.046171,0.044817,0.043144,0.041231,0.039147,0.036950,0.034690,0.032406,0.030133,0.027898,0.025722,0.023624,0.021614,0.019704,0.017900,0.016205,0.014621,0.013148,0.011785,0.010530,0.009377,0.008324,0.007366,0.006497,0.005712,0.005005,0.004372,0.003806,0.003302,0.002855,0.002460,0.002112,0.001806,0.001539,0.001307,0.001105,0.000931,0.000781,0.000652,0.000542,0.000449,0.000370,0.000303,0.000247,0.000201,0.000162,0.000130,0.000104,0.000082,0.000065,0.000051,0.000039,0.000030,0.000023,0.000018,0.000013,0.000010,0.000007,0.000005,0.000004,0.000003,0.000002,0.000001,0.000001,0.000001,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000};
+    float betaDist3[100] = {0.006715,0.012509,0.017463,0.021655,0.025155,0.028031,0.030344,0.032151,0.033506,0.034458,0.035052,0.035331,0.035332,0.035092,0.034643,0.034015,0.033234,0.032327,0.031314,0.030217,0.029054,0.027841,0.026592,0.025322,0.024042,0.022761,0.021489,0.020234,0.019002,0.017799,0.016630,0.015499,0.014409,0.013362,0.012361,0.011407,0.010500,0.009641,0.008830,0.008067,0.007351,0.006681,0.006056,0.005475,0.004936,0.004437,0.003978,0.003555,0.003168,0.002814,0.002492,0.002199,0.001934,0.001695,0.001481,0.001288,0.001116,0.000963,0.000828,0.000708,0.000603,0.000511,0.000431,0.000361,0.000301,0.000250,0.000206,0.000168,0.000137,0.000110,0.000088,0.000070,0.000055,0.000043,0.000033,0.000025,0.000019,0.000014,0.000010,0.000007,0.000005,0.000004,0.000002,0.000002,0.000001,0.000001,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000};
+    float betaDist4[100] = {0.003996,0.007596,0.010824,0.013703,0.016255,0.018501,0.020460,0.022153,0.023597,0.024809,0.025807,0.026607,0.027223,0.027671,0.027963,0.028114,0.028135,0.028038,0.027834,0.027535,0.027149,0.026687,0.026157,0.025567,0.024926,0.024240,0.023517,0.022763,0.021983,0.021184,0.020371,0.019548,0.018719,0.017890,0.017062,0.016241,0.015428,0.014627,0.013839,0.013068,0.012315,0.011582,0.010870,0.010181,0.009515,0.008874,0.008258,0.007668,0.007103,0.006565,0.006053,0.005567,0.005107,0.004673,0.004264,0.003880,0.003521,0.003185,0.002872,0.002581,0.002312,0.002064,0.001835,0.001626,0.001434,0.001260,0.001102,0.000959,0.000830,0.000715,0.000612,0.000521,0.000440,0.000369,0.000308,0.000254,0.000208,0.000169,0.000136,0.000108,0.000084,0.000065,0.000050,0.000037,0.000027,0.000019,0.000014,0.000009,0.000006,0.000004,0.000002,0.000001,0.000001,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000};
+    // float betaDist4[100] = {0.0000119,0.0000470,0.0001048,0.0001843,0.0002850,0.0004061,0.0005469,0.0007066,0.0008846,0.0010801,0.0012924,0.0015208,0.0017645,0.0020229,0.0022952,0.0025807,0.0028787,0.0031885,0.0035093,0.0038404,0.0041811,0.0045307,0.0048884,0.0052536,0.0056256,0.0060035,0.0063867,0.0067744,0.0071660,0.0075608,0.0079579,0.0083567,0.0087564,0.0091564,0.0095560,0.0099543,0.0103507,0.0107444,0.0111348,0.0115212,0.0119027,0.0122787,0.0126484,0.0130112,0.0133663,0.0137131,0.0140506,0.0143784,0.0146956,0.0150015,0.0152954,0.0155766,0.0158443,0.0160979,0.0163366,0.0165597,0.0167665,0.0169563,0.0171282,0.0172817,0.0174160,0.0175304,0.0176241,0.0176965,0.0177468,0.0177743,0.0177782,0.0177579,0.0177127,0.0176418,0.0175444,0.0174200,0.0172677,0.0170868,0.0168767,0.0166365,0.0163657,0.0160634,0.0157289,0.0153615,0.0149606,0.0145253,0.0140550,0.0135489,0.0130063,0.0124265,0.0118088,0.0111525,0.0104568,0.0097210,0.0089444,0.0081263,0.0072659,0.0063626,0.0054155,0.0044241,0.0033876,0.0023052,0.0011762,0.0000000};
+    float single10[100] = {0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,1.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000};
+    float single15[100] = {0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,1.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000};
+    float single20[100] = {0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,1.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000};
+    
+    // float betaDist4[100] = {0.00001,0.00005,0.00010,0.00018,0.00029,0.00041,0.00055,0.00071,0.00088,0.00108,0.00129,0.00152,0.00176,0.00202,0.00230,0.00258,0.00288,0.00319,0.00351,0.00384,0.00418,0.00453,0.00489,0.00525,0.00563,0.00600,0.00639,0.00677,0.00717,0.00756,0.00796,0.00836,0.00876,0.00916,0.00956,0.00995,0.01035,0.01074,0.01113,0.01152,0.01190,0.01228,0.01265,0.01301,0.01337,0.01371,0.01405,0.01438,0.01470,0.01500,0.01530,0.01558,0.01584,0.01610,0.01634,0.01656,0.01677,0.01696,0.01713,0.01728,0.01742,0.01753,0.01762,0.01770,0.01775,0.01777,0.01778,0.01776,0.01771,0.01764,0.01754,0.01742,0.01727,0.01709,0.01688,0.01664,0.01637,0.01606,0.01573,0.01536,0.01496,0.01453,0.01405,0.01355,0.01301,0.01243,0.01181,0.01115,0.01046,0.00972,0.00894,0.00813,0.00727,0.00636,0.00542,0.00442,0.00339,0.00231,0.00118,0.00000};
+    // float uniformDist[40] = {0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500};
+    // float betaDist1[40] = {0.00632,0.02052,0.03731,0.05327,0.06644,0.07587,0.08133,0.08305,0.08153,0.07743,0.07144,0.06421,0.05633,0.04831,0.04052,0.03326,0.02671,0.02098,0.01611,0.01209,0.00884,0.00629,0.00436,0.00292,0.00189,0.00118,0.00070,0.00040,0.00021,0.00011,0.00005,0.00002,0.00001,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000};
+    // float betaDist2[40] = {0.00148,0.00535,0.01081,0.01722,0.02404,0.03083,0.03724,0.04301,0.04794,0.05191,0.05485,0.05672,0.05756,0.05740,0.05633,0.05443,0.05183,0.04864,0.04499,0.04102,0.03683,0.03256,0.02832,0.02419,0.02028,0.01664,0.01334,0.01042,0.00789,0.00577,0.00404,0.00269,0.00168,0.00096,0.00049,0.00021,0.00007,0.00001,0.00000,0.00000};
+    // float betaDist3[40] = {0.00045,0.00169,0.00361,0.00608,0.00897,0.01219,0.01563,0.01920,0.02280,0.02637,0.02981,0.03308,0.03609,0.03882,0.04120,0.04320,0.04479,0.04594,0.04664,0.04688,0.04664,0.04594,0.04479,0.04320,0.04120,0.03882,0.03609,0.03308,0.02981,0.02637,0.02280,0.01920,0.01563,0.01219,0.00897,0.00608,0.00361,0.00169,0.00045,0.00000};
+    // float betaDist4[40] = {0.00018,0.00071,0.00156,0.00270,0.00410,0.00574,0.00758,0.00961,0.01178,0.01407,0.01646,0.01891,0.02140,0.02390,0.02638,0.02882,0.03118,0.03343,0.03556,0.03752,0.03930,0.04086,0.04218,0.04323,0.04397,0.04439,0.04445,0.04413,0.04339,0.04221,0.04057,0.03842,0.03576,0.03253,0.02873,0.02432,0.01926,0.01355,0.00713,0.00000};
+    // float single10[40] = {0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,1.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000};
+    // float single15[40] = {0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,1.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000};
+    // float single20[40] = {0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,1.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000};
+    
+    size_t nThreshold = 100;
+    int nThresholdInt = nThreshold;
+    
+    for (int i = 0; i < nThresholdInt; ++i)
+    {
+        switch (prior) {
+            case 0:
+                distribution.push_back(uniformDist[i]);
+                break;
+            case 1:
+                distribution.push_back(betaDist1[i]);
+                break;
+            case 2:
+                distribution.push_back(betaDist2[i]);
+                break;
+            case 3:
+                distribution.push_back(betaDist3[i]);
+                break;
+            case 4:
+                distribution.push_back(betaDist4[i]);
+                break;
+            case 5:
+                distribution.push_back(single10[i]);
+                break;
+            case 6:
+                distribution.push_back(single15[i]);
+                break;
+            case 7:
+                distribution.push_back(single20[i]);
+                break;
+            default:
+                distribution.push_back(uniformDist[i]);
+        }
+        thresholds.push_back(0.01 + i*0.01);
+    }
+    
+    // double minYin = 2936;
+    // for (size_t i = 2; i < yinBufferSize; ++i)
+    // {
+    //     if (yinBuffer[i] < minYin)
+    //     {
+    //         minYin = yinBuffer[i];
+    //     }
+    // }
+    // if (minYin < 0.01) std::cerr << "min Yin buffer element: " << minYin << std::endl;
+    
+    
+    int currThreshInd = nThreshold-1;
+    tau = 2;
+    
+    // double factor = 1.0 / (0.25 * (nThresholdInt+1) * (nThresholdInt + 1)); // factor to scale down triangular weight
+    size_t minInd = 0;
+    float minVal = 42.f;
+    while (currThreshInd != -1 && tau < yinBufferSize)
+    {
+        if (yinBuffer[tau] < thresholds[currThreshInd])
+        {
+            while (tau + 1 < yinBufferSize && yinBuffer[tau+1] < yinBuffer[tau])
+            {
+                tau++;
+            }
+            // tau is now local minimum
+            // std::cerr << tau << " " << currThreshInd << " "<< thresholds[currThreshInd] << " " << distribution[currThreshInd] << std::endl;
+            if (yinBuffer[tau] < minVal && tau > 2){
+                minVal = yinBuffer[tau];
+                minInd = tau;
+            }
+            peakProb[tau] += distribution[currThreshInd];
+            currThreshInd--;
+        } else {
+            tau++;
+        }
+    }
+    double nonPeakProb = 1;
+    for (size_t i = 0; i < yinBufferSize; ++i)
+    {
+        nonPeakProb -= peakProb[i];
+    }
+    // std::cerr << nonPeakProb << std::endl;
+    if (minInd > 0)
+    {
+        // std::cerr << "min set " << minVal << " " << minInd << " " << nonPeakProb << std::endl; 
+        peakProb[minInd] += nonPeakProb * minWeight;
+    }
+    
+    return peakProb;
+}
+
+double
+YinUtil::parabolicInterpolation(const double *yinBuffer, const size_t tau, const size_t yinBufferSize) 
+{
+    // this is taken almost literally from Joren Six's Java implementation
+    if (tau == yinBufferSize) // not valid anyway.
+    {
+        return static_cast<double>(tau);
+    }
+    
+    double betterTau = 0.0;
+    size_t x0;
+    size_t x2;
+
+    if (tau < 1) 
+    {
+        x0 = tau;
+    } else {
+        x0 = tau - 1;
+    }
+    
+    if (tau + 1 < yinBufferSize) 
+    {
+        x2 = tau + 1;
+    } else {
+        x2 = tau;
+    }
+    
+    if (x0 == tau) 
+    {
+        if (yinBuffer[tau] <= yinBuffer[x2]) 
+        {
+            betterTau = tau;
+        } else {
+            betterTau = x2;
+        }
+    } 
+    else if (x2 == tau) 
+    {
+        if (yinBuffer[tau] <= yinBuffer[x0]) 
+        {
+            betterTau = tau;
+        } 
+        else 
+        {
+            betterTau = x0;
+        }
+    } 
+    else 
+    {
+        float s0, s1, s2;
+        s0 = yinBuffer[x0];
+        s1 = yinBuffer[tau];
+        s2 = yinBuffer[x2];
+        // fixed AUBIO implementation, thanks to Karl Helgason:
+        // (2.0f * s1 - s2 - s0) was incorrectly multiplied with -1
+        betterTau = tau + (s2 - s0) / (2 * (2 * s1 - s2 - s0));
+        
+        // std::cerr << tau << " --> " << betterTau << std::endl;
+        
+    }
+    return betterTau;
+}
+
+double 
+YinUtil::sumSquare(const double *in, const size_t start, const size_t end)
+{
+    double out = 0;
+    for (size_t i = start; i < end; ++i)
+    {
+        out += in[i] * in[i];
+    }
+    return out;
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/YinUtil.h	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,52 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    This file is Copyright (c) 2012 Chris Cannam
+  
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
+    ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+
+#ifndef _YINUTIL_H_
+#define _YINUTIL_H_
+
+#include "vamp-sdk/FFT.h"
+#include "MeanFilter.h"
+
+#include <cmath>
+
+#include <iostream>
+#include <vector>
+#include <exception>
+
+using std::vector;
+
+class YinUtil
+{
+public:
+    static double sumSquare(const double *in, const size_t startInd, const size_t endInd);
+    static void difference(const double *in, double *yinBuffer, const size_t yinBufferSize);
+    static void fastDifference(const double *in, double *yinBuffer, const size_t yinBufferSize);
+    static void cumulativeDifference(double *yinBuffer, const size_t yinBufferSize);
+    static int absoluteThreshold(const double *yinBuffer, const size_t yinBufferSize, const double thresh);
+    static vector<double> yinProb(const double *yinBuffer, const size_t prior, const size_t yinBufferSize);
+    static double parabolicInterpolation(const double *yinBuffer, const size_t tau,
+                                         const size_t yinBufferSize);
+};
+
+#endif

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libmain.cpp	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,46 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    This file is Copyright (c) 2012 Chris Cannam
+  
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
+    ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+
+#include <vamp/vamp.h>
+#include <vamp-sdk/PluginAdapter.h>
+
+#include "PYIN.h"
+#include "VampYin.h"
+
+static Vamp::PluginAdapter<PYIN> yintonyPluginAdapter;
+static Vamp::PluginAdapter<VampYin> vampyinPluginAdapter;
+
+const VampPluginDescriptor *
+vampGetPluginDescriptor(unsigned int version, unsigned int index)
+{
+    if (version < 1) return 0;
+
+    switch (index) {
+    case  0: return yintonyPluginAdapter.getDescriptor();
+    case  1: return vampyinPluginAdapter.getDescriptor();
+    default: return 0;
+    }
+}
+
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/TestFFT.cpp	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,177 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+  This file is Copyright (c) 2012 Chris Cannam
+  
+  Permission is hereby granted, free of charge, to any person
+  obtaining a copy of this software and associated documentation
+  files (the "Software"), to deal in the Software without
+  restriction, including without limitation the rights to use, copy,
+  modify, merge, publish, distribute, sublicense, and/or sell copies
+  of the Software, and to permit persons to whom the Software is
+  furnished to do so, subject to the following conditions:
+
+  The above copyright notice and this permission notice shall be
+  included in all copies or substantial portions of the Software.
+
+  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+  NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
+  ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+  CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+  WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+
+/*
+  This unit test suite for the Vamp SDK FFT implementation is included
+  here mostly for illustrative purposes!
+*/
+
+#include "vamp-sdk/FFT.h"
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MAIN
+
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_SUITE(TestFFT)
+
+#define COMPARE_CONST(a, n) \
+    for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \
+        BOOST_CHECK_SMALL(a[cmp_i] - n, 1e-14);				\
+    }
+
+#define COMPARE_ARRAY(a, b)						\
+    for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \
+        BOOST_CHECK_SMALL(a[cmp_i] - b[cmp_i], 1e-14);			\
+    }
+
+BOOST_AUTO_TEST_CASE(dc)
+{
+    // DC-only signal. The DC bin is purely real
+    double in[] = { 1, 1, 1, 1 };
+    double re[4], im[4];
+    Vamp::FFT::forward(4, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 4.0);
+    BOOST_CHECK_EQUAL(re[1], 0.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    double backim[4];
+    Vamp::FFT::inverse(4, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+}
+
+BOOST_AUTO_TEST_CASE(sine)
+{
+    // Sine. Output is purely imaginary
+    double in[] = { 0, 1, 0, -1 };
+    double re[4], im[4];
+    Vamp::FFT::forward(4, in, 0, re, im);
+    COMPARE_CONST(re, 0.0);
+    BOOST_CHECK_EQUAL(im[0], 0.0);
+    BOOST_CHECK_EQUAL(im[1], -2.0);
+    BOOST_CHECK_EQUAL(im[2], 0.0);
+    double back[4];
+    double backim[4];
+    Vamp::FFT::inverse(4, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+}
+
+BOOST_AUTO_TEST_CASE(cosine)
+{
+    // Cosine. Output is purely real
+    double in[] = { 1, 0, -1, 0 };
+    double re[4], im[4];
+    Vamp::FFT::forward(4, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_EQUAL(re[1], 2.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    double backim[4];
+    Vamp::FFT::inverse(4, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+}
+	
+BOOST_AUTO_TEST_CASE(sineCosine)
+{
+    // Sine and cosine mixed
+    double in[] = { 0.5, 1, -0.5, -1 };
+    double re[4], im[4];
+    Vamp::FFT::forward(4, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_EQUAL(im[0], 0.0);
+    BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12);
+    BOOST_CHECK_EQUAL(im[2], 0.0);
+    double back[4];
+    double backim[4];
+    Vamp::FFT::inverse(4, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+}
+
+BOOST_AUTO_TEST_CASE(nyquist)
+{
+    double in[] = { 1, -1, 1, -1 };
+    double re[4], im[4];
+    Vamp::FFT::forward(4, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_EQUAL(re[1], 0.0);
+    BOOST_CHECK_EQUAL(re[2], 4.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    double backim[4];
+    Vamp::FFT::inverse(4, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+}
+
+BOOST_AUTO_TEST_CASE(dirac)
+{
+    double in[] = { 1, 0, 0, 0 };
+    double re[4], im[4];
+    Vamp::FFT::forward(4, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 1.0);
+    BOOST_CHECK_EQUAL(re[1], 1.0);
+    BOOST_CHECK_EQUAL(re[2], 1.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    double backim[4];
+    Vamp::FFT::inverse(4, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+}
+
+BOOST_AUTO_TEST_CASE(forwardArrayBounds)
+{
+    // initialise bins to something recognisable, so we can tell
+    // if they haven't been written
+    double in[] = { 1, 1, -1, -1 };
+    double re[] = { 999, 999, 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999, 999, 999 };
+    Vamp::FFT::forward(4, in, 0, re+1, im+1);
+    // And check we haven't overrun the arrays
+    BOOST_CHECK_EQUAL(re[0], 999.0);
+    BOOST_CHECK_EQUAL(im[0], 999.0);
+    BOOST_CHECK_EQUAL(re[5], 999.0);
+    BOOST_CHECK_EQUAL(im[5], 999.0);
+}
+
+BOOST_AUTO_TEST_CASE(inverseArrayBounds)
+{
+    // initialise bins to something recognisable, so we can tell
+    // if they haven't been written
+    double re[] = { 0, 1, 0 };
+    double im[] = { 0, -2, 0 };
+    double outre[] = { 999, 999, 999, 999, 999, 999 };
+    double outim[] = { 999, 999, 999, 999, 999, 999 };
+    Vamp::FFT::forward(4, re, im, outre+1, outim+1);
+    // And check we haven't overrun the arrays
+    BOOST_CHECK_EQUAL(outre[0], 999.0);
+    BOOST_CHECK_EQUAL(outim[0], 999.0);
+    BOOST_CHECK_EQUAL(outre[5], 999.0);
+    BOOST_CHECK_EQUAL(outim[5], 999.0);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/TestMeanFilter.cpp	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,91 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+/*
+    This file is Copyright (c) 2012 Chris Cannam
+  
+    Permission is hereby granted, free of charge, to any person
+    obtaining a copy of this software and associated documentation
+    files (the "Software"), to deal in the Software without
+    restriction, including without limitation the rights to use, copy,
+    modify, merge, publish, distribute, sublicense, and/or sell copies
+    of the Software, and to permit persons to whom the Software is
+    furnished to do so, subject to the following conditions:
+
+    The above copyright notice and this permission notice shall be
+    included in all copies or substantial portions of the Software.
+
+    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
+    ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
+    CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+    WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+
+#include "MeanFilter.h"
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MAIN
+
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_SUITE(TestMeanFilter)
+
+BOOST_AUTO_TEST_CASE(simpleFilter)
+{
+    double in[] = { 1.0, 2.0, 3.0 };
+    double out[3];
+    MeanFilter(3).filter(in, out, 3);
+    BOOST_CHECK_EQUAL(out[0], 1.5);
+    BOOST_CHECK_EQUAL(out[1], 2.0);
+    BOOST_CHECK_EQUAL(out[2], 2.5);
+}
+
+BOOST_AUTO_TEST_CASE(simpleSubset)
+{
+    double in[] = { 1.0, 2.0, 3.0, 4.0, 5.0 };
+    double out[3];
+    MeanFilter(3).filterSubsequence(in, out, 5, 3, 1);
+    BOOST_CHECK_EQUAL(out[0], 2.0);
+    BOOST_CHECK_EQUAL(out[1], 3.0);
+    BOOST_CHECK_EQUAL(out[2], 4.0);
+}
+
+BOOST_AUTO_TEST_CASE(flenExceedsArraySize)
+{
+    double in[] = { 1.0, 2.0, 3.0 };
+    double out[3];
+    MeanFilter(5).filter(in, out, 3);
+    BOOST_CHECK_EQUAL(out[0], 2.0);
+    BOOST_CHECK_EQUAL(out[1], 2.0);
+    BOOST_CHECK_EQUAL(out[2], 2.0);
+}
+
+BOOST_AUTO_TEST_CASE(flenIs1)
+{
+    double in[] = { 1.0, 2.0, 3.0 };
+    double out[3];
+    MeanFilter(1).filter(in, out, 3);
+    BOOST_CHECK_EQUAL(out[0], in[0]);
+    BOOST_CHECK_EQUAL(out[1], in[1]);
+    BOOST_CHECK_EQUAL(out[2], in[2]);
+}
+
+BOOST_AUTO_TEST_CASE(arraySizeIs1)
+{
+    double in[] = { 1.0 };
+    double out[1];
+    MeanFilter(3).filter(in, out, 1);
+    BOOST_CHECK_EQUAL(out[0], in[0]);
+}
+
+BOOST_AUTO_TEST_CASE(subsequenceLengthIs1)
+{
+    double in[] = { 1.0, 2.0, 3.0 };
+    double out[1];
+    MeanFilter(3).filterSubsequence(in, out, 3, 1, 2);
+    BOOST_CHECK_EQUAL(out[0], 2.5);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/TestMonoNote.cpp	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,43 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#include "MonoNote.h"
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MAIN
+
+#include <boost/test/unit_test.hpp>
+
+using std::vector;
+using std::pair;
+
+BOOST_AUTO_TEST_SUITE(TestMonoNote)
+
+BOOST_AUTO_TEST_CASE(instantiate)
+{
+    MonoNote mn;
+    vector<vector<pair<double, double> > > pitchProb;
+    size_t n = 8;
+
+    vector<double> pitch;
+    pitch.push_back(48);
+    pitch.push_back(51);
+    pitch.push_back(50);
+    pitch.push_back(51);
+    pitch.push_back(60);
+    pitch.push_back(71);
+    pitch.push_back(51);
+    pitch.push_back(62);
+    
+    for (size_t i = 0; i < n; ++i)
+    {
+        vector<pair<double, double> > temp;
+        temp.push_back(pair<double, double>(pitch[i], 1.0));
+        pitchProb.push_back(temp);
+    }
+
+    mn.process(pitchProb);
+}
+
+
+BOOST_AUTO_TEST_SUITE_END()
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/TestPredominoMethods.cpp	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,132 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#include "PredominoMethods.h"
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MAIN
+
+#include <boost/test/unit_test.hpp>
+#include <boost/test/floating_point_comparison.hpp>
+
+BOOST_AUTO_TEST_SUITE(TestPredominoMethods)
+
+BOOST_AUTO_TEST_CASE(pickNPeaks)
+{
+    vector<float> invec;
+    invec.push_back(0);
+    invec.push_back(1);
+    invec.push_back(2);
+    invec.push_back(1);
+    invec.push_back(3);
+    invec.push_back(3);
+    invec.push_back(4);
+    invec.push_back(0);
+    invec.push_back(1);
+    invec.push_back(0);
+    
+    vector<size_t> outvec = PredominoMethods::pickNPeaks(invec, 2);
+    BOOST_CHECK_EQUAL(outvec[0], 2);
+    BOOST_CHECK_EQUAL(outvec[1], 6);
+    
+}
+
+BOOST_AUTO_TEST_CASE(zeroPad)
+{
+    size_t nIn = 4;
+    size_t nPadded = 8;
+    float *in = new float[nIn];
+    in[0] = 1; in[1] = 2; in[2] = 3; in[3] = 4;
+    float *out = new float[nPadded];
+    PredominoMethods::zeroPad(in, out, nIn, nPadded);
+    
+    BOOST_CHECK_EQUAL(out[0], 3);
+    BOOST_CHECK_EQUAL(out[1], 4);
+    for (size_t i = 2; i < nPadded - nIn/2; ++i) {
+        BOOST_CHECK_EQUAL(out[i], 0);
+    }
+    BOOST_CHECK_EQUAL(out[6], 1);
+    BOOST_CHECK_EQUAL(out[7], 2);
+    
+    delete [] in;
+    delete [] out;
+}
+
+BOOST_AUTO_TEST_CASE(applyWindow)
+{
+    size_t n = 4;
+    vector<float> in;
+    vector<float> win;
+    for (size_t i = 0; i < n; ++i) {
+        in.push_back(i);
+        win.push_back(i);
+    }
+    PredominoMethods::applyWindow(&in, win);
+    for (size_t i = 0; i < n; ++i) {
+        BOOST_CHECK_EQUAL(in[i], i*i);
+    }
+}
+
+BOOST_AUTO_TEST_CASE(hannWindow4)
+{
+    size_t n = 4;
+    vector<float> hw = PredominoMethods::hannWindow(n);
+    // for (size_t i = 0; i < n; ++i) {
+    //     std::cerr << hw[i] << std::endl;
+    // }
+
+    BOOST_CHECK_CLOSE_FRACTION(hw[0], 0.1382, .0001);
+    BOOST_CHECK_CLOSE_FRACTION(hw[1], 0.3618, .0001);
+    BOOST_CHECK_CLOSE_FRACTION(hw[2], 0.3618, .0001);
+    BOOST_CHECK_CLOSE_FRACTION(hw[3], 0.1382, .0001);
+}
+
+BOOST_AUTO_TEST_CASE(parabolicInterpolation)
+{
+    vector<float> mag1;
+    mag1.push_back(1); mag1.push_back(2); mag1.push_back(1);
+    vector<size_t> peakIndices1;
+    peakIndices1.push_back(1);
+    vector<float> interpolated1 = PredominoMethods::parabolicInterpolation(mag1, peakIndices1);
+    BOOST_CHECK_CLOSE_FRACTION(interpolated1[0], 0, 0.001);
+
+    vector<float> mag2;
+    mag2.push_back(1); mag2.push_back(1); mag1.push_back(0);
+    vector<size_t> peakIndices2;
+    peakIndices2.push_back(1);
+    vector<float> interpolated2 = PredominoMethods::parabolicInterpolation(mag2, peakIndices2);
+    BOOST_CHECK_CLOSE_FRACTION(interpolated2[0], -0.5, 0.001);
+}
+
+BOOST_AUTO_TEST_CASE(convertMagnitude2dB)
+{
+    vector<float> mag;
+    mag.push_back(1.01); 
+    vector<float> magdb = PredominoMethods::convertMagnitude2dB(mag);
+    BOOST_CHECK_CLOSE_FRACTION(magdb[0], 0.0864, 0.001);
+}
+
+BOOST_AUTO_TEST_CASE(generateOutputFrequencies)
+{
+    vector<float> f = PredominoMethods::generateOutputFrequencies(110, 2, 50);
+    BOOST_CHECK_CLOSE_FRACTION(f[0], 110, 0.001);
+    BOOST_CHECK_CLOSE_FRACTION(f[1], 113.2232, 0.001);
+}
+
+BOOST_AUTO_TEST_CASE(timedomain2magnitude)
+{
+    float *in = new float[4];
+    in[0] = 1;
+    in[1] = 2;
+    in[2] = 3;
+    in[3] = 2;
+    vector<float> mag = PredominoMethods::timedomain2magnitude(in, 4, 0);
+    // for (size_t i = 0; i < 4; ++i) std::cerr << mag[i] << std::endl;
+    delete [] in;
+    BOOST_CHECK_CLOSE_FRACTION(mag[0], 8, 0.001);
+    BOOST_CHECK_CLOSE_FRACTION(mag[1], 2, 0.001);
+    BOOST_CHECK_CLOSE_FRACTION(mag[2], 0, 0.001);
+    BOOST_CHECK_CLOSE_FRACTION(mag[3], 2, 0.001);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/TestYin.cpp	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,41 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#include "Yin.h"
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MAIN
+
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_SUITE(TestYin)
+
+BOOST_AUTO_TEST_CASE(implicitYinOutput)
+{
+    // this test is just to make sure a YinOutput initialises
+    // -- the actual reason for me to write this is to have written 
+    // a test, so let's plough on...
+    Yin::YinOutput out;
+    
+    BOOST_CHECK_EQUAL(out.f0, 0);
+    BOOST_CHECK_EQUAL(out.periodicity, 0);
+    BOOST_CHECK_EQUAL(out.rms, 0);
+}
+
+BOOST_AUTO_TEST_CASE(sine128)
+{
+    
+    // a very short frame (8) with maximum frequency (128),
+    // assuming a sampling rate of 256 1/sec.
+    // Yin should get it approximately right... 
+    // but because of post-processing we want to be tolerant
+
+    double in[] = { 1, -1, 1, -1, 1, -1, 1, -1 };
+    Yin y(8, 256);
+    Yin::YinOutput yo = y.process(in);
+    
+    BOOST_CHECK(yo.f0 > 114);
+    BOOST_CHECK(yo.f0 < 142);
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/TestYinUtil.cpp	Wed Nov 27 11:59:49 2013 +0000
@@ -0,0 +1,94 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#include "YinUtil.h"
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MAIN
+
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_SUITE(TestYin)
+
+BOOST_AUTO_TEST_CASE(difference1)
+{
+    // difference always needs to be non-negative
+    double in[] = { 1, -1, 1, -1, 1, -1, 1, -1 };
+    double yinBuffer[4];
+    
+    YinUtil::difference(in, yinBuffer, 4);
+    
+    BOOST_CHECK(yinBuffer[0] >= 0);
+    BOOST_CHECK(yinBuffer[1] >= 0);
+    BOOST_CHECK(yinBuffer[2] >= 0);
+    BOOST_CHECK(yinBuffer[3] >= 0);
+}
+
+// BOOST_AUTO_TEST_CASE(difference2)
+// {
+//     // difference for fast sinusoid
+//     double in[] = { 1, -1, 1, -1, 1, -1, 1, -1 };
+//     double yinBuffer[4];
+//     
+//     Yin y(8, 256);
+//     
+//     y.difference(in, yinBuffer);
+//     
+//     BOOST_CHECK_EQUAL(yinBuffer[0], 0);
+//     BOOST_CHECK_EQUAL(yinBuffer[1], 16);
+//     BOOST_CHECK_EQUAL(yinBuffer[2], 0);
+//     BOOST_CHECK_EQUAL(yinBuffer[3], 16);
+// }
+// 
+// BOOST_AUTO_TEST_CASE(difference3)
+// {
+//     // difference for contaminated fast sinusoid
+//     double in[] = { 2, -1, 1, -1, 1, -1, 1, -1 };
+//     double yinBuffer[4];
+//     
+//     Yin y(8, 256);
+//     
+//     y.difference(in, yinBuffer);
+//     
+//     BOOST_CHECK_EQUAL(yinBuffer[0], 0);
+//     BOOST_CHECK_EQUAL(yinBuffer[1], 21);
+//     BOOST_CHECK_EQUAL(yinBuffer[2], 1);
+//     BOOST_CHECK_EQUAL(yinBuffer[3], 21);
+// }
+// 
+// BOOST_AUTO_TEST_CASE(difference4)
+// {
+//     // 
+//     double in[2048];
+//     double yinBuffer[1024];
+//     for (size_t i = 0; i < 2048; i = i + 2) {
+//         in[i] = 1;
+//         in[i+1] = -1;
+//     }
+//     Yin y(2048, 44100);
+//     
+//     y.difference(in, yinBuffer);
+//     
+//     BOOST_CHECK(yinBuffer[0] >= 0);
+//     BOOST_CHECK(yinBuffer[1] >= 0);
+// }
+// 
+// BOOST_AUTO_TEST_CASE(cumulativeDifference1)
+// {
+//     // test against matlab implementation
+//     double yinBuffer[] = {0, 21, 1, 21};
+//     
+//     Yin y(8, 256);
+//     
+//     y.cumulativeDifference(yinBuffer);
+//     
+//     BOOST_CHECK_EQUAL(yinBuffer[0], 1);
+//     BOOST_CHECK_EQUAL(yinBuffer[1], 1); // this is just chance!
+//     BOOST_CHECK(yinBuffer[2] >= 0.0909 && yinBuffer[2] < 0.091);
+//     BOOST_CHECK(yinBuffer[3] >= 1.4651 && yinBuffer[3] < 1.466);
+// 
+// }
+
+
+
+BOOST_AUTO_TEST_SUITE_END()
+