121 the modified sinc smoother implementation (#213)

* Initial modified sinc smoother implemented

* Base FIRfilter class added

* Savitzky_golay_filter updated with FIRfilter

* modified sinc smoother implemented using the FIR class

* ran linting and formatting

* Added authorship to files

* Minor fixes to implementation of modified sinc smoother

* fixed and changed the comments to add more explanation to the approach

* Linting performed

* Small corrections made before creating PR

* added test for kappa corrections

* Fixed type return from "Self" to "_BaseFIRFilter"

* Added changes suggested during PR review

Author: NusretSalli

chemotools/smooth/__init__.py

@@ -1,6 +1,13 @@
 from ._mean_filter import MeanFilter
 from ._median_filter import MedianFilter
+from ._modified_sinc_smoother import ModifiedSincFilter
 from ._savitzky_golay_filter import SavitzkyGolayFilter
 from ._whittaker_smooth import WhittakerSmooth
 
-__all__ = ["MeanFilter", "MedianFilter", "SavitzkyGolayFilter", "WhittakerSmooth"]
+__all__ = [
+    "MeanFilter",
+    "MedianFilter",
+    "ModifiedSincFilter",
+    "SavitzkyGolayFilter",
+    "WhittakerSmooth",
+]

chemotools/smooth/_base.py

@@ -1,5 +1,5 @@
 # _base.py
-# Authors: Niklas Zell <nik.zoe@web.de>, Pau Cabaneros
+# Authors: Niklas Zell <nik.zoe@web.de>, Nusret Emirhan Salli <nusret.emirhan.salli@gmail.com>, Pau Cabaneros
 # License: MIT
 
 from __future__ import annotations
@@ -94,3 +94,109 @@ def _solve_whittaker(
             logger.debug("Primary solver failed (%s); fallback to sparse LU.", e)
             DtD = compute_DtD_sparse(len(x))
             return whittaker_smooth_sparse(x, w, self.lam, DtD)
+
+
+class _BaseFIRFilter(TransformerMixin, OneToOneFeatureMixin, BaseEstimator, ABC):
+    """
+    Base class for linear-phase FIR smoothers.
+
+    Subclasses must implement `_compute_kernel(self) -> np.ndarray`
+    returning a 1D symmetric kernel of odd length whose sum is 1.0.
+
+    Parameters
+    ----------
+    window_size : int, odd >= 3
+        Number of taps in the FIR kernel.
+    mode : {"mirror","constant","nearest","wrap","interp"}, default="interp"
+        Boundary handling. "interp" = linear extrapolation (recommended for MS).  # Schmid et al.
+    axis : int, default=1
+        Axis along which to smooth for 2D inputs (rows × features). Use 1 to
+        smooth along feature axis for each row.
+    """
+
+    def __init__(
+        self,
+        window_size: int = 21,
+        mode: Literal["mirror", "constant", "nearest", "wrap", "interp"] = "interp",
+        axis: int = 1,
+    ) -> None:
+        self.window_size = window_size
+        self.mode = mode
+        self.axis = axis
+
+    # sklearn API
+    def fit(self, X: np.ndarray, y: Optional[np.ndarray] = None) -> "_BaseFIRFilter":
+        X = validate_data(
+            self, X, y="no_validation", ensure_2d=True, reset=True, dtype=np.float64
+        )
+
+        if self.window_size < 3 or self.window_size % 2 == 0:
+            raise ValueError("window_size must be an odd integer >= 3.")
+        self.kernel_ = self._compute_kernel().astype(np.float64, copy=False)
+        if self.kernel_.ndim != 1 or self.kernel_.size != self.window_size:
+            raise ValueError("kernel must be 1D with length equal to window_size.")
+        if not np.allclose(self.kernel_.sum(), 1.0, atol=1e-12):
+            raise ValueError("kernel must be DC-preserving (sum == 1).")
+        self._half_ = (self.window_size - 1) // 2
+        return self
+
+    def transform(self, X: np.ndarray, y: Optional[np.ndarray] = None) -> np.ndarray:
+        check_is_fitted(self, "kernel_")
+        X_ = validate_data(
+            self,
+            X,
+            y="no_validation",
+            ensure_2d=True,
+            copy=True,
+            reset=False,
+            dtype=np.float64,
+        )
+
+        # move smoothing axis to last, convolve row-wise, then move back
+        ax = self.axis if self.axis >= 0 else X_.ndim + self.axis
+        X_sw = np.moveaxis(X_, ax, -1)
+        lead = int(np.prod(X_sw.shape[:-1])) or 1
+        L = X_sw.shape[-1]
+        Z = X_sw.reshape(lead, L)
+        for i in range(lead):
+            Z[i] = self._apply_filter_1d(Z[i])
+        out = X_sw.reshape(*X_sw.shape)
+        return np.moveaxis(out, -1, ax)
+
+    @abstractmethod
+    def _compute_kernel(self) -> np.ndarray:
+        """
+        Subclasses must implement this method to compute the convolution kernel.
+        """
+        raise NotImplementedError
+
+    # --- shared convolution/padding ---
+    def _apply_filter_1d(self, x: np.ndarray) -> np.ndarray:
+        m = self._half_
+        xp = self._pad_1d(x, m)
+        return np.convolve(xp, self.kernel_, mode="valid")  # same length as x
+
+    def _pad_1d(self, x: np.ndarray, m: int) -> np.ndarray:
+        if m == 0:
+            return x.copy()
+        mode = self.mode
+        if mode == "interp":
+            # Linear extrapolation using boundary slopes (paper’s recommendation).
+            if x.size < 2:
+                left = np.repeat(x[0], m)
+                right = np.repeat(x[-1], m)
+            else:
+                ls = x[1] - x[0]
+                rs = x[-1] - x[-2]
+                left = x[0] - ls * np.arange(m, 0, -1, dtype=np.float64)
+                right = x[-1] + rs * np.arange(1, m + 1, dtype=np.float64)
+            return np.concatenate([left, x, right], axis=0)
+        if mode == "nearest":
+            return np.pad(x, (m, m), mode="edge")
+        if mode == "mirror":
+            return np.pad(x, (m, m), mode="reflect")  # mirror without repeating edge
+        if mode == "wrap":
+            return np.pad(x, (m, m), mode="wrap")
+        if mode == "constant":
+            return np.pad(x, (m, m), mode="constant", constant_values=(x[0], x[-1]))
+        raise ValueError(f"Unknown mode='{mode}'")

chemotools/smooth/_modified_sinc_smoother.py

@@ -0,0 +1,203 @@
+"""
+The :mod:chemotools.smooth._modified_sinc_smoother module implements the Modified Sinc Filter (SGF) transformation.
+"""
+
+# Authors: Nusret Emirhan Salli <nusret.emirhan.salli@gmail.com>
+# License: MIT
+
+from __future__ import annotations
+from typing import Literal
+import numpy as np
+
+from ._base import _BaseFIRFilter
+
+
+class ModifiedSincFilter(_BaseFIRFilter):
+    """
+    Modified Sinc smoothing (MS) for denoising while preserving peak positions based on the paper "Why and How Savitzky–Golay Filters Should Be Replaced."
+
+    The Modified Sinc smoother is a linear-phase FIR filter: the signal is
+    convolved with a fixed, symmetric kernel. The kernel is built from:
+      1) a sinc core with argument ((n + 4) / 2) * pi * x (Eq. 3, p. 187),
+      2) a special Gaussian-like window w(x) whose value and slope vanish at
+         the window ends (Eq. 4, p. 187), and
+      3) small optional correction terms that flatten the passband so low
+         frequencies are almost unattenuated (Eqs. 7–8 and Table 1, pp. 187–188).
+
+    Parameters
+    ----------
+    window_size : int, default=21
+        Odd number of taps (2*m + 1). Larger values give stronger smoothing.
+
+    n : int, default=6
+        Even integer >= 2. Controls how many inner zeros the sinc has within
+        the window. The paper discusses n = 2, 4, 6, 8, 10 (Fig. 2).
+
+    alpha : float, default=4.0
+        Positive window width parameter. Larger alpha reduces side lobes more
+        aggressively (steeper roll-off).
+
+    use_corrections : bool, default=True
+        If True and n in {6, 8, 10} and the window is large enough, add the
+        small passband-flattening terms from Eqs. 7–8 (coefficients from Table 1).
+
+    mode : {"mirror", "constant", "nearest", "wrap", "interp"}, default="interp"
+        Boundary strategy, passed to the base FIR class. "interp" performs
+        linear extrapolation (recommended in the paper).
+
+    axis : int, default=1
+        Axis along which to smooth for 2D inputs (rows x features). Use 1 to
+        smooth within each row.
+
+    Methods
+    -------
+    fit(X, y=None)
+        Inherited from the base class. Validates input and builds the kernel.
+
+    transform(X, y=None)
+        Inherited from the base class. Pads, convolves, and returns the smoothed data.
+
+    References
+    ----------
+    [1] Schmid, M.; Rath, D.; Diebold, U. "Why and How Savitzky–Golay Filters Should Be Replaced."
+    ACS measurement science Au 2022, 2 (2), 185-196.
+
+    Examples
+    --------
+    >>> from chemotools.smooth import ModifiedSincFilter
+    >>> import numpy as np
+    >>> X = np.array([[0.0, 1.0, 2.0, 1.0, 0.0]], dtype=float)
+    >>> ms = ModifiedSincFilter(window_size= 9, n=6, alpha=4.0, mode="interp")
+    >>> X_smooth = ms.fit_transform(X)
+    """
+
+    def __init__(
+        self,
+        window_size: int = 21,
+        n: int = 6,
+        alpha: float = 4.0,
+        use_corrections: bool = True,
+        mode: Literal["mirror", "constant", "nearest", "wrap", "interp"] = "interp",
+        axis: int = 1,
+    ) -> None:
+        super().__init__(window_size=window_size, mode=mode, axis=axis)
+        self.n = n
+        self.alpha = alpha
+        self.use_corrections = use_corrections
+
+    def _compute_kernel(self) -> np.ndarray:
+        """
+        Build the Modified Sinc kernel h[i] for i = -m ... m.
+
+        Implementation map to the paper:
+          1) Map index to x = i / (m + 1)  (Eq. 5).
+          2) Base sinc: np.sinc(0.5 * (n + 4) * x)  (Eq. 3).
+          3) Solve a 3x3 system for the window coefficients so that
+             w(0)=1, w(1)=0, w'(1)=0; then we form w(x)  (Eq. 4).
+          4) h = sinc * w.
+          5) Optional corrections from based on Eqs. 7–8 (Table 1).
+          6) Symmetrize and normalize so sum(h)=1  (Eq. 6).
+        """
+
+        # ---- parameter checks ----
+        if self.n % 2 != 0 or self.n < 2:
+            raise ValueError("n must be an even integer >= 2.")
+        if self.alpha <= 0:
+            raise ValueError("alpha must be positive.")
+
+        # ---- 1) Eq. 5: index -> normalized coordinate in [-1, 1] ----
+        m = (self.window_size - 1) // 2
+        i = np.arange(-m, m + 1, dtype=np.float64)
+        x = i / (m + 1) if m >= 0 else np.array([0.0])
+
+        # ---- 2) Eq. 3: modified sinc core (note that numpy's sinc uses sin(pi*u)/(pi*u)) ----
+        core = np.sinc(0.5 * (self.n + 4) * x)
+
+        # ---- 3) Eq. 4: window with w(0)=1, w(1)=0, w'(1)=0 ----
+        # Precompute the exponentials appearing in those constraints.
+        E1 = np.exp(
+            -self.alpha * 1.0
+        )  # exp(-alpha * 1^2)   at x = 1 (central Gaussian)
+        Ep = np.exp(-self.alpha * 1.0)  # exp(-alpha * (1-2)^2) = exp(-alpha) for (x-2)
+        Em = np.exp(
+            -self.alpha * 9.0
+        )  # exp(-alpha * (1+2)^2) = exp(-9*alpha) for (x+2)
+        e4 = np.exp(-self.alpha * 4.0)  # exp(-alpha * (±2)^2) at x = 0
+
+        # Linear system rows correspond to: w(0)=1, w(1)=0, w'(1)=0
+        M = np.array(
+            [
+                [1.0, 2.0 * e4, 1.0],
+                [E1, (Ep + Em), 1.0],
+                [-2 * self.alpha * E1, 2 * self.alpha * (Ep - 3 * Em), 0.0],
+            ],
+            dtype=np.float64,
+        )
+        # solve the system using linear algebra.
+        rhs = np.array([1.0, 0.0, 0.0], dtype=np.float64)
+
+        Acoef, Bcoef, Ccoef = np.linalg.solve(M, rhs)
+
+        # Window values calculated
+        window = (
+            Acoef * np.exp(-self.alpha * x**2)
+            + Bcoef
+            * (
+                np.exp(-self.alpha * (x - 2.0) ** 2)
+                + np.exp(-self.alpha * (x + 2.0) ** 2)
+            )
+            + Ccoef
+        )
+
+        # ---- 4) base kernel = windowed sinc ----
+        h = core * window
+
+        # ---- 5) Eqs. 7–8 + Table 1: optional passband-flattening corrections ----
+        if (
+            self.use_corrections
+            and self._has_kappa_table(self.n)
+            and (m >= self.n // 2 + 2)
+        ):
+            # nu = 1 for n/2 odd (n=6,10); nu = 2 for n=8 (Eq. 7)
+            nu = 1 if ((self.n // 2) % 2 == 1) else 2
+            kappas = self._kappa_coeffs(
+                self.n, m
+            )  # kappa = a + b / (c - m)^3  (Eq. 8; Table 1)
+            corr = 0.0
+            # add the correction term according to Eq. 7
+            for j, kappa in enumerate(kappas):
+                corr += kappa * window * x * np.sin((2 * j + nu) * np.pi * x)
+            h = h + corr
+
+        # ---- 6) Eq. 6: Enforce symmetry and Direct Current = 1 ----
+        h = 0.5 * (h + h[::-1])  # make it symmetric
+        s = h.sum()
+        if not np.isfinite(s) or abs(s) < 1e-15:
+            raise FloatingPointError(
+                "Kernel normalization failed; try different parameters."
+            )
+        h = h / s  # Normalize so sum = 1
+        return h
+
+    # ====== Table 1 (p. 188): kappa fit coefficients for Eq. 8 ======
+    @staticmethod
+    def _has_kappa_table(n: int) -> bool:
+        # The paper provides corrections for n = 6, 8, 10
+        return n in (6, 8, 10)
+
+    @staticmethod
+    def _kappa_coeffs(n: int, m: int) -> np.ndarray:
+        """
+        Return [kappa_0] for n=6; [kappa_0, kappa_1] for n=8 or 10.
+        Formula: kappa_j^(n)(m) = a_j^(n) + b_j^(n) / (c_j^(n) - m)^3
+        (Eq. 8, p. 188; coefficients from Table 1).
+        """
+        if n == 6:
+            ABC = [(0.00172, 0.02437, 1.64375)]
+        elif n == 8:
+            ABC = [(0.00440, 0.08821, 2.35938), (0.00615, 0.02472, 3.63594)]
+        elif n == 10:
+            ABC = [(0.00118, 0.04219, 2.74688), (0.00367, 0.12780, 2.77031)]
+        else:
+            return np.zeros(0, dtype=np.float64)
+        return np.asarray([a + b / ((c - m) ** 3) for a, b, c in ABC], dtype=np.float64)