Median filter bias

Nielsen2011Median.py

# $Id: Nielsen2011Median.py,v 1.1 2011/05/04 12:21:30 fn Exp $

from pylab import *
from scipy.signal import medfilt 
from numpy import median

def g():
    """Generator for data points"""
    if rand() < 0.90: return 12 + randn()
    else: return 12 + 10 * randn()


def medfilt1(data, kernel_size, endcorrection='original'):
    """One-dimensional median filter"""
    halfkernel = int(kernel_size/2)
    data = asarray(data)
    if endcorrection == 'original':
        return medfilt(data, kernel_size)
    if endcorrection == 'zeropad':
        return medfilt(concatenate((asarray([0]).repeat(halfkernel),
                                    data,
                                    asarray([0]).repeat(halfkernel))), 
                       kernel_size)[halfkernel:-halfkernel]
    elif endcorrection == 'extend':
        return medfilt(concatenate((data[0].repeat(halfkernel),
                                    data,
                                    data[-1].repeat(halfkernel))), 
                       kernel_size)[halfkernel:-halfkernel]
    elif endcorrection == 'shrinkkernel':
        filtered_data = zeros(data.shape)
        for n in range(len(data)):
            i1 = max(0, n-halfkernel)
            i2 = min(len(data), n+halfkernel+1)
            filtered_data[n] = median(data[i1:i2])
        return filtered_data

seed(4)   # Overdrivelse fremmer forståelsen
data = [ g() for n in range(50) ]

plot([12]*50, color=(0.8,0.8,0.8), linewidth=12)
plot(medfilt1(data, 51), color=(0.5,1,0.5), linewidth=5)
plot(medfilt1(data, 51, endcorrection='zeropad'), color=(1,0.1,0.5), linestyle="dashed", linewidth=2)
plot(medfilt1(data, 51, endcorrection='extend'), color=(1,0.5,0.5), linewidth=3)
plot(medfilt1(data, 51, endcorrection='shrinkkernel'), color=(0,0,1), linewidth=3)
plot(data, color=(0,0,0))
legend(('"Truth"', 'Original filter (medfilt)', 'Zeropad', 
        'Extend', 'Shrink kernel', 'Original data'))
axis((0, 50, 0, 25))
title('Median filters')
xlabel('Samples')
ylabel('Values')

show()