Mercurial > hg > svcore

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/base/MovingMedian.h	Wed Nov 14 14:21:53 2018 +0000
@@ -0,0 +1,235 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Sonic Visualiser
+    An audio file viewer and annotation editor.
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2007-2015 Particular Programs Ltd,
+    copyright 2018 Queen Mary University of London.
+
+    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 SV_MOVING_MEDIAN_H
+#define SV_MOVING_MEDIAN_H
+
+#include <bqvec/Allocators.h>
+#include <bqvec/VectorOps.h>
+
+#include <algorithm>
+#include <iostream>
+#include <stdexcept>
+
+/**
+ * Obtain the median (or other percentile) of a moving window across a
+ * time series. Construct the MovingMedian object, then push() each
+ * new value in the time series and get() the median of the most
+ * recent window. The size of the window, and the percentile
+ * calculated, can both be changed after construction.
+ *
+ * Note that for even-sized windows, the "median" is taken to be the
+ * value at the start of the second half when sorted, e.g. for size 4,
+ * the element at index 2 (zero-based) in the sorted window.
+ *
+ * Not thread-safe.
+ */
+template <typename T>
+class MovingMedian
+{
+public:
+    MovingMedian(int size, double percentile = 50.f) :
+        m_size(size),
+        m_percentile(percentile) {
+        if (size < 1) throw std::logic_error("size must be >= 1");
+        m_frame = breakfastquay::allocate_and_zero<T>(size);
+	m_sorted = breakfastquay::allocate_and_zero<T>(size);
+        calculateIndex();
+    }
+
+    ~MovingMedian() {
+        breakfastquay::deallocate(m_frame);
+        breakfastquay::deallocate(m_sorted);
+    }
+
+    MovingMedian(const MovingMedian &) =delete;
+    MovingMedian &operator=(const MovingMedian &) =delete;
+
+    void setPercentile(double p) {
+        m_percentile = p;
+        calculateIndex();
+    }
+
+    void push(T value) {
+        if (value != value) {
+            std::cerr << "WARNING: MovingMedian: NaN encountered" << std::endl;
+            value = T();
+        }
+	drop(m_frame[0]);
+        breakfastquay::v_move(m_frame, m_frame+1, m_size-1);
+	m_frame[m_size-1] = value;
+	put(value);
+    }
+
+    T get() const {
+	return m_sorted[m_index];
+    }
+
+    int size() const {
+        return m_size;
+    }
+
+    void reset() {
+        breakfastquay::v_zero(m_frame, m_size);
+        breakfastquay::v_zero(m_sorted, m_size);
+    }
+
+    void resize(int target) {
+        if (target == m_size) return;
+        int diff = std::abs(target - m_size);
+        if (target > m_size) { // grow
+            // we don't want to change the median, so fill spaces with it
+            T fillValue = get();
+            m_frame = breakfastquay::reallocate(m_frame, m_size, target);
+            breakfastquay::v_move(m_frame + diff, m_frame, m_size);
+            breakfastquay::v_set(m_frame, fillValue, diff);
+            m_sorted = breakfastquay::reallocate(m_sorted, m_size, target);
+            for (int sz = m_size + 1; sz <= target; ++sz) {
+                put(m_sorted, sz, fillValue);
+            }
+        } else { // shrink
+            for (int i = 0; i < diff; ++i) {
+                drop(m_sorted, m_size - i, m_frame[i]);
+            }
+            m_sorted = breakfastquay::reallocate(m_sorted, m_size, target);
+            breakfastquay::v_move(m_frame, m_frame + diff, target);
+            m_frame = breakfastquay::reallocate(m_frame, m_size, target);
+        }
+
+        m_size = target;
+        calculateIndex();
+    }
+
+    void checkIntegrity() const {
+        check();
+    }
+
+private:
+    int m_size;
+    double m_percentile;
+    int m_index;
+    T *m_frame;
+    T *m_sorted;
+
+    void calculateIndex() {
+        m_index = int((m_size * m_percentile) / 100.f);
+        if (m_index >= m_size) m_index = m_size-1;
+        if (m_index < 0) m_index = 0;
+    }
+
+    void put(T value) {
+        put(m_sorted, m_size, value);
+    }
+
+    static void put(T *const sorted, int size, T value) {
+
+        // precondition: sorted points to size-1 sorted values,
+        // followed by an unused slot (i.e. only the first size-1
+        // values of sorted are actually sorted)
+        //
+        // postcondition: sorted points to size sorted values
+
+	T *ptr = std::lower_bound(sorted, sorted + size - 1, value);
+        breakfastquay::v_move(ptr + 1, ptr, int(sorted + size - ptr) - 1);
+	*ptr = value;
+    }
+
+    void drop(T value) {
+        drop(m_sorted, m_size, value);
+    }
+
+    static void drop(T *const sorted, int size, T value) {
+
+        // precondition: sorted points to size sorted values, one of
+        // which is value
+        //
+        // postcondition: sorted points to size-1 sorted values,
+        // followed by a slot that has been reset to default value
+        // (i.e. only the first size-1 values of sorted are actually
+        // sorted)
+
+	T *ptr = std::lower_bound(sorted, sorted + size, value);
+	if (*ptr != value) {
+            throw std::logic_error
+                ("MovingMedian::drop: value being dropped is not in array");
+        }
+        breakfastquay::v_move(ptr, ptr + 1, int(sorted + size - ptr) - 1);
+        sorted[size-1] = T();
+    }
+
+    void check() const {
+        bool good = true;
+        for (int i = 1; i < m_size; ++i) {
+            if (m_sorted[i] < m_sorted[i-1]) {
+                std::cerr << "ERROR: MovingMedian::checkIntegrity: "
+                          << "mis-ordered elements in sorted array starting "
+                          << "at index " << i << std::endl;
+                good = false;
+                break;
+            }
+        }
+        for (int i = 0; i < m_size; ++i) {
+            bool found = false;
+            for (int j = 0; j < m_size; ++j) {
+                if (m_sorted[j] == m_frame[i]) {
+                    found = true;
+                    break;
+                }
+            }
+            if (!found) {
+                std::cerr << "ERROR: MovingMedian::checkIntegrity: "
+                          << "element in frame at index " << i
+                          << " not found in sorted array" << std::endl;
+                good = false;
+                break;
+            }
+        }
+        for (int i = 0; i < m_size; ++i) {
+            bool found = false;
+            for (int j = 0; j < m_size; ++j) {
+                if (m_sorted[i] == m_frame[j]) {
+                    found = true;
+                    break;
+                }
+            }
+            if (!found) {
+                std::cerr << "ERROR: MovingMedian::checkIntegrity: "
+                          << "element in sorted array at index " << i
+                          << " not found in source frame" << std::endl;
+                good = false;
+                break;
+            }
+        }
+        if (!good) {
+            std::cerr << "Frame contains:" << std::endl;
+            std::cerr << "[ ";
+            for (int j = 0; j < m_size; ++j) {
+                std::cerr << m_frame[j] << " ";
+            }
+            std::cerr << "]" << std::endl;
+            std::cerr << "Sorted array contains:" << std::endl;
+            std::cerr << "[ ";
+            for (int j = 0; j < m_size; ++j) {
+                std::cerr << m_sorted[j] << " ";
+            }
+            std::cerr << "]" << std::endl;
+            throw std::logic_error("MovingMedian failed integrity check");
+        }
+    }
+};
+
+#endif
+
--- a/base/RealTimeSV.cpp	Wed Nov 07 15:46:36 2018 +0000
+++ b/base/RealTimeSV.cpp	Wed Nov 14 14:21:53 2018 +0000
@@ -100,8 +100,7 @@
     int year = 0, month = 0, day = 0, hour = 0, minute = 0;
     double second = 0.0;

-    char *loc = setlocale(LC_NUMERIC, 0);
-    (void)setlocale(LC_NUMERIC, "C"); // avoid strtod expecting ,-separator in DE
+    char *formerLoc = setlocale(LC_NUMERIC, "C"); // avoid strtod expecting ,-separator in DE

     int i = 0;

@@ -169,7 +168,7 @@

     t = t + fromSeconds(second);

-    setlocale(LC_NUMERIC, loc);
+    setlocale(LC_NUMERIC, formerLoc);

     if (negative) {
         return -t;
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/base/test/TestMovingMedian.h	Wed Nov 14 14:21:53 2018 +0000
@@ -0,0 +1,170 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    Sonic Visualiser
+    An audio file viewer and annotation editor.
+    Centre for Digital Music, Queen Mary, University of London.
+
+    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 TEST_MOVING_MEDIAN_H
+#define TEST_MOVING_MEDIAN_H
+
+#include "../MovingMedian.h"
+
+#include <QObject>
+#include <QtTest>
+#include <QDir>
+
+#include <iostream>
+
+using namespace std;
+
+class TestMovingMedian : public QObject
+{
+    Q_OBJECT
+
+    template <typename T>
+    void checkExpected(const vector<T> &output,
+                       const vector<T> &expected) {
+        if (output.size() != expected.size()) {
+            std::cerr << "ERROR: output array size " << output.size()
+                      << " differs from expected size " << expected.size()
+                      << std::endl;
+        }
+        for (int i = 0; i < int(output.size()); ++i) {
+            if (output[i] != expected[i]) {
+                std::cerr << "ERROR: Value at index " << i
+                          << " in output array differs from expected"
+                          << std::endl;
+                std::cerr << "Output:   ";
+                for (auto v: output) std::cerr << v << " ";
+                std::cerr << "\nExpected: ";
+                for (auto v: expected) std::cerr << v << " ";
+                std::cerr << std::endl;
+                break;
+            }
+        }
+        QCOMPARE(output, expected);
+    }
+
+    template <typename T>
+    void testFixed(int n,
+                   const vector<T> &input,
+                   const vector<T> &expected,
+                   double percentile = 50.0) {
+        vector<T> output;
+        MovingMedian<T> mm(n, percentile);
+        for (auto v: input) {
+            mm.push(v);
+            mm.checkIntegrity();
+            output.push_back(mm.get());
+        }
+        mm.checkIntegrity();
+        checkExpected<T>(output, expected);
+    }
+
+private slots:
+
+    void empty() {
+        MovingMedian<double> mm(3);
+        QCOMPARE(mm.get(), 0.0);
+    }
+
+    void zeros() {
+        vector<double> input { 0.0, 0.0, 0.0, 0.0, 0.0 };
+        vector<double> expected { 0.0, 0.0, 0.0, 0.0, 0.0 };
+        testFixed<double>(3, input, expected);
+    }
+
+    void ascending() {
+        vector<double> input { 1.0, 2.0, 3.0, 4.0, 5.0 };
+        vector<double> expected { 0.0, 1.0, 2.0, 3.0, 4.0 };
+        testFixed<double>(3, input, expected);
+    }
+
+    void ascendingInt() {
+        vector<int> input { 1, 2, 3, 4, 5 };
+        vector<int> expected { 0, 1, 2, 3, 4 };
+        testFixed<int>(3, input, expected);
+    }
+
+    void descending() {
+        vector<double> input { 5.0, 4.0, 3.0, 2.0, 1.0 };
+        vector<double> expected { 0.0, 4.0, 4.0, 3.0, 2.0 };
+        testFixed<double>(3, input, expected);
+    }
+
+    void descendingInt() {
+        vector<int> input { 5, 4, 3, 2, 1 };
+        vector<int> expected { 0, 4, 4, 3, 2 };
+        testFixed<int>(3, input, expected);
+    }
+
+    void duplicates() {
+        vector<double> input { 2.0, 2.0, 3.0, 4.0, 3.0 };
+        vector<double> expected { 0.0, 2.0, 2.0, 3.0, 3.0 };
+        testFixed<double>(3, input, expected);
+    }
+
+    void percentile10() {
+        vector<double> input { 1.0, 2.0, 3.0, 4.0, 5.0 };
+        vector<double> expected { 0.0, 0.0, 1.0, 2.0, 3.0 };
+        testFixed<double>(3, input, expected, 10);
+    }
+
+    void percentile90() {
+        vector<double> input { 1.0, 2.0, 3.0, 4.0, 5.0 };
+        vector<double> expected { 1.0, 2.0, 3.0, 4.0, 5.0 };
+        testFixed<double>(3, input, expected, 90);
+    }
+
+    void even() {
+        vector<double> input { 5.0, 4.0, 3.0, 2.0, 1.0 };
+        vector<double> expected { 0.0, 4.0, 4.0, 4.0, 3.0 };
+        testFixed<double>(4, input, expected);
+    }
+
+    void growing() {
+        vector<double> input { 2.0, 4.0, 3.0, 2.5, 2.5, 3.0, 1.0, 2.0, 1.0, 0.0 };
+        vector<double> expected { 2.0, 4.0, 4.0, 3.0, 2.5, 2.5, 2.5, 2.5, 2.0, 1.0 };
+        vector<double> output;
+        MovingMedian<double> mm(1);
+        for (int i = 0; i < int(input.size()); ++i) {
+            // sizes 1, 1, 2, 2, 3, 3, 4, 4, 5, 5
+            int sz = i/2 + 1;
+            mm.resize(sz);
+            QCOMPARE(mm.size(), sz);
+            mm.push(input[i]);
+            mm.checkIntegrity();
+            output.push_back(mm.get());
+        }
+        mm.checkIntegrity();
+        checkExpected<double>(output, expected);
+    }
+
+    void shrinking() {
+        vector<double> input { 2.0, 4.0, 3.0, 2.5, 2.5, 3.0, 1.0, 2.0, 1.0, 0.0 };
+        vector<double> expected { 0.0, 0.0, 3.0, 3.0, 2.5, 2.5, 3.0, 2.0, 1.0, 0.0 };
+        vector<double> output;
+        MovingMedian<double> mm(99);
+        for (int i = 0; i < int(input.size()); ++i) {
+            // sizes 5, 5, 4, 4, 3, 3, 2, 2, 1, 1
+            int sz = 5 - i/2;
+            mm.resize(sz);
+            QCOMPARE(mm.size(), sz);
+            mm.push(input[i]);
+            mm.checkIntegrity();
+            output.push_back(mm.get());
+        }
+        mm.checkIntegrity();
+        checkExpected<double>(output, expected);
+    }
+};
+
+#endif
--- a/base/test/files.pri	Wed Nov 07 15:46:36 2018 +0000
+++ b/base/test/files.pri	Wed Nov 14 14:21:53 2018 +0000
@@ -1,9 +1,10 @@
 TEST_HEADERS = \
 	     TestColumnOp.h \
 	     TestLogRange.h \
-	     TestRangeMapper.h \
+	     TestMovingMedian.h \
 	     TestOurRealTime.h \
 	     TestPitch.h \
+	     TestRangeMapper.h \
 	     TestScaleTickIntervals.h \
 	     TestStringBits.h \
 	     TestVampRealTime.h
--- a/base/test/svcore-base-test.cpp	Wed Nov 07 15:46:36 2018 +0000
+++ b/base/test/svcore-base-test.cpp	Wed Nov 14 14:21:53 2018 +0000
@@ -19,6 +19,7 @@
 #include "TestOurRealTime.h"
 #include "TestVampRealTime.h"
 #include "TestColumnOp.h"
+#include "TestMovingMedian.h"

 #include <QtTest>

@@ -72,6 +73,11 @@
         if (QTest::qExec(&t, argc, argv) == 0) ++good;
         else ++bad;
     }
+    {
+        TestMovingMedian t;
+        if (QTest::qExec(&t, argc, argv) == 0) ++good;
+        else ++bad;
+    }

     if (bad > 0) {
         SVCERR << "\n********* " << bad << " test suite(s) failed!\n" << endl;
--- a/data/model/FFTModel.cpp	Wed Nov 07 15:46:36 2018 +0000
+++ b/data/model/FFTModel.cpp	Wed Nov 14 14:21:53 2018 +0000
@@ -20,6 +20,7 @@
 #include "base/Pitch.h"
 #include "base/HitCount.h"
 #include "base/Debug.h"
+#include "base/MovingMedian.h"

 #include <algorithm>

@@ -388,7 +389,7 @@
 FFTModel::getPeaks(PeakPickType type, int x, int ymin, int ymax) const
 {
     Profiler profiler("FFTModel::getPeaks");
-
+
     FFTModel::PeakLocationSet peaks;
     if (!isOK()) return peaks;

@@ -429,13 +430,14 @@

     sv_samplerate_t sampleRate = getSampleRate();

-    deque<float> window;
     vector<int> inrange;
-    float dist = 0.5;
+    double dist = 0.5;

     int medianWinSize = getPeakPickWindowSize(type, sampleRate, ymin, dist);
     int halfWin = medianWinSize/2;

+    MovingMedian<float> window(medianWinSize);
+
     int binmin;
     if (ymin > halfWin) binmin = ymin - halfWin;
     else binmin = 0;
@@ -450,26 +452,21 @@

         float value = values[bin];

-        window.push_back(value);
-
         // so-called median will actually be the dist*100'th percentile
         medianWinSize = getPeakPickWindowSize(type, sampleRate, bin, dist);
         halfWin = medianWinSize/2;

-        while ((int)window.size() > medianWinSize) {
-            window.pop_front();
-        }
-
-        int actualSize = int(window.size());
+        int actualSize = std::min(medianWinSize, bin - binmin + 1);
+        window.resize(actualSize);
+        window.setPercentile(dist * 100.0);
+        window.push(value);

         if (type == MajorPitchAdaptivePeaks) {
             if (ymax + halfWin < nv) binmax = ymax + halfWin;
             else binmax = nv - 1;
         }

-        deque<float> sorted(window);
-        sort(sorted.begin(), sorted.end());
-        float median = sorted[int(float(sorted.size()) * dist)];
+        float median = window.get();

         int centrebin = 0;
         if (bin > actualSize/2) centrebin = bin - actualSize/2;
@@ -510,9 +507,9 @@

 int
 FFTModel::getPeakPickWindowSize(PeakPickType type, sv_samplerate_t sampleRate,
-                                int bin, float &percentile) const
+                                int bin, double &dist) const
 {
-    percentile = 0.5;
+    dist = 0.5; // dist is percentile / 100.0
     if (type == MajorPeaks) return 10;
     if (bin == 0) return 3;

@@ -521,10 +518,34 @@

     int hibin = int(lrint((hifreq * m_fftSize) / sampleRate));
     int medianWinSize = hibin - bin;
-    if (medianWinSize < 3) medianWinSize = 3;

-    percentile = 0.5f + float(binfreq / sampleRate);
+    if (medianWinSize < 3) {
+        medianWinSize = 3;
+    }

+    // We want to avoid the median window size changing too often, as
+    // it requires a reallocation. So snap to a nearby round number.
+
+    if (medianWinSize > 20) {
+        medianWinSize = (1 + medianWinSize / 10) * 10;
+    }
+    if (medianWinSize > 200) {
+        medianWinSize = (1 + medianWinSize / 100) * 100;
+    }
+    if (medianWinSize > 2000) {
+        medianWinSize = (1 + medianWinSize / 1000) * 1000;
+    }
+    if (medianWinSize > 20000) {
+        medianWinSize = 20000;
+    }
+
+    if (medianWinSize < 100) {
+        dist = 1.0 - (4.0 / medianWinSize);
+    } else {
+        dist = 1.0 - (8.0 / medianWinSize);
+    }
+    if (dist < 0.5) dist = 0.5;
+
     return medianWinSize;
 }
--- a/data/model/FFTModel.h	Wed Nov 07 15:46:36 2018 +0000
+++ b/data/model/FFTModel.h	Wed Nov 14 14:21:53 2018 +0000
@@ -156,7 +156,7 @@
     mutable breakfastquay::FFT m_fft;

     int getPeakPickWindowSize(PeakPickType type, sv_samplerate_t sampleRate,
-                              int bin, float &percentile) const;
+                              int bin, double &dist) const;

     std::pair<sv_frame_t, sv_frame_t> getSourceSampleRange(int column) const {
         sv_frame_t startFrame = m_windowIncrement * sv_frame_t(column);