"""Core time-series and covariate containers for the nSTAT toolbox.
This module hosts the two foundational classes used everywhere in nSTAT:
- :class:`SignalObj` — multi-dimensional regularly-sampled signal container
(mirrors MATLAB ``SignalObj.m``). All continuous covariates,
stimuli, model rates, etc. inherit from or wrap ``SignalObj``.
- :class:`Covariate` — :class:`SignalObj` subclass that also carries
per-dimension confidence intervals and accepts the convenience aliases
``values`` / ``units`` (mirrors MATLAB ``Covariate.m``).
The ``nspikeTrain`` point-process class (Matlab ``nspikeTrain.m``) was
extracted to :mod:`nstat._spike_train_impl` in v0.3.1 for readability,
but is re-exported here so that ``from nstat.core import nspikeTrain``
continues to work. All time vectors are in **seconds** and sample rates
in **Hz**, matching the MATLAB toolbox.
"""
from __future__ import annotations
from collections.abc import Sequence
from typing import Any
import numpy as np
def _matlab_colon(start: float, step: float, stop: float) -> np.ndarray:
"""Replicate MATLAB ``start:step:stop`` exactly, bit for bit.
Delegates to :func:`_matlab_colon_exact`. An earlier version of this
helper computed ``m = fix((stop-start)/step)`` and returned ``m+1``
samples via ``start + np.arange(m+1) * step``, accumulating float error
element-by-element; against the ``colon_*`` MATLAB gold (487 arrays)
this differed from MATLAB's actual build-from-both-ends algorithm in
385/487 cases, 15 of those in LENGTH too (not just element values).
"""
return _matlab_colon_exact(start, step, stop)
def _matlab_colon_exact(start: float, step: float, stop: float) -> np.ndarray:
"""MATLAB ``start:step:stop`` reproduced bit for bit (double inputs).
``_matlab_colon`` gets the length right but not every element: MATLAB
builds the vector from both ends. With ``n`` intervals (``n`` rounded
from ``(stop-start)/step``, less one if ``start+n*step`` overshoots ``stop``
by more than ``tol = 2*eps*max(|start|,|stop|)``) and end point ``c``
(``stop`` if within ``tol`` of ``start+n*step``, else ``start+n*step``),
element ``k`` is ``start + k*step`` for ``k <= n/2`` and ``c - (n-k)*step``
above; for even ``n`` the middle element is ``start + (c-start)/2`` when
``start`` and ``c`` have the same sign, else ``(start+c)/2``. Verified
bitwise against 4,337 MATLAB R2025b outputs (default history edges
``0:delta:(m+1)*delta`` for 14 deltas and m = 0..80, plus random starts,
signed steps and non-multiple stops). Used where an element-exact edge
vector matters (history windows are binned with ``ceil(t*sampleRate)``).
"""
a, d, b = float(start), float(step), float(stop)
if d == 0.0:
raise ValueError("_matlab_colon_exact: step must be non-zero")
if not (np.isfinite(a) and np.isfinite(d) and np.isfinite(b)):
raise ValueError("_matlab_colon_exact: start, step and stop must be finite")
if (a < b and d < 0.0) or (b < a and d > 0.0):
return np.zeros(0, dtype=float)
tol = 2.0 * np.finfo(float).eps * max(abs(a), abs(b))
sig = 1.0 if d > 0.0 else -1.0
q = (b - a) / d
n = int(np.sign(q) * np.floor(abs(q) + 0.5)) # MATLAB round (half away from zero)
if sig * (a + n * d - b) > tol:
n -= 1
c = a + n * d
if abs(b - c) <= tol:
c = b
out = np.empty(n + 1, dtype=float)
k = np.arange(n // 2 + 1)
out[k] = a + k * d
out[n - k] = c - k * d
if n % 2 == 0:
out[n // 2] = a + (c - a) / 2.0 if np.sign(a) == np.sign(c) else (a + c) / 2.0
return out
def _as_1d_float(values: Sequence[float] | np.ndarray, name: str) -> np.ndarray:
array = np.asarray(values, dtype=float)
if array.ndim == 0:
raise ValueError(f"{name} must be array-like.")
if array.ndim > 2 or (array.ndim == 2 and min(array.shape) != 1):
raise ValueError(f"{name} can only have one dimension.")
return array.reshape(-1)
def _normalize_signal_matrix(data: Sequence[float] | Sequence[Sequence[float]] | np.ndarray, n_time: int) -> np.ndarray:
matrix = np.asarray(data, dtype=float)
if matrix.ndim == 0:
raise ValueError("Data must be array-like.")
if matrix.ndim == 1:
matrix = matrix.reshape(-1, 1)
elif matrix.ndim != 2:
raise ValueError("Data must be one- or two-dimensional.")
if matrix.shape[0] == n_time:
return matrix.astype(float, copy=True)
if matrix.shape[1] == n_time:
return matrix.T.astype(float, copy=True)
raise ValueError("Data dimensions do not match the time vector specified.")
def _coerce_zero_based_indices(values: Sequence[int] | np.ndarray, upper: int) -> list[int]:
out: list[int] = []
for raw in np.asarray(values).reshape(-1):
index = int(raw)
if index < 0 or index >= upper:
raise IndexError("Signal index out of range.")
out.append(index)
return out
def _roundn(values: Sequence[float] | np.ndarray, decimals: int) -> np.ndarray:
return np.round(np.asarray(values, dtype=float), decimals=max(int(decimals), 0))
def _matlab_mode_1d(values: Sequence[float] | np.ndarray) -> float:
array = np.asarray(values, dtype=float).reshape(-1)
if array.size == 0:
return np.nan
unique, counts = np.unique(array, return_counts=True)
best = np.flatnonzero(counts == np.max(counts))
return float(unique[int(best[0])])
def _nearest_sample_matrix(target_time: np.ndarray, source_time: np.ndarray, source_data: np.ndarray) -> np.ndarray:
target = np.asarray(target_time, dtype=float).reshape(-1)
source_t = np.asarray(source_time, dtype=float).reshape(-1)
source = np.asarray(source_data, dtype=float)
if source.ndim == 1:
source = source[:, None]
if source_t.size == 0:
return np.zeros((target.size, source.shape[1]), dtype=float)
right = np.searchsorted(source_t, target, side="left")
right = np.clip(right, 0, source_t.size - 1)
left = np.clip(right - 1, 0, source_t.size - 1)
choose_right = np.abs(source_t[right] - target) <= np.abs(source_t[left] - target)
indices = np.where(choose_right, right, left)
return source[indices]
[docs]
class SignalObj:
"""Multi-dimensional time-series signal object (Matlab ``SignalObj``).
``SignalObj`` is the foundational data container in nSTAT. It stores
one or more signal channels sampled on a common time axis, along with
metadata (name, units, labels) and supports element-wise arithmetic,
resampling, filtering, correlation analysis, and spectral estimation.
Parameters
----------
time : array_like
Monotonically increasing time vector of length *n*.
data : array_like
Signal values. Shape ``(n,)`` for a scalar signal or ``(n, d)``
for a *d*-dimensional signal.
name : str, optional
Human-readable signal name (used as y-axis label in plots).
xlabelval : str, optional
X-axis label string (default ``'time'``).
xunits : str, optional
X-axis unit string (default ``'s'``).
yunits : str, optional
Y-axis unit string.
dataLabels : sequence of str or str, optional
Per-dimension labels. A single string is broadcast to all
dimensions.
plotProps : sequence or str, optional
Per-dimension Matplotlib format strings.
See Also
--------
Covariate : SignalObj subclass with confidence-interval support.
nspikeTrain : Point-process (spike train) companion class.
"""
[docs]
def __init__(
self,
time: Sequence[float],
data: Sequence[float] | Sequence[Sequence[float]] | np.ndarray,
name: str = "",
xlabelval: str = "time",
xunits: str = "s",
yunits: str = "",
dataLabels: Sequence[str] | str | None = None,
plotProps: Sequence[Any] | str | None = None,
**kwargs,
) -> None:
"""Construct a multi-dimensional signal object (Matlab ``SignalObj``).
Parameters
----------
time : array_like, shape (n,)
Monotonically increasing time vector in **seconds**.
data : array_like, shape (n,) or (n, d) or (d, n)
Signal values. Either a 1-D vector of length *n* for a scalar
signal, or a 2-D matrix. Both ``(n, d)`` and ``(d, n)``
orientations are accepted — the layout is inferred from the
time-vector length.
name : str, optional
Human-readable signal name. Used as the y-axis label in
plots. Default ``""``.
xlabelval : str, optional
X-axis label string. Default ``"time"``. Also accepts the
keyword alias ``xlabel=`` for backwards compatibility.
xunits : str, optional
X-axis unit string. Default ``"s"`` (seconds).
yunits : str, optional
Y-axis unit string. Default ``""``.
dataLabels : sequence of str or str or None, optional
Per-dimension labels. A single string is broadcast to all
*d* dimensions. Also accepts the keyword alias
``data_labels=`` for backwards compatibility.
plotProps : sequence or str or None, optional
Per-dimension Matplotlib format strings (one per channel).
**kwargs
Reserved for the ``xlabel`` / ``data_labels`` legacy aliases
above. Any other keyword raises :class:`TypeError`.
Notes
-----
The sample rate is inferred from ``mean(diff(time))``. For
single-sample signals (or non-finite/zero deltas) it falls back
to ``1000.0`` Hz, matching the MATLAB toolbox.
Raises
------
ValueError
If *time* is not a 1-D array-like, or if *data* dimensions do
not match *time* length.
TypeError
If unexpected keyword arguments are passed.
See Also
--------
Covariate : SignalObj subclass with confidence-interval support.
nspikeTrain : Point-process companion class for spike data.
"""
if "xlabel" in kwargs and "xlabelval" not in kwargs:
xlabelval = kwargs.pop("xlabel")
if "data_labels" in kwargs and dataLabels is None:
dataLabels = kwargs.pop("data_labels")
if kwargs:
unexpected = ", ".join(sorted(kwargs))
raise TypeError(f"Unexpected keyword arguments: {unexpected}")
self.time = _as_1d_float(time, "Time vector")
self.data = _normalize_signal_matrix(data, self.time.size)
self.name = str(name)
self.xlabelval = str(xlabelval)
self.xunits = str(xunits)
self.yunits = str(yunits)
self.minTime = float(np.min(self.time)) if self.time.size else 0.0
self.maxTime = float(np.max(self.time)) if self.time.size else 0.0
if self.time.size > 1:
delta_t = float(np.mean(np.diff(self.time)))
else:
delta_t = np.nan
if not np.isfinite(delta_t) or delta_t <= 0:
delta_t = 0.001
self.sampleRate = float(1.0 / delta_t)
self.origSampleRate = float(self.sampleRate)
self.originalTime = self.time.copy()
self.originalData = self.data.copy()
self.dataMask = np.ones(self.dimension, dtype=int)
self.plotProps: list[Any] = []
self.setPlotProps(plotProps)
self.setDataLabels(dataLabels if dataLabels is not None else "")
self.conf_interval: tuple[np.ndarray, np.ndarray] | None = None
@property
def dimension(self) -> int:
"""Number of signal channels (columns in the data matrix)."""
return int(self.data.shape[1])
@property
def values(self) -> np.ndarray:
"""Signal data as a 1-D array (scalar) or 2-D matrix."""
if self.dimension == 1:
return self.data[:, 0]
return self.data
@property
def units(self) -> str:
"""Y-axis unit string (alias for ``yunits``)."""
return self.yunits
@property
def sample_rate(self) -> float:
"""Sampling rate in Hz (alias for ``sampleRate``)."""
return float(self.sampleRate)
def _spawn(
self,
time: np.ndarray,
data: np.ndarray,
*,
data_labels: Sequence[str] | None = None,
plot_props: Sequence[Any] | None = None,
) -> "SignalObj":
labels = list(self.dataLabels) if data_labels is None else list(data_labels)
props = list(self.plotProps) if plot_props is None else list(plot_props)
return self.__class__(
np.asarray(time, dtype=float).copy(),
np.asarray(data, dtype=float).copy(),
self.name,
self.xlabelval,
self.xunits,
self.yunits,
labels,
props,
)
[docs]
def copySignal(self) -> "SignalObj":
"""Return a deep copy of this signal (Matlab ``copySignal``)."""
copied = self._spawn(self.time, self.data)
if self.conf_interval is not None:
copied.conf_interval = (
np.asarray(self.conf_interval[0], dtype=float).copy(),
np.asarray(self.conf_interval[1], dtype=float).copy(),
)
copied.dataMask = np.asarray(self.dataMask, dtype=int).copy()
copied.originalTime = self.originalTime.copy()
copied.originalData = self.originalData.copy()
copied.sampleRate = float(self.sampleRate)
copied.origSampleRate = float(self.origSampleRate)
copied.minTime = float(self.minTime)
copied.maxTime = float(self.maxTime)
return copied
def _binary_operand_matrix(self, other) -> tuple[np.ndarray, list[str]]:
if isinstance(other, SignalObj):
if self.time.shape != other.time.shape or np.max(np.abs(self.time - other.time)) > 1e-9:
raise ValueError("Signals must share an identical time grid for arithmetic operations.")
data = np.asarray(other.data, dtype=float)
if data.shape[1] == 1 and self.dimension > 1:
data = np.repeat(data, self.dimension, axis=1)
labels = list(self.dataLabels)
elif self.dimension == 1 and data.shape[1] > 1:
labels = list(other.dataLabels)
elif data.shape[1] == self.dimension:
labels = list(self.dataLabels)
else:
raise ValueError("Signal dimensions must match for arithmetic operations.")
return data, labels
values = np.asarray(other, dtype=float)
if values.ndim == 0:
data = np.full(self.data.shape, float(values), dtype=float)
return data, list(self.dataLabels)
if values.ndim == 1:
if values.size == self.time.size:
return values.reshape(-1, 1), list(self.dataLabels if self.dimension == 1 else [self.dataLabels[0]])
if values.size == self.dimension:
return np.tile(values.reshape(1, -1), (self.time.size, 1)), list(self.dataLabels)
if values.ndim == 2 and values.shape[0] == self.time.size:
return values, list(self.dataLabels[: values.shape[1]])
raise ValueError("Unsupported arithmetic operand for SignalObj")
def _binary_op(self, other, op) -> "SignalObj":
other_matrix, labels = self._binary_operand_matrix(other)
left = self.data
if left.shape[1] == 1 and other_matrix.shape[1] > 1:
left = np.repeat(left, other_matrix.shape[1], axis=1)
labels = labels if labels else list(self.dataLabels[: other_matrix.shape[1]])
if other_matrix.shape[1] == 1 and left.shape[1] > 1:
other_matrix = np.repeat(other_matrix, left.shape[1], axis=1)
labels = list(self.dataLabels)
result = op(left, other_matrix)
return self._spawn(self.time, result, data_labels=labels)
[docs]
def setName(self, name: str) -> None:
"""Set the signal name (y-axis label)."""
if not isinstance(name, str):
raise TypeError("Name must be a string!")
self.name = name
[docs]
def setXlabel(self, name: str) -> None:
"""Set the x-axis label string."""
self.xlabelval = str(name)
[docs]
def setYLabel(self, name: str) -> None:
"""Set the y-axis label (alias for ``setName``)."""
self.setName(name)
[docs]
def setUnits(self, xUnits: str, yUnits: str | None = None) -> None:
"""Set x-axis and optionally y-axis units."""
if yUnits is not None:
self.setYUnits(yUnits)
self.setXUnits(xUnits)
[docs]
def setXUnits(self, units: str) -> None:
"""Set the x-axis unit string."""
if isinstance(units, str):
self.xunits = units
[docs]
def setYUnits(self, units: str) -> None:
"""Set the y-axis unit string."""
if isinstance(units, str):
self.yunits = units
[docs]
def setSampleRate(self, sampleRate: float) -> None:
"""Set the sample rate, resampling the data if it differs from current."""
requested = float(sampleRate)
current = float(self.sampleRate)
if abs(round(requested, 3) - round(current, 3)) > 0:
self.resampleMe(requested)
[docs]
def setDataLabels(self, dataLabels: Sequence[str] | str | None) -> None:
"""Set per-dimension data labels.
A single string is broadcast to all dimensions. A sequence must
have length equal to ``dimension``.
"""
if dataLabels is None or (isinstance(dataLabels, str) and dataLabels == ""):
self.dataLabels = ["" for _ in range(self.dimension)]
return
if isinstance(dataLabels, str):
self.dataLabels = [dataLabels for _ in range(self.dimension)]
return
labels = [str(label) for label in dataLabels]
if len(labels) != self.dimension:
raise ValueError("Need the number of labels to match the number of dimensions of the SignalObj")
self.dataLabels = labels
[docs]
def setMinTime(self, minTime: float | None = None, holdVals: int = 0) -> None:
"""Extend or trim the signal to start at *minTime*.
If *holdVals* is 1, endpoint values are held when extending;
otherwise the signal is zero-padded.
"""
target = self.time[0] if minTime is None else float(minTime)
timeVec = self.getTime()
if target < float(np.min(timeVec)):
maxTime = float(np.max(timeVec))
dt = 1.0 / self.sampleRate
# MATLAB SignalObj.m:302 uses ``minTime:1/sampleRate:maxTime`` here.
# ``np.arange`` overshoots for stop values just below an integer
# multiple of dt; ``_matlab_colon`` matches MATLAB exactly.
newTime = _matlab_colon(target, dt, maxTime)
numSamples = int(newTime.size - timeVec.size)
if holdVals == 1:
pad = np.tile(self.data[0:1, :], (numSamples, 1))
else:
pad = np.zeros((numSamples, self.dimension), dtype=float)
self.data = np.vstack([pad, self.data])
self.time = newTime
elif target > float(np.min(timeVec)):
startIndex = self.findNearestTimeIndex(target)
self.time = self.time[startIndex:]
self.data = self.data[startIndex:, :]
self.minTime = float(np.min(self.time))
[docs]
def setMaxTime(self, maxTime: float | None = None, holdVals: int = 0) -> None:
"""Extend or trim the signal to end at *maxTime*.
If *holdVals* is 1, endpoint values are held when extending;
otherwise the signal is zero-padded.
"""
target = self.time[-1] if maxTime is None else float(maxTime)
timeVec = self.getTime()
if float(np.max(timeVec)) < target:
minTime = float(np.min(timeVec))
n_samples = int(float(self.sampleRate) * (target - minTime) + 1.0)
n_samples = max(n_samples, timeVec.size)
newTime = np.linspace(minTime, target, n_samples, dtype=float)
numSamples = int(newTime.size - timeVec.size)
if holdVals == 1:
pad = np.tile(self.data[-1:, :], (numSamples, 1))
else:
pad = np.zeros((numSamples, self.dimension), dtype=float)
self.data = np.vstack([self.data, pad])
self.time = newTime
elif float(np.max(timeVec)) > target:
endIndex = self.findNearestTimeIndex(target) + 1
self.time = self.time[:endIndex]
self.data = self.data[:endIndex, :]
self.maxTime = float(np.max(self.time))
[docs]
def setPlotProps(self, plotProps: Sequence[Any] | str | None, index: int | None = None) -> None:
"""Set per-dimension Matplotlib format strings.
When *index* (0-based) is given, only that dimension is updated.
"""
if index is None:
if plotProps is None:
self.plotProps = [None for _ in range(self.dimension)]
elif isinstance(plotProps, str):
self.plotProps = [plotProps for _ in range(self.dimension)]
else:
props = list(plotProps)
if len(props) == 0:
self.plotProps = [None for _ in range(self.dimension)]
return
if len(props) == 1 and self.dimension > 1:
props = props * self.dimension
if len(props) != self.dimension:
raise ValueError("plotProps length must match signal dimension.")
self.plotProps = props
return
indices = _coerce_zero_based_indices([index], self.dimension)
target = indices[0] - 1
if not self.plotProps:
self.plotProps = [None for _ in range(self.dimension)]
if isinstance(plotProps, Sequence) and not isinstance(plotProps, str):
props = list(plotProps)
self.plotProps[target] = props[0] if props else None
else:
self.plotProps[target] = plotProps
[docs]
def setDataMask(self, dataMask: Sequence[int] | np.ndarray) -> None:
"""Set binary data mask (1 = visible, 0 = hidden) for each dimension."""
mask = np.asarray(dataMask, dtype=int).reshape(-1)
if mask.size != self.dimension:
raise ValueError("dataMask must match the number of signal dimensions.")
if np.any((mask != 0) & (mask != 1)):
raise ValueError("dataMask must be binary.")
self.dataMask = mask
[docs]
def setMaskByInd(self, index: Sequence[int] | np.ndarray) -> None:
"""Enable only the dimensions at the given 0-based indices."""
selected = _coerce_zero_based_indices(index, self.dimension)
mask = np.zeros(self.dimension, dtype=int)
mask[np.asarray(selected, dtype=int)] = 1
self.setDataMask(mask)
[docs]
def setMaskByLabels(self, labels: Sequence[str] | str) -> None:
"""Enable only the dimensions whose data labels match *labels*."""
indices = self.getIndicesFromLabels(labels)
if isinstance(indices, list) and indices and isinstance(indices[0], list):
flat = [item for sub in indices for item in sub]
elif isinstance(indices, list):
flat = indices
else:
flat = [indices]
self.setMaskByInd(flat)
[docs]
def setMask(self, mask: Sequence[int] | Sequence[str] | np.ndarray | None = None) -> None:
"""Flexible mask setter accepting indices, labels, or a binary vector.
``None`` clears the mask (all hidden). A binary vector of length
``dimension`` is used directly. A list of labels or 0-based indices
enables only those dimensions.
"""
if mask is None:
self.setDataMask(np.zeros(self.dimension, dtype=int))
return
if isinstance(mask, str):
self.setMaskByLabels(mask)
return
values = list(mask)
if not values:
self.setDataMask(np.zeros(self.dimension, dtype=int))
return
first = values[0]
if isinstance(first, str):
self.setMaskByLabels(values)
return
arr = np.asarray(values)
if arr.size == self.dimension and np.all(np.isin(arr, [0, 1])):
self.setDataMask(arr.astype(int))
return
self.setMaskByInd(arr.astype(int))
[docs]
def getTime(self) -> np.ndarray:
"""Return a copy of the time vector."""
return self.time.copy()
[docs]
def getData(self) -> np.ndarray:
"""Return signal data as a matrix (alias for ``dataToMatrix()``)."""
return self.dataToMatrix()
[docs]
def getOriginalData(self) -> tuple[np.ndarray, np.ndarray]:
"""Return ``(originalTime, originalData)`` copies."""
return self.originalTime.copy(), self.originalData.copy()
[docs]
def getOrigDataSig(self) -> "SignalObj":
"""Return the original (pre-resample) data as a new ``SignalObj``."""
return self._spawn(self.originalTime, self.originalData)
[docs]
def getValueAt(self, x: Sequence[float] | float) -> np.ndarray:
"""Return signal value(s) at time(s) *x* via nearest-neighbour lookup."""
query = np.asarray(x, dtype=float).reshape(-1)
out = np.zeros((query.size, self.dimension), dtype=float)
valid = (query >= self.minTime) & (query <= self.maxTime)
if np.any(valid):
q_valid = query[valid]
right = np.searchsorted(self.time, q_valid, side="left")
right = np.clip(right, 0, self.time.size - 1)
left = np.clip(right - 1, 0, self.time.size - 1)
choose_right = np.abs(self.time[right] - q_valid) <= np.abs(self.time[left] - q_valid)
indices = np.where(choose_right, right, left)
out[valid] = self.data[indices]
return out[0] if np.isscalar(x) else out
[docs]
def getPlotProps(self, index: int) -> Any:
"""Return the plot property for dimension *index* (0-based)."""
idx = _coerce_zero_based_indices([index], self.dimension)[0]
return self.plotProps[idx]
[docs]
def getIndexFromLabel(self, label: str) -> list[int]:
"""Return 0-based indices of dimensions whose label equals *label*."""
matches = [i for i, value in enumerate(self.dataLabels) if value == label]
if not matches:
raise ValueError("Label does not exist!")
return matches
[docs]
def getIndicesFromLabels(self, label: Sequence[str] | str):
"""Return 0-based index(es) for one or more data-label strings."""
if isinstance(label, str):
matches = self.getIndexFromLabel(label)
return matches[0] if len(matches) == 1 else matches
out = [self.getIndexFromLabel(str(item)) for item in label]
counts = [len(item) for item in out]
if counts and max(counts) == 1:
return [item[0] for item in out]
return out
[docs]
def areDataLabelsEmpty(self) -> bool:
"""Return ``True`` if all data labels are empty strings.
Matches Matlab ``SignalObj.areDataLabelsEmpty()``.
"""
return all(not str(label) for label in self.dataLabels)
[docs]
def isLabelPresent(self, label: str) -> bool:
"""Return ``True`` if *label* matches any data label or equals ``'all'``.
Matches Matlab ``SignalObj.isLabelPresent()``.
"""
if str(label).lower() == "all":
return True
try:
self.getIndexFromLabel(label)
return True
except ValueError:
return False
[docs]
def convertNamesToIndices(self, selectorArray) -> list[int] | np.ndarray:
"""Convert label names (or mixed) to 0-based indices.
Matches Matlab ``SignalObj.convertNamesToIndices()`` (the
Python port is 0-based throughout).
"""
if isinstance(selectorArray, str):
if selectorArray == "all":
return list(range(self.dimension))
if self.isLabelPresent(selectorArray):
return self.getIndexFromLabel(selectorArray)
raise ValueError(f"Specified label '{selectorArray}' does not match data label")
if isinstance(selectorArray, (int, float, np.integer)):
return [int(selectorArray)]
if isinstance(selectorArray, np.ndarray):
return selectorArray.astype(int).ravel().tolist()
if isinstance(selectorArray, (list, tuple)):
result: list[int] = []
for item in selectorArray:
if isinstance(item, str):
if self.isLabelPresent(item):
result.extend(self.getIndexFromLabel(item))
else:
result.append(int(item))
return result
return list(range(self.dimension))
def _selector_to_zero_based(self, selectorArray: Sequence[int] | np.ndarray | None) -> np.ndarray:
"""Resolve a selector specification to a 0-based numpy index array.
The name is preserved for historical continuity, but inputs are
now expected to be 0-based; this method validates and normalises
(no ``- 1`` translation).
"""
if selectorArray is None:
if self.isMaskSet():
selected = self.findIndFromDataMask()
else:
selected = list(range(self.dimension))
else:
if isinstance(selectorArray, str):
selected = self.getIndicesFromLabels(selectorArray)
else:
selected = selectorArray
indices = np.asarray(selected, dtype=int).reshape(-1)
if indices.size == 0:
return np.array([], dtype=int)
if np.min(indices) < 0 or np.max(indices) >= self.dimension:
raise IndexError("Signal index out of range.")
return indices
[docs]
def dataToMatrix(self, selectorArray: Sequence[int] | np.ndarray | None = None) -> np.ndarray:
"""Return signal data as an ``(n, d)`` matrix.
*selectorArray* is an optional sequence of 0-based dimension
indices. When ``None``, the data mask selects visible dimensions.
"""
indices = self._selector_to_zero_based(selectorArray)
if indices.size == 0:
return np.zeros((self.time.size, 0), dtype=float)
return self.data[:, indices]
def _labels_for_indices(self, zero_based: np.ndarray) -> list[str]:
return [self.dataLabels[int(i)] for i in zero_based]
def _plot_props_for_indices(self, zero_based: np.ndarray) -> list[Any]:
if not self.plotProps:
return [None for _ in zero_based]
return [self.plotProps[int(i)] for i in zero_based]
[docs]
def getSubSignalFromInd(self, selectorArray: Sequence[int] | np.ndarray) -> "SignalObj":
"""Return a new ``SignalObj`` with only the selected dimensions (0-based)."""
indices = self._selector_to_zero_based(selectorArray)
return self._spawn(
self.time,
self.data[:, indices],
data_labels=self._labels_for_indices(indices),
plot_props=self._plot_props_for_indices(indices),
)
[docs]
def getSubSignalFromNames(self, labels: Sequence[str] | str) -> "SignalObj":
"""Return a sub-signal selected by data-label name(s)."""
indices = self.getIndicesFromLabels(labels)
return self.getSubSignalFromInd(indices if isinstance(indices, list) else [indices])
[docs]
def getSubSignal(self, identifier) -> "SignalObj":
"""Return a sub-signal selected by labels, indices, or mixed."""
if isinstance(identifier, str):
return self.getSubSignalFromNames(identifier)
if isinstance(identifier, np.ndarray):
values = identifier.reshape(-1).tolist()
elif isinstance(identifier, Sequence):
values = list(identifier)
else:
values = [identifier]
if values and isinstance(values[0], str):
return self.getSubSignalFromNames(values)
return self.getSubSignalFromInd(values)
[docs]
def findNearestTimeIndex(self, time: float) -> int:
"""Return the 0-based index of the sample nearest to *time*."""
value = float(time)
if value < self.minTime:
return 0
if value > self.maxTime:
return self.time.size - 1
right = int(np.searchsorted(self.time, value, side="left"))
if right <= 0:
return 0
if right >= self.time.size:
return self.time.size - 1
left = right - 1
if abs(self.time[right] - value) <= abs(self.time[left] - value):
return right
return left
[docs]
def findNearestTimeIndices(self, times: Sequence[float] | np.ndarray) -> np.ndarray:
"""Return 0-based indices of the samples nearest to each time in *times*."""
return np.asarray([self.findNearestTimeIndex(value) for value in np.asarray(times, dtype=float).reshape(-1)], dtype=int)
[docs]
def merge(self, other: "SignalObj", holdVals: int = 0) -> "SignalObj":
"""Merge *other* signal columns into *self*.
Matlab calls ``makeCompatible`` first so that signals with
different time grids are reconciled automatically. The Python
version now does the same.
Parameters
----------
other : SignalObj
Signal whose data columns will be appended.
holdVals : int, optional
Passed to ``makeCompatible`` – ``1`` holds endpoint values
when the time range is extended; ``0`` (default) pads with
zeros.
"""
s1c, s2c = self.makeCompatible(other, holdVals)
merged = s1c._spawn(
s1c.time,
np.column_stack([s1c.data, s2c.data]),
data_labels=[*s1c.dataLabels, *list(s2c.dataLabels)],
plot_props=[*s1c.plotProps, *getattr(s2c, "plotProps", [None for _ in range(s2c.dimension)])],
)
return merged
def __add__(self, other) -> "SignalObj":
return self._binary_op(other, np.add)
def __radd__(self, other) -> "SignalObj":
return self + other
def __sub__(self, other) -> "SignalObj":
return self._binary_op(other, np.subtract)
def __rsub__(self, other) -> "SignalObj":
other_matrix, labels = self._binary_operand_matrix(other)
left = other_matrix
right = self.data
if left.shape[1] == 1 and right.shape[1] > 1:
left = np.repeat(left, right.shape[1], axis=1)
labels = list(self.dataLabels)
if right.shape[1] == 1 and left.shape[1] > 1:
right = np.repeat(right, left.shape[1], axis=1)
return self._spawn(self.time, np.subtract(left, right), data_labels=labels)
def __pos__(self) -> "SignalObj":
return self.copySignal()
def __neg__(self) -> "SignalObj":
return self._spawn(self.time, -self.data, data_labels=list(self.dataLabels))
[docs]
def power(self, exponent: float) -> "SignalObj":
"""Element-wise power ``data ** exponent`` (Matlab ``power``)."""
return self.__class__(
self.time.copy(),
np.power(self.data, float(exponent)),
f"{self.name}^{exponent}",
self.xlabelval,
self.xunits,
self.yunits,
list(self.dataLabels),
list(self.plotProps),
)
[docs]
def sqrt(self) -> "SignalObj":
"""Element-wise square root (Matlab ``sqrt``)."""
return self.power(0.5)
def __mul__(self, other) -> "SignalObj":
return self._binary_op(other, np.multiply)
def __rmul__(self, other) -> "SignalObj":
return self * other
def __matmul__(self, other) -> "SignalObj":
"""Matrix multiply (``@`` operator). Matches Matlab ``mtimes``."""
if isinstance(other, SignalObj):
return self._spawn(self.time, self.data * other.data, data_labels=list(self.dataLabels))
other_arr = np.asarray(other, dtype=float)
result = (self.data.T @ other_arr).T if other_arr.ndim <= 1 else self.data @ other_arr
return self._spawn(self.time[:result.shape[0]] if result.ndim == 2 else self.time, result)
def __truediv__(self, other) -> "SignalObj":
return self._binary_op(other, np.divide)
def __rtruediv__(self, other) -> "SignalObj":
other_matrix, labels = self._binary_operand_matrix(other)
left = other_matrix
right = self.data
if left.shape[1] == 1 and right.shape[1] > 1:
left = np.repeat(left, right.shape[1], axis=1)
labels = list(self.dataLabels)
if right.shape[1] == 1 and left.shape[1] > 1:
right = np.repeat(right, left.shape[1], axis=1)
return self._spawn(self.time, np.divide(left, right), data_labels=labels)
[docs]
def ldivide(self, other) -> "SignalObj":
r"""Element-wise left division (Matlab ``.\``): ``other ./ self``.
Matches Matlab ``SignalObj.ldivide()``.
"""
return self._binary_op(other, lambda a, b: np.divide(b, a))
def __invert__(self) -> "SignalObj":
"""Conjugate transpose (Matlab ``ctranspose`` / ``'``).
For real-valued signals this is equivalent to :attr:`T`.
"""
new_data = self.data.T
new_time = self.time[:new_data.shape[0]] if new_data.shape[0] != self.time.size else self.time
return self._spawn(new_time, new_data)
@property
def T(self) -> "SignalObj":
"""Transpose the data matrix. Matches Matlab ``ctranspose`` / ``transpose``."""
new_data = self.data.T
new_time = self.time[:new_data.shape[0]] if new_data.shape[0] != self.time.size else self.time
return self._spawn(new_time, new_data)
[docs]
def clearPlotProps(self, index=None) -> None:
"""Clear plot properties. Matches Matlab ``clearPlotProps``."""
if index is None:
index = list(range(self.dimension))
else:
index = [i - 1 for i in np.atleast_1d(index)]
for i in index:
if i < len(self.plotProps):
self.plotProps[i] = None
[docs]
def plotPropsSet(self) -> bool:
"""Return ``True`` if any plot property is non-empty.
Matches Matlab ``SignalObj.plotPropsSet()``.
"""
for prop in self.plotProps:
if prop is not None and str(prop) != "":
return True
return False
[docs]
def getSigInTimeWindow(
self,
wMin: Sequence[float] | float | None = None,
wMax: Sequence[float] | float | None = None,
holdVals: int = 0,
) -> "SignalObj":
"""Extract signal within ``[wMin, wMax]``.
Multiple windows can be specified by passing equal-length sequences
for *wMin* and *wMax*; the extracted segments are concatenated as
additional dimensions (Matlab ``getSigInTimeWindow``).
"""
if wMax is None:
wMax = self.maxTime
if wMin is None:
wMin = self.minTime
min_values = np.asarray([wMin] if np.isscalar(wMin) else wMin, dtype=float).reshape(-1)
max_values = np.asarray([wMax] if np.isscalar(wMax) else wMax, dtype=float).reshape(-1)
if min_values.size != max_values.size:
raise ValueError("Window minTimes must contain the same number of elements as window maxTimes")
if min_values.size == 1 and self.minTime == float(min_values[0]) and self.maxTime == float(max_values[0]):
return self.copySignal()
windowed: SignalObj | None = None
for idx, (left, right) in enumerate(zip(min_values, max_values), start=1):
current = self.copySignal()
if left < current.minTime:
current.setMinTime(left, holdVals)
if right > current.maxTime:
current.setMaxTime(right, holdVals)
start = current.findNearestTimeIndex(left)
stop = current.findNearestTimeIndex(right) + 1
current.time = current.time[start:stop]
current.data = current.data[start:stop, :]
labels = list(current.dataLabels)
if min_values.size > 1:
labels = [f"{label}_{{{idx}}}" for label in labels]
current.setDataLabels(labels)
current.setMinTime()
current.setMaxTime()
windowed = current if windowed is None else windowed.merge(current)
return windowed if windowed is not None else self.copySignal()
[docs]
def restoreToOriginal(self, rMask: int = 0) -> None:
"""Restore time, data, and sample rate to their original values.
If *rMask* is 1, the data mask is also reset (all visible).
"""
self.time = self.originalTime.copy()
self.data = self.originalData.copy()
self.minTime = float(np.min(self.time))
self.maxTime = float(np.max(self.time))
self.sampleRate = float(self.origSampleRate)
if rMask == 1:
self.resetMask()
[docs]
def resetMask(self) -> None:
"""Reset the data mask so all dimensions are visible."""
self.dataMask = np.ones(self.dimension, dtype=int)
[docs]
def findIndFromDataMask(self) -> list[int]:
"""Return 0-based indices of dimensions currently visible (mask == 1)."""
return [int(index) for index in np.flatnonzero(self.dataMask == 1)]
[docs]
def isMaskSet(self) -> bool:
"""Return ``True`` if any dimension is currently masked out."""
return bool(np.any(self.dataMask == 0))
[docs]
def derivative(self) -> "SignalObj":
"""Return the per-sample forward-difference derivative.
Matches MATLAB ``SignalObj.m:761`` which is a method (called as
``s.derivative()``), not a property. The previous Python
implementation was a ``@property`` so MATLAB-ported code calling
``sig.derivative()`` silently failed with ``TypeError: 'SignalObj'
object is not callable``. This is now a regular method; access
``sig.derivative()`` (with parentheses) to compute the derivative.
"""
deriv = np.zeros_like(self.data, dtype=float)
if self.data.shape[0] > 1:
deriv[1:, :] = np.diff(self.data, axis=0) * float(self.sampleRate)
deriv[~np.isfinite(deriv)] = 0.0
labels = [f"d_{label}" if label else "" for label in self.dataLabels]
return self._spawn(self.time, deriv, data_labels=labels)
[docs]
def derivativeAt(self, x0: Sequence[float] | float):
"""Return the derivative value(s) at time(s) *x0*."""
deriv = self.derivative()
values = deriv.getValueAt(x0)
return values
[docs]
def integral(self, t0: float | None = None, tf: float | None = None) -> "SignalObj":
"""Cumulative integral of the signal from *t0* to *tf*.
Computed via a causal IIR accumulator:
``y[n] = y[n-1] + x[n] * deltaT``. If *t0* / *tf* are not
specified, ``minTime`` / ``maxTime`` are used.
"""
start = self.minTime if t0 is None else float(t0)
stop = self.maxTime if tf is None else float(tf)
integrated = self.getSigInTimeWindow(start, stop)
dt = 1.0 / max(float(integrated.sampleRate), 1e-12)
integrated = integrated.filter([dt], [1.0, -1.0])
if integrated.yunits and integrated.xunits:
integrated.setYUnits(f"{integrated.yunits}*{integrated.xunits}")
elif integrated.xunits:
integrated.setYUnits(integrated.xunits)
dtstr = " d\\tau"
integrated.setName(f"\\int_{integrated.minTime:g}^{integrated.xlabelval[:1]}\\!\\!{{{integrated.name}{dtstr}}}")
labels_empty = all(not str(label) for label in integrated.dataLabels)
if not labels_empty:
updated_labels: list[str] = []
for label in self.dataLabels:
if label:
updated_labels.append(f"\\int_{integrated.minTime:g}^{integrated.xlabelval[:1]}\\!\\!{{{label}{dtstr}}}")
else:
updated_labels.append("")
integrated.setDataLabels(updated_labels)
return integrated
[docs]
def filter(self, B, A=1) -> "SignalObj":
"""Apply a causal IIR/FIR filter ``(B, A)`` to each dimension.
Equivalent to ``scipy.signal.lfilter(B, A, data)``.
"""
try:
from scipy.signal import lfilter
except Exception as exc: # pragma: no cover
raise ImportError("scipy is required for SignalObj.filter") from exc
b = np.asarray(B, dtype=float).reshape(-1)
a = np.asarray(A, dtype=float).reshape(-1)
filtered = np.column_stack([lfilter(b, a, self.data[:, index]) for index in range(self.dimension)])
return self._spawn(self.time, filtered, data_labels=list(self.dataLabels))
[docs]
def filtfilt(self, B, A=1) -> "SignalObj":
"""Apply a zero-phase IIR/FIR filter ``(B, A)`` to each dimension.
Equivalent to ``scipy.signal.filtfilt(B, A, data)``.
"""
try:
from scipy.signal import filtfilt
except Exception as exc: # pragma: no cover
raise ImportError("scipy is required for SignalObj.filtfilt") from exc
b = np.asarray(B, dtype=float).reshape(-1)
a = np.asarray(A, dtype=float).reshape(-1)
filtered = np.column_stack([filtfilt(b, a, self.data[:, index]) for index in range(self.dimension)])
return self._spawn(self.time, filtered, data_labels=list(self.dataLabels))
[docs]
def makeCompatible(self, other: "SignalObj", holdVals: int = 0) -> tuple["SignalObj", "SignalObj"]:
if (
self.minTime == other.minTime
and self.maxTime == other.maxTime
and round(float(self.sampleRate), 9) == round(float(other.sampleRate), 9)
and self.time.shape == other.time.shape
and np.max(np.abs(self.time - other.time)) <= 1e-9
):
return self, other
s1c = self.copySignal()
s2c = other.copySignal()
min_time = min(s1c.minTime, s2c.minTime)
max_time = max(s1c.maxTime, s2c.maxTime)
sample_rate = max(float(s1c.sampleRate), float(s2c.sampleRate))
s1c.setSampleRate(sample_rate)
s2c.setSampleRate(sample_rate)
s1c.setMinTime(min_time, holdVals)
s2c.setMinTime(min_time, holdVals)
s1c.setMaxTime(max_time, holdVals)
s2c.setMaxTime(max_time, holdVals)
s2c.data = _nearest_sample_matrix(s1c.time, s2c.time, s2c.data)
s2c.time = s1c.time.copy()
s2c.minTime = float(np.min(s2c.time))
s2c.maxTime = float(np.max(s2c.time))
return s1c, s2c
[docs]
def abs(self) -> "SignalObj":
"""Element-wise absolute value (Matlab ``abs``)."""
labels = [f"|{label}|" if label else "" for label in self.dataLabels]
return self._spawn(self.time, np.abs(self.data), data_labels=labels).with_metadata(
name=f"|{self.name}|",
yunits=self.yunits,
)
def __abs__(self) -> "SignalObj":
return self.abs()
[docs]
def log(self) -> "SignalObj":
"""Element-wise natural logarithm (Matlab ``log``)."""
labels = [f"ln({label})" if label else "" for label in self.dataLabels]
yunits = f"ln({self.yunits})" if self.yunits else ""
return self._spawn(self.time, np.log(self.data), data_labels=labels).with_metadata(
name=f"ln({self.name})",
yunits=yunits,
)
[docs]
def mode(self, axis: int | None = None) -> "SignalObj":
"""Column-wise mode of signal data (Matlab ``mode``)."""
axis_arg = 0 if axis is None else axis
if axis_arg == 0:
mode_data = np.asarray([_matlab_mode_1d(self.data[:, i]) for i in range(self.dimension)], dtype=float)
elif axis_arg == 1:
mode_data = np.asarray([_matlab_mode_1d(row) for row in self.data], dtype=float)
else:
raise ValueError("axis must be 0, 1, or None")
array = np.asarray(mode_data, dtype=float)
if array.ndim == 1 and array.size == self.dimension:
labels = [f"mode({label})" if label else "" for label in self.dataLabels]
return self._spawn(
np.asarray([self.time[0], self.time[-1]], dtype=float),
np.vstack([array, array]),
data_labels=labels,
).with_metadata(name=f"mode({self.name})")
reshaped = array.reshape(-1, 1)
return self._spawn(self.time, reshaped, data_labels=[f"mode({self.name})"]).with_metadata(name=f"mode({self.name})")
[docs]
def mean(self, axis: int | None = None) -> "SignalObj":
"""Column-wise mean (default) or row-wise mean of signal data.
``mean()`` or ``mean(0)`` computes the mean of each component
across time. ``mean(1)`` computes the mean value at each time
point across dimensions.
"""
axis_arg = 0 if axis is None else axis
mean_data = np.mean(self.data, axis=axis_arg)
array = np.asarray(mean_data, dtype=float)
if array.ndim == 1 and array.size == self.dimension:
labels = [f"\\mu({label})" if label else "" for label in self.dataLabels]
return self._spawn(
np.asarray([self.time[0], self.time[-1]], dtype=float),
np.vstack([array, array]),
data_labels=labels,
)
reshaped = array.reshape(-1, 1)
return self._spawn(self.time, reshaped, data_labels=[f"\\mu({self.name})"])
[docs]
def std(self, axis: int | None = None) -> "SignalObj":
"""Column-wise standard deviation (sample, ddof=1) of signal data.
``std()`` or ``std(0)`` computes std of each component across
time. ``std(1)`` computes std at each time point across dimensions.
"""
axis_arg = 0 if axis is None else axis
std_data = np.std(self.data, axis=axis_arg, ddof=1)
array = np.asarray(std_data, dtype=float)
if array.ndim == 1 and array.size == self.dimension:
labels = [f"\\sigma({label})" if label else "" for label in self.dataLabels]
return self._spawn(
np.asarray([self.time[0], self.time[-1]], dtype=float),
np.vstack([array, array]),
data_labels=labels,
)
reshaped = array.reshape(-1, 1)
return self._spawn(self.time, reshaped, data_labels=[f"\\sigma({self.name})"])
[docs]
def max(self, axis: int | None = None):
"""Return ``(values, indices, times)`` of column-wise maxima."""
axis_arg = 0 if axis is None else axis
values = np.max(self.data, axis=axis_arg)
indices = np.argmax(self.data, axis=axis_arg)
time = self.time[np.asarray(indices, dtype=int)]
return values, indices, time
[docs]
def min(self, axis: int | None = None):
"""Return ``(values, indices, times)`` of column-wise minima."""
axis_arg = 0 if axis is None else axis
values = np.min(self.data, axis=axis_arg)
indices = np.argmin(self.data, axis=axis_arg)
time = self.time[np.asarray(indices, dtype=int)]
return values, indices, time
[docs]
def resample(self, sample_rate: float) -> "SignalObj":
"""Return a resampled copy at *sample_rate* Hz."""
copied = self.copySignal()
copied.resampleMe(sample_rate)
return copied
[docs]
def resampleMe(self, newSampleRate: float) -> None:
"""Resample data in-place to *newSampleRate* Hz via cubic interpolation."""
try:
from scipy.interpolate import interp1d
except Exception as exc: # pragma: no cover
raise ImportError("scipy is required for SignalObj.resampleMe") from exc
rate = float(newSampleRate)
if rate <= 0:
raise ValueError("sampleRate must be > 0.")
if self.sampleRate == rate:
return
self.restoreToOriginal()
dt = 1.0 / rate
newTime = np.arange(self.time[0], self.time[-1] + 0.5 * dt, dt, dtype=float)
if self.data.shape[0] > 1:
columns = []
if self.time.size >= 4:
interp_kind = "cubic"
elif self.time.size == 3:
interp_kind = "quadratic"
else:
interp_kind = "linear"
for index in range(self.dimension):
interpolator = interp1d(
self.time,
self.data[:, index],
kind=interp_kind,
bounds_error=False,
fill_value=0.0,
)
columns.append(np.asarray(interpolator(newTime), dtype=float))
newData = np.column_stack(columns)
else:
newData = np.asarray(self.data, dtype=float).copy()
self.time = newTime
self.data = newData
self.sampleRate = rate
self.minTime = float(np.min(newTime))
self.maxTime = float(np.max(newTime))
[docs]
def autocorrelation(self) -> "SignalObj":
"""Normalized auto-correlation for each signal dimension.
Returns a new ``SignalObj`` whose time axis is lag (in the original
x-units) and whose data are the correlation coefficients normalised
to unity at lag zero (Matlab ``autocorrelation``).
"""
centered = self.data - np.mean(self.data, axis=0, keepdims=True)
columns: list[np.ndarray] = []
lags: np.ndarray | None = None
for index in range(self.dimension):
series = centered[:, index]
denom = float(np.dot(series, series))
corr = np.correlate(series, series, mode="full")
if denom > 0:
corr = corr / denom
else:
corr = np.zeros_like(corr, dtype=float)
if lags is None:
lags = np.arange(-series.size + 1, series.size, dtype=float) / max(float(self.sampleRate), 1e-12)
columns.append(np.asarray(corr, dtype=float))
data = np.column_stack(columns) if columns else np.zeros((0, 0), dtype=float)
return self.__class__(
lags if lags is not None else np.array([], dtype=float),
data,
f"ACF({self.name})",
"Lag",
self.xunits,
f"{self.yunits}^2" if self.yunits else "",
list(self.dataLabels),
list(self.plotProps),
)
[docs]
def crosscorrelation(self, other: "SignalObj") -> "SignalObj":
"""Normalized cross-correlation between two scalar signals.
Both signals must be one-dimensional. The result is normalised
so that the peak equals the Pearson correlation coefficient
(Matlab ``crosscorrelation``).
"""
if self.dimension != 1 or other.dimension != 1:
raise ValueError("crosscorrelation only supports one-dimensional signals")
s1c, s2c = self.makeCompatible(other)
x = s1c.data[:, 0] - float(np.mean(s1c.data[:, 0]))
y = s2c.data[:, 0] - float(np.mean(s2c.data[:, 0]))
denom = float(np.sqrt(np.dot(x, x) * np.dot(y, y)))
corr = np.correlate(x, y, mode="full")
if denom > 0:
corr = corr / denom
else:
corr = np.zeros_like(corr, dtype=float)
lags = np.arange(-x.size + 1, x.size, dtype=float) / max(float(s1c.sampleRate), 1e-12)
return self.__class__(
lags,
corr,
f"XCORF({self.name})",
"Lag",
self.xunits,
f"{self.yunits}^2" if self.yunits else "",
list(self.dataLabels[:1]),
list(self.plotProps[:1]),
)
[docs]
def xcorr(self, other: "SignalObj" | None = None, maxlag: int | None = None) -> "SignalObj":
"""Raw (un-normalised) cross-correlation (Matlab ``xcorr``).
Computes pairwise cross-correlation for all dimension pairs.
When *other* is ``None`` (auto-correlation), only non-negative
lags are returned. *maxlag* truncates to ``|lag| ≤ maxlag``
samples.
"""
s2 = self if other is None else other
s1c, s2c = self.makeCompatible(s2)
data_columns: list[np.ndarray] = []
data_labels: list[str] = []
lag_index: np.ndarray | None = None
for left_index in range(s1c.dimension):
for right_index in range(s2c.dimension):
corr = np.correlate(s1c.data[:, left_index], s2c.data[:, right_index], mode="full")
lags = np.arange(-s1c.data.shape[0] + 1, s1c.data.shape[0], dtype=int)
if maxlag is not None:
keep = np.abs(lags) <= int(maxlag)
corr = corr[keep]
lags = lags[keep]
if other is None:
keep = lags >= 0
corr = corr[keep]
lags = lags[keep]
if lag_index is None:
lag_index = lags.astype(float) / max(float(s1c.sampleRate), 1e-12)
data_columns.append(np.asarray(corr, dtype=float))
left_label = s1c.dataLabels[left_index] if left_index < len(s1c.dataLabels) else str(left_index + 1)
right_label = s2c.dataLabels[right_index] if right_index < len(s2c.dataLabels) else str(right_index + 1)
data_labels.append(f"corr({left_label},{right_label})")
data = np.column_stack(data_columns) if data_columns else np.zeros((0, 0), dtype=float)
name = f"corr({self.name},{s2.name})"
return self.__class__(
lag_index if lag_index is not None else np.array([], dtype=float),
data,
name,
"\\Delta \\tau",
self.xunits,
f"{self.yunits}^2" if self.yunits else "",
data_labels,
)
# ------------------------------------------------------------------
# Time-shift / alignment helpers (match Matlab SignalObj)
# ------------------------------------------------------------------
[docs]
def shift(self, deltaT: float, updateLabels: bool = False) -> "SignalObj":
"""Return a copy with time shifted by *deltaT* seconds."""
new_time = self.time + float(deltaT)
out = self.__class__(
new_time,
self.data.copy(),
self.name,
self.xlabelval,
self.xunits,
self.yunits,
list(self.dataLabels),
list(self.plotProps),
)
if updateLabels:
out.name = f"{self.name} shifted by {deltaT}"
return out
[docs]
def shiftMe(self, deltaT: float, updateLabels: bool = False) -> None:
"""In-place time shift by *deltaT* seconds (Matlab ``shiftMe``)."""
self.time = self.time + float(deltaT)
self.minTime = float(np.min(self.time)) if self.time.size else 0.0
self.maxTime = float(np.max(self.time)) if self.time.size else 0.0
if updateLabels:
self.name = f"{self.name} shifted by {deltaT}"
[docs]
def alignTime(self, timeMarker: float, newTime: float = 0.0) -> None:
"""Shift so that *timeMarker* becomes *newTime* (Matlab ``alignTime``).
Only shifts if *timeMarker* falls within ``[minTime, maxTime]``,
matching the Matlab implementation's bounds check.
"""
if self.minTime <= float(timeMarker) <= self.maxTime:
self.shiftMe(float(newTime) - float(timeMarker))
# ------------------------------------------------------------------
# ------------------------------------------------------------------
# Peak-finding helpers (match Matlab SignalObj)
# ------------------------------------------------------------------
[docs]
def findPeaks(
self,
peak_type: str = "maxima",
minDistance: int | None = None,
) -> tuple[list[np.ndarray], list[np.ndarray]]:
"""Find local peaks in each signal dimension.
Parameters
----------
peak_type : ``'maxima'`` or ``'minima'``
minDistance : minimum sample distance between peaks (default:
``sampleRate * duration / 10``).
Returns
-------
indices, values : lists of arrays (one per dimension).
Note: The Matlab original has a bug where the ``'minima'`` branch
does not negate the data before calling ``findpeaks``. This Python
port fixes that.
"""
from scipy.signal import find_peaks as _find_peaks
data = np.atleast_2d(self.data)
if data.shape[0] == 1:
data = data.T
N = data.shape[0]
if minDistance is None:
duration = float(self.maxTime - self.minTime)
minDistance = max(1, int(round(self.sampleRate * duration / 10)))
all_indices: list[np.ndarray] = []
all_values: list[np.ndarray] = []
for col in range(data.shape[1]):
sig = data[:, col]
if peak_type == "minima":
sig = -sig
idx, _ = _find_peaks(sig, distance=minDistance)
all_indices.append(idx)
all_values.append(data[idx, col]) # always return actual values
return all_indices, all_values
[docs]
def findMaxima(self) -> tuple[list[np.ndarray], list[np.ndarray]]:
"""Convenience wrapper: ``findPeaks('maxima')``."""
return self.findPeaks("maxima")
[docs]
def findMinima(self) -> tuple[list[np.ndarray], list[np.ndarray]]:
"""Convenience wrapper: ``findPeaks('minima')``."""
return self.findPeaks("minima")
[docs]
def findGlobalPeak(
self, peak_type: str = "maxima"
) -> tuple[np.ndarray, np.ndarray]:
"""Find the global max or min across each dimension.
Returns
-------
times : 1-D array of times at which the peak occurs (one per dim).
values : 1-D array of peak values (one per dim).
Note: The Matlab original has a typo (``sOBj`` instead of ``sObj``)
in the minima branch. This Python port fixes that.
"""
data = np.atleast_2d(self.data)
if data.shape[0] == 1:
data = data.T
if peak_type == "maxima":
idx = np.argmax(data, axis=0)
else:
idx = np.argmin(data, axis=0)
times = self.time[idx]
values = data[idx, np.arange(data.shape[1])]
return np.atleast_1d(times), np.atleast_1d(values)
# ------------------------------------------------------------------
# Alignment / windowing (match Matlab SignalObj)
# ------------------------------------------------------------------
[docs]
def alignToMax(self) -> tuple["SignalObj", float]:
"""Align all dimensions so their peaks coincide at the mean peak time.
Returns ``(aligned_signal, mean_peak_time)``.
Matches Matlab ``SignalObj.alignToMax()``.
"""
peak_times, _ = self.findGlobalPeak("maxima")
mean_time = float(np.mean(peak_times))
delta_t = -(peak_times - mean_time)
aligned = self.getSubSignal(0).shift(float(delta_t[0]))
for i in range(1, self.dimension):
aligned = aligned.merge(self.getSubSignal(i).shift(float(delta_t[i])))
return aligned, mean_time
[docs]
def windowedSignal(self, windowTimes) -> "SignalObj":
"""Extract and concatenate windowed segments.
Matches Matlab ``SignalObj.windowedSignal()``.
"""
windowTimes = np.asarray(windowTimes, dtype=float).ravel()
result = None
for i in range(len(windowTimes) - 1):
seg = self.getSigInTimeWindow(float(windowTimes[i]), float(windowTimes[i + 1]))
if i == 0:
result = seg
else:
seg = seg.shift(-float(windowTimes[i]))
result = result.merge(seg)
return result if result is not None else self.copySignal()
[docs]
def normWindowedSignal(
self,
windowTimes,
numPoints: int = 100,
lbound: float | None = None,
ubound: float | None = None,
) -> "SignalObj":
"""Normalize windowed signal segments to a common time axis.
Matches Matlab ``SignalObj.normWindowedSignal()``.
"""
windowTimes = np.asarray(windowTimes, dtype=float).ravel()
columns: list[np.ndarray] = []
for i in range(len(windowTimes) - 1):
minT = float(windowTimes[i])
maxT = float(windowTimes[i + 1])
dur = abs(maxT - minT)
if lbound is not None and ubound is not None:
if dur > ubound or dur < lbound:
continue
seg = self.getSigInTimeWindow(minT, maxT)
norm_time = np.linspace(minT, maxT, numPoints)
# Matlab uses interp1(..., 'nearest', 0) — nearest-neighbor with 0-fill
from scipy.interpolate import interp1d as _interp1d
_ifn = _interp1d(seg.time, seg.data[:, 0], kind="nearest",
bounds_error=False, fill_value=0.0)
interp_data = _ifn(norm_time)
columns.append(interp_data)
if not columns:
return self.copySignal()
data = np.column_stack(columns)
act_time = np.arange(numPoints, dtype=float) / float(numPoints)
labels = list(self.dataLabels[:1]) * data.shape[1]
return self.__class__(act_time, data, self.name, self.xlabelval, "%", self.yunits, labels)
[docs]
def getSubSignalsWithinNStd(self, nStd: float = 2.0) -> tuple["SignalObj", np.ndarray]:
"""Return sub-signals within *nStd* standard deviations of the mean.
Returns ``(filtered_signal, selected_indices)``.
Matches Matlab ``SignalObj.getSubSignalsWithinNStd()``.
"""
mean_sig = np.mean(self.data, axis=1)
std_sig = np.std(self.data, axis=1, ddof=1)
min_val = mean_sig - nStd * std_sig
max_val = mean_sig + nStd * std_sig
# A column passes if ALL rows are within [minVal, maxVal]
above_min = np.all(self.data >= min_val[:, None], axis=0)
below_max = np.all(self.data <= max_val[:, None], axis=0)
sig_index = np.flatnonzero(above_min & below_max)
if sig_index.size == 0:
return self.copySignal(), sig_index
# 0-based indices for getSubSignal
return self.getSubSignal((sig_index + 1).tolist()), sig_index
# ------------------------------------------------------------------
# Variability plots (match Matlab SignalObj)
# ------------------------------------------------------------------
[docs]
def plotAllVariability(
self,
faceColor=None,
linewidth: float = 3.0,
ciUpper: float | np.ndarray = 1.96,
ciLower: float | np.ndarray | None = None,
ax=None,
):
"""Plot mean ± CI shaded area. Matches Matlab ``plotAllVariability``.
Parameters
----------
faceColor : color, optional
Fill colour (default: tab:blue).
linewidth : float
Width of mean line.
ciUpper, ciLower : float or array
Number of std-devs (scalar) or explicit bounds (array).
ax : matplotlib Axes, optional
"""
import matplotlib.pyplot as plt
if faceColor is None:
faceColor = "tab:blue"
if ciLower is None:
ciLower = ciUpper
if ax is None:
ax = plt.gca()
mean_sig = np.mean(self.data, axis=1)
std_sig = np.std(self.data, axis=1, ddof=1)
ciUpper_arr = np.atleast_1d(ciUpper).ravel()
ciLower_arr = np.atleast_1d(ciLower).ravel()
if ciUpper_arr.size == 1:
ci_top = mean_sig + float(ciUpper_arr[0]) * std_sig
else:
ci_top = mean_sig + ciUpper_arr[:len(mean_sig)]
if ciLower_arr.size == 1:
ci_bottom = mean_sig - float(ciLower_arr[0]) * std_sig
else:
ci_bottom = mean_sig - ciLower_arr[:len(mean_sig)]
ax.fill_between(self.time, ci_bottom, ci_top, color=faceColor, alpha=0.5, edgecolor="none")
(h,) = ax.plot(self.time, mean_sig, "k-", linewidth=linewidth)
return h
[docs]
def plotVariability(self, selectorArray=None, ax=None):
"""Plot mean ± CI for each label group. Matches Matlab ``plotVariability``.
Parameters
----------
selectorArray : list of list[int] or list[int], optional
ax : matplotlib Axes, optional
"""
import matplotlib.pyplot as plt
if ax is None:
ax = plt.gca()
if selectorArray is None:
if not self.areDataLabelsEmpty():
unique_labels = list(dict.fromkeys(self.dataLabels))
selectorArray = [self.getIndexFromLabel(lbl) for lbl in unique_labels]
else:
selectorArray = list(range(self.dimension))
# MATLAB R2014b+ "lines" palette — mirrors the default prop_cycle that
# apply_matlab_style() installs, so axes built via prop_cycle and axes
# built here look consistent.
_TAB_COLORS = [
"#0072BD", "#D95319", "#EDB120", "#7E2F8E",
"#77AC30", "#4DBEEE", "#A2142F", "#A0A0A0",
]
handles = []
if isinstance(selectorArray, list) and selectorArray and isinstance(selectorArray[0], (list, tuple, np.ndarray)):
for i, sel in enumerate(selectorArray):
h = self.getSubSignal(sel).plotAllVariability(
faceColor=_TAB_COLORS[i % len(_TAB_COLORS)], ax=ax
)
handles.append(h)
else:
h = self.getSubSignal(selectorArray).plotAllVariability(ax=ax)
handles.append(h)
return handles
# ------------------------------------------------------------------
# Cross-covariance (match Matlab SignalObj.xcov)
# ------------------------------------------------------------------
[docs]
def xcov(self, other: "SignalObj | None" = None, maxlag: int | None = None,
scaleOpt: str = "none") -> "SignalObj":
"""Cross-covariance (mean-removed xcorr). Matches Matlab ``xcov``.
When called with no *other* argument (auto-covariance), only
non-negative lags are returned — matching Matlab behaviour where
``data=tempC(M-1:end,index)`` and ``lags=tempLags(M-1:end)``.
"""
auto = other is None
s1 = self
s2 = self if auto else other
s1c, s2c = s1.makeCompatible(s2)
data_columns: list[np.ndarray] = []
data_labels: list[str] = []
lag_index: np.ndarray | None = None
for li in range(s1c.dimension):
for ri in range(s2c.dimension):
x = s1c.data[:, li] - np.mean(s1c.data[:, li])
y = s2c.data[:, ri] - np.mean(s2c.data[:, ri])
corr = np.correlate(x, y, mode="full")
N = len(x)
lags = np.arange(-N + 1, N, dtype=int)
# scale
if scaleOpt == "biased":
corr = corr / N
elif scaleOpt == "unbiased":
denom = N - np.abs(lags)
denom[denom <= 0] = 1
corr = corr / denom
elif scaleOpt == "coeff":
corr = corr / corr[N - 1] if corr[N - 1] != 0 else corr
if maxlag is not None:
keep = np.abs(lags) <= int(maxlag)
corr = corr[keep]
lags = lags[keep]
# Matlab returns only non-negative lags for auto-covariance
if auto:
nonneg = lags >= 0
corr = corr[nonneg]
lags = lags[nonneg]
if lag_index is None:
lag_index = lags.astype(float) / max(float(s1c.sampleRate), 1e-12)
data_columns.append(np.asarray(corr, dtype=float))
ll = s1c.dataLabels[li] if li < len(s1c.dataLabels) else str(li + 1)
rl = s2c.dataLabels[ri] if ri < len(s2c.dataLabels) else str(ri + 1)
data_labels.append(f"xcov({ll},{rl})")
data = np.column_stack(data_columns) if data_columns else np.zeros((0, 0), dtype=float)
return self.__class__(
lag_index if lag_index is not None else np.array([], dtype=float),
data,
f"xcov({self.name},{s2.name})",
"\\Delta \\tau",
self.xunits,
f"{self.yunits}^2" if self.yunits else "",
data_labels,
)
# ------------------------------------------------------------------
# Spectral methods (match Matlab SignalObj)
# ------------------------------------------------------------------
[docs]
def periodogram(self, NFFT: int | None = None) -> tuple[np.ndarray, np.ndarray]:
"""Power spectral density via periodogram (Matlab ``periodogram``).
Loops over all signal dimensions like the Matlab implementation.
When *NFFT* is ``None`` we match MATLAB's default of
``max(256, 2**nextpow2(N))`` rather than scipy's default of
``N``, producing the denser frequency grid (and characteristic
zero-notch comb) seen in the MATLAB reference figures.
Returns ``(frequencies, psd)`` where *psd* has shape
``(nfreqs,)`` for 1-D signals or ``(nfreqs, dimension)`` for
multi-dimensional signals.
"""
from scipy.signal import periodogram as _periodogram
fs = float(self.sampleRate)
N = self.data.shape[0]
if NFFT is None:
NFFT = max(256, int(2 ** np.ceil(np.log2(N))))
psd_cols: list[np.ndarray] = []
f_out: np.ndarray | None = None
ndim = self.dimension
for i in range(ndim):
x = self.data[:, i] if self.data.ndim == 2 else self.data
f, Pxx = _periodogram(x, fs=fs, nfft=NFFT, window="boxcar",
scaling="density")
if f_out is None:
f_out = f
psd_cols.append(Pxx)
if ndim == 1:
return f_out, psd_cols[0]
return f_out, np.column_stack(psd_cols)
[docs]
def MTMspectrum(self, NW: float = 4.0, NFFT: int | None = None,
Pval: float = 0.95,
Kmax: int | None = None,
) -> tuple[np.ndarray, np.ndarray, np.ndarray | None]:
"""Multi-taper spectral estimate (Matlab ``MTMspectrum``).
Uses discrete prolate spheroidal sequences (DPSS / Slepian tapers).
Loops over all signal dimensions like the Matlab implementation.
Parameters
----------
NW : float
Time-bandwidth product (default 4).
NFFT : int, optional
FFT length (default next power of 2 >= N).
Pval : float, optional
Confidence level for the chi-squared confidence interval
(default 0.95). Set to ``None`` to skip CI computation.
Kmax : int, optional
Number of tapers (default ``2*NW - 1``).
Returns
-------
frequencies : ndarray
psd : ndarray
Shape ``(nfreqs,)`` for 1-D or ``(nfreqs, dimension)``.
psd_ci : ndarray or None
Shape ``(nfreqs, 2)`` for 1-D or ``(nfreqs, 2*dimension)``
containing ``[lower, upper]`` columns per dimension.
``None`` when *Pval* is ``None``.
"""
from scipy.signal.windows import dpss
N = self.data.shape[0]
fs = float(self.sampleRate)
if Kmax is None:
Kmax = int(2 * NW - 1)
if NFFT is None:
# MATLAB's pmtm pads to max(256, 2**nextpow2(2*N)) for a denser
# frequency grid -- a plain 2**nextpow2(N) misses the harmonic
# peaks of v2 = sin(v1^2) visible in SignalObjExamples_14.png.
NFFT = max(256, int(2 ** np.ceil(np.log2(2 * N))))
tapers, eigenvalues = dpss(N, NW, Kmax, return_ratios=True)
frequencies = np.fft.rfftfreq(NFFT, d=1.0 / fs)
nfreqs = len(frequencies)
# chi-squared CI bounds (degrees of freedom = 2*Kmax)
ci_lo_factor = ci_hi_factor = None
if Pval is not None:
from scipy.stats import chi2
dof = 2 * Kmax
alpha = 1.0 - Pval
ci_lo_factor = dof / chi2.ppf(1.0 - alpha / 2.0, dof)
ci_hi_factor = dof / chi2.ppf(alpha / 2.0, dof)
ndim = self.dimension
psd_cols: list[np.ndarray] = []
ci_cols: list[np.ndarray] = []
for di in range(ndim):
x = self.data[:, di] if self.data.ndim == 2 else self.data
Sk = np.zeros((Kmax, nfreqs))
for k in range(Kmax):
xw = x * tapers[k]
Xf = np.fft.rfft(xw, n=NFFT)
Sk[k] = np.abs(Xf) ** 2
weights = eigenvalues / eigenvalues.sum()
psd = np.dot(weights, Sk) * (2.0 / fs)
psd[0] /= 2.0
if NFFT % 2 == 0:
psd[-1] /= 2.0
psd_cols.append(psd)
if Pval is not None:
ci_cols.append(psd * ci_lo_factor)
ci_cols.append(psd * ci_hi_factor)
if ndim == 1:
psd_out = psd_cols[0]
ci_out = np.column_stack(ci_cols) if ci_cols else None
else:
psd_out = np.column_stack(psd_cols)
ci_out = np.column_stack(ci_cols) if ci_cols else None
return frequencies, psd_out, ci_out
[docs]
def spectrogram(self, nperseg: int = 256, noverlap: int | None = None,
NFFT: int | None = None,
window: str = "hann") -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Short-time Fourier transform spectrogram (Matlab ``spectrogram``).
Returns ``(frequencies, times, Sxx)``.
"""
from scipy.signal import spectrogram as _spectrogram
x = self.data[:, 0] if self.data.ndim == 2 else self.data
fs = float(self.sampleRate)
if noverlap is None:
noverlap = nperseg // 2
if NFFT is None:
NFFT = nperseg
f, t, Sxx = _spectrogram(x, fs=fs, window=window,
nperseg=nperseg, noverlap=noverlap,
nfft=NFFT)
# offset times to match signal start
t = t + self.minTime
return f, t, Sxx
[docs]
def setConfInterval(self, bounds: tuple[np.ndarray, np.ndarray]) -> None:
"""Attach ``(lower, upper)`` confidence bounds aligned with time."""
low, high = bounds
low_arr = np.asarray(low, dtype=float)
high_arr = np.asarray(high, dtype=float)
if low_arr.shape[0] != self.time.shape[0] or high_arr.shape[0] != self.time.shape[0]:
raise ValueError("confidence interval bounds must align with time.")
self.conf_interval = (low_arr, high_arr)
[docs]
def dataToStructure(self, selectorArray: Sequence[int] | np.ndarray | None = None) -> dict[str, Any]:
"""Serialize signal data to a plain dict (Matlab ``dataToStructure``)."""
data = self.dataToMatrix(selectorArray)
plot_props = list(self.plotProps)
if all(prop is None for prop in plot_props):
plot_props = []
return {
"time": self.time.tolist(),
"data": data.tolist(),
"name": self.name,
"xlabelval": self.xlabelval,
"xunits": self.xunits,
"yunits": self.yunits,
"dataLabels": list(self.dataLabels),
"plotProps": plot_props,
}
[docs]
def toStructure(self) -> dict[str, Any]:
"""Serialize the full signal to a plain dict (Matlab ``toStructure``)."""
return self.dataToStructure()
[docs]
@staticmethod
def signalFromStruct(structure: dict[str, Any]) -> "SignalObj":
"""Reconstruct a ``SignalObj`` from a dict (Matlab ``signalFromStruct``)."""
return SignalObj(
structure["time"],
structure["data"],
structure.get("name", ""),
structure.get("xlabelval", "time"),
structure.get("xunits", "s"),
structure.get("yunits", ""),
structure.get("dataLabels"),
structure.get("plotProps"),
)
[docs]
def plot(self, selectorArray=None, plotPropsIn=None, handle=None):
"""Plot selected signal dimensions (Matlab ``plot``).
Parameters
----------
selectorArray : optional
Dimension selector (labels, 0-based indices, or ``None`` for all
visible dimensions).
plotPropsIn : optional
Override Matplotlib format strings for each dimension.
handle : matplotlib Axes, optional
Axes to draw into; defaults to ``plt.gca()``.
Returns
-------
list of Line2D
"""
import matplotlib.pyplot as plt
from .confidence_interval import MATLAB_COLOR_ORDER
ax = plt.gca() if handle is None else handle
signal = self.getSubSignal(selectorArray) if selectorArray is not None else self.getSubSignal(self.findIndFromDataMask() or list(range(self.dimension)))
props = signal.plotProps if plotPropsIn is None else list(plotPropsIn)
if len(props) == 1 and signal.dimension > 1:
props = props * signal.dimension
if not props:
props = [None for _ in range(signal.dimension)]
lines = []
for index in range(signal.dimension):
kwargs = {}
prop = props[index]
if isinstance(prop, str) and prop:
kwargs["fmt"] = prop
elif prop is None:
kwargs["color"] = MATLAB_COLOR_ORDER[index % MATLAB_COLOR_ORDER.shape[0]]
if "fmt" in kwargs:
fmt = kwargs.pop("fmt")
line = ax.plot(signal.time, signal.data[:, index], fmt, **kwargs)
else:
line = ax.plot(signal.time, signal.data[:, index], **kwargs)
lines.extend(line)
if handle is None:
xunits = f" [{signal.xunits}]" if signal.xunits else ""
yunits = f" [{signal.yunits}]" if signal.yunits else ""
ax.set_xlabel(f"{signal.xlabelval}{xunits}")
ax.set_ylabel(f"{signal.name}{yunits}")
return lines
class Covariate(SignalObj):
"""Signal with per-dimension confidence intervals (Matlab ``Covariate``).
``Covariate`` extends :class:`SignalObj` with a list of
:class:`~nstat.confidence_interval.ConfidenceInterval` objects (one per
dimension) and propagates those intervals through ``+`` and ``-``
arithmetic. It also provides ``'zero-mean'`` and ``'standard'``
signal representations used by the GLM design-matrix builder.
Parameters
----------
*args, **kwargs
Forwarded to :class:`SignalObj`. The keyword aliases ``values``
(→ ``data``) and ``units`` (→ ``yunits``) are accepted for
convenience.
See Also
--------
SignalObj : Base time-series container.
ConfidenceInterval : CI storage class used by ``ci``.
"""
[docs]
def __init__(self, *args, **kwargs) -> None:
"""Construct a covariate signal (Matlab ``Covariate``).
Accepts the same positional / keyword arguments as
:class:`SignalObj`, with two convenience aliases for users coming
from earlier Python releases:
- ``values=`` is forwarded to the :class:`SignalObj` ``data=``
parameter.
- ``units=`` is forwarded to the :class:`SignalObj` ``yunits=``
parameter.
Parameters
----------
*args, **kwargs
Forwarded to :class:`SignalObj.__init__`. Time is in
**seconds** and the inferred sample rate is in **Hz**.
Notes
-----
After construction the ``ci`` attribute is ``None``. Attach
confidence intervals via :meth:`setConfInterval` or
:meth:`computeMeanPlusCI`.
See Also
--------
SignalObj : Base time-series container.
nstat.confidence_interval.ConfidenceInterval :
CI storage class used by ``ci``.
"""
if "values" in kwargs and "data" not in kwargs:
kwargs["data"] = kwargs.pop("values")
if "units" in kwargs and "yunits" not in kwargs:
kwargs["yunits"] = kwargs.pop("units")
super().__init__(*args, **kwargs)
self.ci: list[Any] | None = None
[docs]
def computeMeanPlusCI(self, alphaVal: float = 0.05) -> "Covariate":
"""Compute row-wise mean with empirical confidence intervals.
Treats each column as a replicate. Returns a scalar ``Covariate``
whose CI bounds are the *alphaVal*/2 and 1−*alphaVal*/2 quantiles
of the empirical CDF across replicates (Matlab ``computeMeanPlusCI``).
"""
from .confidence_interval import ConfidenceInterval
sorted_data = np.sort(self.data, axis=1)
n_rep = sorted_data.shape[1]
if n_rep == 0:
raise ValueError("Covariate must contain at least one column to compute confidence intervals.")
ecdf = np.arange(1, n_rep + 1, dtype=float) / float(n_rep)
lower = np.empty(sorted_data.shape[0], dtype=float)
upper = np.empty(sorted_data.shape[0], dtype=float)
for row_idx in range(sorted_data.shape[0]):
row = sorted_data[row_idx]
lower_idx = np.flatnonzero(ecdf < (alphaVal / 2.0))
upper_idx = np.flatnonzero(ecdf > (1.0 - alphaVal / 2.0))
lower[row_idx] = row[int(lower_idx[-1])] if lower_idx.size else row[0]
upper[row_idx] = row[int(upper_idx[0])] if upper_idx.size else row[-1]
confInt = ConfidenceInterval(self.time, np.column_stack([lower, upper]))
mean_signal = np.mean(self.data, axis=1)
newCov = Covariate(
self.time.copy(),
mean_signal,
self.name,
self.xlabelval,
self.xunits,
self.yunits,
[f"\\mu({self.name})" if self.name else "\\mu"],
)
newCov.setConfInterval(confInt)
return newCov
[docs]
def plot(self, selectorArray=None, plotPropsIn=None, handle=None):
"""Plot signal dimensions with shaded confidence intervals."""
lines = super().plot(selectorArray, plotPropsIn, handle)
if self.isConfIntervalSet():
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
ax = plt.gca() if handle is None else handle
selectors = self.findIndFromDataMask() if selectorArray is None else (
self.getIndicesFromLabels(selectorArray) if isinstance(selectorArray, str) else list(np.asarray(selectorArray).reshape(-1))
)
if not isinstance(selectors, list):
selectors = [selectors]
if selectors and isinstance(selectors[0], list):
selectors = [item[0] for item in selectors]
for line_index, selector in enumerate(selectors):
color = getattr(lines[line_index], "get_color", lambda: "b")()
if isinstance(color, (str, bytes)):
color = mcolors.to_rgb(color)
self.ci[selector].plot(color, ax=ax)
return lines
[docs]
def getSubSignal(self, identifier) -> "Covariate":
"""Return a sub-covariate preserving matching CIs."""
sub = super().getSubSignal(identifier)
cov = Covariate(
sub.time,
sub.data,
sub.name,
sub.xlabelval,
sub.xunits,
sub.yunits,
list(sub.dataLabels),
list(sub.plotProps),
)
if self.isConfIntervalSet():
selected: list[int] = []
for label in cov.dataLabels:
if label:
match = next((i for i, original in enumerate(self.dataLabels) if original == label), None)
if match is None:
raise ValueError("Unable to align Covariate confidence interval with sub-signal labels.")
selected.append(match)
else:
selected.append(len(selected))
cov.setConfInterval([self.ci[index] for index in selected])
return cov
[docs]
def getSigRep(self, repType: str = "standard") -> "Covariate":
"""Return a signal representation of this covariate.
Parameters
----------
repType : str
``'standard'`` returns ``self`` unchanged.
``'zero-mean'`` returns ``self - mean(self)`` with confidence
intervals propagated (Matlab parity: uses operator overload so
CIs shift by the same constant).
"""
rep = str(repType).strip().lower()
if rep == "standard":
return self
if rep == "zero-mean":
# Build a constant Covariate holding the per-column mean so that
# the CI-propagating __sub__ is invoked (Matlab: ``self - self.mu``).
mu_vals = np.mean(self.data, axis=0, keepdims=True)
mu_broadcast = np.repeat(mu_vals, len(self.time), axis=0)
mu_cov = Covariate(
self.time.copy(),
mu_broadcast,
self.name,
self.xlabelval,
self.xunits,
self.yunits,
list(self.dataLabels),
list(self.plotProps),
)
return self - mu_cov
raise ValueError("repType must be either 'zero-mean' or 'standard'")
[docs]
def get(self, identifier):
"""MATLAB ``get`` / subscript alias for :meth:`__getitem__`."""
return self.__getitem__(identifier)
[docs]
def filtfilt(self, B, A=1) -> "Covariate":
"""Apply a zero-phase filter (Matlab ``filtfilt`` override).
Thin wrapper around :meth:`SignalObj.filtfilt` that returns a
``Covariate`` so confidence-interval propagation continues to work
on the result.
"""
filtered = super().filtfilt(B, A)
cov = Covariate(
filtered.time,
filtered.data,
filtered.name,
filtered.xlabelval,
filtered.xunits,
filtered.yunits,
list(filtered.dataLabels),
list(filtered.plotProps),
)
return cov
[docs]
def toStructure(self) -> dict[str, Any]:
"""Serialize to a dict, including CI payload if present."""
structure = super().toStructure()
if self.isConfIntervalSet():
ci_payload: list[dict[str, Any]] = []
for item in self.ci or []:
if hasattr(item, "dataToStructure"):
ci_payload.append(item.dataToStructure())
if ci_payload:
structure["ci"] = ci_payload[0] if len(ci_payload) == 1 else ci_payload
return structure
[docs]
def isConfIntervalSet(self) -> bool:
"""Return ``True`` if at least one dimension has a CI attached."""
return bool(self.ci)
[docs]
def setConfInterval(self, ciObj) -> None:
"""Attach one or more ``ConfidenceInterval`` objects to this covariate."""
if isinstance(ciObj, list):
self.ci = list(ciObj)
else:
self.ci = [ciObj]
[docs]
def copySignal(self) -> "Covariate":
"""Deep-copy including confidence intervals (Matlab ``copySignal``)."""
copied = Covariate(
self.time.copy(),
self.data.copy(),
self.name,
self.xlabelval,
self.xunits,
self.yunits,
list(self.dataLabels),
list(self.plotProps),
)
copied.dataMask = np.asarray(self.dataMask, dtype=int).copy()
copied.originalTime = self.originalTime.copy()
copied.originalData = self.originalData.copy()
copied.sampleRate = float(self.sampleRate)
copied.origSampleRate = float(self.origSampleRate)
copied.minTime = float(self.minTime)
copied.maxTime = float(self.maxTime)
copied.ci = None if not self.ci else list(self.ci)
if self.conf_interval is not None:
copied.conf_interval = (
np.asarray(self.conf_interval[0], dtype=float).copy(),
np.asarray(self.conf_interval[1], dtype=float).copy(),
)
return copied
def __add__(self, other):
"""Add two covariates, propagating confidence intervals (Matlab ``plus``)."""
covOut = super().__add__(other)
if isinstance(other, Covariate):
if self.isConfIntervalSet() and not other.isConfIntervalSet():
covOut.setConfInterval([self.ci[index] + other.getSubSignal(index + 1) for index in range(self.dimension)])
elif self.isConfIntervalSet() and other.isConfIntervalSet():
covOut.setConfInterval([self.ci[index] + other.ci[index] for index in range(self.dimension)])
elif (not self.isConfIntervalSet()) and other.isConfIntervalSet():
covOut.setConfInterval([other.ci[index] + self.getSubSignal(index + 1) for index in range(other.dimension)])
return covOut
def __sub__(self, other):
"""Subtract two covariates, propagating confidence intervals (Matlab ``minus``)."""
covOut = super().__sub__(other)
if isinstance(other, Covariate):
if self.isConfIntervalSet() and not other.isConfIntervalSet():
covOut.setConfInterval([self.ci[index] - other.getSubSignal(index + 1) for index in range(self.dimension)])
elif self.isConfIntervalSet() and other.isConfIntervalSet():
covOut.setConfInterval([self.ci[index] - other.ci[index] for index in range(self.dimension)])
elif (not self.isConfIntervalSet()) and other.isConfIntervalSet():
covOut.setConfInterval([self.getSubSignal(index + 1) - other.ci[index] for index in range(other.dimension)])
return covOut
[docs]
def dataToStructure(self, selectorArray: Sequence[int] | np.ndarray | None = None) -> dict[str, Any]:
"""Serialize data payload (Matlab ``dataToStructure`` override).
Thin wrapper around :meth:`SignalObj.dataToStructure` kept here so
the MATLAB method-order parity audit pairs this Covariate-level
declaration with the MATLAB ``dataToStructure`` slot just before
``fromStructure``.
"""
return super().dataToStructure(selectorArray)
[docs]
@staticmethod
def fromStructure(structure: dict[str, Any]) -> "Covariate":
"""Reconstruct a ``Covariate`` (with optional CIs) from a dict."""
from .confidence_interval import ConfidenceInterval
cov = Covariate(
structure["time"],
structure["data"],
structure.get("name", ""),
structure.get("xlabelval", "time"),
structure.get("xunits", "s"),
structure.get("yunits", ""),
structure.get("dataLabels"),
structure.get("plotProps"),
)
ci_payload = structure.get("ci")
if ci_payload is None:
return cov
if isinstance(ci_payload, list):
cov.setConfInterval([ConfidenceInterval.fromStructure(item) for item in ci_payload])
elif isinstance(ci_payload, tuple):
cov.setConfInterval([ConfidenceInterval.fromStructure(item) for item in ci_payload])
else:
cov.setConfInterval(ConfidenceInterval.fromStructure(ci_payload))
return cov
# ------------------------------------------------------------------
# Python-only extensions (parity: python-extension).
# Declared after the MATLAB-mirrored block so the class_method_parity
# scanner sees the MATLAB-named methods in MATLAB source order.
# ------------------------------------------------------------------
@property
def mu(self) -> SignalObj: # parity: python-extension (Matlab dependent property)
"""Column-wise mean as a ``SignalObj`` (Matlab ``mu`` property)."""
return self.mean()
@property
def sigma(self) -> SignalObj: # parity: python-extension (Matlab dependent property)
"""Column-wise standard deviation as a ``SignalObj`` (Matlab ``sigma``)."""
return self.std()
def __getitem__(self, identifier): # parity: python-extension (subsref shim)
"""Subscript a covariate (Matlab ``get`` / ``subsref`` accessor)."""
return self.getSubSignal(identifier)
# MATLAB-name aliases so callers can use ``cov.plus(other)`` /
# ``cov.minus(other)`` directly (mirroring the MATLAB method names);
# the canonical operator overloads (``__add__`` / ``__sub__``) above
# carry the parity contract.
[docs]
def plus(self, other): # parity: python-extension (alias for __add__)
"""MATLAB ``plus`` alias for :meth:`__add__`."""
return self.__add__(other)
[docs]
def minus(self, other): # parity: python-extension (alias for __sub__)
"""MATLAB ``minus`` alias for :meth:`__sub__`."""
return self.__sub__(other)
# nspikeTrain was extracted to ``nstat._spike_train_impl`` as part of the
# v0.3.1 module split. We re-export it here so all legacy import paths
# (``from nstat.core import nspikeTrain``, ``from nstat.nspikeTrain import
# nspikeTrain``, ``from nstat import nspikeTrain``) continue to work.
from ._spike_train_impl import nspikeTrain # noqa: E402
# Backward-compatible alias used by earlier Python scaffolding.
SpikeTrain = nspikeTrain