Coverage for physiodsp / activity / step_count.py: 88%
311 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-07-14 19:36 +0000
1from __future__ import annotations
3from dataclasses import dataclass, field
4from typing import List, Tuple
6import numpy as np
7import pandas as pd
8from pydantic import BaseModel, Field, PositiveFloat, PositiveInt
9from scipy.signal import butter, sosfiltfilt, welch
11from physiodsp.base import BaseAlgorithm
12from physiodsp.sensors.imu.accelerometer import AccelerometerData
15@dataclass
16class StepBout:
17 """A continuous bout of detected locomotion."""
18 start_s: float
19 end_s: float
20 step_count: int
21 mean_cadence_spm: float
22 activity: str # "walk" | "run"
23 step_timestamps: List[float] = field(default_factory=list)
26class StepCountSettings(BaseModel):
27 """Configuration for the StepCount algorithm.
29 All thresholds are expressed in SI-consistent units (g for acceleration,
30 Hz for frequency, seconds for time).
31 """
33 # --- Stage 1: Pre-processing ---
34 hp_cutoff_hz: PositiveFloat = Field(
35 default=0.5,
36 description="High-pass Butterworth cutoff in Hz for gravity removal.",
37 )
39 # --- Stage 2: Motion gate ---
40 sma_sleep_threshold: PositiveFloat = Field(
41 default=0.02,
42 description="SMA below this value (g) -> window classified as sleep/rest.",
43 )
44 sma_adl_threshold: PositiveFloat = Field(
45 default=0.05,
46 description="SMA above this value (g) is treated as clear locomotion energy.",
47 )
49 # --- Stage 3: Periodicity detector ---
50 periodicity_threshold: PositiveFloat = Field(
51 default=0.50,
52 description="Base autocorrelation peak threshold for locomotion classification.",
53 )
54 spectral_purity_threshold: PositiveFloat = Field(
55 default=8.0,
56 description="PSD peak-to-mean ratio used for secondary confirmation.",
57 )
59 # --- Stage 4: Peak detector ---
60 threshold_factor: PositiveFloat = Field(
61 default=0.30,
62 description="Step peak threshold = threshold_factor * local RMS amplitude.",
63 )
64 analysis_window_s: PositiveFloat = Field(
65 default=2.0,
66 description="Window length in seconds for processing.",
67 )
69 # --- Stage 5: Post-processing ---
70 min_bout_steps: PositiveInt = Field(
71 default=3,
72 description="Minimum steps required for a bout to be retained.",
73 )
74 bout_merge_gap_s: float = Field(
75 default=3.0,
76 ge=0.0,
77 description="Merge consecutive bouts whose inter-bout gap is shorter than this.",
78 )
79 interpolation_max_void_s: PositiveFloat = Field(
80 default=2.0,
81 description="Maximum gap duration (seconds) for missing step interpolation.",
82 )
84 # --- Output / aggregation ---
85 aggregation_window: PositiveInt = Field(
86 default=60,
87 description="Aggregation window length in seconds for aggregate().",
88 )
91class StepCount(BaseAlgorithm):
92 """Wrist-accelerometer step counting algorithm.
94 Five-stage pipeline:
95 1. Pre-processing - VM, high-pass filter
96 2. Motion gate - SMA threshold (sleep / rest rejection)
97 3. Periodicity - autocorrelation + PSD to confirm locomotion
98 4. Peak detection - adaptive-threshold peak detector on BP signal
99 5. Validation - bout merging, short-bout filter, step interpolation
100 """
102 _algorithm_name = "StepCount"
103 _version = "0.1.0"
105 def __init__(self, settings: StepCountSettings = StepCountSettings()) -> None:
106 self.settings = settings
108 def run(self, accelerometer: AccelerometerData) -> "StepCount":
109 """Run the full step counting pipeline on tri-axial accelerometer data."""
110 cfg = self.settings
111 fs = float(accelerometer.fs)
113 # ----------------------------------------------------------------
114 # Stage 1 — Pre-processing
115 # ----------------------------------------------------------------
116 # Vector magnitude (includes gravity)
117 vm = np.sqrt(accelerometer.x**2 + accelerometer.y**2 + accelerometer.z**2)
119 # Remove gravity — high-pass at hp_cutoff_hz
120 vm_dynamic = highpass_filter(vm, fs=fs, cutoff_hz=cfg.hp_cutoff_hz)
122 t = np.asarray(accelerometer.timestamps, dtype=float)
123 total_duration_s = t[-1] - t[0] if len(t) > 1 else 0.0
125 # ----------------------------------------------------------------
126 # Stage 2–4 — Window-by-window processing
127 # ----------------------------------------------------------------
128 win_samples = int(cfg.analysis_window_s * fs)
129 n_windows = len(vm_dynamic) // win_samples
131 all_step_times: List[float] = []
132 all_step_cadences: List[float] = []
133 all_step_activities: List[str] = []
134 window_activity: List[str] = []
135 window_timestamps: List[float] = []
137 # --- Stage 2–3: classify every window (motion gate + periodicity) ---
138 # Detection is deferred to Stage 4 so that contiguous locomotion windows
139 # can be filtered and peak-picked as one continuous signal. Detecting
140 # each 2 s window in isolation drops steps at window boundaries
141 # (filter edge transients + per-window detector reset), which
142 # systematically underestimates the count during sustained walking.
143 window_is_loco: List[bool] = []
144 window_cadence_hz: List[float] = []
145 prev_cadence_hz: float = 0.0
147 for w in range(n_windows):
148 s = w * win_samples
149 e = s + win_samples
150 seg = vm_dynamic[s:e]
151 window_timestamps.append(t[s])
153 # --- Stage 2: SMA gate ---
154 sma = float(np.mean(np.abs(seg)))
156 # --- Stage 3: Periodicity / locomotion classifier ---
157 is_loco, cadence_hz, activity_state = is_locomotion_window(
158 dynamic_vm=seg,
159 fs=fs,
160 sma=sma,
161 sma_sleep_thr=cfg.sma_sleep_threshold,
162 sma_adl_thr=cfg.sma_adl_threshold,
163 periodicity_threshold=cfg.periodicity_threshold,
164 spectral_purity_threshold=cfg.spectral_purity_threshold,
165 )
167 if is_loco:
168 if cadence_hz <= 0.0:
169 cadence_hz = prev_cadence_hz * 0.9
170 if cadence_hz <= 0.0:
171 # No usable cadence estimate: demote to rest.
172 is_loco = False
173 activity_state = "rest"
174 else:
175 if prev_cadence_hz > 0.0:
176 cadence_hz = 0.7 * prev_cadence_hz + 0.3 * cadence_hz
177 prev_cadence_hz = cadence_hz
179 if not is_loco:
180 prev_cadence_hz = 0.0
182 window_activity.append(activity_state)
183 window_is_loco.append(is_loco)
184 window_cadence_hz.append(cadence_hz if is_loco else 0.0)
186 # --- Stage 4: continuous detection over contiguous locomotion runs ---
187 for run_start, run_end, run_cadence_hz in group_locomotion_runs(
188 window_is_loco=window_is_loco,
189 window_cadence_hz=window_cadence_hz,
190 win_samples=win_samples,
191 n_samples=len(vm_dynamic),
192 ):
193 seg = vm_dynamic[run_start:run_end]
194 bp_seg = adaptive_bandpass_filter(
195 seg, center_freq_hz=run_cadence_hz, fs=fs
196 )
197 peak_indices = detect_steps_in_window(
198 bp_signal=bp_seg,
199 fs=fs,
200 cadence_hz=run_cadence_hz,
201 threshold_factor=cfg.threshold_factor,
202 )
204 cadence_spm = run_cadence_hz * 60.0
205 for idx in peak_indices:
206 abs_idx = run_start + idx
207 win_idx = min(abs_idx // win_samples, n_windows - 1)
208 all_step_times.append(float(t[abs_idx]))
209 all_step_cadences.append(cadence_spm)
210 all_step_activities.append(window_activity[win_idx])
212 # ----------------------------------------------------------------
213 # Stage 5 — Post-processing
214 # ----------------------------------------------------------------
215 if len(all_step_times) == 0:
216 self._set_empty_outputs(t, total_duration_s, window_timestamps,
217 window_activity)
218 return self
220 step_times = np.array(all_step_times)
221 step_cadences = np.array(all_step_cadences)
222 step_activities = np.array(all_step_activities)
224 mean_cadence_hz = float(np.mean(step_cadences)) / 60.0
225 step_times_interp = interpolate_missed_steps(
226 step_times_s=step_times.tolist(),
227 cadence_hz=mean_cadence_hz,
228 max_void_s=cfg.interpolation_max_void_s,
229 )
230 step_times = np.array(step_times_interp)
231 step_cadences = np.full(len(step_times), np.mean(step_cadences))
232 step_activities = np.full(len(step_times), "walk", dtype=object)
234 bouts = group_steps_into_bouts(
235 step_times_s=step_times,
236 cadence_per_step=step_cadences,
237 activities=step_activities,
238 max_gap_s=cfg.bout_merge_gap_s if cfg.bout_merge_gap_s > 0 else 1e9,
239 )
241 bouts = merge_close_bouts(bouts, merge_gap_s=cfg.bout_merge_gap_s)
242 bouts = filter_short_bouts(bouts, min_steps=cfg.min_bout_steps)
244 surviving_steps = sorted(
245 t for b in bouts for t in b.step_timestamps
246 )
248 self.bouts: List[StepBout] = bouts
249 self.total_steps: int = len(surviving_steps)
250 self.step_timestamps: np.ndarray = np.array(surviving_steps)
252 t0 = float(t[0])
253 relative_steps = self.step_timestamps - t0
254 bin_ts, cadence_spm = steps_to_cadence_series(
255 step_times_s=relative_steps,
256 total_duration_s=total_duration_s,
257 bin_s=1.0,
258 )
260 activity_per_bin = self._map_activity_to_bins(
261 window_timestamps=window_timestamps,
262 window_activity=window_activity,
263 bin_timestamps=bin_ts + t0,
264 )
266 step_count_per_bin = np.histogram(
267 relative_steps,
268 bins=np.append(bin_ts, bin_ts[-1] + 1.0),
269 )[0]
271 self.biomarker = pd.DataFrame({
272 "timestamps": bin_ts + t0,
273 "step_count": step_count_per_bin,
274 "cadence_spm": cadence_spm,
275 "activity_state": activity_per_bin,
276 })
278 self.timestamps = self.biomarker["timestamps"].values
279 self.values = self.biomarker["cadence_spm"].values
281 return self
283 def aggregate(self, method: str = "mean") -> "StepCount":
284 """Aggregate cadence_spm over aggregation_window-second bins."""
285 super().aggregate(self.timestamps, self.values, method)
286 return self
288 def _set_empty_outputs(
289 self,
290 t: np.ndarray,
291 total_duration_s: float,
292 window_timestamps: List[float],
293 window_activity: List[str],
294 ) -> None:
295 self.bouts = []
296 self.total_steps = 0
297 self.step_timestamps = np.array([])
298 t0 = float(t[0]) if len(t) > 0 else 0.0
299 n_bins = max(1, int(total_duration_s))
300 bin_ts = np.arange(n_bins, dtype=float) + t0
301 activity_per_bin = self._map_activity_to_bins(
302 window_timestamps, window_activity, bin_ts
303 )
304 self.biomarker = pd.DataFrame({
305 "timestamps": bin_ts,
306 "step_count": np.zeros(n_bins, dtype=int),
307 "cadence_spm": np.zeros(n_bins, dtype=float),
308 "activity_state": activity_per_bin,
309 })
310 self.timestamps = self.biomarker["timestamps"].values
311 self.values = self.biomarker["cadence_spm"].values
313 @staticmethod
314 def _map_activity_to_bins(
315 window_timestamps: List[float],
316 window_activity: List[str],
317 bin_timestamps: np.ndarray,
318 ) -> np.ndarray:
319 result = np.full(len(bin_timestamps), "rest", dtype=object)
320 if not window_timestamps:
321 return result
322 wt = np.array(window_timestamps)
323 for i, bt in enumerate(bin_timestamps):
324 idx = int(np.argmin(np.abs(wt - bt)))
325 result[i] = window_activity[idx]
326 return result
329# ---------------------------------------------------------------------------
330# Utilities
331# ---------------------------------------------------------------------------
333def highpass_filter(signal: np.ndarray, fs: float, cutoff_hz: float, order: int = 4) -> np.ndarray:
334 nyq = fs / 2.0
335 sos = butter(order, cutoff_hz / nyq, btype="high", output="sos")
336 return sosfiltfilt(sos, signal)
339def adaptive_bandpass_filter(
340 signal: np.ndarray,
341 center_freq_hz: float,
342 fs: float,
343 order: int = 4,
344 width_factor: float = 0.5,
345) -> np.ndarray:
346 nyq = fs / 2.0
347 low = max(0.5, center_freq_hz * (1.0 - width_factor))
348 high = min(nyq - 0.1, center_freq_hz * (1.0 + width_factor))
349 if low >= high:
350 low, high = 0.5, min(nyq - 0.1, 3.5)
351 sos = butter(order, [low / nyq, high / nyq], btype="band", output="sos")
352 return sosfiltfilt(sos, signal)
355def autocorrelation_peak(
356 signal: np.ndarray,
357 fs: float,
358 freq_min_hz: float = 0.5,
359 freq_max_hz: float = 3.5,
360) -> Tuple[float, float]:
361 sig = signal - signal.mean()
362 if np.std(sig) < 1e-9:
363 return 0.0, 0.0
365 acf = np.correlate(sig, sig, mode="full")
366 acf = acf[len(acf) // 2:]
367 acf = acf / (acf[0] + 1e-12)
369 lag_min = max(1, int(fs / freq_max_hz))
370 lag_max = min(len(acf) - 1, int(fs / freq_min_hz))
372 if lag_min >= lag_max:
373 return 0.0, 0.0
375 sub = acf[lag_min:lag_max + 1]
376 best_rel = int(np.argmax(sub))
377 best_lag = lag_min + best_rel
378 peak_val = float(acf[best_lag])
379 dominant_freq = fs / best_lag if best_lag > 0 else 0.0
380 return peak_val, dominant_freq
383def psd_dominant_freq(
384 signal: np.ndarray,
385 fs: float,
386 freq_min_hz: float = 0.5,
387 freq_max_hz: float = 4.0,
388) -> Tuple[float, float]:
389 nperseg = min(len(signal), 256)
390 freqs, psd = welch(signal, fs=fs, nperseg=nperseg)
392 mask = (freqs >= freq_min_hz) & (freqs <= freq_max_hz)
393 if not np.any(mask):
394 return 0.0, 0.0
396 sub_freqs = freqs[mask]
397 sub_psd = psd[mask]
398 peak_idx = int(np.argmax(sub_psd))
399 dominant_freq = float(sub_freqs[peak_idx])
400 spectral_purity = float(sub_psd[peak_idx] / (np.mean(sub_psd) + 1e-12))
401 return dominant_freq, spectral_purity
404def is_locomotion_window(
405 dynamic_vm: np.ndarray,
406 fs: float,
407 sma: float,
408 sma_sleep_thr: float,
409 sma_adl_thr: float,
410 periodicity_threshold: float,
411 spectral_purity_threshold: float,
412) -> Tuple[bool, float, str]:
413 if sma < sma_sleep_thr:
414 return False, 0.0, "sleep"
416 if sma > sma_adl_thr:
417 eff_threshold = periodicity_threshold * 0.85
418 else:
419 eff_threshold = periodicity_threshold * 1.15
421 ac_peak, dominant_freq = autocorrelation_peak(dynamic_vm, fs=fs)
423 if ac_peak >= eff_threshold and dominant_freq > 0.0:
424 activity = "run" if dominant_freq > 2.2 else "walk"
425 return True, dominant_freq, activity
427 ambiguous_low = eff_threshold * 0.70
428 if ambiguous_low <= ac_peak < eff_threshold:
429 psd_freq, purity = psd_dominant_freq(dynamic_vm, fs=fs)
430 if purity >= spectral_purity_threshold and psd_freq > 0.0:
431 activity = "run" if psd_freq > 2.2 else "walk"
432 return True, psd_freq, activity
434 return False, 0.0, "rest" if sma > sma_sleep_thr else "sleep"
437def group_locomotion_runs(
438 window_is_loco: List[bool],
439 window_cadence_hz: List[float],
440 win_samples: int,
441 n_samples: int,
442) -> List[Tuple[int, int, float]]:
443 """Group consecutive locomotion windows into continuous sample-index runs.
445 Returns a list of ``(start_sample, end_sample, mean_cadence_hz)`` tuples,
446 one per contiguous block of locomotion windows. Detecting steps over these
447 runs (rather than per fixed window) avoids losing steps to filter edge
448 effects and detector resets at every 2 s boundary.
450 The final run is extended to ``n_samples`` so the trailing remainder
451 samples that the integer window count drops are still searched for steps.
452 """
453 runs: List[Tuple[int, int, float]] = []
454 n_windows = len(window_is_loco)
455 i = 0
456 while i < n_windows:
457 if not window_is_loco[i]:
458 i += 1
459 continue
460 j = i
461 cadences: List[float] = []
462 while j < n_windows and window_is_loco[j]:
463 cadences.append(window_cadence_hz[j])
464 j += 1
465 start = i * win_samples
466 end = n_samples if j == n_windows else j * win_samples
467 runs.append((start, end, float(np.mean(cadences))))
468 i = j
469 return runs
472def detect_steps_in_window(
473 bp_signal: np.ndarray,
474 fs: float,
475 cadence_hz: float,
476 threshold_factor: float = 0.30,
477 cadence_tolerance: float = 0.60,
478) -> List[int]:
479 if len(bp_signal) == 0 or cadence_hz <= 0:
480 return []
482 rms_win = max(1, int(fs * 2.0))
483 sq = bp_signal ** 2
484 kernel = np.ones(rms_win) / rms_win
485 rms_smooth = np.sqrt(np.convolve(sq, kernel, mode="same") + 1e-12)
486 threshold = threshold_factor * rms_smooth
488 min_gap = max(1, int((cadence_tolerance / cadence_hz) * fs))
489 max_gap = int((2.0 / (cadence_hz * cadence_tolerance)) * fs)
491 steps: List[int] = []
492 last_step_idx = -max_gap
493 in_peak = False
494 peak_val = -np.inf
495 peak_idx = 0
497 for i, (val, thr) in enumerate(zip(bp_signal, threshold)):
498 if val > thr:
499 if not in_peak:
500 in_peak = True
501 peak_val = val
502 peak_idx = i
503 elif val > peak_val:
504 peak_val = val
505 peak_idx = i
506 else:
507 if in_peak:
508 in_peak = False
509 gap = peak_idx - last_step_idx
510 if not steps:
511 # Always accept the first peak to start the sequence
512 steps.append(peak_idx)
513 last_step_idx = peak_idx
514 elif min_gap <= gap <= max_gap:
515 steps.append(peak_idx)
516 last_step_idx = peak_idx
517 elif gap > max_gap:
518 # Possible missed step(s) - still accept this peak
519 steps.append(peak_idx)
520 last_step_idx = peak_idx
522 if in_peak:
523 gap = peak_idx - last_step_idx
524 if min_gap <= gap:
525 steps.append(peak_idx)
527 return steps
530def group_steps_into_bouts(
531 step_times_s: np.ndarray,
532 cadence_per_step: np.ndarray,
533 activities: np.ndarray,
534 max_gap_s: float = 3.0,
535) -> List[StepBout]:
536 if len(step_times_s) == 0:
537 return []
539 bouts: List[StepBout] = []
540 bout_steps = [step_times_s[0]]
541 bout_cadences = [cadence_per_step[0]]
542 bout_activities = [activities[0]]
544 for i in range(1, len(step_times_s)):
545 gap = step_times_s[i] - step_times_s[i - 1]
546 if gap <= max_gap_s:
547 bout_steps.append(step_times_s[i])
548 bout_cadences.append(cadence_per_step[i])
549 bout_activities.append(activities[i])
550 else:
551 bouts.append(_make_bout(bout_steps, bout_cadences, bout_activities))
552 bout_steps = [step_times_s[i]]
553 bout_cadences = [cadence_per_step[i]]
554 bout_activities = [activities[i]]
556 bouts.append(_make_bout(bout_steps, bout_cadences, bout_activities))
557 return bouts
560def _make_bout(
561 step_times: List[float],
562 cadences: List[float],
563 activities: List[str],
564) -> StepBout:
565 dominant_activity = max(set(activities), key=activities.count)
566 return StepBout(
567 start_s=step_times[0],
568 end_s=step_times[-1],
569 step_count=len(step_times),
570 mean_cadence_spm=float(np.mean(cadences)),
571 activity=dominant_activity,
572 step_timestamps=list(step_times),
573 )
576def filter_short_bouts(bouts: List[StepBout], min_steps: int = 3) -> List[StepBout]:
577 return [b for b in bouts if b.step_count >= min_steps]
580def merge_close_bouts(bouts: List[StepBout], merge_gap_s: float = 3.0) -> List[StepBout]:
581 if not bouts or merge_gap_s <= 0:
582 return bouts
584 merged: List[StepBout] = [bouts[0]]
585 for current in bouts[1:]:
586 prev = merged[-1]
587 gap = current.start_s - prev.end_s
588 if gap <= merge_gap_s:
589 all_steps = prev.step_timestamps + current.step_timestamps
590 all_cadences = [prev.mean_cadence_spm] * prev.step_count + [
591 current.mean_cadence_spm] * current.step_count
592 all_acts = [prev.activity] * prev.step_count + \
593 [current.activity] * current.step_count
594 merged[-1] = _make_bout(all_steps, all_cadences, all_acts)
595 else:
596 merged.append(current)
597 return merged
600def interpolate_missed_steps(
601 step_times_s: List[float],
602 cadence_hz: float,
603 max_void_s: float = 2.0,
604) -> List[float]:
605 if len(step_times_s) < 2 or cadence_hz <= 0:
606 return list(step_times_s)
608 expected_period_s = 1.0 / cadence_hz
609 result = [step_times_s[0]]
611 for i in range(1, len(step_times_s)):
612 gap = step_times_s[i] - step_times_s[i - 1]
613 if gap > expected_period_s * 1.3 and gap <= max_void_s:
614 n_missing = round(gap / expected_period_s) - 1
615 if n_missing >= 1:
616 for k in range(1, n_missing + 1):
617 interp_t = step_times_s[i - 1] + k * expected_period_s
618 result.append(interp_t)
619 result.append(step_times_s[i])
621 return sorted(result)
624def steps_to_cadence_series(
625 step_times_s: np.ndarray,
626 total_duration_s: float,
627 bin_s: float = 1.0,
628) -> Tuple[np.ndarray, np.ndarray]:
629 n_bins = max(1, int(np.ceil(total_duration_s / bin_s)))
630 bin_edges = np.arange(0, n_bins + 1) * bin_s
631 counts, _ = np.histogram(step_times_s, bins=bin_edges)
632 cadence_spm = counts * (60.0 / bin_s)
633 bin_timestamps = bin_edges[:-1]
634 return bin_timestamps, cadence_spm