Local Median Approach

Mathematical formulation and iterative procedure for robust EMD

Mathematical Setup

The local median approach replaces the standard weighted mean with a weighted median, providing robustness against outliers while preserving the adaptive nature of kernel smoothing.

Core Optimization Problem

At each point \(t\), we find the local median as the solution to a weighted L1 minimization:

\[ \tilde{X}_t^{(1)} = \underset{m}{\text{argmin}} \sum_s |X_s - m| \cdot K_h(s-t) \]

This is exactly the weighted median of the observations \(\{X_s\}\) with weights \(\{K_h(s-t)\}\).

Kernel Function

We use the Epanechnikov kernel with bandwidth \(h\):

\[ K_h(u) = \frac{1}{h} K\left(\frac{u}{h}\right) \]
\[ K(u) = \frac{3}{4}(1 - u^2) \cdot \mathbf{1}_{|u| \leq 1} \]

Why Weighted Median?

Property Weighted Mean Weighted Median
Optimization \(\min_m \sum w_i(y_i - m)^2\) \(\min_m \sum w_i|y_i - m|\)
Breakdown Point 0% 50%
Single Outlier Effect Unbounded influence Bounded influence
Computation Simple average Requires sorting

Interactive Demonstration

Explore the local median approach with these interactive visualizations.

Click anywhere on the chart to see the weighted median calculation at that point

0.80
0.50
MSE (Local Median)
-
MSE (Local Median)
-
MSE (Local Mean)
-

Iterative Sifting Procedure

The local median smoothing is applied iteratively with decreasing bandwidth to extract multiple IMFs.

Sifting Step k

\[ \tilde{X}_t^{(k)} = \underset{m}{\text{argmin}} \sum_s |X_s^{(k)} - m| \cdot K_{h_k}(s-t) \]

Valid for \(t \in [h_k, 1-h_k]\) (boundary handling).

Residual Update

\[ X_t^{(k+1)} = X_t^{(k)} - \tilde{X}_t^{(k)} \]

The residual from iteration \(k\) becomes the input for iteration \(k+1\).

Algorithm Steps

  1. Initialize: Set \(X^{(1)} = X\) (original signal) and choose initial bandwidth \(h_1\)
  2. Smooth: Compute local median \(\tilde{X}^{(k)}\) using bandwidth \(h_k\)
  3. Extract: Store \(\tilde{X}^{(k)}\) as the k-th IMF component
  4. Residual: Compute \(X^{(k+1)} = X^{(k)} - \tilde{X}^{(k)}\)
  5. Halve: Update bandwidth \(h_{k+1} = h_k / \sqrt{2}\)
  6. Check: If \(h_{k+1} < h_{\min}\), stop. Otherwise, repeat from step 2

Bandwidth Selection

The bandwidth sequence controls the scale of features extracted at each iteration.

Geometric Decay Rule

\[ h_{k+1} = \frac{h_k}{\sqrt{2}}, \quad k = 1, 2, \ldots \]

Equivalently: \(h_k = h_1 \cdot 2^{-(k-1)/2}\)

Initialization

Rule of thumb: Start with \(h_1 \approx 0.1\) to smooth out noise while capturing gross features of the signal.

Stopping Criterion

Stop when: \(h_{k+1} < h_{\min}\), where \(h_{\min}\) is the minimum meaningful bandwidth (typically around 0.02-0.03 depending on sample size).

Bandwidth Sequence Visualization

h1
0.1000
h2
0.0707
h3
0.0500
h4
0.0354
h5
0.0250
Iteration Bandwidth Formula Captures
1 0.1 \(h_1\) Gross trend / low frequency
2 0.0707 \(h_1/\sqrt{2}\) Medium-scale features
3 0.05 \(h_1/2\) Finer oscillations
4 0.0354 \(h_1/(2\sqrt{2})\) High frequency details
5 0.025 \(h_1/4\) Finest resolvable scale

Output: IMF Collection

The sifting process produces a sequence of Intrinsic Mode Functions.

IMF Sequence

\[ \text{IMFs} = \left\{ \tilde{X}_t^{(1)}, \tilde{X}_t^{(2)}, \ldots, \tilde{X}_t^{(K)} \right\} \]

where \(K\) is determined by the stopping criterion.

Signal Reconstruction

\[ \bar{X}_t = \sum_{k=1}^{K} \tilde{X}_t^{(k)} \]

The sum of all IMFs approximates the original signal (minus the final residual).

Properties of the Decomposition

Pseudocode

Local Median EMD Sifting

function emd_local_median(X, h_init, h_min):
IMFs = []
residual = X
h = h_init
while h >= h_min:
# Compute local median with bandwidth h
trend = local_median_smooth(residual, h)
# Store as IMF
IMFs.append(trend)
# Update residual
residual = residual - trend
# Halve bandwidth
h = h / sqrt(2)
return IMFs, residual

Next: See Results