"""Per-window statistical features of a time-domain signal. The cheapest feature front end there is: instead of transforming a window of samples, you reduce it to a handful of summary statistics computed straight from the raw values. No FFT, no filterbank, just sums and a few powers. These are the descriptors a microcontroller can afford to compute faster than the data arrives, and they are the classic front end for accelerometer activity recognition [Bao & Intille 2004] and for the energy / zero-crossing gate of a simple voice-activity detector [Rabiner & Juang 1993]. Each function takes a window ``x`` and reduces along the **last axis**, so a single window of shape ``(n,)`` yields a scalar and a stack of windows of shape ``(n_windows, n)`` yields one value per window of shape ``(n_windows,)``. This is the natural layout for a sliding window over a 1-D signal. The features fall into four families: - **Level**: ``mean`` (the DC / gravity component). - **Spread / power**: ``variance``, ``std``, ``rms``, ``energy``, ``peak_to_peak``. How big the swings are. These carry most of the discriminative power for "is anything happening". - **Shape**: ``skewness`` (asymmetry) and ``kurtosis`` (tailedness / peakiness), the third and fourth standardised moments. They separate distributions that share a mean and variance. - **Oscillation**: ``zero_crossing_rate``, a coarse, amplitude-free proxy for dominant frequency, the time-domain cousin of the spectral centroid. See the [zero-crossing](../../basics/zero-crossing.qmd) chapter for ZCR as a frequency estimator; here it is used as one feature among several. Conventions, chosen to match the float arithmetic of the on-device [`features.h`](embedded.qmd) library and SciPy's defaults: - ``variance`` and the moments are **population** statistics (divide by ``n``, i.e. ``ddof=0``), not sample statistics (``n-1``). On a fixed window the difference is a constant factor a classifier absorbs; the population form is what the integer/embedded code computes. - ``kurtosis`` is **excess** kurtosis (Fisher): 0 for a Gaussian, positive for heavy-tailed / peaky, negative for flat / bounded. Both moments use the biased (population) estimator, matching ``scipy.stats`` defaults. - ``energy`` is the sum of squares; ``rms`` is ``sqrt(energy / n)``. """ import numpy as np _EPS = 1e-12 def mean(x, axis=-1): """Arithmetic mean of the window: its DC level. For an accelerometer axis at rest this is the projection of gravity; for an audio frame it is the (usually near-zero) DC offset. On its own it rarely discriminates *activity*, but it anchors the central moments below. """ return np.mean(x, axis=axis) def variance(x, axis=-1): """Population variance (divide by ``n``): the mean squared deviation. The single most useful activity feature. A still sensor has tiny variance; a moving one has large variance, regardless of orientation. This is the feature the embedded decision-tree classifier on the [hardware page](embedded.qmd) keys on. """ return np.var(x, axis=axis) # numpy default ddof=0 (population) def std(x, axis=-1): """Standard deviation: ``sqrt(variance)``, back in the signal's own units.""" return np.std(x, axis=axis) def rms(x, axis=-1): """Root-mean-square level: ``sqrt(mean(x**2))``. The square root of the average power. For a zero-mean window it equals the standard deviation; for a window with a DC offset it also counts that offset, so ``rms**2 = mean**2 + variance``. A robust loudness / intensity proxy. """ x = np.asarray(x, dtype=float) return np.sqrt(np.mean(x * x, axis=axis)) def energy(x, axis=-1): """Sum of squares over the window: its total energy. Differs from ``rms`` only by the window length (``energy = n * rms**2``), so it is window-length dependent. Useful as a raw "how much signal is here" gate when every window is the same length. """ x = np.asarray(x, dtype=float) return np.sum(x * x, axis=axis) def peak_to_peak(x, axis=-1): """Range of the window: ``max - min``. A one-glance amplitude measure that, unlike variance, is dominated by the single largest excursion. Cheap, but sensitive to a lone spike or a clipped sample, which is exactly when it disagrees with variance. """ return np.ptp(x, axis=axis) def skewness(x, axis=-1): """Skewness: the third standardised moment (asymmetry of the distribution). Zero for a symmetric window, positive when the long tail is on the high side, negative when it is on the low side. Distinguishes, for instance, a one-sided impulse train from a symmetric oscillation that shares its variance. Returns 0 where the variance is ~0 (a constant window has no defined shape). """ x = np.asarray(x, dtype=float) mu = np.mean(x, axis=axis, keepdims=True) d = x - mu m2 = np.mean(d * d, axis=axis) m3 = np.mean(d ** 3, axis=axis) denom = m2 ** 1.5 return np.where(denom > _EPS, m3 / np.where(denom > _EPS, denom, 1.0), 0.0) def kurtosis(x, axis=-1): """Excess kurtosis: the fourth standardised moment minus 3 (Fisher). Zero for a Gaussian, **positive** for a heavy-tailed or peaky window (rare large excursions: impulsive noise, spikes, transients), **negative** for a flat-topped or bounded one (a sine wave sits at -1.5, uniform noise at -1.2). The classic discriminator for "is this clean background, or does it contain outliers". Returns 0 where the variance is ~0. Beware the mixture trap (see the page): adding a *wide* secondary component to a narrow one makes the mixture heavy-tailed (positive excess), even if that secondary component is platykurtic in isolation. Reason about the mixture, not the parts. """ x = np.asarray(x, dtype=float) mu = np.mean(x, axis=axis, keepdims=True) d = x - mu m2 = np.mean(d * d, axis=axis) m4 = np.mean(d ** 4, axis=axis) denom = m2 * m2 return np.where(denom > _EPS, m4 / np.where(denom > _EPS, denom, 1.0) - 3.0, 0.0) def zero_crossing_rate(x, axis=-1): """Fraction of adjacent sample pairs that change sign. A number in [0, 1]: the count of strict sign changes divided by the ``n-1`` adjacent pairs in the window. High for noisy or high-frequency content, low for a slow or DC-dominated window, so it acts as a coarse, amplitude-free pitch / brightness proxy and the classic voiced/unvoiced cue in speech. Exact zeros are not counted as crossings (a sign change needs ``x[i]*x[i+1] < 0``). This is the same quantity the [zero-crossing](../../basics/zero-crossing.qmd) chapter uses to estimate frequency (there as a count; here normalised to a rate so windows of different length compare). To recover that chapter's frequency estimate, multiply by ``(n-1)`` to get the count, then by ``fs / (2 * (n-1))``. """ x = np.asarray(x, dtype=float) crossings = (x[..., :-1] * x[..., 1:]) < 0.0 return np.mean(crossings, axis=axis) def frame_features(x, axis=-1): """Bundle every per-window feature into a dict. Returns ``mean``, ``variance``, ``rms``, ``energy``, ``peak_to_peak``, ``skewness``, ``kurtosis`` and ``zero_crossing_rate``, each a scalar for one window or a length-``n_windows`` array for a stack. ``std`` is omitted as the redundant ``sqrt(variance)``. """ return { "mean": mean(x, axis=axis), "variance": variance(x, axis=axis), "rms": rms(x, axis=axis), "energy": energy(x, axis=axis), "peak_to_peak": peak_to_peak(x, axis=axis), "skewness": skewness(x, axis=axis), "kurtosis": kurtosis(x, axis=axis), "zero_crossing_rate": zero_crossing_rate(x, axis=axis), }