Chris@3: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@3: /*
Chris@3:     Permission is hereby granted, free of charge, to any person
Chris@3:     obtaining a copy of this software and associated documentation
Chris@3:     files (the "Software"), to deal in the Software without
Chris@3:     restriction, including without limitation the rights to use, copy,
Chris@3:     modify, merge, publish, distribute, sublicense, and/or sell copies
Chris@3:     of the Software, and to permit persons to whom the Software is
Chris@3:     furnished to do so, subject to the following conditions:
Chris@3: 
Chris@3:     The above copyright notice and this permission notice shall be
Chris@3:     included in all copies or substantial portions of the Software.
Chris@3: 
Chris@3:     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
Chris@3:     EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
Chris@3:     MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
Chris@3:     NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
Chris@3:     ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
Chris@3:     CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
Chris@3:     WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
Chris@3: */
Chris@3: 
Chris@3: #include "CepstrumPitchTracker.h"
Chris@3: 
Chris@3: #include <vector>
Chris@3: #include <algorithm>
Chris@3: 
Chris@3: #include <cstdio>
Chris@3: #include <cmath>
Chris@3: #include <complex>
Chris@3: 
Chris@3: using std::string;
Chris@7: using std::vector;
Chris@16: using Vamp::RealTime;
Chris@7: 
Chris@8: CepstrumPitchTracker::Hypothesis::Hypothesis()
Chris@7: {
Chris@8:     m_state = New;
Chris@7:     m_age = 0;
Chris@7: }
Chris@7: 
Chris@11: CepstrumPitchTracker::Hypothesis::~Hypothesis()
Chris@11: {
Chris@11: }
Chris@11: 
Chris@7: bool
Chris@7: CepstrumPitchTracker::Hypothesis::isWithinTolerance(Estimate s)
Chris@7: {
Chris@7:     if (m_pending.empty()) {
Chris@7:         return true;
Chris@7:     }
Chris@16: 
Chris@16:     // check we are within a relatively close tolerance of the last
Chris@16:     // candidate
Chris@17:     Estimate last = m_pending[m_pending.size()-1];
Chris@7:     double r = s.freq / last.freq;
Chris@7:     int cents = lrint(1200.0 * (log(r) / log(2.0)));
Chris@17:     if (cents < -60 || cents > 60) return false;
Chris@16: 
Chris@17:     // and within a slightly bigger tolerance of the current mean
Chris@17:     double meanFreq = getMeanFrequency();
Chris@17:     r = s.freq / meanFreq;
Chris@16:     cents = lrint(1200.0 * (log(r) / log(2.0)));
Chris@16:     if (cents < -80 || cents > 80) return false;
Chris@16:     
Chris@16:     return true;
Chris@7: }
Chris@7: 
Chris@7: bool 
Chris@7: CepstrumPitchTracker::Hypothesis::isSatisfied()
Chris@7: {
Chris@15:     if (m_pending.empty()) return false;
Chris@15:     
Chris@15:     double meanConfidence = 0.0;
Chris@15:     for (int i = 0; i < m_pending.size(); ++i) {
Chris@15:         meanConfidence += m_pending[i].confidence;
Chris@15:     }
Chris@15:     meanConfidence /= m_pending.size();
Chris@15: 
Chris@15:     int lengthRequired = int(2.0 / meanConfidence + 0.5);
Chris@15:     std::cerr << "meanConfidence = " << meanConfidence << ", lengthRequired = " << lengthRequired << ", length = " << m_pending.size() << std::endl;
Chris@15: 
Chris@15:     return (m_pending.size() > lengthRequired);
Chris@7: }
Chris@7: 
Chris@8: void
Chris@8: CepstrumPitchTracker::Hypothesis::advanceTime()
Chris@8: {
Chris@8:     ++m_age;
Chris@8: }
Chris@8: 
Chris@7: bool
Chris@7: CepstrumPitchTracker::Hypothesis::test(Estimate s)
Chris@7: {
Chris@8:     bool accept = false;
Chris@8: 
Chris@8:     switch (m_state) {
Chris@8: 
Chris@8:     case New:
Chris@8:         m_state = Provisional;
Chris@8:         accept = true;
Chris@8:         break;
Chris@8: 
Chris@8:     case Provisional:
Chris@8:         if (m_age > 3) {
Chris@8:             m_state = Rejected;
Chris@8:         } else if (isWithinTolerance(s)) {
Chris@8:             accept = true;
Chris@8:         }
Chris@8:         break;
Chris@8:         
Chris@8:     case Satisfied:
Chris@8:         if (m_age > 3) {
Chris@8:             m_state = Expired;
Chris@8:         } else if (isWithinTolerance(s)) {
Chris@8:             accept = true;
Chris@8:         }
Chris@8:         break;
Chris@8: 
Chris@8:     case Rejected:
Chris@8:         break;
Chris@8: 
Chris@8:     case Expired:
Chris@8:         break;
Chris@7:     }
Chris@7: 
Chris@8:     if (accept) {
Chris@8:         m_pending.push_back(s);
Chris@8:         m_age = 0;
Chris@8:         if (m_state == Provisional && isSatisfied()) {
Chris@8:             m_state = Satisfied;
Chris@7:         }
Chris@7:     }
Chris@7: 
Chris@14:     return accept && (m_state == Satisfied);
Chris@8: }        
Chris@7: 
Chris@7: CepstrumPitchTracker::Hypothesis::State
Chris@7: CepstrumPitchTracker::Hypothesis::getState()
Chris@7: {
Chris@7:     return m_state;
Chris@7: }
Chris@7: 
Chris@12: int
Chris@12: CepstrumPitchTracker::Hypothesis::getPendingLength()
Chris@12: {
Chris@12:     return m_pending.size();
Chris@12: }
Chris@12: 
Chris@7: CepstrumPitchTracker::Hypothesis::Estimates
Chris@7: CepstrumPitchTracker::Hypothesis::getAcceptedEstimates()
Chris@7: {
Chris@7:     if (m_state == Satisfied || m_state == Expired) {
Chris@7:         return m_pending;
Chris@7:     } else {
Chris@7:         return Estimates();
Chris@7:     }
Chris@7: }
Chris@7: 
Chris@17: double
Chris@17: CepstrumPitchTracker::Hypothesis::getMeanFrequency()
Chris@17: {
Chris@17:     double acc = 0.0;
Chris@17:     for (int i = 0; i < m_pending.size(); ++i) {
Chris@17:         acc += m_pending[i].freq;
Chris@17:     }
Chris@17:     acc /= m_pending.size();
Chris@17:     return acc;
Chris@17: }
Chris@17: 
Chris@16: CepstrumPitchTracker::Hypothesis::Note
Chris@16: CepstrumPitchTracker::Hypothesis::getAveragedNote()
Chris@16: {
Chris@16:     Note n;
Chris@16: 
Chris@16:     if (!(m_state == Satisfied || m_state == Expired)) {
Chris@16:         n.freq = 0.0;
Chris@16:         n.time = RealTime::zeroTime;
Chris@16:         n.duration = RealTime::zeroTime;
Chris@16:         return n;
Chris@16:     }
Chris@16: 
Chris@16:     n.time = m_pending.begin()->time;
Chris@16: 
Chris@16:     Estimates::iterator i = m_pending.end();
Chris@16:     --i;
Chris@16:     n.duration = i->time - n.time;
Chris@16: 
Chris@17:     // just mean frequency for now, but this isn't at all right perceptually
Chris@17:     n.freq = getMeanFrequency();
Chris@16:     
Chris@16:     return n;
Chris@16: }
Chris@16: 
Chris@11: void
Chris@16: CepstrumPitchTracker::Hypothesis::addFeatures(FeatureSet &fs)
Chris@11: {
Chris@11:     for (int i = 0; i < m_pending.size(); ++i) {
Chris@11: 	Feature f;
Chris@11: 	f.hasTimestamp = true;
Chris@11: 	f.timestamp = m_pending[i].time;
Chris@11: 	f.values.push_back(m_pending[i].freq);
Chris@16: 	fs[0].push_back(f);
Chris@11:     }
Chris@16: 
Chris@16:     Feature nf;
Chris@16:     nf.hasTimestamp = true;
Chris@16:     nf.hasDuration = true;
Chris@16:     Note n = getAveragedNote();
Chris@16:     nf.timestamp = n.time;
Chris@16:     nf.duration = n.duration;
Chris@16:     nf.values.push_back(n.freq);
Chris@16:     fs[1].push_back(nf);
Chris@11: }
Chris@3: 
Chris@3: CepstrumPitchTracker::CepstrumPitchTracker(float inputSampleRate) :
Chris@3:     Plugin(inputSampleRate),
Chris@3:     m_channels(0),
Chris@3:     m_stepSize(256),
Chris@3:     m_blockSize(1024),
Chris@3:     m_fmin(50),
Chris@3:     m_fmax(1000),
Chris@18:     m_vflen(1),
Chris@3:     m_binFrom(0),
Chris@3:     m_binTo(0),
Chris@15:     m_bins(0)
Chris@3: {
Chris@3: }
Chris@3: 
Chris@3: CepstrumPitchTracker::~CepstrumPitchTracker()
Chris@3: {
Chris@3: }
Chris@3: 
Chris@3: string
Chris@3: CepstrumPitchTracker::getIdentifier() const
Chris@3: {
Chris@3:     return "cepstrum-pitch";
Chris@3: }
Chris@3: 
Chris@3: string
Chris@3: CepstrumPitchTracker::getName() const
Chris@3: {
Chris@3:     return "Cepstrum Pitch Tracker";
Chris@3: }
Chris@3: 
Chris@3: string
Chris@3: CepstrumPitchTracker::getDescription() const
Chris@3: {
Chris@3:     return "Estimate f0 of monophonic material using a cepstrum method.";
Chris@3: }
Chris@3: 
Chris@3: string
Chris@3: CepstrumPitchTracker::getMaker() const
Chris@3: {
Chris@3:     return "Chris Cannam";
Chris@3: }
Chris@3: 
Chris@3: int
Chris@3: CepstrumPitchTracker::getPluginVersion() const
Chris@3: {
Chris@3:     // Increment this each time you release a version that behaves
Chris@3:     // differently from the previous one
Chris@3:     return 1;
Chris@3: }
Chris@3: 
Chris@3: string
Chris@3: CepstrumPitchTracker::getCopyright() const
Chris@3: {
Chris@3:     return "Freely redistributable (BSD license)";
Chris@3: }
Chris@3: 
Chris@3: CepstrumPitchTracker::InputDomain
Chris@3: CepstrumPitchTracker::getInputDomain() const
Chris@3: {
Chris@3:     return FrequencyDomain;
Chris@3: }
Chris@3: 
Chris@3: size_t
Chris@3: CepstrumPitchTracker::getPreferredBlockSize() const
Chris@3: {
Chris@3:     return 1024;
Chris@3: }
Chris@3: 
Chris@3: size_t 
Chris@3: CepstrumPitchTracker::getPreferredStepSize() const
Chris@3: {
Chris@3:     return 256;
Chris@3: }
Chris@3: 
Chris@3: size_t
Chris@3: CepstrumPitchTracker::getMinChannelCount() const
Chris@3: {
Chris@3:     return 1;
Chris@3: }
Chris@3: 
Chris@3: size_t
Chris@3: CepstrumPitchTracker::getMaxChannelCount() const
Chris@3: {
Chris@3:     return 1;
Chris@3: }
Chris@3: 
Chris@3: CepstrumPitchTracker::ParameterList
Chris@3: CepstrumPitchTracker::getParameterDescriptors() const
Chris@3: {
Chris@3:     ParameterList list;
Chris@3:     return list;
Chris@3: }
Chris@3: 
Chris@3: float
Chris@3: CepstrumPitchTracker::getParameter(string identifier) const
Chris@3: {
Chris@3:     return 0.f;
Chris@3: }
Chris@3: 
Chris@3: void
Chris@3: CepstrumPitchTracker::setParameter(string identifier, float value) 
Chris@3: {
Chris@3: }
Chris@3: 
Chris@3: CepstrumPitchTracker::ProgramList
Chris@3: CepstrumPitchTracker::getPrograms() const
Chris@3: {
Chris@3:     ProgramList list;
Chris@3:     return list;
Chris@3: }
Chris@3: 
Chris@3: string
Chris@3: CepstrumPitchTracker::getCurrentProgram() const
Chris@3: {
Chris@3:     return ""; // no programs
Chris@3: }
Chris@3: 
Chris@3: void
Chris@3: CepstrumPitchTracker::selectProgram(string name)
Chris@3: {
Chris@3: }
Chris@3: 
Chris@3: CepstrumPitchTracker::OutputList
Chris@3: CepstrumPitchTracker::getOutputDescriptors() const
Chris@3: {
Chris@3:     OutputList outputs;
Chris@3: 
Chris@3:     int n = 0;
Chris@3: 
Chris@3:     OutputDescriptor d;
Chris@3: 
Chris@3:     d.identifier = "f0";
Chris@3:     d.name = "Estimated f0";
Chris@3:     d.description = "Estimated fundamental frequency";
Chris@3:     d.unit = "Hz";
Chris@3:     d.hasFixedBinCount = true;
Chris@3:     d.binCount = 1;
Chris@3:     d.hasKnownExtents = true;
Chris@3:     d.minValue = m_fmin;
Chris@3:     d.maxValue = m_fmax;
Chris@3:     d.isQuantized = false;
Chris@3:     d.sampleType = OutputDescriptor::FixedSampleRate;
Chris@3:     d.sampleRate = (m_inputSampleRate / m_stepSize);
Chris@3:     d.hasDuration = false;
Chris@3:     outputs.push_back(d);
Chris@3: 
Chris@16:     d.identifier = "notes";
Chris@16:     d.name = "Notes";
Chris@16:     d.description = "Derived fixed-pitch note frequencies";
Chris@16:     d.unit = "Hz";
Chris@16:     d.hasFixedBinCount = true;
Chris@16:     d.binCount = 1;
Chris@16:     d.hasKnownExtents = true;
Chris@16:     d.minValue = m_fmin;
Chris@16:     d.maxValue = m_fmax;
Chris@16:     d.isQuantized = false;
Chris@16:     d.sampleType = OutputDescriptor::FixedSampleRate;
Chris@16:     d.sampleRate = (m_inputSampleRate / m_stepSize);
Chris@16:     d.hasDuration = true;
Chris@16:     outputs.push_back(d);
Chris@16: 
Chris@3:     return outputs;
Chris@3: }
Chris@3: 
Chris@3: bool
Chris@3: CepstrumPitchTracker::initialise(size_t channels, size_t stepSize, size_t blockSize)
Chris@3: {
Chris@3:     if (channels < getMinChannelCount() ||
Chris@3: 	channels > getMaxChannelCount()) return false;
Chris@3: 
Chris@3: //    std::cerr << "CepstrumPitchTracker::initialise: channels = " << channels
Chris@3: //	      << ", stepSize = " << stepSize << ", blockSize = " << blockSize
Chris@3: //	      << std::endl;
Chris@3: 
Chris@3:     m_channels = channels;
Chris@3:     m_stepSize = stepSize;
Chris@3:     m_blockSize = blockSize;
Chris@3: 
Chris@3:     m_binFrom = int(m_inputSampleRate / m_fmax);
Chris@3:     m_binTo = int(m_inputSampleRate / m_fmin); 
Chris@3: 
Chris@3:     if (m_binTo >= (int)m_blockSize / 2) {
Chris@3:         m_binTo = m_blockSize / 2 - 1;
Chris@3:     }
Chris@3: 
Chris@3:     m_bins = (m_binTo - m_binFrom) + 1;
Chris@3: 
Chris@3:     reset();
Chris@3: 
Chris@3:     return true;
Chris@3: }
Chris@3: 
Chris@3: void
Chris@3: CepstrumPitchTracker::reset()
Chris@3: {
Chris@3: }
Chris@3: 
Chris@3: void
Chris@15: CepstrumPitchTracker::filter(const double *cep, double *data)
Chris@3: {
Chris@3:     for (int i = 0; i < m_bins; ++i) {
Chris@5:         double v = 0;
Chris@5:         int n = 0;
Chris@5:         // average according to the vertical filter length
Chris@5:         for (int j = -m_vflen/2; j <= m_vflen/2; ++j) {
Chris@5:             int ix = i + m_binFrom + j;
Chris@5:             if (ix >= 0 && ix < m_blockSize) {
Chris@5:                 v += cep[ix];
Chris@5:                 ++n;
Chris@5:             }
Chris@5:         }
Chris@15:         data[i] = v / n;
Chris@3:     }
Chris@6: }
Chris@6: 
Chris@18: double
Chris@18: CepstrumPitchTracker::cubicInterpolate(const double y[4], double x)
Chris@18: {
Chris@18:     double a0 = y[3] - y[2] - y[0] + y[1];
Chris@18:     double a1 = y[0] - y[1] - a0;
Chris@18:     double a2 = y[2] - y[0];
Chris@18:     double a3 = y[1];
Chris@18:     return
Chris@18:         a0 * x * x * x +
Chris@18:         a1 * x * x +
Chris@18:         a2 * x +
Chris@18:         a3;
Chris@18: }
Chris@18: 
Chris@18: double
Chris@18: CepstrumPitchTracker::findInterpolatedPeak(const double *in, int maxbin)
Chris@18: {
Chris@18:     if (maxbin < 2 || maxbin > m_bins - 3) {
Chris@18:         return maxbin;
Chris@18:     }
Chris@18: 
Chris@18:     double maxval = 0.0;
Chris@18:     double maxidx = maxbin;
Chris@18: 
Chris@18:     const int divisions = 10;
Chris@18:     double y[4];
Chris@18: 
Chris@18:     y[0] = in[maxbin-1];
Chris@18:     y[1] = in[maxbin];
Chris@18:     y[2] = in[maxbin+1];
Chris@18:     y[3] = in[maxbin+2];
Chris@18:     for (int i = 0; i < divisions; ++i) {
Chris@18:         double probe = double(i) / double(divisions);
Chris@18:         double value = cubicInterpolate(y, probe);
Chris@18:         if (value > maxval) {
Chris@18:             maxval = value; 
Chris@18:             maxidx = maxbin + probe;
Chris@18:         }
Chris@18:     }
Chris@18: 
Chris@18:     y[3] = y[2];
Chris@18:     y[2] = y[1];
Chris@18:     y[1] = y[0];
Chris@18:     y[0] = in[maxbin-2];
Chris@18:     for (int i = 0; i < divisions; ++i) {
Chris@18:         double probe = double(i) / double(divisions);
Chris@18:         double value = cubicInterpolate(y, probe);
Chris@18:         if (value > maxval) {
Chris@18:             maxval = value; 
Chris@18:             maxidx = maxbin - 1 + probe;
Chris@18:         }
Chris@18:     }
Chris@18: 
Chris@18: /*
Chris@18:     std::cerr << "centre = " << maxbin << ": ["
Chris@18:               << in[maxbin-2] << ","
Chris@18:               << in[maxbin-1] << ","
Chris@18:               << in[maxbin] << ","
Chris@18:               << in[maxbin+1] << ","
Chris@18:               << in[maxbin+2] << "] -> " << maxidx << std::endl;
Chris@18: */
Chris@18: 
Chris@18:     return maxidx;
Chris@18: }
Chris@18: 
Chris@3: CepstrumPitchTracker::FeatureSet
Chris@16: CepstrumPitchTracker::process(const float *const *inputBuffers, RealTime timestamp)
Chris@3: {
Chris@3:     FeatureSet fs;
Chris@3: 
Chris@3:     int bs = m_blockSize;
Chris@3:     int hs = m_blockSize/2 + 1;
Chris@3: 
Chris@3:     double *rawcep = new double[bs];
Chris@3:     double *io = new double[bs];
Chris@3:     double *logmag = new double[bs];
Chris@3: 
Chris@4:     // The "inverse symmetric" method. Seems to be the most reliable
Chris@3:         
Chris@3:     for (int i = 0; i < hs; ++i) {
Chris@3: 
Chris@3: 	double power =
Chris@3: 	    inputBuffers[0][i*2  ] * inputBuffers[0][i*2  ] +
Chris@3: 	    inputBuffers[0][i*2+1] * inputBuffers[0][i*2+1];
Chris@3: 	double mag = sqrt(power);
Chris@3: 	
Chris@3: 	double lm = log(mag + 0.00000001);
Chris@3: 	
Chris@4: 	logmag[i] = lm;
Chris@4: 	if (i > 0) logmag[bs - i] = lm;
Chris@3:     }
Chris@3: 
Chris@4:     fft(bs, true, logmag, 0, rawcep, io);
Chris@3:     
Chris@3:     delete[] logmag;
Chris@3:     delete[] io;
Chris@3: 
Chris@22:     double cep1 = rawcep[1];
Chris@22: 
Chris@3:     int n = m_bins;
Chris@3:     double *data = new double[n];
Chris@3:     filter(rawcep, data);
Chris@3:     delete[] rawcep;
Chris@3: 
Chris@3:     double maxval = 0.0;
Chris@6:     int maxbin = -1;
Chris@3: 
Chris@3:     for (int i = 0; i < n; ++i) {
Chris@3:         if (data[i] > maxval) {
Chris@3:             maxval = data[i];
Chris@3:             maxbin = i;
Chris@3:         }
Chris@3:     }
Chris@3: 
Chris@15:     if (maxbin < 0) {
Chris@15:         delete[] data;
Chris@15:         return fs;
Chris@15:     }
Chris@15: 
Chris@15:     double nextPeakVal = 0.0;
Chris@15:     for (int i = 1; i+1 < n; ++i) {
Chris@15:         if (data[i] > data[i-1] &&
Chris@15:             data[i] > data[i+1] &&
Chris@15:             i != maxbin &&
Chris@15:             data[i] > nextPeakVal) {
Chris@15:             nextPeakVal = data[i];
Chris@15:         }
Chris@15:     }
Chris@8: 
Chris@18:     double cimax = findInterpolatedPeak(data, maxbin);
Chris@18:     double peakfreq = m_inputSampleRate / (cimax + m_binFrom);
Chris@15: 
Chris@15:     double confidence = 0.0;
Chris@15:     if (nextPeakVal != 0.0) {
Chris@22:         std::cerr << "maxval = " << maxval << ", cep1 = " << cep1 << std::endl;
Chris@22:         double conf0 = (maxval - nextPeakVal) / 80.0;
Chris@22:         double conf1 = (cep1 / bs) / 2;
Chris@22:         if (conf0 > 1.0) conf0 = 1.0;
Chris@22:         if (conf1 > 1.0) conf1 = 1.0;
Chris@22:         confidence = conf0 * conf1;
Chris@22:         std::cerr << "conf0 = " << conf0 << ", conf1 = " << conf1 << ", confidence = " << confidence << std::endl;
Chris@15:     }
Chris@15: 
Chris@8:     Hypothesis::Estimate e;
Chris@8:     e.freq = peakfreq;
Chris@8:     e.time = timestamp;
Chris@15:     e.confidence = confidence;
Chris@8: 
Chris@8:     m_accepted.advanceTime();
Chris@11: 
Chris@8:     for (int i = 0; i < m_possible.size(); ++i) {
Chris@8:         m_possible[i].advanceTime();
Chris@8:     }
Chris@8: 
Chris@11:     if (!m_accepted.test(e)) {
Chris@13: 
Chris@11:         int candidate = -1;
Chris@13:         bool accepted = false;
Chris@13: 
Chris@11:         for (int i = 0; i < m_possible.size(); ++i) {
Chris@11:             if (m_possible[i].test(e)) {
Chris@13:                 accepted = true;
Chris@11:                 if (m_possible[i].getState() == Hypothesis::Satisfied) {
Chris@11:                     candidate = i;
Chris@11:                 }
Chris@11:                 break;
Chris@11:             }
Chris@11:         }
Chris@12: 
Chris@13:         if (!accepted) {
Chris@13:             Hypothesis h;
Chris@13:             h.test(e); //!!! must succeed as h is new, so perhaps there should be a ctor for this
Chris@13:             m_possible.push_back(h);
Chris@13:         }
Chris@13: 
Chris@11:         if (m_accepted.getState() == Hypothesis::Expired) {
Chris@16:             m_accepted.addFeatures(fs);
Chris@12:         }
Chris@12:         
Chris@12:         if (m_accepted.getState() == Hypothesis::Expired ||
Chris@12:             m_accepted.getState() == Hypothesis::Rejected) {
Chris@11:             if (candidate >= 0) {
Chris@11:                 m_accepted = m_possible[candidate];
Chris@11:             } else {
Chris@11:                 m_accepted = Hypothesis();
Chris@11:             }
Chris@11:         }
Chris@8: 
Chris@14:         // reap rejected/expired hypotheses from possible list
Chris@14:         Hypotheses toReap = m_possible;
Chris@14:         m_possible.clear();
Chris@14:         for (int i = 0; i < toReap.size(); ++i) {
Chris@14:             Hypothesis h = toReap[i];
Chris@14:             if (h.getState() != Hypothesis::Rejected && 
Chris@14:                 h.getState() != Hypothesis::Expired) {
Chris@14:                 m_possible.push_back(h);
Chris@14:             }
Chris@14:         }
Chris@14:     }  
Chris@14: 
Chris@15:     std::cerr << "accepted length = " << m_accepted.getPendingLength()
Chris@15:               << ", state = " << m_accepted.getState()
Chris@15:               << ", hypothesis count = " << m_possible.size() << std::endl;
Chris@12: 
Chris@3:     delete[] data;
Chris@3:     return fs;
Chris@3: }
Chris@3: 
Chris@3: CepstrumPitchTracker::FeatureSet
Chris@3: CepstrumPitchTracker::getRemainingFeatures()
Chris@3: {
Chris@3:     FeatureSet fs;
Chris@15:     if (m_accepted.getState() == Hypothesis::Satisfied) {
Chris@16:         m_accepted.addFeatures(fs);
Chris@11:     }
Chris@3:     return fs;
Chris@3: }
Chris@3: 
Chris@3: void
Chris@3: CepstrumPitchTracker::fft(unsigned int n, bool inverse,
Chris@3:                     double *ri, double *ii, double *ro, double *io)
Chris@3: {
Chris@3:     if (!ri || !ro || !io) return;
Chris@3: 
Chris@3:     unsigned int bits;
Chris@3:     unsigned int i, j, k, m;
Chris@3:     unsigned int blockSize, blockEnd;
Chris@3: 
Chris@3:     double tr, ti;
Chris@3: 
Chris@3:     if (n < 2) return;
Chris@3:     if (n & (n-1)) return;
Chris@3: 
Chris@3:     double angle = 2.0 * M_PI;
Chris@3:     if (inverse) angle = -angle;
Chris@3: 
Chris@3:     for (i = 0; ; ++i) {
Chris@3: 	if (n & (1 << i)) {
Chris@3: 	    bits = i;
Chris@3: 	    break;
Chris@3: 	}
Chris@3:     }
Chris@3: 
Chris@23:     int *table = new int[n];
Chris@3: 
Chris@23:     for (i = 0; i < n; ++i) {
Chris@23: 	
Chris@23:         m = i;
Chris@3: 
Chris@23:         for (j = k = 0; j < bits; ++j) {
Chris@23:             k = (k << 1) | (m & 1);
Chris@23:             m >>= 1;
Chris@23:         }
Chris@3: 
Chris@23:         table[i] = k;
Chris@3:     }
Chris@3: 
Chris@3:     if (ii) {
Chris@3: 	for (i = 0; i < n; ++i) {
Chris@3: 	    ro[table[i]] = ri[i];
Chris@3: 	    io[table[i]] = ii[i];
Chris@3: 	}
Chris@3:     } else {
Chris@3: 	for (i = 0; i < n; ++i) {
Chris@3: 	    ro[table[i]] = ri[i];
Chris@3: 	    io[table[i]] = 0.0;
Chris@3: 	}
Chris@3:     }
Chris@3: 
Chris@3:     blockEnd = 1;
Chris@3: 
Chris@3:     for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
Chris@3: 
Chris@3: 	double delta = angle / (double)blockSize;
Chris@3: 	double sm2 = -sin(-2 * delta);
Chris@3: 	double sm1 = -sin(-delta);
Chris@3: 	double cm2 = cos(-2 * delta);
Chris@3: 	double cm1 = cos(-delta);
Chris@3: 	double w = 2 * cm1;
Chris@3: 	double ar[3], ai[3];
Chris@3: 
Chris@3: 	for (i = 0; i < n; i += blockSize) {
Chris@3: 
Chris@3: 	    ar[2] = cm2;
Chris@3: 	    ar[1] = cm1;
Chris@3: 
Chris@3: 	    ai[2] = sm2;
Chris@3: 	    ai[1] = sm1;
Chris@3: 
Chris@3: 	    for (j = i, m = 0; m < blockEnd; j++, m++) {
Chris@3: 
Chris@3: 		ar[0] = w * ar[1] - ar[2];
Chris@3: 		ar[2] = ar[1];
Chris@3: 		ar[1] = ar[0];
Chris@3: 
Chris@3: 		ai[0] = w * ai[1] - ai[2];
Chris@3: 		ai[2] = ai[1];
Chris@3: 		ai[1] = ai[0];
Chris@3: 
Chris@3: 		k = j + blockEnd;
Chris@3: 		tr = ar[0] * ro[k] - ai[0] * io[k];
Chris@3: 		ti = ar[0] * io[k] + ai[0] * ro[k];
Chris@3: 
Chris@3: 		ro[k] = ro[j] - tr;
Chris@3: 		io[k] = io[j] - ti;
Chris@3: 
Chris@3: 		ro[j] += tr;
Chris@3: 		io[j] += ti;
Chris@3: 	    }
Chris@3: 	}
Chris@3: 
Chris@3: 	blockEnd = blockSize;
Chris@3:     }
Chris@3: }
Chris@3: 
Chris@3: