Intrinsic multiscale filtering

Install & quickstart

pip install pimf

import numpy as np
from pimf import IMF
from pimf.contrasts import SmoothAbs
from pimf.kernels import epanechnikov

t = np.arange(1000) / 1000
rng = np.random.default_rng(0)
y = np.sin(2 * np.pi * t) + rng.normal(0, 0.1, len(t))

method = IMF(contrast=SmoothAbs(H=0.2), kernel=epanechnikov)
result = method.decompose(y, h1=0.25, a=np.sqrt(2), k_max=8)

result.imfs      # eight components, each with 1,000 samples
result.residual  # what remains after extracting them

IMF() uses a weighted mean with the squared-triangular kernel. Reuse an object with a different signal or schedule:

linear = IMF().decompose(y)
second_result = method.decompose(y, h1=0.15, k_max=5)

Reading the result

Sum selected components to keep particular scales. Adding every component and the residual reconstructs the original signal, including its noise.

unfinished = [stage.stage for stage in result.stages if not stage.converged]
coarse_signal = result.imfs[:2].sum(axis=0)

Window schedule

A bandwidth h sets how far a window reaches on each side, as a fraction of the signal length. With 1,000 samples, h1=0.25 reaches 250 samples each way: a 501-sample window.

Stop after k_max components, before the next bandwidth falls below h_min, or after a window shrinks to one sample. The result may therefore contain fewer than k_max components.

# A fixed maximum number of scales.
four_scales = method.decompose(y, h1=0.2, a=2, k_max=4)

# Bandwidths: 0.25, 0.125, 0.0625.
to_minimum = method.decompose(y, a=2, k_max=None, h_min=0.0625)

# Specify window sizes in samples.
explicit = method.decompose(y, window_sizes=[151, 75, 31])

window_sizes replaces the generated schedule. Use positive odd integers and leave the other schedule arguments at their defaults.

Kernels

The kernel controls how much each sample in a window contributes to the fit. Weights sum to one.

from pimf.kernels import epanechnikov, squared_triangle, triangle, uniform

weighted_mean = IMF(kernel=epanechnikov).decompose(y)
weights = epanechnikov.weights(31)

Contrasts

The contrast controls how the fit responds to differences in sample values, including outliers.

from pimf.contrasts import Quadratic, SmoothAbs

linear_method = IMF(contrast=Quadratic())
robust_method = IMF(contrast=SmoothAbs(H=0.2))

H must be finite and positive, in the same units as your sample values. For noise with standard deviation 0.1, H=0.2 is a starting point. Unlike h1, it changes outlier sensitivity, not window size.

Boundaries & convergence

Near either end, a window extends past the available samples. boundary determines which values fill that missing part.

For [1, 2, 3, 4], adding two samples on each side gives:

wrap:     3 4 | 1 2 3 4 | 1 2
reflect:  3 2 | 1 2 3 4 | 3 2
edge:     1 1 | 1 2 3 4 | 4 4

These choices first affect fits near the ends. Those differences can propagate into later components.

Solver accuracy

Robust fits use an iterative solver. max_iter limits iterations per stage; tol sets the stopping tolerance. Smaller tol requires smaller updates before stopping and may need more iterations.

careful = IMF(
    contrast=SmoothAbs(H=0.2),
    kernel=epanechnikov,
    boundary="reflect",
    max_iter=200,
    tol=1e-8,
)
careful_result = careful.decompose(y, k_max=6)

Check stage.converged: False means the estimate was returned before it met the tolerance. Increase max_iter to allow more work, or increase tol to accept a larger final update. This is separate from k_max, which limits the number of components.

Python API

IMF

IMF(contrast=None, kernel=None, *, boundary="wrap", max_iter=60, tol=1e-6)

IMF.decompose(y, *, h1=0.25, a=sqrt(2), k_max=8,
              h_min=None, window_sizes=None) -> IMFResult

None selects Quadratic() and squared_triangle. y must be nonempty, finite and one-dimensional. It is converted to floating point without modifying the input. Invalid numerical settings raise ValueError.

Kernel

kernel.profile(u)
kernel.weights(window_size, *, bandwidth=None)
kernel(window_size, *, bandwidth=None)

profile(u) returns finite nonnegative values with the same shape as u. Weights outside [−1, 1] are zero; the remaining weights are normalized and must have a positive sum.

Here bandwidth is a radius in samples, unlike the normalized half-widths passed to decompose. Omitting it uses the integer window radius. A one-sample window has weight one.

Contrast

contrast(r)                      # loss rho(r)
contrast.psi(r)                  # score rho'(r)
contrast.curvature()             # positive global upper bound on rho''
contrast.solve(windows, weights) # optional direct minimizer

Loss and score methods must accept NumPy arrays. A direct solver receives rows of local observations and a normalized weight vector, and returns one estimate per row. Without a direct solver, supply a finite positive curvature bound for a convex differentiable contrast.

IMFResult and StageInfo

See Reading the result for their fields. A stage converges when final_max_delta <= tol * (1 + max(abs(component))). Direct solvers report iterations=1, final_max_delta=nan and converged=True.

Custom kernels & contrasts

Subclass Kernel and define a profile. The base class handles sampling, compact support and normalization.

from pimf import Kernel

class Cosine(Kernel):
    def profile(self, u):
        return np.cos(np.pi * u / 2)

cosine_result = IMF(kernel=Cosine()).decompose(y)

This Huber contrast uses a quadratic loss for small differences and a linear loss beyond delta:

from pimf import Contrast

class Huber(Contrast):
    def __init__(self, delta):
        if not np.isfinite(delta) or delta <= 0:
            raise ValueError("delta must be finite and positive")
        self.delta = delta

    def __call__(self, r):
        r = np.asarray(r, dtype=float)
        magnitude = np.abs(r)
        return np.where(
            magnitude <= self.delta,
            0.5 * r**2,
            self.delta * (magnitude - 0.5 * self.delta),
        )

    def psi(self, r):
        return np.clip(r, -self.delta, self.delta)

    def curvature(self):
        return 1.0

huber_result = IMF(contrast=Huber(0.3)).decompose(y)