C中的滚动中值算法

Posted 2023-02-16

【中文标题】C中的滚动中值算法

【英文标题】：Rolling median algorithm in C

【发布时间】：2010-11-21 12:14:09

【问题描述】：

我目前正在研究一种在 C 中实现滚动中值滤波器（类似于滚动均值滤波器）的算法。根据我对文献的搜索，似乎有两种相当有效的方法可以做到这一点。首先是对值的初始窗口进行排序，然后执行二分搜索以插入新值并在每次迭代时删除现有值。

第二个（来自 Hardle 和 Steiger，1995，JRSS-C，算法 296）构建了一个双端堆结构，一端是最大堆，另一端是最小堆，中间是中间值。这产生了一种线性时间算法，而不是 O(n log n)。

这是我的问题：实现前者是可行的，但我需要在数百万个时间序列上运行它，所以效率很重要。事实证明，后者非常难以实施。我在 R 的 stats 包的代码的 Trunmed.c 文件中找到了代码，但它相当难以理解。

有人知道线性时间滚动中值算法的编写良好的 C 实现吗？

编辑：链接到 Trunmed.c 代码http://google.com/codesearch/p?hl=en&sa=N&cd=1&ct=rc#mYw3h_Lb_e0/R-2.2.0/src/library/stats/src/Trunmed.c

【问题讨论】：

刚刚实现了移动均值...移动中值有点棘手。尝试使用谷歌搜索移动中位数。

试过谷歌和谷歌代码搜索。它出现了 Trunmed.c 代码和另一种语言的实现，用于 Trunmed 代码的 SGI 端口（据我所知）。此外，我引用的 JRSS 算法显然是该期刊系列中唯一没有原始代码存档的算法。

每个时间序列有多少个数字？即使有一百万个，如果你只有几千个数字，运行时间可能不会超过一两分钟（如果你的代码编写得很好的话）。

两堆解如何线性？它是 O(n log k)，其中 k 是窗口大小，因为堆的删除是 O(log k)。

一些实现和比较：github.com/suomela/median-filter

【参考方案1】：

我已经看过 R 的 src/library/stats/src/Trunmed.c 几次，因为我也希望在独立的 C++ 类/C 子例程中也有类似的东西。请注意，这实际上是两个实现合而为一，请参阅src/library/stats/man/runmed.Rd（帮助文件的来源），其中说

\details
  Apart from the end values, the result \codey = runmed(x, k) simply has
  \codey[j] = median(x[(j-k2):(j+k2)]) (k = 2*k2+1), computed very
  efficiently.

  The two algorithms are internally entirely different:
  \describe
    \item"Turlach"is the Härdle-Steiger
      algorithm (see Ref.) as implemented by Berwin Turlach.
      A tree algorithm is used, ensuring performance \eqnO(n \log
        k)O(n * log(k)) where \coden <- length(x) which is
      asymptotically optimal.
    \item"Stuetzle"is the (older) Stuetzle-Friedman implementation
      which makes use of median \emphupdating when one observation
      enters and one leaves the smoothing window.  While this performs as
      \eqnO(n \times k)O(n * k) which is slower asymptotically, it is
      considerably faster for small \eqnk or \eqnn.
  

很高兴看到它以更独立的方式重新使用。你在做志愿者吗？我可以帮助一些 R 位。

编辑 1：除了上面旧版 Trunmed.c 的链接，这里是当前的 SVN 副本

Srunmed.c（适用于 Stuetzle 版本） Trunmed.c（适用于 Turlach 版本） runmed.R 用于调用这些的 R 函数

编辑 2：Ryan Tibshirani 在 fast median binning 上有一些 C 和 Fortran 代码，这可能是窗口化方法的合适起点。

【讨论】：

谢谢德克。一旦我得到一个干净的解决方案，我计划在 GPL 下发布它。我也有兴趣设置 R 和 Python 接口。

@AWB 这个想法最终发生了什么？您是否将您的解决方案整合到一个包中？

【参考方案2】：

我找不到具有 order-statistic 的 c++ 数据结构的现代实现，因此最终在 MAK 建议的 top coders 链接中实现了这两个想法（Match Editorial：向下滚动到 FloatingMedian）。

两个多组

第一个想法将数据划分为两个数据结构（堆、多集等），每次插入/删除 O(ln N) 不允许在没有大成本的情况下动态更改分位数。 IE。我们可以有滚动的中位数，或者滚动的 75%，但不能同时使用两者。

段树

第二个想法使用 O(ln N) 的段树进行插入/删除/查询，但更灵活。最好的“N”是数据范围的大小。因此，如果您的滚动中位数具有一百万个项目的窗口，但您的数据与 1..65536 不同，那么每次移动 100 万个滚动窗口只需要 16 次操作！！

c++ 代码与上面 Denis 发布的类似（“这是一个简单的量化数据算法”）

GNU 顺序统计树

就在放弃之前，发现stdlibc++中包含顺序统计树！！！

这些有两个关键操作：

iter = tree.find_by_order(value)
order = tree.order_of_key(value)

参见libstdc++ manual policy_based_data_structures_test（搜索“拆分并加入”）。

我已经包装了树，以便在支持 c++0x/c++11 样式部分类型定义的编译器的便利标头中使用：

#if !defined(GNU_ORDER_STATISTIC_SET_H)
#define GNU_ORDER_STATISTIC_SET_H
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

// A red-black tree table storing ints and their order
// statistics. Note that since the tree uses
// tree_order_statistics_node_update as its update policy, then it
// includes its methods by_order and order_of_key.
template <typename T>
using t_order_statistic_set = __gnu_pbds::tree<
                                  T,
                                  __gnu_pbds::null_type,
                                  std::less<T>,
                                  __gnu_pbds::rb_tree_tag,
                                  // This policy updates nodes'  metadata for order statistics.
                                  __gnu_pbds::tree_order_statistics_node_update>;

#endif //GNU_ORDER_STATISTIC_SET_H

【讨论】：

实际上，libstdc++ 扩展容器不允许允许多个值！按设计！正如我上面的名字 (t_order_statistic_set) 所建议的那样，多个值被合并。所以，为了我们的目的，他们需要做更多的工作:-(

我们需要 1) 制作要计数的值映射（而不是集合） 2) 分支大小应反映键的计数 (libstdc++-v3/include/ext/pb_ds/detail/ tree_policy/order_statistics_imp.hpp) 从树继承，以及 3) 重载 insert() 以增加计数/如果该值已经存在则调用 update_to_top() 4) 重载 erase() 以减少计数/调用 update_to_top() 如果值不是唯一的（参见 libstdc++-v3/include/ext/pb_ds/detail/rb_tree_map_/rb_tree_.hpp）有志愿者吗？？

【参考方案3】：

我已经完成了C implementation here。这个问题有更多细节：Rolling median in C - Turlach implementation。

示例用法：

int main(int argc, char* argv[])

   int i, v;
   Mediator* m = MediatorNew(15);
 
   for (i=0; i<30; i++) 
      v = rand() & 127;
      printf("Inserting %3d \n", v);
      MediatorInsert(m, v);
      v = MediatorMedian(m);
      printf("Median = %3d.\n\n", v);
      ShowTree(m);
   

【讨论】：

基于 min-median-max 堆的出色、快速和清晰的实现。非常好的工作。

如何找到这个解决方案的 Java 版本？

【参考方案4】：

我使用这个增量中值估计器：

median += eta * sgn(sample - median)

与更常见的均值估计器具有相同的形式：

mean += eta * (sample - mean)

这里 eta 是一个小的学习率参数（例如0.001），sgn() 是返回-1, 0, 1 之一的符号函数。 （如果数据是非固定的并且您想跟踪随时间的变化，请使用像这样的常量eta；否则，对于固定源使用类似eta = 1 / n 的东西来收敛，其中n 是看到的样本数远。）

另外，我修改了中位数估计器，使其适用于任意分位数。通常，quantile function 告诉您将数据分成两个部分的值：p 和 1 - p。下面逐步估计这个值：

quantile += eta * (sgn(sample - quantile) + 2.0 * p - 1.0)

p 的值应在[0, 1] 范围内。这实质上将sgn() 函数的对称输出-1, 0, 1 向一侧倾斜，将数据样本分成两个大小不等的箱（数据的分数p 和1 - p 小于/大于分位数估计， 分别）。请注意，对于 p = 0.5，这会减少到中值估计量。

【讨论】：

酷，这是一个根据运行均值调整“eta”的修改...（均值用作中位数的粗略估计，因此它以与收敛相同的速率收敛于大值微小的值）。即 eta 自动调整。 ***.com/questions/11482529/…

对于类似的技术，请参阅这篇关于节俭流的论文：arxiv.org/pdf/1407.1121v1.pdf 它可以估计任何四分位数并适应均值的变化。它要求您只存储两个值：上次估计和上次调整的方向（+1 或 -1）。该算法易于实现。我发现大约 97% 的时间误差在 5% 以内。

【参考方案5】：

这是一个简单的量化数据算法（几个月后）：

""" median1.py: moving median 1d for quantized, e.g. 8-bit data

Method: cache the median, so that wider windows are faster.
    The code is simple -- no heaps, no trees.

Keywords: median filter, moving median, running median, numpy, scipy

See Perreault + Hebert, Median Filtering in Constant Time, 2007,
    http://nomis80.org/ctmf.html: nice 6-page paper and C code,
    mainly for 2d images

Example:
    y = medians( x, window=window, nlevel=nlevel )
    uses:
    med = Median1( nlevel, window, counts=np.bincount( x[0:window] ))
    med.addsub( +, - )  -- see the picture in Perreault
    m = med.median()  -- using cached m, summ

How it works:
    picture nlevel=8, window=3 -- 3 1s in an array of 8 counters:
        counts: . 1 . . 1 . 1 .
        sums:   0 1 1 1 2 2 3 3
                        ^ sums[3] < 2 <= sums[4] <=> median 4
        addsub( 0, 1 )  m, summ stay the same
        addsub( 5, 1 )  slide right
        addsub( 5, 6 )  slide left

Updating `counts` in an `addsub` is trivial, updating `sums` is not.
But we can cache the previous median `m` and the sum to m `summ`.
The less often the median changes, the faster;
so fewer levels or *wider* windows are faster.
(Like any cache, run time varies a lot, depending on the input.)

See also:
    scipy.signal.medfilt -- runtime roughly ~ window size
    http://***.com/questions/1309263/rolling-median-algorithm-in-c

"""

from __future__ import division
import numpy as np  # bincount, pad0

__date__ = "2009-10-27 oct"
__author_email__ = "denis-bz-py at t-online dot de"


#...............................................................................
class Median1:
    """ moving median 1d for quantized, e.g. 8-bit data """

    def __init__( s, nlevel, window, counts ):
        s.nlevel = nlevel  # >= len(counts)
        s.window = window  # == sum(counts)
        s.half = (window // 2) + 1  # odd or even
        s.setcounts( counts )

    def median( s ):
        """ step up or down until sum cnt to m-1 < half <= sum to m """
        if s.summ - s.cnt[s.m] < s.half <= s.summ:
            return s.m
        j, sumj = s.m, s.summ
        if sumj <= s.half:
            while j < s.nlevel - 1:
                j += 1
                sumj += s.cnt[j]
                # print "j sumj:", j, sumj
                if sumj - s.cnt[j] < s.half <= sumj:  break
        else:
            while j > 0:
                sumj -= s.cnt[j]
                j -= 1
                # print "j sumj:", j, sumj
                if sumj - s.cnt[j] < s.half <= sumj:  break
        s.m, s.summ = j, sumj
        return s.m

    def addsub( s, add, sub ):
        s.cnt[add] += 1
        s.cnt[sub] -= 1
        assert s.cnt[sub] >= 0, (add, sub)
        if add <= s.m:
            s.summ += 1
        if sub <= s.m:
            s.summ -= 1

    def setcounts( s, counts ):
        assert len(counts) <= s.nlevel, (len(counts), s.nlevel)
        if len(counts) < s.nlevel:
            counts = pad0__( counts, s.nlevel )  # numpy array / list
        sumcounts = sum(counts)
        assert sumcounts == s.window, (sumcounts, s.window)
        s.cnt = counts
        s.slowmedian()

    def slowmedian( s ):
        j, sumj = -1, 0
        while sumj < s.half:
            j += 1
            sumj += s.cnt[j]
        s.m, s.summ = j, sumj

    def __str__( s ):
        return ("median %d: " % s.m) + \
            "".join([ (" ." if c == 0 else "%2d" % c) for c in s.cnt ])

#...............................................................................
def medianfilter( x, window, nlevel=256 ):
    """ moving medians, y[j] = median( x[j:j+window] )
        -> a shorter list, len(y) = len(x) - window + 1
    """
    assert len(x) >= window, (len(x), window)
    # np.clip( x, 0, nlevel-1, out=x )
        # cf http://scipy.org/Cookbook/Rebinning
    cnt = np.bincount( x[0:window] )
    med = Median1( nlevel=nlevel, window=window, counts=cnt )
    y = (len(x) - window + 1) * [0]
    y[0] = med.median()
    for j in xrange( len(x) - window ):
        med.addsub( x[j+window], x[j] )
        y[j+1] = med.median()
    return y  # list
    # return np.array( y )

def pad0__( x, tolen ):
    """ pad x with 0 s, numpy array or list """
    n = tolen - len(x)
    if n > 0:
        try:
            x = np.r_[ x, np.zeros( n, dtype=x[0].dtype )]
        except NameError:
            x += n * [0]
    return x

#...............................................................................
if __name__ == "__main__":
    Len = 10000
    window = 3
    nlevel = 256
    period = 100

    np.set_printoptions( 2, threshold=100, edgeitems=10 )
    # print medians( np.arange(3), 3 )

    sinwave = (np.sin( 2 * np.pi * np.arange(Len) / period )
        + 1) * (nlevel-1) / 2
    x = np.asarray( sinwave, int )
    print "x:", x
    for window in ( 3, 31, 63, 127, 255 ):
        if window > Len:  continue
        print "medianfilter: Len=%d window=%d nlevel=%d:" % (Len, window, nlevel)
            y = medianfilter( x, window=window, nlevel=nlevel )
        print np.array( y )

# end median1.py

【参考方案6】：

可以通过维护两个数字分区来找到滚动中位数。

要维护分区，请使用 Min Heap 和 Max Heap。

最大堆将包含小于等于中位数的数字。

Min Heap 将包含大于等于中位数的数字。

平衡约束： 如果元素总数是偶数，那么两个堆应该有相等的元素。

如果元素总数为奇数，则 Max Heap 将比 Min Heap 多一个元素。

中值元素：如果两个分区的元素数量相等，则中值将是第一个分区的最大元素和第二个分区的最小元素之和的一半。

否则中位数将是第一个分区的最大元素。

算法- 1-取两个堆（1个最小堆和1个最大堆） Max Heap 将包含前半数的元素 最小堆将包含后半数的元素 2-将流中的新数字与最大堆顶部进行比较， 如果它更小或相等，则将该数字添加到最大堆中。 否则在最小堆中添加数字。 3- 如果 min Heap 比 Max Heap 有更多的元素 然后删除最小堆的顶部元素并添加最大堆。 如果 max Heap 比 Min Heap 有多个元素 然后删除 Max Heap 的顶部元素并添加 Min Heap。 4-如果两个堆都有相同数量的元素，那么 中位数将是最大堆中的最大元素和最小堆中的最小元素之和的一半。 否则，中位数将是第一个分区的最大元素。

public class Solution 

    public static void main(String[] args) 
        Scanner in = new Scanner(System.in);
        RunningMedianHeaps s = new RunningMedianHeaps();
        int n = in.nextInt();
        for(int a_i=0; a_i < n; a_i++)
            printMedian(s,in.nextInt());
        
        in.close();       
    

    public static void printMedian(RunningMedianHeaps s, int nextNum)
            s.addNumberInHeap(nextNum);
            System.out.printf("%.1f\n",s.getMedian());
    


class RunningMedianHeaps
    PriorityQueue<Integer> minHeap = new PriorityQueue<Integer>();
    PriorityQueue<Integer> maxHeap = new PriorityQueue<Integer>(Comparator.reverseOrder());

    public double getMedian() 

        int size = minHeap.size() + maxHeap.size();     
        if(size % 2 == 0)
            return (maxHeap.peek()+minHeap.peek())/2.0;
        return maxHeap.peek()*1.0;
    

    private void balanceHeaps() 
        if(maxHeap.size() < minHeap.size())
        
            maxHeap.add(minHeap.poll());
           
        else if(maxHeap.size() > 1+minHeap.size())
        
            minHeap.add(maxHeap.poll());
        
    

    public void addNumberInHeap(int num) 
        if(maxHeap.size()==0 || num <= maxHeap.peek())
        
            maxHeap.add(num);
        
        else
        
            minHeap.add(num);
        
        balanceHeaps();
    

【讨论】：

我不清楚第三个 Java 答案为 C 问题提供了多少好处。您应该提出一个新问题，然后在该问题中提供您的 Java 答案。

逻辑在阅读此“然后删除最小堆的顶部元素并添加最小堆”后死亡。 .至少在发帖前有礼貌地阅读算法

这个算法不是针对滚动中位数，而是针对越来越多的元素的中位数。对于滚动中位数，还必须从堆中删除一个元素，这需要首先找到。

【参考方案7】：

可能值得指出的是，有一种特殊情况有一个简单的精确解：当流中的所有值都是（相对）小的定义范围内的整数时。例如，假设它们都必须位于 0 到 1023 之间。在这种情况下，只需定义一个包含 1024 个元素和一个计数的数组，然后清除所有这些值。对于流中的每个值，递增相应的 bin 和计数。在流结束后，找到包含计数/2 最高值的 bin - 通过添加从 0 开始的连续 bin 轻松完成。使用相同的方法可以找到任意排名顺序的值。 （如果需要在运行期间检测 bin 饱和并将存储 bin 的大小“升级”为更大的类型，则会有一个小麻烦。）

这种特殊情况看似人为，但在实践中却很常见。如果实数在一个范围内并且已知“足够好”的精度水平，它也可以用作实数的近似值。这几乎适用于一组“真实世界”对象的任何一组测量。例如，一群人的身高或体重。集数不够大？它同样适用于地球上所有（单个）细菌的长度或重量——假设有人可以提供数据！

看起来我误读了原文——它似乎想要一个滑动窗口中位数，而不是一个很长的流的中位数。这种方法仍然适用。为初始窗口加载前 N 个流值，然后为第 N+1 个流值增加相应的 bin，同时减少与第 0 个流值对应的 bin。在这种情况下，有必要保留最后 N 个值以允许递减，这可以通过循环寻址大小为 N 的数组来有效地完成。由于中位数的位置只能改变 -2,-1,0,1 ,2 在滑动窗口的每一步上，没有必要在每一步将所有 bin 相加到中值，只需根据修改了哪些边 bin 来调整“中值指针”。例如，如果新值和被删除的值都低于当前中值，则它不会改变（偏移量 = 0）。当 N 变得太大而不能方便地保存在内存中时，该方法就会失效。

【参考方案8】：

如果您有能力将值作为时间点的函数进行引用，则可以对值进行替换采样，应用 bootstrapping 在置信区间内生成自举中值。与不断将传入值排序到数据结构中相比，这可以让您更高效地计算近似中位数。

【参考方案9】：

对于那些需要在 Java 中运行中位数的人...PriorityQueue 是您的朋友。 O(log N) 插入，O(1) 当前中位数和 O(N) 删除。如果您知道数据的分布，您可以做得比这更好。

public class RunningMedian 
  // Two priority queues, one of reversed order.
  PriorityQueue<Integer> lower = new PriorityQueue<Integer>(10,
          new Comparator<Integer>() 
              public int compare(Integer arg0, Integer arg1) 
                  return (arg0 < arg1) ? 1 : arg0 == arg1 ? 0 : -1;
              
          ), higher = new PriorityQueue<Integer>();

  public void insert(Integer n) 
      if (lower.isEmpty() && higher.isEmpty())
          lower.add(n);
      else 
          if (n <= lower.peek())
              lower.add(n);
          else
              higher.add(n);
          rebalance();
      
  

  void rebalance() 
      if (lower.size() < higher.size() - 1)
          lower.add(higher.remove());
      else if (higher.size() < lower.size() - 1)
          higher.add(lower.remove());
  

  public Integer getMedian() 
      if (lower.isEmpty() && higher.isEmpty())
          return null;
      else if (lower.size() == higher.size())
          return (lower.peek() + higher.peek()) / 2;
      else
          return (lower.size() < higher.size()) ? higher.peek() : lower
                  .peek();
  

  public void remove(Integer n) 
      if (lower.remove(n) || higher.remove(n))
          rebalance();
  

【讨论】：

c++ 在标准库的扩展中具有来自 gnu 的顺序统计树。请参阅下面的帖子。

我认为您的代码没有正确放在这里。那里有一些不完整的部分，例如：), higher = new PriorityQueue<Integer>(); 或 new PriorityQueue<Integer>(10,。我无法运行代码。

@Hengameh Java 以分号结束语句——换行根本不重要。你一定是复制错了。

您应该提出一个新问题，然后在该问题中提供您的 Java 答案。

【参考方案10】：

当精确输出不重要时（用于显示目的等），可以使用此方法 您需要 totalcount 和 lastmedian，再加上 newvalue。

totalcount++;
newmedian=lastmedian+(newvalue>lastmedian?1:-1)*(lastmedian==0?newvalue: lastmedian/totalcount*2);

为诸如 page_display_time 之类的内容生成非常精确的结果。

规则：输入流需要按页面显示时间顺序平滑，数量大（>30 等），并且中位数不为零。

示例： 页面加载时间，800 条，10ms...3000ms，平均 90ms，实际中位数：11ms

30 次输入后，中位数误差一般

另一个有类似解决方案的思想家在这里：Median Filter Super efficient implementation

【参考方案11】：

这里是java实现

package MedianOfIntegerStream;

import java.util.Comparator;
import java.util.HashSet;
import java.util.Iterator;
import java.util.Set;
import java.util.TreeSet;


public class MedianOfIntegerStream 

    public Set<Integer> rightMinSet;
    public Set<Integer> leftMaxSet;
    public int numOfElements;

    public MedianOfIntegerStream() 
        rightMinSet = new TreeSet<Integer>();
        leftMaxSet = new TreeSet<Integer>(new DescendingComparator());
        numOfElements = 0;
    

    public void addNumberToStream(Integer num) 
        leftMaxSet.add(num);

        Iterator<Integer> iterMax = leftMaxSet.iterator();
        Iterator<Integer> iterMin = rightMinSet.iterator();
        int maxEl = iterMax.next();
        int minEl = 0;
        if (iterMin.hasNext()) 
            minEl = iterMin.next();
        

        if (numOfElements % 2 == 0) 
            if (numOfElements == 0) 
                numOfElements++;
                return;
             else if (maxEl > minEl) 
                iterMax.remove();

                if (minEl != 0) 
                    iterMin.remove();
                
                leftMaxSet.add(minEl);
                rightMinSet.add(maxEl);
            
         else 

            if (maxEl != 0) 
                iterMax.remove();
            

            rightMinSet.add(maxEl);
        
        numOfElements++;
    

    public Double getMedian() 
        if (numOfElements % 2 != 0)
            return new Double(leftMaxSet.iterator().next());
        else
            return (leftMaxSet.iterator().next() + rightMinSet.iterator().next()) / 2.0;
    

    private class DescendingComparator implements Comparator<Integer> 
        @Override
        public int compare(Integer o1, Integer o2) 
            return o2 - o1;
        
    

    public static void main(String[] args) 
        MedianOfIntegerStream streamMedian = new MedianOfIntegerStream();

        streamMedian.addNumberToStream(1);
        System.out.println(streamMedian.getMedian()); // should be 1

        streamMedian.addNumberToStream(5);
        streamMedian.addNumberToStream(10);
        streamMedian.addNumberToStream(12);
        streamMedian.addNumberToStream(2);
        System.out.println(streamMedian.getMedian()); // should be 5

        streamMedian.addNumberToStream(3);
        streamMedian.addNumberToStream(8);
        streamMedian.addNumberToStream(9);
        System.out.println(streamMedian.getMedian()); // should be 6.5
    

【讨论】：

您应该提出一个新问题，然后在该问题中提供您的 Java 答案。

【参考方案12】：

基于@mathog 的想法，这是一个 C# 实现，用于在具有已知值范围的字节数组上运行中位数。可以扩展到其他整数类型。

  /// <summary>
  /// Median estimation by histogram, avoids multiple sorting operations for a running median
  /// </summary>
  public class MedianEstimator
  
    private readonly int m_size2;
    private readonly byte[] m_counts;

    /// <summary>
    /// Estimated median, available right after calling <see cref="Init"/> or <see cref="Update"/>.
    /// </summary>
    public byte Median  get; private set; 

    /// <summary>
    /// Ctor
    /// </summary>
    /// <param name="size">Median size in samples</param>
    /// <param name="maxValue">Maximum expected value in input data</param>
    public MedianEstimator(
      int size,
      byte maxValue)
    
      m_size2 = size / 2;
      m_counts = new byte[maxValue + 1];
    

    /// <summary>
    /// Initializes the internal histogram with the passed sample values
    /// </summary>
    /// <param name="values">Array of values, usually the start of the array for a running median</param>
    public void Init(byte[] values)
    
      for (var i = 0; i < values.Length; i++)
        m_counts[values[i]]++;

      UpdateMedian();
    

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private void UpdateMedian()
    
      // The median is the first value up to which counts add to size / 2
      var sum = 0;
      Median = 0;
      for (var i = 0; i < m_counts.Length; i++)
      
        sum += m_counts[i];
        Median = (byte) i;
        if (sum > m_size2) break;
      
    

    /// <summary>
    /// Updates the median estimation by removing <paramref name="last"/> and adding <paramref name="next"/>. These
    /// values must be updated as the running median is applied. If the median length is <i>N</i>, at the sample
    /// <i>i</i>, <paramref name="last"/> is sample at index <i>i</i>-<i>N</i>/2 and <paramref name="next"/> is sample
    /// at index <i>i</i>+<i>N</i>/2+1.
    /// </summary>
    /// <param name="last">Sample at the start of the moving window that is to be removed</param>
    /// <param name="next">Sample at the end of the moving window + 1 that is to be added</param>
    public void Update(byte last, byte next)
    
      m_counts[last]--;
      m_counts[next]++;

      // The conditions below do not change median value so there is no need to update it
      if (last == next ||
          last < Median && next < Median || // both below median
          last > Median && next > Median) // both above median
        return;

      UpdateMedian();
    

针对正在运行的中位数进行测试，有时间：

private void TestMedianEstimator()

  var r = new Random();
  const int SIZE = 15;
  const byte MAX_VAL = 80;
  var values = new byte[100000];
  for (var i = 0; i < values.Length; i++)
    values[i] = (byte) (MAX_VAL * r.NextDouble());

  var timer = Stopwatch.StartNew();
  // Running median
  var window = new byte[2 * SIZE + 1];
  var medians = new byte[values.Length];
  for (var i = SIZE; i < values.Length - SIZE - 1; i++)
  
    for (int j = i - SIZE, k = 0; j <= i + SIZE; j++, k++)
      window[k] = values[j];
    Array.Sort(window);
    medians[i] = window[SIZE];
  
  timer.Stop();
  var elapsed1 = timer.Elapsed;

  timer.Restart();
  var me = new MedianEstimator(2 * SIZE + 1, MAX_VAL);
  me.Init(values.Slice(0, 2 * SIZE + 1));
  var meMedians = new byte[values.Length];
  for (var i = SIZE; i < values.Length - SIZE - 1; i++)
  
    meMedians[i] = me.Median;
    me.Update(values[i - SIZE], values[i + SIZE + 1]);
  
  timer.Stop();
  var elapsed2 = timer.Elapsed;

  WriteLineToLog($"elapsed1.TotalMilliseconds / elapsed2.TotalMilliseconds:0.00");

  var diff = 0;
  for (var i = 0; i < meMedians.Length; i++)
    diff += Math.Abs(meMedians[i] - medians[i]);
  WriteLineToLog($"Diff: diff");

【参考方案13】：

如果您只需要平滑平均值，一种快速/简单的方法是将最新值乘以 x，将平均值乘以 (1-x)，然后将它们相加。这将成为新的平均值。

编辑：不是用户要求的，也不是统计上有效的，但对于很多用途来说已经足够了。 我会把它留在这里（尽管投反对票）以供搜索！

【讨论】：

计算平均值。他想要中位数。此外，他正在计算一个滑动窗口的中值，而不是整个集合的中值。

这会计算一个窗口值的运行平均值，其衰减常数取决于 X - 在性能很重要且您不必费心做卡尔曼滤波器的情况下，它非常有用。我把它放进去以便搜索可以找到它。

这也是我立即想到的，已经实现了这样一个滤波器，作为音频应用程序的一个非常基本且便宜的低通滤波器。