这是从我上面引用的参考资料中复制的“滚动中位数”的代码。我相信算法的计算成本（每个样本）是。$O(\log_2(N))$

//Copyright (c) 2011 ashelly.myopenid.com under <http://w...content-available-to-author-only...e.org/licenses/mit-license>

#include <stdlib.h>

//Customize for your data Item type
typedef int Item;
#define ItemLess(a,b)  ((a)<(b))
#define ItemMean(a,b)  (((a)+(b))/2)

typedef struct Mediator_t
{
   Item* data;  //circular queue of values
   int*  pos;   //index into `heap` for each value
   int*  heap;  //max/median/min heap holding indexes into `data`.
   int   N;     //allocated size.
   int   idx;   //position in circular queue
   int   ct;    //count of items in queue
} Mediator;

/*--- Helper Functions ---*/

#define minCt(m) (((m)->ct-1)/2) //count of items in minheap
#define maxCt(m) (((m)->ct)/2)   //count of items in maxheap 

//returns 1 if heap[i] < heap[j]
int mmless(Mediator* m, int i, int j)
{
   return ItemLess(m->data[m->heap[i]],m->data[m->heap[j]]);
}

//swaps items i&j in heap, maintains indexes
int mmexchange(Mediator* m, int i, int j)
{
   int t = m->heap[i];
   m->heap[i]=m->heap[j];
   m->heap[j]=t;
   m->pos[m->heap[i]]=i;
   m->pos[m->heap[j]]=j;
   return 1;
}

//swaps items i&j if i<j;  returns true if swapped
int mmCmpExch(Mediator* m, int i, int j)
{
   return (mmless(m,i,j) && mmexchange(m,i,j));
}

//maintains minheap property for all items below i/2.
void minSortDown(Mediator* m, int i)
{
   for (; i <= minCt(m); i*=2)
   {  if (i>1 && i < minCt(m) && mmless(m, i+1, i)) { ++i; }
      if (!mmCmpExch(m,i,i/2)) { break; }
   }
}

//maintains maxheap property for all items below i/2. (negative indexes)
void maxSortDown(Mediator* m, int i)
{
   for (; i >= -maxCt(m); i*=2)
   {  if (i<-1 && i > -maxCt(m) && mmless(m, i, i-1)) { --i; }
      if (!mmCmpExch(m,i/2,i)) { break; }
   }
}

//maintains minheap property for all items above i, including median
//returns true if median changed
int minSortUp(Mediator* m, int i)
{
   while (i>0 && mmCmpExch(m,i,i/2)) i/=2;
   return (i==0);
}

//maintains maxheap property for all items above i, including median
//returns true if median changed
int maxSortUp(Mediator* m, int i)
{
   while (i<0 && mmCmpExch(m,i/2,i))  i/=2;
   return (i==0);
}

/*--- Public Interface ---*/


//creates new Mediator: to calculate `nItems` running median. 
//mallocs single block of memory, caller must free.
Mediator* MediatorNew(int nItems)
{
   int size = sizeof(Mediator)+nItems*(sizeof(Item)+sizeof(int)*2);
   Mediator* m=  malloc(size);
   m->data= (Item*)(m+1);
   m->pos = (int*) (m->data+nItems);
   m->heap = m->pos+nItems + (nItems/2); //points to middle of storage.
   m->N=nItems;
   m->ct = m->idx = 0;
   while (nItems--)  //set up initial heap fill pattern: median,max,min,max,...
   {  m->pos[nItems]= ((nItems+1)/2) * ((nItems&1)?-1:1);
      m->heap[m->pos[nItems]]=nItems;
   }
   return m;
}


//Inserts item, maintains median in O(lg nItems)
void MediatorInsert(Mediator* m, Item v)
{
   int isNew=(m->ct<m->N);
   int p = m->pos[m->idx];
   Item old = m->data[m->idx];
   m->data[m->idx]=v;
   m->idx = (m->idx+1) % m->N;
   m->ct+=isNew;
   if (p>0)         //new item is in minHeap
   {  if (!isNew && ItemLess(old,v)) { minSortDown(m,p*2);  }
      else if (minSortUp(m,p)) { maxSortDown(m,-1); }
   }
   else if (p<0)   //new item is in maxheap
   {  if (!isNew && ItemLess(v,old)) { maxSortDown(m,p*2); }
      else if (maxSortUp(m,p)) { minSortDown(m, 1); }
   }
   else            //new item is at median
   {  if (maxCt(m)) { maxSortDown(m,-1); }
      if (minCt(m)) { minSortDown(m, 1); }
   }
}

//returns median item (or average of 2 when item count is even)
Item MediatorMedian(Mediator* m)
{
   Item v= m->data[m->heap[0]];
   if ((m->ct&1)==0) { v= ItemMean(v,m->data[m->heap[-1]]); }
   return v;
}


/*--- Test Code ---*/
#include <stdio.h>
void PrintMaxHeap(Mediator* m)
{
   int i;
   if(maxCt(m))
      printf("Max: %3d",m->data[m->heap[-1]]);
   for (i=2;i<=maxCt(m);++i)
   {
      printf("|%3d ",m->data[m->heap[-i]]);
      if(++i<=maxCt(m)) printf("%3d",m->data[m->heap[-i]]);
   }
   printf("\n");
}
void PrintMinHeap(Mediator* m)
{
   int i;
   if(minCt(m))
      printf("Min: %3d",m->data[m->heap[1]]);
   for (i=2;i<=minCt(m);++i)
   {
      printf("|%3d ",m->data[m->heap[i]]);
      if(++i<=minCt(m)) printf("%3d",m->data[m->heap[i]]);
   }
   printf("\n");
}

void ShowTree(Mediator* m)
{
   PrintMaxHeap(m);
   printf("Mid: %3d\n",m->data[m->heap[0]]);
   PrintMinHeap(m);
   printf("\n");
}

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);
   }
}

这是C++中的滚动中值算法，每步复杂度为是中值滤波器的长度（仅支持奇数）。通过每一步，您需要使用一个输入值进行过滤并返回一个新的中值，该中值也存储在变量中。该算法将中值保存在一个变量中个输入的历史记录保存在一个循环缓冲区中。在新输入时，算法执行以下操作：$O(N)$$N$update()medianmedian$N$

1. 从历史记录中删除最旧的值（但不要忘记它）并使用新输入更新历史记录。

2. 如果新输入高于median且移除的值等于或低于median：

   2.1。如果上述历史记录中的值的数量median大于，则设置为其中的最低值。$(N-1)/2$median

3. 否则，如果新输入低于median且移除的值等于或高于median：

   3.1。如果历史记录中低于的值的数量median大于，则设置为其中的最高值。$(N-1)/2$median

4. 返回median。

补充步骤 2.1 中的计数。和 3.1。用于处理历史中位数的重复项。如果两个步骤均未输入，则新输入和删除的值均高于或均低于median，因此保持不变。

该算法可以很容易地转换为使用浮点数、双精度数等而不是整数，只要记住用合适的替换INT_MAX和。INT_MIN所有整数（但不是无符号整数）都应更改为新的数据类型。

我从Robert's answer测试了这个和 Ashly 的代码的运行时间。盈亏平衡点似乎在。对于较大的，Ashelly 的代码更快，对于较小的，我的代码更快：$N=35$$N$$O(\log N)$$N$

图 1. 运行时间比较（越少越好）。编译完成，g++ -O3代码在 x64 Intel i7 PC 上运行。

我的代码一开始不会给你一个更短的中值滤波器。相反，它将假定输入为零的历史。

#include <limits.h>

const unsigned int N = 9; // odd please
const unsigned int sideN = (N-1)/2;
int history[N];
unsigned int historyIndex;
int median;

void init() {
  // Assume zeros before actual input. If you need to init with other values,
  // just feed those in using update()
  for (unsigned int i = 0; i < N; i++) {
    history[i] = 0; 
  }  
  historyIndex = 0;
  median = 0;
}

int update(int addValue) {  
  int removeValue = history[historyIndex];
  history[historyIndex] = addValue;
  if (addValue > median && removeValue <= median) {
    int lowestAboveMedian = INT_MAX;
    unsigned int numAboveMedian = 0;
    for (unsigned int i = 0; i < N; i++) {
      if (history[i] > median) {
        numAboveMedian++;
        if (history[i] < lowestAboveMedian) {
          lowestAboveMedian = history[i];
        }
      }
    }
    if (numAboveMedian > sideN) {
      median = lowestAboveMedian;
    }
  } else if (addValue < median && removeValue >= median) {
    int highestBelowMedian = INT_MIN;
    unsigned int numBelowMedian = 0;
    for (unsigned int i = 0; i < N; i++) {
      if (history[i] < median) {
        numBelowMedian++;
        if (history[i] > highestBelowMedian) {
          highestBelowMedian = history[i];
        }
      }
    }
    if (numBelowMedian > sideN) {
      median = highestBelowMedian;
    }
  }
  historyIndex++;
  if (historyIndex >= N) {
    historyIndex = 0;
  }
  return median;
}

// The rest is only for testing of the median algo

#include <stdio.h>
#include <algorithm>
#include <vector> 
#include <stdlib.h>
#include <time.h>

int main() {
  const unsigned int NData = 1000;
  printf ("N = %u, sideN = %u\n", N, sideN);
  srand (time(NULL));
  init();
  for (unsigned int i = 0; i < NData; i++) {
    int input = rand() % 10;
    printf("input = %d, median = %d, sortedHistory =", input, update(input));
    std::vector<int> sortedHistory;
    for (unsigned int j = 0; j < N; j++) {
      sortedHistory.push_back(history[j]);
    }
    std::sort(sortedHistory.begin(), sortedHistory.end());
    for (unsigned int j = 0; j < N; j++) {
      if (j == sideN) {
        printf(" (%d)", sortedHistory[j]);
      } else {
        printf(" %d", sortedHistory[j]);
      }
    }
    printf("\n");
    if (sortedHistory[sideN] != median) {
      printf("Mismatch!\n");
      return -1;
    }
  }
  printf("%u tests OK!\n", NData);
}

我想你的意思是一个有效的滑动中位数？ 我在 SE.SE 也问过同样的问题。

他们指出了一些论文，但我不知道如何执行算法。$O(\log_2(N))$

进行中位数的蛮力方法是对数据缓冲区进行排序（使用快速排序，它是，就像 FFT 一样）并挑选出恰好位于已排序数据缓冲区中间的值。如果缓冲区有偶数个样本，则平均两个中间值。$O(N \log_2(N))$

如果您的缓冲区不长（如 = 5 或 15），一个简单的方法是扫描最大次，每次重新标记您的最大值具有最大可能值的负值（导致下一个最大值成为下一次扫描的新最大值）。这样做次，最后剩下的最大值是中位数。这是简单的代码，但不是特别有效，因此缓冲区长度必须很小。这是破坏性扫描算法，因此您必须将数据复制到每个样本的临时缓冲区。$N$$O(N^2)$$\frac{N+1}{2}$$\frac{N+1}{2}$

让我们分开功能。排序算法在很多地方都被认为是已知的，并且通常可以通过高效的库获得。由于@robert bristow-johnson 提供了排序部分，这里是一个使用 Matlab 的 Matlab 版本sort.m，它适应边界（当没有足够的样本用于完整中值时）。有趣的是，您可以访问许多高效的排序算法，其中一些利用数据类型（例如整数或非整数）或近似，例如：

• 通过连续分箱快速计算中位数
• 恒定时间内的中值滤波

这是代码，希望有有意义的变量名和注释：

lWindow = 3; % Maximum size of the window; preferably odd for uniqueness; would be adapted at the borders
nSample = 25; % Number of data sample
data = randn(nSample,1)+(1:nSample)'; % Your data
dataRunningMedian = zeros(size(data));

if ~mod(lWindow,2)
error('Even median length is not supported right now');
else
lWindowHalf = (lWindow-1)/2;
end

% For each sample of the data
for iSample = 1:nSample
% Choose a frame around the sample, of maximum length 'lWindow'
% max(iSample-lWindowHalf,1) avoids indices below 1
% min(iSample+lWindowHalf,nSample)  avoids indices above nSample
lFrameIndex = (max(iSample-lWindowHalf,1):min(iSample+lWindowHalf,nSample));
% Get the data frame
dataFrame = data(lFrameIndex);
% Here, use any from-the-book sorting algorithm, or the one from @robert bristow-johnson
dataFrameSort = sort(dataFrame);
% Pick the central samples, and their average if they are two (even frame length)
lFrameLengthHalf = length(lFrameIndex)/2;
lFrameCentralIndex = [floor(lFrameLengthHalf),ceil(lFrameLengthHalf)];
dataRunningMedian(iSample) = mean(dataFrameSort(lFrameCentralIndex));
end

plot([data,dataRunningMedian]);axis tight; grid on
xlabel('Sample index')
ylabel('Amplitude (u.a.)')
legend('Data','Run. Med.')