Mercurial > hg > devuvuzelator

changeset 4:d90abfa9585a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
* Adaptive method
author Chris Cannam
date Fri, 11 Jun 2010 15:38:44 +0100
parents 8b79175c9f02
children 45bcfa3d5da7
files devuvuzelator-ladspa.cpp devuvuzelator-vst.cpp devuvuzelator.cpp median.h
diffstat 4 files changed, 261 insertions(+), 128 deletions(-) [+]
line wrap: on
line diff
--- a/devuvuzelator-ladspa.cpp	Fri Jun 11 11:45:41 2010 +0100
+++ b/devuvuzelator-ladspa.cpp	Fri Jun 11 15:38:44 2010 +0100
@@ -5,7 +5,10 @@
 #include <iostream>
 #include <cmath>
 
-#define FFTSIZE 1024
+#include "median.h"
+
+#define FFTSIZE 4096
+#define WINSIZE 1024
 
 class Devuvuzelator
 {
@@ -19,13 +22,12 @@
     enum {
         InputPort     = 0,
 	OutputPort    = 1,
-        LowPort       = 2,
-        HighPort      = 3,
-        FundamentalPort = 4,
-        BandwidthPort = 5,
-        HarmonicsPort = 6,
-        ReductionPort = 7,
-	PortCount     = 8,
+        LatencyPort   = 2,
+        FundamentalPort = 3,
+        BandwidthPort = 4,
+        HarmonicsPort = 5,
+        ReductionPort = 6,
+	PortCount     = 7,
     };
 
     static const char *const portNames[PortCount];
@@ -46,6 +48,7 @@
     void runImpl(unsigned long);
     void processFrame();
     void processSpectralFrame();
+    void updateParameters();
 
     static void fft(unsigned int n, bool inverse,
                     double *ri, double *ii, double *ro, double *io);
@@ -53,14 +56,20 @@
     int m_sampleRate;
     float *m_input;
     float *m_output;
-    float *m_low;
-    float *m_high;
-    float *m_fundamental;
-    float *m_bandwidth;
-    float *m_harmonics;
-    float *m_reduction;
+
+    float m_fundamental;
+    float m_bandwidth;
+    float m_harmonics;
+    float m_reduction;
+
+    float *m_platency;
+    float *m_pfundamental;
+    float *m_pbandwidth;
+    float *m_pharmonics;
+    float *m_preduction;
 
     const int m_fftsize;
+    const int m_winsize;
     const int m_increment;
     int m_fill;
     int m_read;
@@ -69,6 +78,7 @@
     double *m_real;
     double *m_imag;
     double *m_window;
+    MedianFilter<double> **m_medians;
 };
 
 const char *const
@@ -76,12 +86,11 @@
 {
     "Input",
     "Output",
-    "Low threshold (dB)",
-    "High threshold (dB)",
+    "latency",
     "Fundamental frequency (Hz)",
     "Bandwidth of fundamental (Hz)",
     "Number of partials",
-    "Reduction (dB)",
+    "Reduction (%)",
 };
 
 const LADSPA_PortDescriptor 
@@ -89,8 +98,7 @@
 {
     LADSPA_PORT_INPUT  | LADSPA_PORT_AUDIO,
     LADSPA_PORT_OUTPUT | LADSPA_PORT_AUDIO,
-    LADSPA_PORT_INPUT  | LADSPA_PORT_CONTROL,
-    LADSPA_PORT_INPUT  | LADSPA_PORT_CONTROL,
+    LADSPA_PORT_OUTPUT | LADSPA_PORT_CONTROL,
     LADSPA_PORT_INPUT  | LADSPA_PORT_CONTROL,
     LADSPA_PORT_INPUT  | LADSPA_PORT_CONTROL,
     LADSPA_PORT_INPUT  | LADSPA_PORT_CONTROL,
@@ -102,18 +110,15 @@
 {
     { 0, 0, 0 },
     { 0, 0, 0 },
-    { LADSPA_HINT_DEFAULT_MIDDLE |
-      LADSPA_HINT_BOUNDED_BELOW | LADSPA_HINT_BOUNDED_ABOVE, -80, 0 },
-    { LADSPA_HINT_DEFAULT_HIGH |
-      LADSPA_HINT_BOUNDED_BELOW | LADSPA_HINT_BOUNDED_ABOVE, -80, 0 },
+    { 0, 0, 0 },
     { LADSPA_HINT_DEFAULT_LOW |
-      LADSPA_HINT_BOUNDED_BELOW | LADSPA_HINT_BOUNDED_ABOVE, 110, 550 },
+      LADSPA_HINT_BOUNDED_BELOW | LADSPA_HINT_BOUNDED_ABOVE, 50, 770 },
     { LADSPA_HINT_DEFAULT_MIDDLE |
       LADSPA_HINT_BOUNDED_BELOW | LADSPA_HINT_BOUNDED_ABOVE, 20, 100 },
     { LADSPA_HINT_DEFAULT_MIDDLE | LADSPA_HINT_INTEGER |
       LADSPA_HINT_BOUNDED_BELOW | LADSPA_HINT_BOUNDED_ABOVE, 0, 6 },
     { LADSPA_HINT_DEFAULT_MIDDLE |
-      LADSPA_HINT_BOUNDED_BELOW | LADSPA_HINT_BOUNDED_ABOVE, 0, 20 },
+      LADSPA_HINT_BOUNDED_BELOW | LADSPA_HINT_BOUNDED_ABOVE, 0, 100 },
 };
 
 const LADSPA_Properties
@@ -154,22 +159,37 @@
     m_sampleRate(sampleRate),
     m_input(0),
     m_output(0),
-    m_low(0),
-    m_high(0),
+    m_platency(0),
+    m_pfundamental(0),
+    m_pbandwidth(0),
+    m_pharmonics(0),
+    m_preduction(0),
     m_fftsize(FFTSIZE),
-    m_increment(m_fftsize/4),
+    m_winsize(WINSIZE),
+    m_increment(m_winsize/2),
     m_fill(0),
     m_read(0)
 {
-    m_buffer = new float[m_fftsize];
-    m_outacc = new float[m_fftsize * 2];
+    m_buffer = new float[m_winsize];
+    m_outacc = new float[m_winsize * 2];
+    m_window = new double[m_winsize];
     m_real = new double[m_fftsize];
     m_imag = new double[m_fftsize];
-    m_window = new double[m_fftsize];
+    m_medians = new MedianFilter<double> *[m_fftsize/2+1];
 
-    for (int i = 0; i < m_fftsize; ++i) {
-        m_window[i] = 0.5 - 0.5 * cos(2 * M_PI * i / m_fftsize);
+    for (int i = 0; i < m_winsize; ++i) {
+        m_window[i] = 0.5 - 0.5 * cos(2 * M_PI * i / m_winsize);
     }
+    for (int i = 0; i < m_fftsize/2+1; ++i) {
+        m_medians[i] = 0;
+    }
+
+    m_fundamental = 230;
+    m_bandwidth = 60;
+    m_harmonics = 3;
+    m_reduction = 50;
+
+    std::cerr << "note: latency = " << (float(m_winsize) / m_sampleRate)*1000 << " msec" << std::endl;
 
     reset();
 }
@@ -181,6 +201,10 @@
     delete[] m_real;
     delete[] m_imag;
     delete[] m_window;
+    for (int i = 0; i < m_fftsize/2+1; ++i) {
+        delete m_medians[i];
+    }
+    delete[] m_medians;
 }
     
 LADSPA_Handle
@@ -199,12 +223,11 @@
     float **ports[PortCount] = {
 	&devuvu->m_input,
 	&devuvu->m_output,
-        &devuvu->m_low,
-        &devuvu->m_high,
-        &devuvu->m_fundamental,
-        &devuvu->m_bandwidth,
-        &devuvu->m_harmonics,
-        &devuvu->m_reduction,
+        &devuvu->m_platency,
+        &devuvu->m_pfundamental,
+        &devuvu->m_pbandwidth,
+        &devuvu->m_pharmonics,
+        &devuvu->m_preduction,
     };
 
     *ports[port] = (float *)location;
@@ -239,20 +262,34 @@
 void
 Devuvuzelator::reset()
 {
-    for (int i = 0; i < m_fftsize; ++i) {
+    for (int i = 0; i < m_winsize; ++i) {
         m_buffer[i] = 0.f;
     }
-    for (int i = 0; i < m_fftsize*2; ++i) {
+    for (int i = 0; i < m_winsize*2; ++i) {
         m_outacc[i] = 0.f;
     }
     m_fill = 0;
     m_read = 0;
+    for (int i = 0; i < m_fftsize/2+1; ++i) {
+        if (m_medians[i]) m_medians[i]->reset();
+    }
+}
+
+void
+Devuvuzelator::updateParameters()
+{
+    if (m_platency) *m_platency = m_winsize;
+    if (m_pfundamental) m_fundamental = *m_pfundamental;
+    if (m_pbandwidth) m_bandwidth = *m_pbandwidth;
+    if (m_pharmonics) m_harmonics = *m_pharmonics;
+    if (m_preduction) m_reduction = *m_preduction;
 }
 
 void
 Devuvuzelator::runImpl(unsigned long sampleCount)
 {
     if (!m_input || !m_output) return;
+    updateParameters();
 
     int ii = 0;
     int oi = 0;
@@ -260,21 +297,21 @@
 
     while (ii < sc) {
 
-        m_output[oi++] = m_outacc[m_read++] / 1.5f;
+        m_output[oi++] = m_outacc[m_read++];
 
-        if (m_fill == m_fftsize) {
+        if (m_fill == m_winsize) {
 
             processFrame();
 
-            for (int j = m_increment; j < m_fftsize; ++j) {
+            for (int j = m_increment; j < m_winsize; ++j) {
                 m_buffer[j - m_increment] = m_buffer[j];
             }
 
-            for (int j = m_increment; j < m_fftsize*2; ++j) {
+            for (int j = m_increment; j < m_winsize*2; ++j) {
                 m_outacc[j - m_increment] = m_outacc[j];
             }
 
-            for (int j = m_fftsize*2 - m_increment; j < m_fftsize*2; ++j) {
+            for (int j = m_winsize*2 - m_increment; j < m_winsize*2; ++j) {
                 m_outacc[j] = 0.f;
             }
 
@@ -290,9 +327,14 @@
 Devuvuzelator::processFrame()
 {
     double *frame = (double *)alloca(m_fftsize * sizeof(double));
-    int ix = m_fftsize/2;
     for (int i = 0; i < m_fftsize; ++i) {
-        frame[ix++] = m_buffer[i] * m_window[i];
+        frame[i] = 0.0;
+    }
+
+    int ix = m_fftsize - m_winsize/2;
+    while (ix < 0) ix += m_fftsize;
+    for (int i = 0; i < m_winsize; ++i) {
+        frame[ix++] += m_buffer[i] * m_window[i];
         if (ix == m_fftsize) ix = 0;
     }
 
@@ -308,9 +350,10 @@
     double *spare = (double *)alloca(m_fftsize * sizeof(double));
     fft(m_fftsize, true, m_real, m_imag, frame, spare);
 
-    ix = m_fftsize/2;
-    for (int i = 0; i < m_fftsize; ++i) {
-        m_outacc[m_fftsize + i] += frame[ix++] * m_window[i];
+    ix = m_fftsize - m_winsize/2;
+    while (ix < 0) ix += m_fftsize;
+    for (int i = 0; i < m_winsize; ++i) {
+        m_outacc[m_winsize + i] += frame[ix++];
         if (ix == m_fftsize) ix = 0;
     }
 }
--- a/devuvuzelator-vst.cpp	Fri Jun 11 11:45:41 2010 +0100
+++ b/devuvuzelator-vst.cpp	Fri Jun 11 15:38:44 2010 +0100
@@ -12,18 +12,19 @@
 #define snprintf _snprintf
 #define alloca _alloca
 
-#define FFTSIZE 1024
+#define FFTSIZE 4096
+#define WINSIZE 1024
+
+#include "median.h"
 
 class Devuvuzelator : public AudioEffect
 {
     enum {
-        LowParam       = 0,
-        HighParam      = 1,
-        FundamentalParam = 2,
-        BandwidthParam = 3,
-        HarmonicsParam = 4,
-        ReductionParam = 5,
-        NumParams      = 6
+        FundamentalParam = 1,
+        BandwidthParam = 2,
+        HarmonicsParam = 3,
+        ReductionParam = 4,
+        NumParams      = 5
     };
 
 public:
@@ -78,6 +79,7 @@
     float m_reduction;
 
     const int m_fftsize;
+    const int m_winsize;
     const int m_increment;
     int m_fill;
     int m_read;
@@ -86,6 +88,7 @@
     double *m_real;
     double *m_imag;
     double *m_window;
+    MedianFilter<double> **m_medians;
 };
 
 // VST params 0->1
@@ -94,12 +97,10 @@
 Devuvuzelator::setParameter(VstInt32 index, float value)
 {
     switch (index) {
-    case 0: m_low = -80 + 80 * value; break;
-    case 1: m_high = -80 + 80 * value; break;
-    case 2: m_fundamental = 110 + 440 * value; break;
-    case 3: m_bandwidth = 20 + 80 * value; break;
-    case 4: m_harmonics = int(value * 6 + 0.5); break;
-    case 5: m_reduction = 20 * value; break;
+    case 0: m_fundamental = 50 + 720 * value; break;
+    case 1: m_bandwidth = 20 + 80 * value; break;
+    case 2: m_harmonics = int(value * 6 + 0.5); break;
+    case 3: m_reduction = 100 * value; break;
     }
 }
 
@@ -107,14 +108,12 @@
 Devuvuzelator::getParameter(VstInt32 index)
 {
     switch (index) {
-    case 0: return (m_low + 80) / 80;
-    case 1: return (m_high + 80) / 80;
-    case 2: return (m_fundamental - 110) / 440;
-    case 3: return (m_bandwidth - 20) / 80;
-    case 4: return (m_harmonics / 6.f);
-    case 5: return m_reduction / 20;
+    case 0: return (m_fundamental - 50) / 720;
+    case 1: return (m_bandwidth - 20) / 80;
+    case 2: return (m_harmonics / 6.f);
+    case 3: return m_reduction / 100;
     }
-	return 0;
+    return 0;
 }
 
 // NB! The max name length for VST parameter names, labels
@@ -125,12 +124,10 @@
 Devuvuzelator::getParameterLabel(VstInt32 index, char *label)
 {
     const char *units[NumParams] = {
-        "dB",
-        "dB",
         "Hz",
         "Hz",
         "",
-        "dB",
+        "%",
     };
     
     vst_strncpy(label, units[index], kVstMaxParamStrLen);
@@ -139,22 +136,18 @@
 void
 Devuvuzelator::getParameterDisplay(VstInt32 index, char *label)
 {
-	switch (index) {
-		case 0: snprintf(label, kVstMaxParamStrLen, "%f", m_low); break;
-		case 1: snprintf(label, kVstMaxParamStrLen, "%f", m_high); break;
-		case 2: snprintf(label, kVstMaxParamStrLen, "%f", m_fundamental); break;
-		case 3: snprintf(label, kVstMaxParamStrLen, "%f", m_bandwidth); break;
-		case 4: snprintf(label, kVstMaxParamStrLen, "%d", m_harmonics); break;
-		case 5: snprintf(label, kVstMaxParamStrLen, "%f", m_reduction); break;
-	}
+    switch (index) {
+    case 0: snprintf(label, kVstMaxParamStrLen, "%f", m_fundamental); break;
+    case 1: snprintf(label, kVstMaxParamStrLen, "%f", m_bandwidth); break;
+    case 2: snprintf(label, kVstMaxParamStrLen, "%d", m_harmonics); break;
+    case 3: snprintf(label, kVstMaxParamStrLen, "%f", m_reduction); break;
+    }
 }
 
 void
 Devuvuzelator::getParameterName(VstInt32 index, char *label)
 {
     const char *names[NumParams] = {
-        "Floor",
-        "Ceiling",
         "Pitch",
         "B/W",
         "Partials",
@@ -169,29 +162,29 @@
     m_sampleRate(0),
     m_input(0),
     m_output(0),
-    m_low(0),
-    m_high(0),
     m_fftsize(FFTSIZE),
     m_increment(m_fftsize/4),
     m_fill(0),
     m_read(0)
 {
-    m_buffer = new float[m_fftsize];
-    m_outacc = new float[m_fftsize * 2];
+    m_buffer = new float[m_winsize];
+    m_outacc = new float[m_winsize * 2];
+    m_window = new double[m_winsize];
     m_real = new double[m_fftsize];
     m_imag = new double[m_fftsize];
-    m_window = new double[m_fftsize];
+    m_medians = new MedianFilter<double> *[m_fftsize/2+1];
 
-    for (int i = 0; i < m_fftsize; ++i) {
-        m_window[i] = 0.5 - 0.5 * cos(2 * M_PI * i / m_fftsize);
+    for (int i = 0; i < m_winsize; ++i) {
+        m_window[i] = 0.5 - 0.5 * cos(2 * M_PI * i / m_winsize);
+    }
+    for (int i = 0; i < m_fftsize/2+1; ++i) {
+        m_medians[i] = 0;
     }
 
-    m_low = -40;
-    m_high = -20;
-    m_fundamental = 220;
+    m_fundamental = 230;
     m_bandwidth = 60;
     m_harmonics = 3;
-    m_reduction = 10;
+    m_reduction = 30;
     
     setUniqueID('qmvz');
     setNumInputs(1);
@@ -209,19 +202,26 @@
     delete[] m_real;
     delete[] m_imag;
     delete[] m_window;
+    for (int i = 0; i < m_fftsize/2+1; ++i) {
+        delete m_medians[i];
+    }
+    delete[] m_medians;
 }
 
 void
 Devuvuzelator::reset()
 {
-    for (int i = 0; i < m_fftsize; ++i) {
+    for (int i = 0; i < m_winsize; ++i) {
         m_buffer[i] = 0.f;
     }
-    for (int i = 0; i < m_fftsize*2; ++i) {
+    for (int i = 0; i < m_winsize*2; ++i) {
         m_outacc[i] = 0.f;
     }
     m_fill = 0;
     m_read = 0;
+    for (int i = 0; i < m_fftsize/2+1; ++i) {
+        if (m_medians[i]) m_medians[i]->reset();
+    }
 }
 
 void
@@ -231,25 +231,25 @@
 
     int ii = 0;
     int oi = 0;
-	const int sc = sampleCount;
+    const int sc = sampleCount;
 
     while (ii < sc) {
 
-        m_output[oi++] = m_outacc[m_read++] / 1.5f;
+        m_output[oi++] = m_outacc[m_read++];
 
-        if (m_fill == m_fftsize) {
+        if (m_fill == m_winsize) {
 
             processFrame();
 
-            for (int j = m_increment; j < m_fftsize; ++j) {
+            for (int j = m_increment; j < m_winsize; ++j) {
                 m_buffer[j - m_increment] = m_buffer[j];
             }
 
-            for (int j = m_increment; j < m_fftsize*2; ++j) {
+            for (int j = m_increment; j < m_winsize*2; ++j) {
                 m_outacc[j - m_increment] = m_outacc[j];
             }
 
-            for (int j = m_fftsize*2 - m_increment; j < m_fftsize*2; ++j) {
+            for (int j = m_winsize*2 - m_increment; j < m_winsize*2; ++j) {
                 m_outacc[j] = 0.f;
             }
 
@@ -265,9 +265,14 @@
 Devuvuzelator::processFrame()
 {
     double *frame = (double *)alloca(m_fftsize * sizeof(double));
-    int ix = m_fftsize/2;
     for (int i = 0; i < m_fftsize; ++i) {
-        frame[ix++] = m_buffer[i] * m_window[i];
+        frame[i] = 0.0;
+    }
+
+    int ix = m_fftsize - m_winsize/2;
+    while (ix < 0) ix += m_fftsize;
+    for (int i = 0; i < m_winsize; ++i) {
+        frame[ix++] += m_buffer[i] * m_window[i];
         if (ix == m_fftsize) ix = 0;
     }
 
@@ -283,9 +288,10 @@
     double *spare = (double *)alloca(m_fftsize * sizeof(double));
     fft(m_fftsize, true, m_real, m_imag, frame, spare);
 
-    ix = m_fftsize/2;
-    for (int i = 0; i < m_fftsize; ++i) {
-        m_outacc[m_fftsize + i] += frame[ix++] * m_window[i];
+    ix = m_fftsize - m_winsize/2;
+    while (ix < 0) ix += m_fftsize;
+    for (int i = 0; i < m_winsize; ++i) {
+        m_outacc[m_winsize + i] += frame[ix++];
         if (ix == m_fftsize) ix = 0;
     }
 }
--- a/devuvuzelator.cpp	Fri Jun 11 11:45:41 2010 +0100
+++ b/devuvuzelator.cpp	Fri Jun 11 15:38:44 2010 +0100
@@ -22,29 +22,25 @@
         int lowbin = 0.5 + (m_fftsize * lowfreq) / m_sampleRate;
         int highbin = 0.5 + (m_fftsize * highfreq) / m_sampleRate;
 
-//        std::cerr << "partials " << h << ": freqs " << lowfreq << "->"
-//                  << highfreq << ", bins " << lowbin << "->"  << highbin << std::endl;
+        for (int i = lowbin; i <= highbin; ++i) {
 
-        for (int i = lowbin; i <= highbin; ++i) {
-//            std::cerr << "bin " << i << " freq " << (m_sampleRate * i) / m_fftsize << std::endl;
-            ratios[i] = 1.0;
-            double db = 10 * log10(mags[i]);
-            if (db > m_low && db < m_high) {
-                double r = m_reduction;
-                ratios[i] = pow(10, -r / 10);
+            if (!m_medians[i]) {
+                const float filtersecs = 5.f;
+                int filterlen = (filtersecs * m_sampleRate) / m_increment;
+                m_medians[i] = new MedianFilter<double>(filterlen);
+            }
+
+            m_medians[i]->push(mags[i]);
+            double threshold = m_medians[i]->getAt(m_reduction);
+            if (mags[i] > threshold) {
+                double target = mags[i] - threshold;
+                ratios[i] = (target / mags[i]);
+            } else {
+                ratios[i] = 0.0;
             }
         }
     }
-/*
-    for (int i = 0; i < hs-1; ++i) {
-        if (ratios[i] == 1.0 && ratios[i+1] < 1.0) {
-            ratios[i] = (ratios[i+1] + 1) / 2;
-        } else if (ratios[i] < 1.0 && ratios[i+1] == 1.0) {
-            ratios[i+1] = (ratios[i] + 1) / 2;
-            ++i;
-        }
-    }
-*/
+
     for (int i = 0; i < hs; ++i) {
         m_real[i] *= ratios[i];
         m_imag[i] *= ratios[i];
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/median.h	Fri Jun 11 15:38:44 2010 +0100
@@ -0,0 +1,88 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#ifndef _MEDIAN_H_
+#define _MEDIAN_H_
+
+#include <algorithm>
+#include <cassert>
+
+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) {
+	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];
+    }
+
+    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;
+    }
+
+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);
+    }
+};
+
+#endif
+