我目前正在编写一个用于在 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
这是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
}
}
根据 @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}");
}