MedianFilter

Median of the last N samples: outlier rejection.

A linear filter can only attenuate an outlier; a median discards it. A run of up to (N - 1) / 2 consecutive corrupted samples cannot pass a median of odd length N at all – with the default N = 3, a single outlier is rejected entirely – whereas it would pass an exponential filter scaled by the filter parameter. For this reason a median filter is placed before any smoothing filter in a measurement chain.

The cost is latency: a genuine step change is delayed by (N - 1) / 2 samples. N should be odd; for even N the output is the mean of the two middle samples.

The window is stored twice: taps in arrival order, and window sorted ascending. Each tick deletes the sample leaving the window from window and inserts the new one, which is O(N), and the median is then read from the middle in O(1). This keeps the block usable at large N; recomputing the median of the whole window each tick would cost O(N log N).

The initial condition is selected through the initialization enum. Since the state is the window of past samples, every variant fills the window with one common value:

• InitialOutput(y0=...): the window is filled with y0, so for N >= 3 the first output equals y0.
• SteadyState: the window is filled with the first input sample, so the filter starts in equilibrium.
• InitialState(x0=...): the window is filled with x0.

Usage

DiscreteComponents.MedianFilter(taps_init=_median_taps_init(initialization))

Parameters:

Connectors

• u - This connector represents a real signal as an input to a component (RealInput)
• y - This connector represents a real signal as an output from a component (RealOutput)

Variables

Behavior

9-element Vector{Symbolics.Equation}:
 (taps(t))[1] ~ u(t)
 r_out(t) ~ DiscreteComponents._median_rank_of(Shift(t, -1)(window(t)), Shift(t, -1)((taps(t))[3]))
 r_in(t) ~ DiscreteComponents._median_insert_rank(Shift(t, -1)(window(t)), u(t))
 y(t) ~ DiscreteComponents._median_pick(window(t))
 (taps(t))[2] ~ Shift(t, -1)((taps(t))[1])
 (taps(t))[3] ~ Shift(t, -1)((taps(t))[2])
 (window(t))[1] ~ DiscreteComponents._median_slot(Shift(t, -1)(window(t)), u(t), r_out(t), r_in(t), 1)
 (window(t))[2] ~ DiscreteComponents._median_slot(Shift(t, -1)(window(t)), u(t), r_out(t), r_in(t), 2)
 (window(t))[3] ~ DiscreteComponents._median_slot(Shift(t, -1)(window(t)), u(t), r_out(t), r_in(t), 3)

Source

"""
Median of the last `N` samples: outlier rejection.

A linear filter can only *attenuate* an outlier; a median *discards* it. A run of up to `(N - 1) / 2` consecutive corrupted samples cannot pass a median of odd length `N` at all -- with the default `N = 3`, a single outlier is rejected entirely -- whereas it would pass an exponential filter scaled by the filter parameter. For this reason a median filter is placed *before* any smoothing filter in a measurement chain.

The cost is latency: a genuine step change is delayed by `(N - 1) / 2` samples. `N` should be odd; for even `N` the output is the mean of the two middle samples.

The window is stored twice: `taps` in arrival order, and `window` sorted ascending. Each tick deletes the sample leaving the window from `window` and inserts the new one, which is O(N), and the median is then read from the middle in O(1). This keeps the block usable at large `N`; recomputing the median of the whole window each tick would cost O(N log N).

The initial condition is selected through the `initialization` enum. Since the state is the window of past samples, every variant fills the window with one common value:
- `InitialOutput(y0=...)`: the window is filled with `y0`, so for `N >= 3` the first output equals `y0`.
- `SteadyState`: the window is filled with the first input sample, so the filter starts in equilibrium.
- `InitialState(x0=...)`: the window is filled with `x0`.
"""
component MedianFilter@[input clk extends Discrete]
  "Input signal"
  u = RealInput@[clk]() {
    "Dyad": {
      "placement": {
        "diagram": {"iconName": "default", "x1": -100, "y1": 450, "x2": 0, "y2": 550, "rot": 0}
      },
      "tags": []
    }
  }
  "Median of the last N samples"
  y = RealOutput@[clk]() {
    "Dyad": {
      "placement": {
        "diagram": {"iconName": "default", "x1": 1000, "y1": 450, "x2": 1100, "y2": 550, "rot": 0}
      },
      "tags": []
    }
  }
  "Window length: number of samples the median is taken over (odd)"
  structural parameter N::Integer = 3
  "Initial-condition specification"
  structural parameter initialization::InitialCondition = DiscreteComponents.InitialCondition.InitialOutput(y0 = 0.0)
  "Common initial value of the window (used by the InitialOutput and InitialState variants)"
  final parameter taps_init::Real = _median_taps_init(initialization)
  "The last N input samples, newest first"
  variable taps::Real[N]
  "The same N samples in ascending order"
  variable window::Real[N]
  "Position in `window` of the sample leaving on this tick"
  variable r_out::Real
  "Rank at which this tick's sample enters `window`"
  variable r_in::Real
relations
  taps[1]@clk = u@clk
  for i in 2:N
    taps[i]@clk = taps[i - 1]@(clk-1)
  end
  r_out@clk = _median_rank_of(window@(clk-1), taps[N]@(clk-1))
  r_in@clk = _median_insert_rank(window@(clk-1), u@clk)
  for i in 1:N
    window[i]@clk = _median_slot(window@(clk-1), u@clk, r_out@clk, r_in@clk, i)
  end
  y = _median_pick(window)
  switch initialization
    case InitialOutput
      initial taps@(clk-1) = fill(taps_init, N)
      initial window@(clk-1) = fill(taps_init, N)
    case SteadyState
      initial taps@(clk-1) = fill(u@clk, N)
      initial window@(clk-1) = fill(u@clk, N)
    case InitialState
      initial taps@(clk-1) = fill(taps_init, N)
      initial window@(clk-1) = fill(taps_init, N)
  end
metadata {"Dyad": {"icons": {"default": "dyad://DiscreteComponents/MedianFilter.svg"}}}
end

"""
Median of the last `N` samples: outlier rejection.

A linear filter can only *attenuate* an outlier; a median *discards* it. A run of up to `(N - 1) / 2` consecutive corrupted samples cannot pass a median of odd length `N` at all -- with the default `N = 3`, a single outlier is rejected entirely -- whereas it would pass an exponential filter scaled by the filter parameter. For this reason a median filter is placed *before* any smoothing filter in a measurement chain.

The cost is latency: a genuine step change is delayed by `(N - 1) / 2` samples. `N` should be odd; for even `N` the output is the mean of the two middle samples.

The window is stored twice: `taps` in arrival order, and `window` sorted ascending. Each tick deletes the sample leaving the window from `window` and inserts the new one, which is O(N), and the median is then read from the middle in O(1). This keeps the block usable at large `N`; recomputing the median of the whole window each tick would cost O(N log N).

The initial condition is selected through the `initialization` enum. Since the state is the window of past samples, every variant fills the window with one common value:
- `InitialOutput(y0=...)`: the window is filled with `y0`, so for `N >= 3` the first output equals `y0`.
- `SteadyState`: the window is filled with the first input sample, so the filter starts in equilibrium.
- `InitialState(x0=...)`: the window is filled with `x0`.
"""
component MedianFilter
  "Input signal"
  u = RealInput@[clk]() {
    "Dyad": {
      "placement": {
        "diagram": {"iconName": "default", "x1": -100, "y1": 450, "x2": 0, "y2": 550, "rot": 0}
      },
      "tags": []
    }
  }
  "Median of the last N samples"
  y = RealOutput@[clk]() {
    "Dyad": {
      "placement": {
        "diagram": {"iconName": "default", "x1": 1000, "y1": 450, "x2": 1100, "y2": 550, "rot": 0}
      },
      "tags": []
    }
  }
  "Window length: number of samples the median is taken over (odd)"
  structural parameter N::Integer = 3
  "Initial-condition specification"
  structural parameter initialization::InitialCondition = DiscreteComponents.InitialCondition.InitialOutput(y0 = 0.0)
  "Common initial value of the window (used by the InitialOutput and InitialState variants)"
  final parameter taps_init::Real = _median_taps_init(initialization)
  "The last N input samples, newest first"
  variable taps::Real[N]
  "The same N samples in ascending order"
  variable window::Real[N]
  "Position in `window` of the sample leaving on this tick"
  variable r_out::Real
  "Rank at which this tick's sample enters `window`"
  variable r_in::Real
relations
  taps[1]@clk = u@clk
  for i in 2:N
    taps[i]@clk = taps[i - 1]@(clk-1)
  end
  r_out@clk = _median_rank_of(window@(clk-1), taps[N]@(clk-1))
  r_in@clk = _median_insert_rank(window@(clk-1), u@clk)
  for i in 1:N
    window[i]@clk = _median_slot(window@(clk-1), u@clk, r_out@clk, r_in@clk, i)
  end
  y = _median_pick(window)
  switch initialization
    case InitialOutput
      initial taps@(clk-1) = fill(taps_init, N)
      initial window@(clk-1) = fill(taps_init, N)
    case SteadyState
      initial taps@(clk-1) = fill(u@clk, N)
      initial window@(clk-1) = fill(u@clk, N)
    case InitialState
      initial taps@(clk-1) = fill(taps_init, N)
      initial window@(clk-1) = fill(taps_init, N)
  end
metadata {"Dyad": {"icons": {"default": "dyad://DiscreteComponents/MedianFilter.svg"}}}
end

Test Cases

No test cases defined.

Related

• Examples
• Experiments
• Analyses