Mercurial > hg > match-vamp
changeset 141:d6f22887283e adaptive_diagonals
Merge from branch refactors
| author | Chris Cannam |
|---|---|
| date | Fri, 09 Jan 2015 17:20:52 +0000 |
| parents | 50712e4b8c89 (diff) cfba9aec7569 (current diff) |
| children | |
| files | Makefile.inc Makefile.linux src/FeatureConditioner.h src/Finder.cpp src/MatchPipeline.cpp src/MatchPipeline.h src/MatchVampPlugin.cpp src/MatchVampPlugin.h test/TestFeatureExtractor.cpp test/expected.csv test/regressiontest.sh |
| diffstat | 10 files changed, 504 insertions(+), 72 deletions(-) [+] |
line wrap: on
line diff
--- a/Makefile.inc Thu Jan 08 12:11:27 2015 +0000 +++ b/Makefile.inc Fri Jan 09 17:20:52 2015 +0000 @@ -37,12 +37,13 @@ # DO NOT DELETE src/DistanceMetric.o: src/DistanceMetric.h +src/Path.o: src/Path.h src/FeatureConditioner.o: src/FeatureConditioner.h -src/Path.o: src/Path.h src/MatchFeatureFeeder.o: src/MatchFeatureFeeder.h src/Matcher.h src/MatchFeatureFeeder.o: src/DistanceMetric.h src/Finder.h src/FeatureExtractor.o: src/FeatureExtractor.h src/Finder.o: src/Finder.h src/Matcher.h src/DistanceMetric.h src/Path.h +src/Finder.o: src/MedianFilter.h src/Matcher.o: src/Matcher.h src/DistanceMetric.h src/MatchPipeline.o: src/MatchPipeline.h src/Matcher.h src/DistanceMetric.h src/MatchPipeline.o: src/Finder.h src/FeatureExtractor.h @@ -50,7 +51,7 @@ src/MatchVampPlugin.o: src/MatchVampPlugin.h src/MatchPipeline.h src/MatchVampPlugin.o: src/Matcher.h src/DistanceMetric.h src/Finder.h src/MatchVampPlugin.o: src/FeatureExtractor.h src/FeatureConditioner.h -src/MatchVampPlugin.o: src/MatchFeatureFeeder.h src/Path.h +src/MatchVampPlugin.o: src/MatchFeatureFeeder.h src/MatchFeatureFeeder.o: src/Matcher.h src/DistanceMetric.h src/Finder.h src/Finder.o: src/Matcher.h src/DistanceMetric.h src/Matcher.o: src/DistanceMetric.h @@ -60,4 +61,6 @@ src/MatchVampPlugin.o: src/MatchPipeline.h src/Matcher.h src/DistanceMetric.h src/MatchVampPlugin.o: src/Finder.h src/FeatureExtractor.h src/MatchVampPlugin.o: src/FeatureConditioner.h src/MatchFeatureFeeder.h +test/TestFeatureConditioner.o: src/FeatureConditioner.h +test/TestDistanceMetric.o: src/DistanceMetric.h test/TestFeatureExtractor.o: src/FeatureExtractor.h
--- a/Makefile.linux Thu Jan 08 12:11:27 2015 +0000 +++ b/Makefile.linux Fri Jan 09 17:20:52 2015 +0000 @@ -1,6 +1,6 @@ -#CXXFLAGS += -fPIC -ffast-math -O3 -Wall -Werror -Wfloat-conversion -CXXFLAGS += -fPIC -g -Wall -Werror -DPERFORM_ERROR_CHECKS=1 +#CXXFLAGS += -fPIC -ffast-math -O3 -Wall -Werror +CXXFLAGS += -fPIC -g -Wall -Werror -Wfloat-conversion -DPERFORM_ERROR_CHECKS=1 LDFLAGS += -Wl,-Bstatic -lvamp-sdk -Wl,-Bdynamic PLUGIN_LDFLAGS += -shared -Wl,-Bsymbolic -Wl,-z,defs -lpthread -Wl,--version-script=vamp-plugin.map
--- a/src/Finder.cpp Thu Jan 08 12:11:27 2015 +0000 +++ b/src/Finder.cpp Fri Jan 09 17:20:52 2015 +0000 @@ -17,6 +17,7 @@ #include "Finder.h" #include "Path.h" +#include "MedianFilter.h" #include <algorithm> #include <iomanip> @@ -96,11 +97,19 @@ } void -Finder::recalculatePathCostMatrix(int r1, int c1, int r2, int c2) +Finder::recalculatePathCostMatrix(int r1, int c1, int r2, int c2) { + recalculatePathCostMatrix(r1, c1, r2, c2, m_m->getDiagonalWeight()); +} + +void +Finder::recalculatePathCostMatrix(int r1, int c1, int r2, int c2, + float diagonalWeight) +{ + cerr << "recalculatePathCostMatrix: (" << r1 << "," << c1 << ") -> (" + << r2 << "," << c2 << ")" << endl; + int prevRowStart = 0, prevRowStop = 0; - - float diagonalWeight = m_m->getDiagonalWeight(); for (int r = r1; r <= r2; r++) { @@ -311,7 +320,13 @@ for (int j = -4; j <= 0; ++j) { cerr << setw(w) << j; for (int i = -4; i <= 0; ++i) { - cerr << setw(ww) << m_m->getDistance(err.r + j, err.c + i); + int r = err.r + j; + int c = err.c + i; + float dist = 0.f; + if (r >= 0 && c >= 0) { + dist = m_m->getDistance(r, c); + } + cerr << setw(ww) << dist; } cerr << endl; } @@ -326,7 +341,13 @@ for (int j = -4; j <= 0; ++j) { cerr << setw(w) << j; for (int i = -4; i <= 0; ++i) { - cerr << setw(ww) << m_m->getPathCost(err.r + j, err.c + i); + int r = err.r + j; + int c = err.c + i; + double cost = 0.0; + if (r >= 0 && c >= 0) { + cost = m_m->getPathCost(r, c); + } + cerr << setw(ww) << cost; } cerr << endl; } @@ -335,25 +356,18 @@ } #endif -int -Finder::retrievePath(bool smooth, vector<int> &pathx, vector<int> &pathy) +void +Finder::getEndPoint(int &x, int &y) { - pathx.clear(); - pathy.clear(); - -#ifdef PERFORM_ERROR_CHECKS - checkAndReport(); -#endif - int ex = m_m->getOtherFrameCount() - 1; int ey = m_m->getFrameCount() - 1; if (ex < 0 || ey < 0) { - return 0; + return; } - int x = ex; - int y = ey; + x = ex; + y = ey; #ifdef DEBUG_FINDER cerr << "*** retrievePath: smooth = " << smooth << endl; @@ -377,11 +391,26 @@ y = ey; } - recalculatePathCostMatrix(0, 0, y, x); - #ifdef DEBUG_FINDER cerr << "*** retrievePath: start: x = " << x << ", y = " << y << endl; #endif +} + +void +Finder::retrievePath(vector<int> &pathx, + vector<int> &pathy, + vector<float> &distances) +{ + pathx.clear(); + pathy.clear(); + distances.clear(); + +#ifdef PERFORM_ERROR_CHECKS + checkAndReport(); +#endif + + int x, y; + getEndPoint(x, y); while (m_m->isAvailable(y, x) && (x > 0 || y > 0)) { @@ -389,6 +418,7 @@ pathx.push_back(x); pathy.push_back(y); + distances.push_back(m_m->getDistance(y, x)); switch (m_m->getAdvance(y, x)) { case Matcher::AdvanceThis: @@ -422,13 +452,146 @@ reverse(pathx.begin(), pathx.end()); reverse(pathy.begin(), pathy.end()); - - if (smooth) { - int smoothedLen = Path().smooth(pathx, pathy, pathx.size()); - return smoothedLen; - } else { - return pathx.size(); - } + reverse(distances.begin(), distances.end()); } +void +Finder::smooth(const vector<double> &mag1, const vector<double> &mag2) +{ + if (m_m->getDiagonalWeight() <= 1.0) { + cerr << "Finder::smooth: Diagonal weight is already " + << m_m->getDiagonalWeight() << ", adaptive smoothing will have " + << "no effect, skipping it" << endl; + return; + } + vector<double> confidence; + + vector<int> pathx, pathy; + vector<float> distances; + retrievePath(pathx, pathy, distances); + + int len = pathx.size(); + + for (int i = 0; i < len; ++i) { + + int x = pathx[i]; + int y = pathy[i]; + + double magSum = mag1[x] + mag2[y]; + double distance = distances[i]; + double c = magSum - distance * magSum; + confidence.push_back(c); + + /* + if (x != prevx) { + Feature f; + f.values.push_back(magSum); + returnFeatures[m_magOutNo].push_back(f); + + f.values.clear(); + f.values.push_back(distance); + returnFeatures[m_distOutNo].push_back(f); + } + + prevx = x; + prevy = y; + */ + } + + confidence = MedianFilter<double>::filter(3, confidence); + vector<double> filtered = MedianFilter<double>::filter(50, confidence); + for (int i = 0; i < len; ++i) { + confidence[i] -= filtered[i]; + if (confidence[i] < 0.f) { + confidence[i] = 0.f; + } + } + vector<double> deriv; + deriv.resize(len, 0.f); + for (int i = 1; i < len; ++i) { + deriv[i] = confidence[i] - confidence[i-1]; + } + vector<int> inflections; + for (int i = 1; i < len; ++i) { + if (deriv[i-1] > 0 && deriv[i] < 0) { + inflections.push_back(i); + } + } + + /* + for (int i = 0; i < len; ++i) { + + int x = pathx[i]; + int y = pathy[i]; + + if (x != prevx) { + Feature f; + double c = confidence[i]; + f.values.push_back(c); + returnFeatures[m_confidenceOutNo].push_back(f); + } + + prevx = x; + prevy = y; + } + */ + + map<int, int> pinpoints; + + for (int ii = 0; ii < int(inflections.size()); ++ii) { + + int i = inflections[ii]; + + int x = pathx[i]; + int y = pathy[i]; + + pinpoints[x] = y; + } + + cerr << "Pin points are:" << endl; + for (map<int, int>::const_iterator i = pinpoints.begin(); + i != pinpoints.end(); ++i) { + cerr << "[" << i->first << "," << i->second << "] "; + } + cerr << endl; + + if (pinpoints.size() < 2) return; + + int ex, ey; + getEndPoint(ex, ey); + + pair<int, int> prev(0, 0); + + for (map<int, int>::const_iterator i = pinpoints.begin(); + i != pinpoints.end(); ++i) { + + if (i == pinpoints.begin()) continue; + + map<int, int>::const_iterator j = i; + ++j; + if (j != pinpoints.end() && i->second == j->second) { + // skip this one + continue; + } + + pair<int, int> curr = *i; + + if (curr.first > ex || curr.second > ey) { + curr = pair<int, int>(ex, ey); + } + + recalculatePathCostMatrix(prev.second, prev.first, + curr.second, curr.first, + 1.0); + + prev = curr; + } + + recalculatePathCostMatrix(prev.second, prev.first, + m_m->getFrameCount() - 1, + m_m->getOtherFrameCount() - 1, + 1.0); +} + +
--- a/src/Finder.h Thu Jan 08 12:11:27 2015 +0000 +++ b/src/Finder.h Fri Jan 09 17:20:52 2015 +0000 @@ -18,6 +18,7 @@ #define _FINDER_H_ #include <vector> +#include <map> #include <iostream> #include "Matcher.h" @@ -54,17 +55,41 @@ void recalculatePathCostMatrix(int r1, int c1, int r2, int c2); /** + * Recalculate a rectangle of the path cost matrix using the given + * diagonal weight instead of the weight obtained from the Matcher + * parameters. + */ + void recalculatePathCostMatrix(int r1, int c1, int r2, int c2, + float diagonalWeight); + + /** * Track back after all of the matchers have been fed in order to * obtain the lowest cost path available. Path x and y coordinate * pairs are returned in corresponding elements of pathx and - * pathy. Return value is the length of the returned path: only - * this many elements from pathx and pathy are valid (any - * subsequent ones may be spurious). + * pathy, and the corresponding elements of the returned distance + * vector contain the distance between the pair of features + * selected at each step in the path. + */ + void retrievePath(std::vector<int> &pathx, + std::vector<int> &pathy, + std::vector<float> &distances); + + /** + * Using the provided magnitude arrays for the first and second + * inputs as an indication of salience, smooth the cost matrix by + * recalculating with locally adaptive diagonal weights in order + * to obtain a smoother path that retains the original locations + * for the most salient points. Subsequent calls to retrievePath + * will retrieve the smoothed version. The original can be + * restored with a call to recalculatePathCostMatrix(). * - * @param smooth whether to smooth the path before returning it + * Note that this is quite separate from Path::smooth() which + * smooths a path without trying to maintain locations for salient + * points. That function will achieve a smoother path, but may + * smooth out significant locations such as onsets. */ - int retrievePath(bool smooth, std::vector<int> &pathx, std::vector<int> &pathy); - + void smooth(const vector<double> &mag1, const vector<double> &mag2); + protected: #ifdef PERFORM_ERROR_CHECKS struct ErrorPosition { @@ -83,6 +108,8 @@ ErrorPosition checkPathCostMatrix(); void checkAndReport(); #endif + + void getEndPoint(int &x, int &y); Matcher *m_m; int m_duration1;
--- a/src/MatchPipeline.cpp Thu Jan 08 12:11:27 2015 +0000 +++ b/src/MatchPipeline.cpp Fri Jan 09 17:20:52 2015 +0000 @@ -99,6 +99,23 @@ f2 = m_f2; } +static double +magOf(const vector<double> &f) +{ + double mag = 0.0; + for (int j = 0; j < (int)f.size(); ++j) { + mag += f[j] * f[j]; + } + return sqrt(mag); +} + +void +MatchPipeline::extractFeatureMagnitudes(double &mag1, double &mag2) +{ + mag1 = magOf(m_f1); + mag2 = magOf(m_f2); +} + void MatchPipeline::extractConditionedFeatures(vector<double> &c1, vector<double> &c2) {
--- a/src/MatchPipeline.h Thu Jan 08 12:11:27 2015 +0000 +++ b/src/MatchPipeline.h Fri Jan 09 17:20:52 2015 +0000 @@ -74,11 +74,21 @@ * If a frame was just fed in at the first or second pipeline * stage, it can be retrieved from the second stage here. That is, * if you provided frequency-domain audio, extractFeatures will - * give you back the FeatureExtractor's features. + * give you back the FeatureExtractor's features; if you provided + * features, extractFeatures will simply give you those back. If + * you provided conditioned features, the return from this + * function is undefined. */ void extractFeatures(vector<double> &f1, vector<double> &f2); /** + * If a frame was just fed in at the first or second pipeline + * stage, you can obtain the magnitudes of its features for both + * inputs here. + */ + void extractFeatureMagnitudes(double &mag1, double &mag2); + + /** * Retrieve the conditioned features from the third pipeline stage. */ void extractConditionedFeatures(vector<double> &f1, vector<double> &f2);
--- a/src/MatchVampPlugin.cpp Thu Jan 08 12:11:27 2015 +0000 +++ b/src/MatchVampPlugin.cpp Fri Jan 09 17:20:52 2015 +0000 @@ -19,7 +19,6 @@ #include "Matcher.h" #include "MatchFeatureFeeder.h" #include "FeatureExtractor.h" -#include "Path.h" #include <vamp/vamp.h> #include <vamp-sdk/PluginAdapter.h> @@ -27,6 +26,11 @@ #include <vector> #include <algorithm> +#include <map> + +#include <cstdio> + +using namespace std; //static int extant = 0; @@ -66,10 +70,10 @@ m_defaultFcParams() { if (inputSampleRate < sampleRateMin) { - std::cerr << "MatchVampPlugin::MatchVampPlugin: input sample rate " + cerr << "MatchVampPlugin::MatchVampPlugin: input sample rate " << inputSampleRate << " < min supported rate " << sampleRateMin << ", plugin will refuse to initialise" - << std::endl; + << endl; } if (!m_serialisingMutexInitialised) { @@ -82,12 +86,11 @@ } m_pipeline = 0; -// std::cerr << "MatchVampPlugin::MatchVampPlugin(" << this << "): extant = " << ++extant << std::endl; } MatchVampPlugin::~MatchVampPlugin() { -// std::cerr << "MatchVampPlugin::~MatchVampPlugin(" << this << "): extant = " << --extant << std::endl; +// cerr << "MatchVampPlugin::~MatchVampPlugin(" << this << "): extant = " << --extant << endl; delete m_pipeline; @@ -209,7 +212,7 @@ desc.description = "Total frame energy threshold below which a feature will be regarded as silent"; desc.minValue = 0; desc.maxValue = 1; - desc.defaultValue = m_defaultFcParams.silenceThreshold; + desc.defaultValue = (float)m_defaultFcParams.silenceThreshold; desc.isQuantized = false; list.push_back(desc); @@ -259,7 +262,7 @@ } float -MatchVampPlugin::getParameter(std::string name) const +MatchVampPlugin::getParameter(string name) const { if (name == "serialise") { return m_serialise ? 1.0 : 0.0; @@ -280,14 +283,14 @@ } else if (name == "smooth") { return m_smooth ? 1.0 : 0.0; } else if (name == "silencethreshold") { - return m_fcParams.silenceThreshold; + return (float)m_fcParams.silenceThreshold; } return 0.0; } void -MatchVampPlugin::setParameter(std::string name, float value) +MatchVampPlugin::setParameter(string name, float value) { if (name == "serialise") { m_serialise = (value > 0.5); @@ -337,9 +340,9 @@ MatchVampPlugin::initialise(size_t channels, size_t stepSize, size_t blockSize) { if (m_inputSampleRate < sampleRateMin) { - std::cerr << "MatchVampPlugin::MatchVampPlugin: input sample rate " + cerr << "MatchVampPlugin::MatchVampPlugin: input sample rate " << m_inputSampleRate << " < min supported rate " - << sampleRateMin << std::endl; + << sampleRateMin << endl; return false; } if (channels < getMinChannelCount() || @@ -364,6 +367,10 @@ delete m_pipeline; m_pipeline = 0; m_frameNo = 0; + + m_mag1.clear(); + m_mag2.clear(); + createMatchers(); m_begin = true; m_locked = false; @@ -497,6 +504,59 @@ m_cbFeaturesOutNo = list.size(); list.push_back(desc); + desc.identifier = "feature_distance"; + desc.name = "Feature Distance"; + desc.description = "Value of distance metric between features at each point-in-A along the chosen alignment path"; + desc.unit = ""; + desc.hasFixedBinCount = true; + desc.binCount = 1; + desc.hasKnownExtents = false; + desc.isQuantized = false; + desc.sampleType = OutputDescriptor::FixedSampleRate; + desc.sampleRate = outRate; + m_distOutNo = list.size(); + list.push_back(desc); + + desc.identifier = "feature_mag"; + desc.name = "Feature Magnitudes"; + desc.description = "Sum of magnitudes of feature pairs for each point-in-A along the chosen alignment path"; + desc.unit = ""; + desc.hasFixedBinCount = true; + desc.binCount = 1; + desc.hasKnownExtents = false; + desc.isQuantized = false; + desc.sampleType = OutputDescriptor::FixedSampleRate; + desc.sampleRate = outRate; + m_magOutNo = list.size(); + list.push_back(desc); +/* + desc.identifier = "confidence"; + desc.name = "Confidence"; + desc.description = "Confidence metric for the quality of match at each point-in-A along the chosen alignment path"; + desc.unit = ""; + desc.hasFixedBinCount = true; + desc.binCount = 1; + desc.hasKnownExtents = false; + desc.isQuantized = false; + desc.sampleType = OutputDescriptor::FixedSampleRate; + desc.sampleRate = outRate; + m_confidenceOutNo = list.size(); + list.push_back(desc); + + desc.identifier = "confidence_peaks"; + desc.name = "Confidence Peaks"; + desc.description = "Peak locations for the confidence metric"; + desc.unit = ""; + desc.hasFixedBinCount = true; + desc.binCount = 0; + desc.hasKnownExtents = false; + desc.isQuantized = false; + desc.sampleType = OutputDescriptor::FixedSampleRate; + desc.sampleRate = outRate; + m_confPeakOutNo = list.size(); + list.push_back(desc); +*/ + return list; } @@ -517,10 +577,16 @@ m_begin = false; } -// std::cerr << timestamp.toString(); +// cerr << timestamp.toString(); m_pipeline->feedFrequencyDomainAudio(inputBuffers[0], inputBuffers[1]); + + double m1, m2; + m_pipeline->extractFeatureMagnitudes(m1, m2); + m_mag1.push_back(m1); + m_mag2.push_back(m2); + FeatureSet returnFeatures; vector<double> f1, f2; @@ -572,18 +638,25 @@ FeatureSet returnFeatures; Finder *finder = m_pipeline->getFinder(); - std::vector<int> pathx; - std::vector<int> pathy; - int len = finder->retrievePath(m_smooth, pathx, pathy); + + if (m_smooth) { + finder->smooth(m_mag1, m_mag2); + } + + vector<int> pathx; + vector<int> pathy; + vector<float> distances; + finder->retrievePath(pathx, pathy, distances); int prevx = 0; int prevy = 0; + int len = pathx.size(); for (int i = 0; i < len; ++i) { int x = pathx[i]; int y = pathy[i]; - + Vamp::RealTime xt = Vamp::RealTime::frame2RealTime (x * m_stepSize, lrintf(m_inputSampleRate)); Vamp::RealTime yt = Vamp::RealTime::frame2RealTime @@ -654,31 +727,31 @@ /* for (int i = 0; i < len; ++i) { - std::cerr << i << ": [" << pathx[i] << "," << pathy[i] << "]" << std::endl; + cerr << i << ": [" << pathx[i] << "," << pathy[i] << "]" << endl; } - std::cerr << std::endl; - std::cerr << "File: A" << std::endl; - std::cerr << "Marks: -1" << std::endl; - std::cerr << "FixedPoints: true 0" << std::endl; - std::cerr << "0" << std::endl; - std::cerr << "0" << std::endl; - std::cerr << "0" << std::endl; - std::cerr << "0" << std::endl; - std::cerr << "File: B" << std::endl; - std::cerr << "Marks: 0" << std::endl; - std::cerr << "FixedPoints: true 0" << std::endl; - std::cerr << "0.02" << std::endl; - std::cerr << "0.02" << std::endl; + cerr << endl; + cerr << "File: A" << endl; + cerr << "Marks: -1" << endl; + cerr << "FixedPoints: true 0" << endl; + cerr << "0" << endl; + cerr << "0" << endl; + cerr << "0" << endl; + cerr << "0" << endl; + cerr << "File: B" << endl; + cerr << "Marks: 0" << endl; + cerr << "FixedPoints: true 0" << endl; + cerr << "0.02" << endl; + cerr << "0.02" << endl; - std::cerr << len << std::endl; + cerr << len << endl; for (int i = 0; i < len; ++i) { - std::cerr << pathx[i] << std::endl; + cerr << pathx[i] << endl; } - std::cerr << len << std::endl; + cerr << len << endl; for (int i = 0; i < len; ++i) { - std::cerr << pathy[i] << std::endl; + cerr << pathy[i] << endl; } */ }
--- a/src/MatchVampPlugin.h Thu Jan 08 12:11:27 2015 +0000 +++ b/src/MatchVampPlugin.h Fri Jan 09 17:20:52 2015 +0000 @@ -90,6 +90,9 @@ FeatureConditioner::Parameters m_fcParams; FeatureConditioner::Parameters m_defaultFcParams; + std::vector<double> m_mag1; + std::vector<double> m_mag2; + mutable int m_pathOutNo; mutable int m_abOutNo; mutable int m_baOutNo; @@ -99,6 +102,12 @@ mutable int m_bFeaturesOutNo; mutable int m_caFeaturesOutNo; mutable int m_cbFeaturesOutNo; + mutable int m_magOutNo; + mutable int m_distOutNo; +/* + mutable int m_confidenceOutNo; + mutable int m_confPeakOutNo; +*/ #ifdef _WIN32 static HANDLE m_serialisingMutex;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/MedianFilter.h Fri Jan 09 17:20:52 2015 +0000 @@ -0,0 +1,130 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ + +/* + QM DSP Library + + Centre for Digital Music, Queen Mary, University of London. + This file Copyright 2010 Chris Cannam. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. See the file + COPYING included with this distribution for more information. +*/ + +#ifndef MEDIAN_FILTER_H +#define MEDIAN_FILTER_H + +#include <algorithm> +#include <cassert> +#include <cmath> +#include <iostream> +#include <vector> + +template <typename T> +class MedianFilter +{ +public: + MedianFilter(int size, float percentile = 50.f) : + m_size(size), + m_frame(new T[size]), + m_sorted(new T[size]), + m_sortend(m_sorted + size - 1) { + setPercentile(percentile); + reset(); + } + + ~MedianFilter() { + delete[] m_frame; + delete[] m_sorted; + } + + void setPercentile(float p) { + m_index = int((m_size * p) / 100.f); + if (m_index >= m_size) m_index = m_size-1; + if (m_index < 0) m_index = 0; + } + + void push(T value) { + if (value != value) { + std::cerr << "WARNING: MedianFilter::push: attempt to push NaN" << std::endl; + return; // nan + } + drop(m_frame[0]); + const int sz1 = m_size-1; + for (int i = 0; i < sz1; ++i) m_frame[i] = m_frame[i+1]; + m_frame[m_size-1] = value; + put(value); + } + + T get() const { + return m_sorted[m_index]; + } + + int getSize() const { + return m_size; + } + + T getAt(float percentile) { + int ix = int((m_size * percentile) / 100.f); + if (ix >= m_size) ix = m_size-1; + if (ix < 0) ix = 0; + return m_sorted[ix]; + } + + void reset() { + for (int i = 0; i < m_size; ++i) m_frame[i] = 0; + for (int i = 0; i < m_size; ++i) m_sorted[i] = 0; + } + + static std::vector<T> filter(int size, const std::vector<T> &in) { + std::vector<T> out; + MedianFilter<T> f(size); + for (int i = 0; i < int(in.size()); ++i) { + f.push(in[i]); + T median = f.get(); + if (i >= size/2) out.push_back(median); + } + while (out.size() < in.size()) { + f.push(T()); + out.push_back(f.get()); + } + return out; + } + +private: + const int m_size; + T *const m_frame; + T *const m_sorted; + T *const m_sortend; + int m_index; + + void put(T value) { + // precondition: m_sorted contains m_size-1 values, packed at start + // postcondition: m_sorted contains m_size values, one of which is value + T *point = std::lower_bound(m_sorted, m_sortend, value); + const int n = m_sortend - point; + for (int i = n; i > 0; --i) point[i] = point[i-1]; + *point = value; + } + + void drop(T value) { + // precondition: m_sorted contains m_size values, one of which is value + // postcondition: m_sorted contains m_size-1 values, packed at start + T *point = std::lower_bound(m_sorted, m_sortend + 1, value); + if (*point != value) { + std::cerr << "WARNING: MedianFilter::drop: *point is " << *point + << ", expected " << value << std::endl; + } + const int n = m_sortend - point; + for (int i = 0; i < n; ++i) point[i] = point[i+1]; + *m_sortend = T(0); + } + + MedianFilter(const MedianFilter &); // not provided + MedianFilter &operator=(const MedianFilter &); // not provided +}; + +#endif +
--- a/test/TestFeatureExtractor.cpp Thu Jan 08 12:11:27 2015 +0000 +++ b/test/TestFeatureExtractor.cpp Fri Jan 09 17:20:52 2015 +0000 @@ -14,7 +14,7 @@ static int freq2mid(double freq) { - return round(57.0 + 12.0 * log(freq / 220.) / log(2.)); + return int(round(57.0 + 12.0 * log(freq / 220.) / log(2.))); } static int freq2chroma(double freq)