From 6940ad29d892bd8dfb40f158a01b66d2bf8f7b6b Mon Sep 17 00:00:00 2001 From: Chris Cannam Date: Fri, 10 Jun 2022 12:35:15 +0100 Subject: [PATCH] Introduce a histogram filter for filtering discrete values --- src/common/HistogramFilter.h | 176 ++++++++++++++++++++++++++++ src/common/MovingMedian.h | 5 + src/common/RingBuffer.h | 1 - src/common/SingleThreadRingBuffer.h | 144 +++++++++++++++++++++++ src/finer/BinSegmenter.h | 7 +- src/test/TestSignalBits.cpp | 109 +++++++++++++++++ 6 files changed, 438 insertions(+), 4 deletions(-) create mode 100644 src/common/HistogramFilter.h create mode 100644 src/common/SingleThreadRingBuffer.h diff --git a/src/common/HistogramFilter.h b/src/common/HistogramFilter.h new file mode 100644 index 0000000..50aa137 --- /dev/null +++ b/src/common/HistogramFilter.h @@ -0,0 +1,176 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ + +/* + Rubber Band Library + An audio time-stretching and pitch-shifting library. + Copyright 2007-2022 Particular Programs Ltd. + + 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. + + Alternatively, if you have a valid commercial licence for the + Rubber Band Library obtained by agreement with the copyright + holders, you may redistribute and/or modify it under the terms + described in that licence. + + If you wish to distribute code using the Rubber Band Library + under terms other than those of the GNU General Public License, + you must obtain a valid commercial licence before doing so. +*/ + +#ifndef RUBBERBAND_HISTOGRAM_FILTER_H +#define RUBBERBAND_HISTOGRAM_FILTER_H + +#include "SingleThreadRingBuffer.h" + +#include +#include + +namespace RubberBand { + +/** + * A median or modal filter implemented using a histogram. The values + * must come from a compact set of integers within [0,n) where n is + * specified in the constructor. Pushing a value updates the + * histogram, after which you can get either the median or the mode. + * You can call drop() to drop the oldest value without pushing a new + * one, for example to drain the filter at the tail of the sequence. + */ +class HistogramFilter +{ +public: + HistogramFilter(int nValues, int filterLength) : + m_buffer(filterLength), + m_histogram(nValues, 0) { + } + ~HistogramFilter() { } + + int getFilterLength() const { + return m_buffer.getSize(); + } + + void reset() { + m_buffer.reset(); + for (int i = 0; i < m_histogram.size(); ++i) { + m_histogram[i] = 0; + } + } + + void push(int value) { + if (m_buffer.getWriteSpace() == 0) { + int toDrop = m_buffer.readOne(); + --m_histogram[toDrop]; + } + m_buffer.writeOne(value); + ++m_histogram[value]; + } + + void drop() { + if (m_buffer.getReadSpace() > 0) { + int toDrop = m_buffer.readOne(); + --m_histogram[toDrop]; + } + } + + /** Return the median of the values in the filter currently. If + * the median lies between two values, return the first of them. + */ + int getMedian() const { + int half = (m_buffer.getReadSpace() + 1) / 2; + int acc = 0; + for (int i = 0; i < m_histogram.size(); ++i) { + acc += m_histogram[i]; + if (acc >= half) { + return i; + } + } + return 0; + } + + /** Return the modal value, that is, the value that occurs the + * most often in the filter currently. If multiple values occur + * an equal number of times, return the smallest of them. + */ + int getMode() const { + int max = 0; + int mode = 0; + for (int i = 0; i < m_histogram.size(); ++i) { + int h = m_histogram[i]; + if (i == 0 || h > max) { + max = h; + mode = i; + } + } + return mode; + } + + // Convenience function that filters an array in-place. Filter is + // a median filter unless modal is true. Array has length n. + // Modifies both the filter and the array. + // + static void filter(HistogramFilter &f, int *v, int n, bool modal = false) { + f.reset(); + int flen = f.getFilterLength(); + int i = -(flen / 2); + int j = 0; + while (i != n) { + if (j < n) { + f.push(v[j]); + } else if (j >= flen) { + f.drop(); + } + if (i >= 0) { + int m = modal ? f.getMode() : f.getMedian(); + v[i] = m; + } + ++i; + ++j; + } + } + + // As above but with a vector argument + // + static void filter(HistogramFilter &f, std::vector &v, bool modal = false) { + filter(f, v.data(), v.size(), modal); + } + + // Convenience function that median-filters an array + // in-place. Array has length n. Modifies both the filter and the + // array. + // + static void medianFilter(HistogramFilter &f, int *v, int n) { + filter(f, v, n, false); + } + + // As above but with a vector argument + // + static void medianFilter(HistogramFilter &f, std::vector &v) { + medianFilter(f, v.data(), v.size()); + } + + // Convenience function that modal-filters an array + // in-place. Array has length n. Modifies both the filter and the + // array. + // + static void modalFilter(HistogramFilter &f, int *v, int n) { + filter(f, v, n, true); + } + + // As above but with a vector argument + // + static void modalFilter(HistogramFilter &f, std::vector &v) { + modalFilter(f, v.data(), v.size()); + } + +private: + SingleThreadRingBuffer m_buffer; + std::vector m_histogram; + int m_mode; +}; + +} + +#endif diff --git a/src/common/MovingMedian.h b/src/common/MovingMedian.h index e8d40ed..1dfc3ed 100644 --- a/src/common/MovingMedian.h +++ b/src/common/MovingMedian.h @@ -107,6 +107,11 @@ private: // precondition: sorted contains m_length values, one of which is toDrop // postcondition: sorted contains m_length values, one of which is toPut // (and one instance of toDrop has been removed) + + // This implementation was timed for rather short filters (no + // longer than maybe 16 items). Two binary searches plus a + // memmove should be faster for longer ones. + const int n = m_length; T *sorted = sortedFor(filter); int dropIx; diff --git a/src/common/RingBuffer.h b/src/common/RingBuffer.h index a905a33..2fabd14 100644 --- a/src/common/RingBuffer.h +++ b/src/common/RingBuffer.h @@ -44,7 +44,6 @@ namespace RubberBand { * RingBuffer is thread-safe provided only one thread writes and only * one thread reads. */ - template class RingBuffer { diff --git a/src/common/SingleThreadRingBuffer.h b/src/common/SingleThreadRingBuffer.h new file mode 100644 index 0000000..e3e4925 --- /dev/null +++ b/src/common/SingleThreadRingBuffer.h @@ -0,0 +1,144 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ + +/* + Rubber Band Library + An audio time-stretching and pitch-shifting library. + Copyright 2007-2022 Particular Programs Ltd. + + 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. + + Alternatively, if you have a valid commercial licence for the + Rubber Band Library obtained by agreement with the copyright + holders, you may redistribute and/or modify it under the terms + described in that licence. + + If you wish to distribute code using the Rubber Band Library + under terms other than those of the GNU General Public License, + you must obtain a valid commercial licence before doing so. +*/ + +#ifndef RUBBERBAND_SINGLE_THREAD_RINGBUFFER_H +#define RUBBERBAND_SINGLE_THREAD_RINGBUFFER_H + +#include + +namespace RubberBand { + +/** + * SingleThreadRingBuffer implements a ring buffer to be used to store + * a sample type T, for reading and writing within a single + * thread. SingleThreadRingBuffer is a simple container, not a + * thread-safe lock-free structure: use RingBuffer for the situation + * with reader and writer in separate threads. Currently this + * implementation only supports reading and writing a single sample at + * a time. + */ +template +class SingleThreadRingBuffer +{ +public: + /** + * Create a ring buffer with room to write n samples. + * + * Note that the internal storage size will actually be n+1 + * samples, as one element is unavailable for administrative + * reasons. Since the ring buffer performs best if its size is a + * power of two, this means n should ideally be some power of two + * minus one. + */ + SingleThreadRingBuffer(int n) : + m_buffer(n + 1, T()), + m_writer(0), + m_reader(0), + m_size(n + 1) { } + + virtual ~SingleThreadRingBuffer() { } + + /** + * Return the total capacity of the ring buffer in samples. + * (This is the argument n passed to the constructor.) + */ + int getSize() const { + return m_size - 1; + } + + /** + * Reset read and write pointers, thus emptying the buffer. + */ + void reset() { + m_writer = m_reader; + } + + /** + * Return the amount of data available for reading, in samples. + */ + int getReadSpace() const { + if (m_writer > m_reader) return m_writer - m_reader; + else if (m_writer < m_reader) return (m_writer + m_size) - m_reader; + else return 0; + } + + /** + * Return the amount of space available for writing, in samples. + */ + int getWriteSpace() const { + int space = (m_reader + m_size - m_writer - 1); + if (space >= m_size) space -= m_size; + return space; + } + + /** + * Read one sample from the buffer. If no sample is available, + * silently return zero. + */ + T readOne() { + if (m_writer == m_reader) { + return {}; + } + auto value = m_buffer[m_reader]; + if (++m_reader == m_size) m_reader = 0; + return value; + } + + /** + * Pretend to read one sample from the buffer, without actually + * returning it (i.e. discard the next sample). + */ + void skipOne() { + if (m_writer == m_reader) { + return; + } + if (++m_reader == m_size) m_reader = 0; + } + + /** + * Write one sample to the buffer. If insufficient space is + * available, the sample will not be written. Returns the number + * of samples actually written, i.e. 0 or 1. + */ + int writeOne(const T &value) { + if (getWriteSpace() == 0) return 0; + m_buffer[m_writer] = value; + if (++m_writer == m_size) m_writer = 0; + return 1; + } + +protected: + std::vector m_buffer; + int m_writer; + int m_reader; + int m_size; + +private: + SingleThreadRingBuffer(const SingleThreadRingBuffer &); // not provided + SingleThreadRingBuffer &operator=(const SingleThreadRingBuffer &); // not provided +}; + +} + +#endif + diff --git a/src/finer/BinSegmenter.h b/src/finer/BinSegmenter.h index dc4dee0..5b2c7d3 100644 --- a/src/finer/BinSegmenter.h +++ b/src/finer/BinSegmenter.h @@ -25,6 +25,7 @@ #define RUBBERBAND_BIN_SEGMENTER_H #include "BinClassifier.h" +#include "../common/HistogramFilter.h" #include @@ -54,7 +55,7 @@ public: BinSegmenter(Parameters parameters) : m_parameters(parameters), m_numeric(m_parameters.binCount, 0), - m_classFilter(15) + m_classFilter(3, 15) { } @@ -70,7 +71,7 @@ public: m_numeric[i] = 2; break; } } - MovingMedian::filter(m_classFilter, m_numeric); + HistogramFilter::modalFilter(m_classFilter, m_numeric); /* std::cout << "c:"; for (int i = 0; i < n; ++i) { @@ -125,7 +126,7 @@ public: protected: Parameters m_parameters; std::vector m_numeric; - MovingMedian m_classFilter; + HistogramFilter m_classFilter; //!!! dupes int binForFrequency(double f) const { diff --git a/src/test/TestSignalBits.cpp b/src/test/TestSignalBits.cpp index 73dc467..383fbcd 100644 --- a/src/test/TestSignalBits.cpp +++ b/src/test/TestSignalBits.cpp @@ -25,6 +25,7 @@ #include #include "../common/MovingMedian.h" +#include "../common/HistogramFilter.h" #include "../finer/Peak.h" using namespace RubberBand; @@ -94,6 +95,114 @@ BOOST_AUTO_TEST_CASE(moving_median_n_1) BOOST_TEST(arr == expected, tt::per_element()); } +BOOST_AUTO_TEST_CASE(histogram_median_simple_3) +{ + HistogramFilter hf(5, 3); // nValues, filterLength + vector arr { 1, 2, 3 }; + vector expected { 1, 2, 2 }; + HistogramFilter::medianFilter(hf, arr); + BOOST_TEST(arr == expected, tt::per_element()); +} + +BOOST_AUTO_TEST_CASE(histogram_median_simple_4) +{ + HistogramFilter hf(5, 4); // nValues, filterLength + vector arr { 1, 2, 3, 4 }; + vector expected { 2, 2, 3, 3 }; + HistogramFilter::medianFilter(hf, arr); + BOOST_TEST(arr == expected, tt::per_element()); +} + +BOOST_AUTO_TEST_CASE(histogram_median_simple_3_4) +{ + HistogramFilter hf(5, 3); // nValues, filterLength + vector arr { 3, 1, 0, 2 }; + vector expected { 1, 1, 1, 0 }; + HistogramFilter::medianFilter(hf, arr); + BOOST_TEST(arr == expected, tt::per_element()); +} + +BOOST_AUTO_TEST_CASE(histogram_median_simple_5_4) +{ + HistogramFilter hf(5, 5); + vector arr { 3, 1, 0, 2 }; + vector expected { 1, 1, 1, 1 }; + HistogramFilter::medianFilter(hf, arr); + BOOST_TEST(arr == expected, tt::per_element()); +} + +BOOST_AUTO_TEST_CASE(histogram_median_order_1) +{ + HistogramFilter hf(4, 1); + vector arr { 3, 1, 0, 0 }; + vector expected = arr; + HistogramFilter::medianFilter(hf, arr); + BOOST_TEST(arr == expected, tt::per_element()); +} + +BOOST_AUTO_TEST_CASE(histogram_median_n_1) +{ + HistogramFilter hf(3, 6); + vector arr { 1 }; + vector expected { 1 }; + HistogramFilter::medianFilter(hf, arr); + BOOST_TEST(arr == expected, tt::per_element()); +} + +BOOST_AUTO_TEST_CASE(histogram_mode_simple_3) +{ + HistogramFilter hf(5, 3); // nValues, filterLength + vector arr { 1, 2, 2 }; + vector expected { 1, 2, 2 }; + HistogramFilter::modalFilter(hf, arr); + BOOST_TEST(arr == expected, tt::per_element()); +} + +BOOST_AUTO_TEST_CASE(histogram_mode_simple_4) +{ + HistogramFilter hf(5, 4); // nValues, filterLength + vector arr { 1, 2, 2, 4 }; + vector expected { 2, 2, 2, 2 }; + HistogramFilter::modalFilter(hf, arr); + BOOST_TEST(arr == expected, tt::per_element()); +} + +BOOST_AUTO_TEST_CASE(histogram_mode_simple_3_4) +{ + HistogramFilter hf(5, 3); // nValues, filterLength + vector arr { 3, 1, 0, 0 }; + vector expected { 1, 0, 0, 0 }; + HistogramFilter::modalFilter(hf, arr); + BOOST_TEST(arr == expected, tt::per_element()); +} + +BOOST_AUTO_TEST_CASE(histogram_mode_simple_5_4) +{ + HistogramFilter hf(5, 5); + vector arr { 3, 1, 0, 0 }; + vector expected { 0, 0, 0, 0 }; + HistogramFilter::modalFilter(hf, arr); + BOOST_TEST(arr == expected, tt::per_element()); +} + +BOOST_AUTO_TEST_CASE(histogram_mode_order_1) +{ + HistogramFilter hf(4, 1); + vector arr { 3, 1, 0, 0 }; + vector expected = arr; + HistogramFilter::modalFilter(hf, arr); + BOOST_TEST(arr == expected, tt::per_element()); +} + +BOOST_AUTO_TEST_CASE(histogram_mode_n_1) +{ + HistogramFilter hf(3, 6); + vector arr { 1 }; + vector expected { 1 }; + HistogramFilter::modalFilter(hf, arr); + BOOST_TEST(arr == expected, tt::per_element()); +} + BOOST_AUTO_TEST_CASE(peakpick_nearest_2_1) { Peak pp(1);