diff --git a/CMakeLists.txt b/CMakeLists.txt index 45295d3..212e3e6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,4 +1,4 @@ -cmake_minimum_required(VERSION 3.10) +cmake_minimum_required(VERSION 3.24) project(sensor_processing_native) # Set C++ standard @@ -11,14 +11,39 @@ if(NOT CMAKE_BUILD_TYPE) endif() # Compiler-specific optimizations +# Note: -ffast-math is intentionally omitted — it breaks FINUFFT (IEEE-754 violations) if(CMAKE_CXX_COMPILER_ID STREQUAL "GNU") - set(CMAKE_CXX_FLAGS_RELEASE "-O3 -march=native -mtune=native -ffast-math") + set(CMAKE_CXX_FLAGS_RELEASE "-O3 -march=native -mtune=native") elseif(CMAKE_CXX_COMPILER_ID STREQUAL "Clang") - set(CMAKE_CXX_FLAGS_RELEASE "-O3 -march=native -mtune=native -ffast-math") + set(CMAKE_CXX_FLAGS_RELEASE "-O3 -march=native -mtune=native") elseif(CMAKE_CXX_COMPILER_ID STREQUAL "MSVC") set(CMAKE_CXX_FLAGS_RELEASE "/O2 /arch:AVX2") endif() +# Apple Clang needs explicit help finding Homebrew's libomp. +# Flags must be semicolon-separated CMake lists, not space-separated strings, +# so that each flag is passed as a separate argument to the compiler. +if(APPLE) + execute_process(COMMAND brew --prefix libomp OUTPUT_VARIABLE LIBOMP_PREFIX OUTPUT_STRIP_TRAILING_WHITESPACE) + set(OpenMP_C_FLAGS "-Xpreprocessor;-fopenmp;-I${LIBOMP_PREFIX}/include") + set(OpenMP_CXX_FLAGS "-Xpreprocessor;-fopenmp;-I${LIBOMP_PREFIX}/include") + set(OpenMP_C_LIB_NAMES "omp") + set(OpenMP_CXX_LIB_NAMES "omp") + set(OpenMP_omp_LIBRARY "${LIBOMP_PREFIX}/lib/libomp.dylib") +endif() + +# Download CPM.cmake +file( + DOWNLOAD + https://github.com/cpm-cmake/CPM.cmake/releases/download/v0.40.8/CPM.cmake + ${CMAKE_CURRENT_BINARY_DIR}/cmake/CPM.cmake + EXPECTED_HASH SHA256=78ba32abdf798bc616bab7c73aac32a17bbd7b06ad9e26a6add69de8f3ae4791 +) +include(${CMAKE_CURRENT_BINARY_DIR}/cmake/CPM.cmake) + +# Add finufft dependency +CPMAddPackage("gh:flatironinstitute/finufft@2.5.0") + # Find required packages find_package(PkgConfig) @@ -27,7 +52,7 @@ include_directories(${CMAKE_CURRENT_SOURCE_DIR}/src) # Add the shared library add_library(sensor_processing SHARED - sensor_processing_native.cpp + src/sensor_processing_native.cpp ) # Set library properties @@ -42,11 +67,15 @@ target_include_directories(sensor_processing PRIVATE ${CMAKE_CURRENT_SOURCE_DIR} ) +# Link finufft and enable its include guard in the source +target_link_libraries(sensor_processing PRIVATE finufft) +target_compile_definitions(sensor_processing PRIVATE USE_FINUFFT) + # Platform-specific settings if(ANDROID) # Android-specific settings find_library(log-lib log) - target_link_libraries(sensor_processing ${log-lib}) + target_link_libraries(sensor_processing PRIVATE ${log-lib}) elseif(APPLE) # macOS/iOS-specific settings @@ -61,7 +90,7 @@ elseif(WIN32) elseif(UNIX) # Linux-specific settings - target_link_libraries(sensor_processing m pthread) + target_link_libraries(sensor_processing PRIVATE m pthread) endif() # Install rules diff --git a/finufft_migration_plan.md b/finufft_migration_plan.md new file mode 100644 index 0000000..23f0527 --- /dev/null +++ b/finufft_migration_plan.md @@ -0,0 +1,154 @@ +# Migration Plan: Self-contained finufft dependency in senpy + +## Context +`senpy/setup.py` hardcodes `/tmp/finufft/` paths and silently builds without NUFFT support if those paths don't exist. Getting a working build requires undocumented manual steps (clone finufft, run cmake with specific flags, build). This caused extensive debugging of the Dockerfile. The fix: make `setup.py` auto-fetch and build finufft if it's not present, so `pip install senpy` just works. + +**Constraint**: finufft must link dynamically (shared `.so`). Static linking causes undefined symbol errors at import time. Use `-DFINUFFT_STATIC_LINKING=OFF`. + +## Critical Files +- [senpy/setup.py](../sensor_preprocessing_cpp/senpy/setup.py) — core rewrite +- [senpy/pyproject.toml](../sensor_preprocessing_cpp/senpy/pyproject.toml) — add cmake to build deps +- [senpy/.gitignore](../sensor_preprocessing_cpp/senpy/.gitignore) — add `_deps/` +- [Dockerfile](Dockerfile) — remove manual finufft build block + +--- + +## Changes + +### 1. `senpy/setup.py` — Add auto-fetch logic + +Pin a version and cache inside the repo at `_deps/finufft/` (gitignored): + +```python +FINUFFT_VERSION = "v2.5.0" +FINUFFT_REPO = "https://github.com/flatironinstitute/finufft.git" + +# Cache location: inside the senpy directory (gitignored) +FINUFFT_BUILD_ROOT = os.path.join(ROOT_DIR, "_deps", f"finufft-{FINUFFT_VERSION}") +``` + +Add a `_ensure_finufft()` helper that: +1. Checks if `_deps/finufft-v2.5.0/build/src/libfinufft.so` already exists → skips build +2. Otherwise: shallow-clones the pinned tag, runs cmake with `-DFINUFFT_STATIC_LINKING=OFF -DBUILD_TESTING=OFF`, runs `make -j` +3. Raises `RuntimeError` if the `.so` isn't present after build (no silent fallback) + +Add a `BuildExtWithFinufft(build_ext)` subclass whose `run()`: +1. Calls `_ensure_finufft()` to get the build root +2. Patches `self.extensions[0]` with the correct include dirs, `-DUSE_FINUFFT`, and link args +3. Calls `super().run()` + +Restructure the `Extension` object to be finufft-free at definition time (just pybind11 + numpy + source); `BuildExtWithFinufft.run()` augments it dynamically. + +The `use_finufft` conditional is eliminated entirely — finufft is always fetched and always used. Remove the old module-level path checks. + +```python +def _run(cmd, cwd=None): + import subprocess + print(f"[senpy] Running: {' '.join(cmd)}") + subprocess.check_call(cmd, cwd=cwd) + +def _ensure_finufft(): + so_path = os.path.join(FINUFFT_BUILD_ROOT, "build", "src", "libfinufft.so") + if os.path.isfile(so_path): + print(f"[senpy] finufft already built at {so_path}") + return FINUFFT_BUILD_ROOT + + print(f"[senpy] Fetching finufft {FINUFFT_VERSION} to {FINUFFT_BUILD_ROOT}") + os.makedirs(FINUFFT_BUILD_ROOT, exist_ok=True) + _run(["git", "clone", "--branch", FINUFFT_VERSION, "--depth", "1", + FINUFFT_REPO, FINUFFT_BUILD_ROOT]) + + cmake_build = os.path.join(FINUFFT_BUILD_ROOT, "build") + os.makedirs(cmake_build, exist_ok=True) + _run(["cmake", "-DCMAKE_BUILD_TYPE=Release", + "-DFINUFFT_STATIC_LINKING=OFF", "-DBUILD_TESTING=OFF", ".."], + cwd=cmake_build) + + import multiprocessing + _run(["make", f"-j{multiprocessing.cpu_count()}"], cwd=cmake_build) + + if not os.path.isfile(so_path): + raise RuntimeError(f"[senpy] finufft build failed: {so_path} not found") + return FINUFFT_BUILD_ROOT + + +class BuildExtWithFinufft(build_ext): + def run(self): + root = _ensure_finufft() + lib_dir = os.path.join(root, "build", "src") + for ext in self.extensions: + ext.include_dirs += [ + os.path.join(root, "include"), + os.path.join(root, "build", "_deps", "fftw3-src", "api"), + os.path.join(root, "build", "_deps", "fftw3-build"), + ] + ext.extra_compile_args += ["-DUSE_FINUFFT"] + ext.extra_link_args += [ + os.path.join(lib_dir, "libfinufft.so"), + f"-Wl,-rpath,{lib_dir}", + ] + super().run() +``` + +`setup()` call: change `cmdclass={'build_ext': build_ext}` → `cmdclass={'build_ext': BuildExtWithFinufft}`. + +### 2. `senpy/pyproject.toml` — Add cmake + +```toml +[build-system] +requires = ["setuptools>=42", "wheel", "pybind11>=2.6", "cmake>=3.18"] +build-backend = "setuptools.build_meta" +``` + +The `cmake` PyPI package installs the cmake binary into pip's isolated build environment, removing the system cmake requirement. + +### 3. `senpy/.gitignore` (create if missing) + +``` +_deps/ +*.egg-info/ +build/ +__pycache__/ +*.so +``` + +### 4. `pisces2/Dockerfile` — Remove manual finufft block + +Remove: +- The entire `RUN mkdir -p /tmp/finufft && ...` block (current lines 36–52) +- `ENV LD_LIBRARY_PATH=/tmp/finufft/build/src:/usr/local/lib` — the rpath is baked into `_core.so` by setup.py; no `LD_LIBRARY_PATH` needed +- `ENV PKG_CONFIG_PATH=...` — no longer needed + +Simplify the senpy install step from the multi-line verification + `pip install -v .` to: + +```dockerfile +RUN git clone https://github.com/Arcascope/sensor_preprocessing_cpp.git /opt/sensor_preprocessing_cpp +RUN pip install /opt/sensor_preprocessing_cpp/senpy +``` + +The `apt-get install git cmake build-essential` line stays (git is needed for the clone inside setup.py; cmake is a belt-and-suspenders fallback even though pyproject.toml now provides it). + +--- + +## Developer Workflow After Migration + +```bash +# Fresh clone — just works +pip install -e sensor_preprocessing_cpp/senpy +# → clones finufft v2.5.0 to senpy/_deps/finufft-v2.5.0/ (~3 min first time) +# → subsequent installs skip the build entirely + +# Force finufft rebuild +rm -rf sensor_preprocessing_cpp/senpy/_deps/ +pip install -e sensor_preprocessing_cpp/senpy + +# Upgrade finufft: change FINUFFT_VERSION in setup.py, rebuild +``` + +## Verification + +1. Delete `_deps/` and `pip install -e senpy/` — confirm output shows `[senpy] Fetching finufft v2.5.0...` followed by a successful build +2. Run again — confirm output shows `[senpy] finufft already built at ...` (skip) +3. `python -c "import senpy; print('ok')"` — no ImportError +4. Docker: `docker build --no-cache -t pisces2 .` — finufft is fetched inside the `pip install` step, no pre-steps needed +5. Confirm `ldd` on the installed `_core.so` shows `libfinufft.so` resolved via the rpath diff --git a/senpy/CMakeLists.txt b/senpy/CMakeLists.txt index 2eaa617..cf3da0c 100644 --- a/senpy/CMakeLists.txt +++ b/senpy/CMakeLists.txt @@ -16,8 +16,23 @@ include_directories(${CPP_DIR}) # Find pybind11 find_package(pybind11 REQUIRED) +# Finufft configuration +set(FINUFFT_INCLUDE_DIR "/tmp/finufft/include") +set(FINUFFT_LIB "/tmp/finufft/build/src/libfinufft.a") +set(FFTW3_LIB "/tmp/finufft/build/_deps/fftw3-build/libfftw3.a") +set(FFTW3F_LIB "/tmp/finufft/build/_deps/fftw3f-build/libfftw3f.a") + +if(EXISTS "${FINUFFT_INCLUDE_DIR}" AND EXISTS "${FINUFFT_LIB}") + message(STATUS "Finufft found at ${FINUFFT_LIB}") + set(USE_FINUFFT ON) + add_definitions(-DUSE_FINUFFT) +else() + message(STATUS "Finufft not found, building without NUFFT support") + set(USE_FINUFFT OFF) +endif() + # Create the Python module with the new name -pybind11_add_module(_core +pybind11_add_module(_core ${CPP_DIR}/sensor_processing_native.cpp ) @@ -27,8 +42,17 @@ set_target_properties(_core PROPERTIES OUTPUT_NAME "_core") # expose include directory to the target target_include_directories(_core PRIVATE ${CPP_DIR}) +# Add finufft include directory if available +if(USE_FINUFFT) + target_include_directories(_core PRIVATE ${FINUFFT_INCLUDE_DIR}) +endif() + # Link against necessary libraries -target_link_libraries(_core PRIVATE) +if(USE_FINUFFT) + target_link_libraries(_core PRIVATE ${FINUFFT_LIB} ${FFTW3_LIB} ${FFTW3F_LIB} m gomp) +else() + target_link_libraries(_core PRIVATE) +endif() # Set RPATH for the module set_target_properties(_core PROPERTIES diff --git a/senpy/_core.cpython-312-x86_64-linux-gnu.so b/senpy/_core.cpython-312-x86_64-linux-gnu.so index 00aca3f..0da01a9 100755 Binary files a/senpy/_core.cpython-312-x86_64-linux-gnu.so and b/senpy/_core.cpython-312-x86_64-linux-gnu.so differ diff --git a/senpy/api.py b/senpy/api.py index beca90b..a256734 100644 --- a/senpy/api.py +++ b/senpy/api.py @@ -10,6 +10,14 @@ # Import the C++ module import senpy._core as _senpy +# Try to import finufft; fall back to a slow numpy DFT if unavailable +try: + import finufft + + _HAS_FINUFFT = True +except ImportError: + _HAS_FINUFFT = False + # Add a simple test function to verify imports work def test_import(): @@ -335,6 +343,147 @@ def resample_accelerometer_microseconds( ) +# ── Cubic-spline resampling (C++ backend) ───────────────────────── + + +def resample_accelerometer_cubic( + timestamps: NDArray[np.float64], + x: NDArray[np.float64], + y: NDArray[np.float64], + z: NDArray[np.float64], + target_fs: float, + ts_unit: str = "s", +) -> AccelerometerData: + """Resample accelerometer data using natural cubic spline interpolation. + + Same interface as ``resample_accelerometer`` but uses C3-continuous + cubic splines instead of piecewise-linear interpolation, giving much + better high-frequency rolloff (sinc^4-like vs sinc^2). + + Args: + timestamps: Sample timestamps. + x, y, z: Acceleration components. + target_fs: Target sampling frequency in Hz. + ts_unit: Timestamp unit — ``'s'``, ``'ms'``, or ``'us'``. + + Returns: + AccelerometerData on a uniform grid at *target_fs*. + """ + if not (len(timestamps) == len(x) == len(y) == len(z)): + raise ValueError("All input arrays must have the same length") + + conversion_scalar = {"s": 1e6, "ms": 1e3, "us": 1.0}.get(ts_unit, 1e6) + timestamps_us = (timestamps * conversion_scalar).astype(np.int64) + + result = _senpy.resample_accelerometer_cubic(timestamps_us, x, y, z, target_fs) + return AccelerometerData( + timestamps_us=result["timestamps"], + x=result["x"], + y=result["y"], + z=result["z"], + ) + + +def resample_accelerometer_cubic_microseconds( + timestamps: NDArray[np.int64], + x: NDArray[np.float64], + y: NDArray[np.float64], + z: NDArray[np.float64], + target_fs: float, +) -> AccelerometerData: + """Cubic-spline resampling with microsecond timestamps.""" + if not (len(timestamps) == len(x) == len(y) == len(z)): + raise ValueError("All input arrays must have the same length") + + result = _senpy.resample_accelerometer_cubic(timestamps, x, y, z, target_fs) + return AccelerometerData( + timestamps_us=result["timestamps"], + x=result["x"], + y=result["y"], + z=result["z"], + ) + + +# ── NUFFT-based spectrogram (finufft / numpy-DFT backend) ──────── + + +def _nufft_type1(x_nonuniform: NDArray[np.float64], + c: NDArray[np.complex128], + n_modes: int) -> NDArray[np.complex128]: + """Type-1 NUFFT: non-uniform points → Fourier coefficients. + + Computes f[k] = sum_j c_j * exp(+i * k * x_j) + for k = -n_modes//2 … n_modes//2-1, where x_j ∈ [-pi, pi). + + Uses finufft if available, otherwise falls back to a direct + (slow, O(NM)) DFT so that tests can run without finufft installed. + """ + if _HAS_FINUFFT: + return finufft.nufft1d1(x_nonuniform, c, n_modes) + else: + ks = np.arange(n_modes) - n_modes // 2 + return np.dot(np.exp(1j * np.outer(ks, x_nonuniform)), c) + + +def compute_spectrogram_nufft( + timestamps: NDArray[np.float64], + signal: NDArray[np.float64], + secperseg: float, + secoverlap: float, + ts_unit: str = "s", + target_fs: Optional[float] = None, +) -> SpectrogramResult: + """Compute a spectrogram directly from non-uniformly sampled data + using the Non-Uniform FFT, bypassing time-domain resampling entirely. + + For each STFT window the algorithm: + + 1. Selects the non-uniform samples falling within the window (fixed duration in seconds). + 2. Applies a Hann window evaluated at the true sample times. + 3. Computes the spectrum via type-1 NUFFT using finufft (via C++ backend). + 4. Interpolates to a common frequency grid (if target_fs is specified). + + Because no resampling interpolation is involved, the result is free + of the spectral artifacts that linear or spline resampling introduce + into the 0–10 Hz band. + + Args: + timestamps: Sample timestamps (same length as *signal*). + signal: 1-D signal values (e.g. jerk magnitude). + secperseg: Window duration in seconds. + secoverlap: Overlap duration in seconds. + ts_unit: Timestamp unit — ``'s'``, ``'ms'``, or ``'us'``. + target_fs: If specified, interpolate output to a fixed frequency grid + with this sampling rate. Frequency spacing will be 1/secperseg. + Fmax will be target_fs/2. + + Returns: + SpectrogramResult with *frequencies* (Hz), *times* (s), and + *Sxx* power spectral density array shaped ``(n_times, n_freqs)``. + """ + if len(timestamps) != len(signal): + raise ValueError("timestamps and signal must have the same length") + + # Convert timestamps to seconds if needed + conversion = {"s": 1.0, "ms": 1e-3, "us": 1e-6}.get(ts_unit, 1.0) + t = timestamps.astype(np.float64) * conversion + + # Call C++ backend with optional target_fs for resampling + target_fs_val = target_fs if target_fs is not None else 0.0 + + try: + result_dict = _senpy.compute_spectrogram_nufft( + t, signal.astype(np.float64), secperseg, secoverlap, target_fs_val + ) + return SpectrogramResult( + frequencies=result_dict["freqs"], + times=result_dict["times"], + Sxx=result_dict["Sxx"] + ) + except RuntimeError as e: + raise RuntimeError(f"C++ NUFFT computation failed: {e}") + + def compute_jerk( timestamps: NDArray[np.float64], x: NDArray[np.float64], @@ -756,6 +905,9 @@ class SamplingRates: "ShortTimeFTResult", "MotionFeatures", "resample_accelerometer", + "resample_accelerometer_cubic", + "resample_accelerometer_cubic_microseconds", + "compute_spectrogram_nufft", "compute_jerk", "compute_magnitude", "compute_spectrogram", diff --git a/senpy/plot_spectrogram_comparison.py b/senpy/plot_spectrogram_comparison.py new file mode 100644 index 0000000..c13682c --- /dev/null +++ b/senpy/plot_spectrogram_comparison.py @@ -0,0 +1,317 @@ +"""Spectrogram comparison: linear vs. cubic vs. NUFFT-direct. + +For each dataset in tests/data/, produces a vertical multi-panel figure with one +row per resampling method (Linear, Cubic spline, NUFFT direct), showing the full +time duration. Each row is a pcolormesh of the spectrogram restricted to 0–15 Hz +(where physiologically relevant sleep-scoring content lives), with its own +colorbar for perfect visual alignment. + +Usage (from the senpy/ directory): + python plot_spectrogram_comparison.py + python plot_spectrogram_comparison.py --duration full --datasets SleepAccel +""" + +import argparse +import sys +from pathlib import Path + +import numpy as np +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import matplotlib.colors as mcolors +from matplotlib.ticker import MultipleLocator + +# Make sure senpy is importable when run from the repo root +sys.path.insert(0, str(Path(__file__).parent.parent)) + +from senpy.api import ( + resample_accelerometer, + resample_accelerometer_cubic, + compute_spectrogram, + compute_spectrogram_nufft, + compute_jerk, +) + +DATA_DIR = Path(__file__).parent / "tests" / "data" +OUT_DIR = Path(__file__).parent.parent # workspace root + + +# ── data loading ───────────────────────────────────────────────── + +def load_csv(path: Path): + """Return (t, x, y, z) in seconds, auto-detecting format.""" + with open(path) as f: + first = f.readline() + if first.strip().upper().startswith("TIMESTAMP"): + data = np.loadtxt(path, delimiter=",", skiprows=1) + else: + data = np.loadtxt(path) + return data[:, 0], data[:, 1], data[:, 2], data[:, 3] + + +# ── spectrogram helpers ─────────────────────────────────────────── + +TARGET_FS = 32.0 +NPERSEG = 512 # ~10 s windows @ 50 Hz +NOVERLAP = 384 # 75% overlap → ~2.5 s stride +FREQ_MAX = 15.0 # Hz (0–15 Hz is the range of interest) + + +def _spectrogram_linear(t, x, y, z): + r = resample_accelerometer(t, x, y, z, TARGET_FS) + # compute_jerk on the uniform grid — timestamps drive the C++ difference + j = compute_jerk(r.timestamps_s, r.x, r.y, r.z) + return compute_spectrogram(j.jerk, TARGET_FS, NPERSEG, NOVERLAP) + + +def _spectrogram_cubic(t, x, y, z): + r = resample_accelerometer_cubic(t, x, y, z, TARGET_FS) + j = compute_jerk(r.timestamps_s, r.x, r.y, r.z) + return compute_spectrogram(j.jerk, TARGET_FS, NPERSEG, NOVERLAP) + + +def _spectrogram_nufft(t, x, y, z): + # Compute jerk directly on the raw non-uniform timestamps — no resampling. + # The C++ jerk computes |Δacc| between consecutive samples; with non-uniform + # timestamps this is the L2 displacement in acceleration space, which is the + # same signal the ML pipeline uses (compute_jerk does not normalise by dt). + # + # np.ascontiguousarray is essential: column slices from a 2-D array + # (e.g. data[:,1]) have a non-unit stride, but the C++ wrapper reads the + # raw pointer assuming contiguous memory. t is safe because compute_jerk + # copies it into a new int64 array internally. + j = compute_jerk( + t, + np.ascontiguousarray(x), + np.ascontiguousarray(y), + np.ascontiguousarray(z), + ) + # jerk timestamps come back as microseconds; convert to seconds for NUFFT + t_jerk_s = j.timestamps_us.astype(np.float64) / 1e6 + return compute_spectrogram_nufft(t_jerk_s, j.jerk, NPERSEG, NOVERLAP) + + +METHODS = { + "Linear\n(baseline)": _spectrogram_linear, + "Cubic spline\n(C++)": _spectrogram_cubic, + "NUFFT direct\n(no resampling)": _spectrogram_nufft, +} + + +# ── colour-scale helper ─────────────────────────────────────────── + +def _percentile_norm(Sxx, plo=2, phi=98): + """Robust log-scale colour normalisation.""" + flat = Sxx[Sxx > 0].ravel() + if len(flat) == 0: + return mcolors.LogNorm(vmin=1e-10, vmax=1.0) + vmin = np.percentile(flat, plo) + vmax = np.percentile(flat, phi) + vmin = max(vmin, vmax * 1e-5) # clamp dynamic range + return mcolors.LogNorm(vmin=vmin, vmax=vmax) + + +# ── figure builder ──────────────────────────────────────────────── + +CMAP = "inferno" + +METHOD_COLORS = { + "Linear\n(baseline)": "#e07040", + "Cubic spline\n(C++)": "#4db8e8", + "NUFFT direct\n(no resampling)": "#7ddf74", +} + + +def make_comparison_figure(name: str, t, x, y, z, duration_s: float): + n_methods = len(METHODS) + + # 3 rows (one per method) × 1 column layout + fig, axes = plt.subplots( + n_methods, 1, + figsize=(12, 3.5 * n_methods), + constrained_layout=False, + ) + fig.patch.set_facecolor("#0d0d0d") + + # Ensure axes is always a list (for single-method case compatibility) + if n_methods == 1: + axes = [axes] + + # Compute spectrograms for all methods + specs = {} + for label, fn in METHODS.items(): + specs[label] = fn(t, x, y, z) + + # ── per-method rows ────────────────────────────────────────── + # Each row gets its own percentile-based colour normalisation and colorbar + # so that the structure in every method is independently legible and rows + # are perfectly aligned in time. + for row_idx, (ax, (label, fn)) in enumerate(zip(axes, METHODS.items())): + spec = specs[label] + fmask = spec.frequencies <= FREQ_MAX + freqs = spec.frequencies[fmask] + times = spec.times / 60.0 # minutes + Sxx = spec.Sxx[:, fmask] + + # Guard against all-zero or near-zero Sxx + Sxx = np.maximum(Sxx, 1e-30) + + # Per-panel colour normalisation — computed from this method's own Sxx + norm = _percentile_norm(Sxx, plo=1, phi=99) + + # pcolormesh expects (n_times+1,) and (n_freqs+1,) edge arrays + if len(times) > 1: + dt = times[1] - times[0] + t_edges = np.append(times - dt / 2, times[-1] + dt / 2) + else: + t_edges = np.array([0.0, duration_s / 60.0]) + + if len(freqs) > 1: + df = freqs[1] - freqs[0] + f_edges = np.append(freqs - df / 2, freqs[-1] + df / 2) + else: + f_edges = np.array([0.0, FREQ_MAX]) + + pcm = ax.pcolormesh( + t_edges, f_edges, Sxx.T, + norm=norm, + cmap=CMAP, + rasterized=True, + shading="flat", + ) + + # Annotate leakage: mean power in 0–10 Hz band (text box) + band_mask = (freqs >= 0.5) & (freqs <= 10.0) + mean_pw = np.mean(Sxx[:, band_mask]) if np.any(band_mask) else 0.0 + + color = METHOD_COLORS[label] + + ax.set_facecolor("#0d0d0d") + for spine in ax.spines.values(): + spine.set_edgecolor("#444444") + + ax.tick_params(colors="#cccccc", labelsize=9) + ax.xaxis.label.set_color("#cccccc") + ax.yaxis.label.set_color("#cccccc") + + ax.set_xlabel("Time (min)", fontsize=10, color="#cccccc") + ax.set_ylabel("Frequency (Hz)", fontsize=10, color="#cccccc") + + ax.set_ylim(0, FREQ_MAX) + ax.yaxis.set_major_locator(MultipleLocator(2)) + ax.yaxis.set_minor_locator(MultipleLocator(1)) + ax.tick_params(which="minor", length=2, color="#444444") + + + # Title with colour accent bar + ax.set_title( + f"{label}", + fontsize=11, + color=color, + fontweight="bold", + pad=8, + ) + + # Inset text: mean 0–10 Hz power + ax.text( + 0.97, 0.97, + f"mean 0–10 Hz\n{mean_pw:.2e}", + transform=ax.transAxes, + ha="right", va="top", + fontsize=7.5, + color="#dddddd", + bbox=dict(boxstyle="round,pad=0.3", facecolor="#1a1a1a", + edgecolor="#555555", alpha=0.85), + ) + + # Horizontal reference lines for physio bands + ax.axhline(9/60, color="#4488ff", lw=0.7, alpha=0.5, ls="--") # BR min + ax.axhline(25/60, color="#4488ff", lw=0.7, alpha=0.5, ls="--") # BR max + ax.axhline(0.5, color="#ff8844", lw=0.7, alpha=0.5, ls="--") # HR min + ax.axhline(2.0, color="#ff8844", lw=0.7, alpha=0.5, ls="--") # HR max + + # ── individual colourbar per row ────────────────────────── + # Add colourbar to the right of each row, ensuring rows stay aligned + cbar_ax = fig.add_axes([0.92, 0.12 + (n_methods - row_idx - 1) * (0.78 / n_methods), + 0.015, (0.78 / n_methods) * 0.95]) + cbar = fig.colorbar(pcm, cax=cbar_ax) + cbar.set_label("PSD", color="#cccccc", fontsize=8) + cbar.ax.tick_params(colors="#cccccc", labelsize=7) + cbar.outline.set_edgecolor("#444444") + + # ── layout adjustment ──────────────────────────────────────── + fig.subplots_adjust(left=0.10, right=0.90, top=0.94, bottom=0.10, + hspace=0.35) + + # ── super-title ────────────────────────────────────────────── + jitter_std_ms = np.std(np.diff(t)) * 1000 + nominal_fs = 1.0 / np.median(np.diff(t)) + fig.suptitle( + f"{name} · {duration_s/60:.1f} min · " + f"nominal {nominal_fs:.1f} Hz · jitter σ = {jitter_std_ms:.1f} ms\n" + f"Magnitude spectrogram — dashed: BR band (blue), HR band (orange)", + color="#e0e0e0", + fontsize=10, + y=0.98, + ) + + return fig + + +# ── main ────────────────────────────────────────────────────────── + +def main(): + parser = argparse.ArgumentParser( + description="Plot resampling method spectrogram comparison") + parser.add_argument( + "--duration", type=str, default="full", + help="Seconds of data to use per dataset, or 'full' for entire dataset (default: full)") + parser.add_argument( + "--datasets", nargs="+", + default=["SleepAccel", "Weaver", "Dreamt"], + help="Dataset stem names to process") + args = parser.parse_args() + + saved = [] + for stem in args.datasets: + csv_path = DATA_DIR / f"{stem}.csv" + if not csv_path.exists(): + print(f" [skip] {csv_path} not found") + continue + + print(f"\nLoading {stem}.csv …", flush=True) + t_all, x_all, y_all, z_all = load_csv(csv_path) + + # Determine duration + if args.duration.lower() == "full": + t, x, y, z = t_all, x_all, y_all, z_all + else: + duration_s = float(args.duration) + mask = t_all - t_all[0] < duration_s + t, x, y, z = t_all[mask], x_all[mask], y_all[mask], z_all[mask] + + if len(t) < 500: + print(f" [skip] only {len(t)} samples after trimming") + continue + + actual_dur = t[-1] - t[0] + print(f" {len(t)} samples, {actual_dur:.1f} s, " + f"jitter σ = {np.std(np.diff(t))*1000:.2f} ms") + + print(" Computing spectrograms …", flush=True) + fig = make_comparison_figure(stem, t, x, y, z, actual_dur) + + out_path = OUT_DIR / f"spectrogram_comparison_{stem}.png" + fig.savefig(out_path, dpi=150, bbox_inches="tight", + facecolor=fig.get_facecolor()) + plt.close(fig) + print(f" Saved → {out_path}") + saved.append(out_path) + + print(f"\nDone. {len(saved)} figure(s) written.") + return saved + + +if __name__ == "__main__": + main() diff --git a/senpy/pyproject.toml b/senpy/pyproject.toml index 9eac6e1..91d56d9 100644 --- a/senpy/pyproject.toml +++ b/senpy/pyproject.toml @@ -1,3 +1,3 @@ [build-system] -requires = ["setuptools>=42", "wheel", "pybind11>=2.6"] -build-backend = "setuptools.build_meta" \ No newline at end of file +requires = ["setuptools>=42", "wheel", "pybind11>=2.6", "numpy>=1.19.0"] +build-backend = "setuptools.build_meta" diff --git a/senpy/setup.py b/senpy/setup.py index d949e1b..127626d 100644 --- a/senpy/setup.py +++ b/senpy/setup.py @@ -1,9 +1,12 @@ from setuptools import setup, Extension, find_packages from setuptools.command.build_ext import build_ext +import subprocess +import shutil import sys import setuptools import os import glob +import platform class get_pybind_include(object): """Helper class to determine the pybind11 include path (imported lazily)""" @@ -16,7 +19,7 @@ def __str__(self): CPP_DIR = os.path.abspath(os.path.join(ROOT_DIR, '..', 'src')) # explicit binding file in this directory -PROCESSING_SRC = os.path.join(CPP_DIR, 'sensor_processing_native.cpp') +PROCESSING_SRC = os.path.join('..', 'src', 'sensor_processing_native.cpp') # ensure we compile only the processing TU which now contains the pybind11 module sources = [PROCESSING_SRC] @@ -31,6 +34,69 @@ def __str__(self): # numpy may be installed in the isolated build env; it's okay if not available here pass +# Finufft built by the top-level CMake into build/_deps/ +BUILD_DIR = os.path.abspath(os.path.join(ROOT_DIR, '..', 'build')) +finufft_include_dir = os.path.join(BUILD_DIR, "_deps/finufft-src/include") +finufft_lib_dir = os.path.join(BUILD_DIR, "_deps/finufft-build/src") +_lib_ext = "dylib" if platform.system() == "Darwin" else "so" +finufft_lib_so = os.path.join(finufft_lib_dir, f"libfinufft.{_lib_ext}") +fftw3_lib = os.path.join(BUILD_DIR, "_deps/fftw3-build/libfftw3.a") +fftw3f_lib = os.path.join(BUILD_DIR, "_deps/fftw3f-build/libfftw3f.a") + +def required_artifacts(): + return [finufft_include_dir, finufft_lib_so, fftw3_lib, fftw3f_lib] + +def artifacts_ready(): + return all(os.path.exists(path) for path in required_artifacts()) + +def ensure_native_artifacts(): + if artifacts_ready(): + print("senpy: reusing existing FINUFFT build artifacts", flush=True) + return + + repo_root = os.path.abspath(os.path.join(ROOT_DIR, '..')) + print(f"senpy: configuring native build in {BUILD_DIR}", flush=True) + subprocess.run( + [ + "cmake", + "-S", repo_root, + "-B", BUILD_DIR, + "-DCMAKE_BUILD_TYPE=Release", + "-DFINUFFT_STATIC_LINKING=OFF", + ], + check=True, + ) + print("senpy: building native artifacts with CMake", flush=True) + subprocess.run( + ["cmake", "--build", BUILD_DIR, "--parallel"], + check=True, + ) + print("senpy: native artifact build complete", flush=True) + + missing = [path for path in required_artifacts() if not os.path.exists(path)] + if missing: + raise RuntimeError( + "Required finufft build artifact not found after CMake build: " + + ", ".join(missing) + ) + +class SenpyBuildExt(build_ext): + def run(self): + ensure_native_artifacts() + super().run() + self._copy_runtime_libraries() + + def _copy_runtime_libraries(self): + for ext in self.extensions: + ext_path = self.get_ext_fullpath(ext.name) + ext_dir = os.path.dirname(ext_path) + os.makedirs(ext_dir, exist_ok=True) + target_lib = os.path.join(ext_dir, os.path.basename(finufft_lib_so)) + shutil.copy2(finufft_lib_so, target_lib) + print(f"senpy: copied runtime library to {target_lib}", flush=True) + +include_dirs.append(finufft_include_dir) + ext_modules = [ Extension( 'senpy._core', # Full module path - creates senpy/_core.so @@ -38,6 +104,13 @@ def __str__(self): include_dirs=include_dirs, language='c++', extra_compile_args=['-std=c++17', '-O3', '-march=native', '-DPYTHON'], + extra_link_args=[ + f'-L{finufft_lib_dir}', + '-lfinufft', + fftw3_lib, fftw3f_lib, + '-Wl,-rpath,@loader_path' if platform.system() == 'Darwin' else '-Wl,-rpath,$ORIGIN', + '-lomp' if platform.system() == 'Darwin' else '-lgomp', + ], ), ] @@ -50,10 +123,12 @@ def __str__(self): ext_modules=ext_modules, packages=['senpy'], # Define senpy as a package package_dir={'senpy': '.'}, # The senpy package is in the current directory + package_data={'senpy': ['*.so']}, + include_package_data=True, py_modules=['senpy', 'senpy.api', "senpy._core"], install_requires=['pybind11>=2.6.0', 'numpy>=1.19.0'], setup_requires=['pybind11>=2.6.0'], - cmdclass={'build_ext': build_ext}, + cmdclass={'build_ext': SenpyBuildExt}, zip_safe=False, python_requires='>=3.11', -) \ No newline at end of file +) diff --git a/senpy/tests/__init__.py b/senpy/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/senpy/tests/conftest.py b/senpy/tests/conftest.py new file mode 100644 index 0000000..10b47ca --- /dev/null +++ b/senpy/tests/conftest.py @@ -0,0 +1,40 @@ +"""Shared fixtures for resampling artifact tests.""" + +import os +import pytest +import numpy as np +from pathlib import Path + +DATA_DIR = Path(__file__).parent / "data" + + +# ── helpers ──────────────────────────────────────────────────────── + + +def load_csv(name: str): + """Load one of the test CSVs, returning (timestamps_s, x, y, z). + + Handles the three formats found in the data directory: + - SleepAccel: space-delimited, no header, columns = t x y z (seconds) + - Weaver/Dreamt: comma-delimited with TIMESTAMP,ACC_X,ACC_Y,ACC_Z header + """ + path = DATA_DIR / name + # Sniff first line + with open(path) as f: + first = f.readline() + + if first.strip().startswith("TIMESTAMP"): + data = np.loadtxt(path, delimiter=",", skiprows=1) + else: + data = np.loadtxt(path) + + t, x, y, z = data[:, 0], data[:, 1], data[:, 2], data[:, 3] + return t, x, y, z + + +@pytest.fixture(params=["SleepAccel.csv", "Weaver.csv", "Dreamt.csv"]) +def real_accel(request): + """Parametrised fixture yielding (name, t, x, y, z) for each test CSV.""" + name = request.param + t, x, y, z = load_csv(name) + return name, t, x, y, z diff --git a/senpy/tests/test_resampling_artifacts.py b/senpy/tests/test_resampling_artifacts.py new file mode 100644 index 0000000..e89906f --- /dev/null +++ b/senpy/tests/test_resampling_artifacts.py @@ -0,0 +1,444 @@ +"""Tests for resampling artifact characterisation. + +Compares three approaches: + - **linear**: existing piecewise-linear resampling (baseline) + - **cubic**: natural cubic-spline resampling (C++ backend) + - **nufft_spec**: NUFFT-based spectrogram computed directly from + non-uniform samples (bypasses resampling entirely) + +Each test computes a concrete numerical metric — spurious spectral +power in the 0–10 Hz band, SNR, reconstruction RMSE, etc. — and +prints the results so you get a side-by-side comparison in the +pytest output. + +Run with: + pytest tests/test_resampling_artifacts.py -v -s +""" + +import numpy as np +import pytest +from numpy.typing import NDArray +from senpy.api import ( + resample_accelerometer, + resample_accelerometer_cubic, + resample_accelerometer_nufft, + compute_spectrogram, + compute_spectrogram_nufft, + compute_magnitude, + AccelerometerData, +) + +TARGET_FS = 32.0 # Hz — standard accelerometer target rate + +# ── resampling dispatcher ───────────────────────────────────────── + +RESAMPLE_METHODS = { + "linear": resample_accelerometer, + "cubic": resample_accelerometer_cubic, +} + + +def _resample(method: str, t, x, y, z, fs, ts_unit="s") -> AccelerometerData: + return RESAMPLE_METHODS[method](t, x, y, z, fs, ts_unit=ts_unit) + + +# ── spectral helpers ────────────────────────────────────────────── + + +def band_power(signal: NDArray, fs: float, fmin: float, fmax: float) -> float: + """Integrated PSD in [fmin, fmax] Hz via periodogram.""" + n = len(signal) + freqs = np.fft.rfftfreq(n, 1.0 / fs) + spectrum = np.abs(np.fft.rfft(signal)) ** 2 / n + mask = (freqs >= fmin) & (freqs <= fmax) + return float(np.sum(spectrum[mask])) + + +def spectrogram_band_power(spec, fmin: float, fmax: float) -> float: + """Mean power in [fmin, fmax] across all time windows.""" + mask = (spec.frequencies >= fmin) & (spec.frequencies <= fmax) + return float(np.mean(spec.Sxx[:, mask])) + + +# ══════════════════════════════════════════════════════════════════ +# 1. SYNTHETIC TONE LEAKAGE (resampling comparison) +# ══════════════════════════════════════════════════════════════════ + + +class TestSyntheticToneLeakage: + """Inject a pure tone above the 0–10 Hz band, resample from + jittered timestamps, and measure how much power leaks into 0–10 Hz. + """ + + @staticmethod + def _make_jittered_tone(freq_hz=15.0, duration_s=60.0, + nominal_fs=55.0, jitter_std_s=0.002, + seed=42): + rng = np.random.default_rng(seed) + n = int(duration_s * nominal_fs) + dt = 1.0 / nominal_fs + t = np.arange(n) * dt + rng.normal(0, jitter_std_s, n) + t.sort() + x = np.sin(2 * np.pi * freq_hz * t) + y = np.zeros_like(x) + z = np.zeros_like(x) + return t, x, y, z + + @pytest.mark.parametrize("method", ["linear", "cubic"]) + def test_tone_leakage_into_0_10hz(self, method): + """Power in 0–10 Hz should be small relative to the tone.""" + t, x, y, z = self._make_jittered_tone() + resampled = _resample(method, t, x, y, z, TARGET_FS) + + tone_power = band_power(resampled.x, TARGET_FS, 14.0, 16.0) + leak_power = band_power(resampled.x, TARGET_FS, 0.5, 10.0) + ratio = leak_power / (tone_power + 1e-30) + + print(f"\n [{method}] tone_leakage_ratio = {ratio:.6e}" + f" (leak={leak_power:.4e}, tone={tone_power:.4e})") + assert ratio < 1.0, f"{method}: leakage ratio {ratio} >= 1.0" + + @pytest.mark.parametrize("method", ["linear", "cubic"]) + @pytest.mark.parametrize("jitter_ms", [0.5, 1.0, 2.0, 5.0]) + def test_jitter_sensitivity_sweep(self, method, jitter_ms): + """Leakage as a function of jitter magnitude.""" + t, x, y, z = self._make_jittered_tone( + jitter_std_s=jitter_ms / 1000.0) + resampled = _resample(method, t, x, y, z, TARGET_FS) + + leak = band_power(resampled.x, TARGET_FS, 0.5, 10.0) + tone = band_power(resampled.x, TARGET_FS, 14.0, 16.0) + ratio = leak / (tone + 1e-30) + print(f"\n [{method}, jitter={jitter_ms}ms] ratio={ratio:.6e}") + assert np.isfinite(ratio) + + +# ══════════════════════════════════════════════════════════════════ +# 2. WHITE-NOISE FLOOR +# ══════════════════════════════════════════════════════════════════ + + +class TestWhiteNoiseFloor: + + @staticmethod + def _make_jittered_white_noise(duration_s=60.0, nominal_fs=55.0, + jitter_std_s=0.002, seed=99): + rng = np.random.default_rng(seed) + n = int(duration_s * nominal_fs) + t = np.arange(n) / nominal_fs + rng.normal(0, jitter_std_s, n) + t.sort() + x = rng.standard_normal(n) + y = rng.standard_normal(n) + z = rng.standard_normal(n) + return t, x, y, z + + @pytest.mark.parametrize("method", ["linear", "cubic"]) + def test_flat_spectrum_preserved(self, method): + """Low and high bands should have roughly equal power.""" + t, x, y, z = self._make_jittered_white_noise() + resampled = _resample(method, t, x, y, z, TARGET_FS) + + p_low = band_power(resampled.x, TARGET_FS, 1.0, 10.0) + p_high = band_power(resampled.x, TARGET_FS, 10.0, 20.0) + ratio = p_low / (p_high + 1e-30) + + print(f"\n [{method}] low/high = {ratio:.4f}") + assert 0.3 < ratio < 3.0 + + +# ══════════════════════════════════════════════════════════════════ +# 3. KNOWN-EMPTY-BAND TEST +# ══════════════════════════════════════════════════════════════════ + + +class TestKnownEmptyBand: + """Signal has energy only above 12 Hz.""" + + @staticmethod + def _make_high_freq(duration_s=60.0, nominal_fs=55.0, + jitter_std_s=0.002, seed=77): + rng = np.random.default_rng(seed) + n = int(duration_s * nominal_fs) + t = np.arange(n) / nominal_fs + rng.normal(0, jitter_std_s, n) + t.sort() + x = (np.sin(2 * np.pi * 13 * t) + + 0.8 * np.sin(2 * np.pi * 15 * t) + + 0.6 * np.sin(2 * np.pi * 18 * t)) + y = np.zeros_like(x) + z = np.zeros_like(x) + return t, x, y, z + + @pytest.mark.parametrize("method", ["linear", "cubic"]) + def test_no_energy_below_10hz(self, method): + t, x, y, z = self._make_high_freq() + resampled = _resample(method, t, x, y, z, TARGET_FS) + + p_signal = band_power(resampled.x, TARGET_FS, 12.0, 20.0) + p_artifact = band_power(resampled.x, TARGET_FS, 0.5, 10.0) + ratio = p_artifact / (p_signal + 1e-30) + + print(f"\n [{method}] artifact/signal = {ratio:.6e}") + assert ratio < 1.0 + + +# ══════════════════════════════════════════════════════════════════ +# 4. SINUSOID RECONSTRUCTION ACCURACY +# ══════════════════════════════════════════════════════════════════ + + +class TestSinusoidReconstruction: + + @pytest.mark.parametrize("method", ["linear", "cubic"]) + @pytest.mark.parametrize("freq", [1.0, 5.0, 10.0]) + def test_reconstruction_rmse(self, method, freq): + duration = 10.0 + n_in = int(duration * 55.0) + t_in = np.linspace(0, duration, n_in, endpoint=False) + x_in = np.sin(2 * np.pi * freq * t_in) + y_in = np.zeros_like(x_in) + z_in = np.zeros_like(x_in) + + resampled = _resample(method, t_in, x_in, y_in, z_in, TARGET_FS) + t_out = resampled.timestamps_s + x_true = np.sin(2 * np.pi * freq * t_out) + rmse = np.sqrt(np.mean((resampled.x - x_true) ** 2)) + + print(f"\n [{method}, freq={freq}Hz] RMSE = {rmse:.6e}") + assert rmse < 0.1 + + @pytest.mark.parametrize("method", ["linear", "cubic"]) + def test_cubic_beats_linear(self, method): + """Cubic should always have lower RMSE than linear at 5 Hz.""" + freq = 5.0 + duration = 10.0 + n_in = int(duration * 55.0) + t_in = np.linspace(0, duration, n_in, endpoint=False) + x_in = np.sin(2 * np.pi * freq * t_in) + y_in = z_in = np.zeros_like(x_in) + + resampled = _resample(method, t_in, x_in, y_in, z_in, TARGET_FS) + x_true = np.sin(2 * np.pi * freq * resampled.timestamps_s) + rmse = np.sqrt(np.mean((resampled.x - x_true) ** 2)) + + # Store for cross-method comparison + if not hasattr(TestSinusoidReconstruction, "_rmses"): + TestSinusoidReconstruction._rmses = {} + TestSinusoidReconstruction._rmses[method] = rmse + + if len(TestSinusoidReconstruction._rmses) == 2: + r = TestSinusoidReconstruction._rmses + print(f"\n linear={r['linear']:.6e} cubic={r['cubic']:.6e}" + f" improvement={r['linear']/r['cubic']:.1f}x") + assert r["cubic"] <= r["linear"] + + +# ══════════════════════════════════════════════════════════════════ +# 5. NUFFT SPECTROGRAM vs. RESAMPLE-THEN-SPECTROGRAM +# ══════════════════════════════════════════════════════════════════ + + +class TestNufftSpectrogram: + """Compare the NUFFT spectrogram (no resampling) against the + traditional resample-then-FFT pipeline. + """ + + @staticmethod + def _make_known_spectrum(duration_s=60.0, nominal_fs=55.0, + jitter_std_s=0.002, seed=42): + """Signal with known spectral content: tones at 3 and 7 Hz.""" + rng = np.random.default_rng(seed) + n = int(duration_s * nominal_fs) + t = np.arange(n) / nominal_fs + rng.normal(0, jitter_std_s, n) + t.sort() + signal = (np.sin(2 * np.pi * 3.0 * t) + + 0.5 * np.sin(2 * np.pi * 7.0 * t)) + return t, signal + + def test_nufft_spectrogram_peaks_at_known_freqs(self): + """NUFFT spectrogram should show peaks at 3 Hz and 7 Hz.""" + t, signal = self._make_known_spectrum() + nperseg = 256 + noverlap = 128 + + spec = compute_spectrogram_nufft(t, signal, nperseg, noverlap) + mean_psd = np.mean(spec.Sxx, axis=0) + + # Find the peak frequency within each expected band + mask_3 = (spec.frequencies >= 2.0) & (spec.frequencies <= 4.0) + mask_7 = (spec.frequencies >= 6.0) & (spec.frequencies <= 8.0) + + peak_3 = spec.frequencies[mask_3][np.argmax(mean_psd[mask_3])] + peak_7 = spec.frequencies[mask_7][np.argmax(mean_psd[mask_7])] + + # Also check that the 3 Hz band has more power (signal is 1.0 vs 0.5 amplitude) + power_3 = np.max(mean_psd[mask_3]) + power_7 = np.max(mean_psd[mask_7]) + + print(f"\n NUFFT peaks: {peak_3:.2f} Hz (power={power_3:.4e}), " + f"{peak_7:.2f} Hz (power={power_7:.4e})") + + assert abs(peak_3 - 3.0) < 0.5, f"3 Hz peak at {peak_3}" + assert abs(peak_7 - 7.0) < 0.5, f"7 Hz peak at {peak_7}" + assert power_3 > power_7, "3 Hz should be stronger than 7 Hz" + + def test_nufft_vs_resample_leakage(self): + """For a high-freq tone, NUFFT spectrogram should show less + artifact power in 0–10 Hz than resample-then-spectrogram. + """ + rng = np.random.default_rng(42) + n = int(60 * 55) + t = np.arange(n) / 55.0 + rng.normal(0, 0.002, n) + t.sort() + signal = np.sin(2 * np.pi * 15.0 * t) + nperseg = 256 + noverlap = 128 + + # NUFFT spectrogram: directly from non-uniform samples + spec_nufft = compute_spectrogram_nufft(t, signal, nperseg, noverlap) + nufft_leak = spectrogram_band_power(spec_nufft, 0.5, 10.0) + nufft_tone = spectrogram_band_power(spec_nufft, 14.0, 16.0) + + # Resample-then-spectrogram (linear) + from senpy.api import resample_accelerometer + resampled = resample_accelerometer( + t, signal, np.zeros_like(signal), np.zeros_like(signal), + TARGET_FS) + spec_linear = compute_spectrogram( + resampled.x, TARGET_FS, nperseg, noverlap) + linear_leak = spectrogram_band_power(spec_linear, 0.5, 10.0) + linear_tone = spectrogram_band_power(spec_linear, 14.0, 16.0) + + # Resample-then-spectrogram (cubic) + from senpy.api import resample_accelerometer_cubic + resampled_c = resample_accelerometer_cubic( + t, signal, np.zeros_like(signal), np.zeros_like(signal), + TARGET_FS) + spec_cubic = compute_spectrogram( + resampled_c.x, TARGET_FS, nperseg, noverlap) + cubic_leak = spectrogram_band_power(spec_cubic, 0.5, 10.0) + cubic_tone = spectrogram_band_power(spec_cubic, 14.0, 16.0) + + print(f"\n Spectrogram leakage comparison (0–10 Hz / 14–16 Hz):") + print(f" linear resample: {linear_leak/(linear_tone+1e-30):.6e}") + print(f" cubic resample: {cubic_leak/(cubic_tone+1e-30):.6e}") + print(f" NUFFT direct: {nufft_leak/(nufft_tone+1e-30):.6e}") + + # All should have finite values + assert np.isfinite(nufft_leak) + assert np.isfinite(linear_leak) + + @pytest.mark.parametrize("jitter_ms", [0.5, 1.0, 2.0, 5.0]) + def test_nufft_spectrogram_jitter_robustness(self, jitter_ms): + """NUFFT spectrogram should resolve 3 Hz and 7 Hz tones + even under heavy jitter. + """ + rng = np.random.default_rng(42) + n = int(60 * 55) + t = np.arange(n) / 55.0 + rng.normal(0, jitter_ms / 1000, n) + t.sort() + signal = (np.sin(2 * np.pi * 3.0 * t) + + 0.5 * np.sin(2 * np.pi * 7.0 * t)) + + spec = compute_spectrogram_nufft(t, signal, 256, 128) + mean_psd = np.mean(spec.Sxx, axis=0) + + mask_3 = (spec.frequencies >= 2.0) & (spec.frequencies <= 4.0) + mask_7 = (spec.frequencies >= 6.0) & (spec.frequencies <= 8.0) + peak_3 = spec.frequencies[mask_3][np.argmax(mean_psd[mask_3])] + peak_7 = spec.frequencies[mask_7][np.argmax(mean_psd[mask_7])] + + print(f"\n [jitter={jitter_ms}ms] peaks: 3Hz→{peak_3:.2f}, 7Hz→{peak_7:.2f}") + assert abs(peak_3 - 3.0) < 1.0 + assert abs(peak_7 - 7.0) < 1.0 + + +# ══════════════════════════════════════════════════════════════════ +# 6. REAL DATA: SPECTROGRAM COMPARISON +# ══════════════════════════════════════════════════════════════════ + + +class TestRealDataSpectrogram: + """Compare spectrograms across methods on real accelerometer data. + Uses only the first 60 seconds to keep memory reasonable. + """ + + @pytest.mark.parametrize("method", ["linear", "cubic"]) + def test_spectrogram_smoothness(self, real_accel, method): + name, t, x, y, z = real_accel + + # Use first 60 seconds only + mask = t - t[0] < 60 + t, x, y, z = t[mask], x[mask], y[mask], z[mask] + if len(t) < 500: + pytest.skip(f"{name}: not enough data") + + resampled = _resample(method, t, x, y, z, TARGET_FS) + mag = compute_magnitude(resampled.x, resampled.y, resampled.z) + + nperseg = min(256, len(mag) // 4) + noverlap = nperseg // 2 + if len(mag) < nperseg * 2: + pytest.skip(f"{name}: signal too short") + + spec = compute_spectrogram(mag, TARGET_FS, nperseg, noverlap) + mask_f = (spec.frequencies >= 0.5) & (spec.frequencies <= 10.0) + band = spec.Sxx[:, mask_f] + + if band.shape[0] > 1: + temporal_diff = np.mean(np.abs(np.diff(band, axis=0))) + flicker = temporal_diff / (np.mean(band) + 1e-30) + else: + flicker = 0.0 + + print(f"\n [{method}, {name}] flicker = {flicker:.6f}" + f" mean_power = {np.mean(band):.4e}") + assert np.isfinite(flicker) + + def test_nufft_spectrogram_on_real_data(self, real_accel): + """NUFFT spectrogram on real data should produce finite output.""" + name, t, x, y, z = real_accel + mask = t - t[0] < 60 + t, x, y, z = t[mask], x[mask], y[mask], z[mask] + if len(t) < 500: + pytest.skip(f"{name}: not enough data") + + mag = compute_magnitude(x, y, z) + spec = compute_spectrogram_nufft(t, mag, 256, 128) + + print(f"\n [nufft, {name}] shape={spec.Sxx.shape}" + f" freq_range=[{spec.frequencies[0]:.2f}, {spec.frequencies[-1]:.2f}]" + f" mean_power={np.mean(spec.Sxx):.4e}") + + assert spec.Sxx.shape[0] > 0 + assert spec.Sxx.shape[1] > 0 + assert np.all(np.isfinite(spec.Sxx)) + + +# ══════════════════════════════════════════════════════════════════ +# 7. TIMESTAMP JITTER CHARACTERISATION (real data) +# ══════════════════════════════════════════════════════════════════ + + +class TestTimestampJitter: + + def test_jitter_statistics(self, real_accel): + name, t, x, y, z = real_accel + dt = np.diff(t) + median_dt = np.median(dt) + nominal_fs = 1.0 / median_dt if median_dt > 0 else 0 + + jitter = dt - median_dt + jitter_std = np.std(jitter) + jitter_max = np.max(np.abs(jitter)) + jitter_pct99 = np.percentile(np.abs(jitter), 99) + n_gaps = np.sum(dt > 3 * median_dt) + + print(f"\n [{name}]" + f"\n nominal_fs = {nominal_fs:.2f} Hz" + f"\n median_dt = {median_dt * 1000:.3f} ms" + f"\n jitter_std = {jitter_std * 1000:.3f} ms" + f"\n jitter_max = {jitter_max * 1000:.3f} ms" + f"\n jitter_p99 = {jitter_pct99 * 1000:.3f} ms" + f"\n n_gaps(>3x) = {n_gaps}" + f"\n n_samples = {len(t)}") + assert nominal_fs > 0 diff --git a/src/sensor_processing_native.cpp b/src/sensor_processing_native.cpp index aa0a1f9..738fd0b 100644 --- a/src/sensor_processing_native.cpp +++ b/src/sensor_processing_native.cpp @@ -7,6 +7,8 @@ #include #include +#include + // Logging macros for Android // #define LOG_TAG "SensorProcessing" // #define LOGD(...) __android_log_print(ANDROID_LOG_DEBUG, LOG_TAG, __VA_ARGS__) @@ -289,6 +291,208 @@ class SensorProcessor return result; } + // Compute spectrogram from non-uniformly sampled data using NUFFT + // Bypasses resampling to avoid spectral artifacts in the 0-10 Hz band + // secperseg and secoverlap are in seconds (not samples) + // target_fs: if > 0, resample output to fixed freq grid (fmax = target_fs/2, spacing = 1/secperseg) + static SpectrogramResult computeSpectrogramNUFFT( + const std::vector ×tamps, + const std::vector &signal, + double secperseg, + double secoverlap, + double target_fs = 0.0) + { + if (timestamps.size() != signal.size()) + throw std::runtime_error("timestamps and signal must have the same length"); + + // Estimate median sampling rate + std::vector diffs; + for (size_t i = 1; i < timestamps.size(); ++i) + diffs.push_back(timestamps[i] - timestamps[i-1]); + + std::sort(diffs.begin(), diffs.end()); + double dt_median = diffs[diffs.size() / 2]; + double median_fs = 1.0 / dt_median; + + // Window duration and hop in seconds (input is already in seconds) + double win_dur = secperseg; + double hop_dur = secperseg - secoverlap; + + // Compute nfft based on actual window duration and median sampling rate + int nfft = (int)(secperseg * median_fs); + // Round up to next power of 2 for efficiency + int nfft_padded = nfft; + { + int p = 1; + while (p < nfft) p <<= 1; + nfft_padded = p; + } + + // Frequency axis for positive frequencies + int n_pos_freqs = nfft_padded / 2; + std::vector freqs(n_pos_freqs); + for (int i = 0; i < n_pos_freqs; ++i) + freqs[i] = i / win_dur; + + // Prepare output + std::vector window_centres; + std::vector> spectra; + + double t_start = timestamps[0]; + double t_end = timestamps[timestamps.size() - 1]; + + // Sliding window loop + double win_start = t_start; + while (win_start + win_dur <= t_end + dt_median) + { + double win_end = win_start + win_dur; + + // Select samples in window + std::vector t_win, s_win; + for (size_t i = 0; i < timestamps.size(); ++i) + { + if (timestamps[i] >= win_start && timestamps[i] < win_end) + { + t_win.push_back(timestamps[i]); + s_win.push_back(signal[i]); + } + } + + if (t_win.size() < 4) + { + // Skip windows with too few samples + win_start += hop_dur; + continue; + } + + // Compute tau and Hann window + std::vector tau(t_win.size()); + std::vector hann(t_win.size()); + double s_mean = 0.0; + for (size_t i = 0; i < s_win.size(); ++i) + s_mean += s_win[i]; + s_mean /= s_win.size(); + + // Apply Hann window and detrend + for (size_t i = 0; i < t_win.size(); ++i) + { + tau[i] = (t_win[i] - win_start) / win_dur; + hann[i] = 0.5 * (1.0 - std::cos(2.0 * M_PI * tau[i])); + s_win[i] = (s_win[i] - s_mean) * hann[i]; + } + + // Map to NUFFT domain [-pi, pi) + std::vector x(t_win.size()); + for (size_t i = 0; i < tau.size(); ++i) + x[i] = 2.0 * M_PI * tau[i] - M_PI; + + // Compute NUFFT using finufft + // Allocate arrays for finufft (uses std::complex) + std::vector> c_complex(t_win.size()); + std::vector> fhat_complex(nfft_padded); + + // Fill input array with signal values as complex + for (size_t i = 0; i < t_win.size(); ++i) + c_complex[i] = std::complex(s_win[i], 0.0); + + // Make plan and execute + int ier = 0; + finufft_plan plan = nullptr; + finufft_opts opts; + finufft_default_opts(&opts); + + int64_t nfft_long = nfft_padded; + ier = finufft_makeplan(1, 1, &nfft_long, +1, 1, 1e-14, &plan, &opts); + if (ier != 0) + throw std::runtime_error("finufft_makeplan failed"); + + ier = finufft_setpts(plan, (int64_t)t_win.size(), x.data(), nullptr, nullptr, 0, nullptr, nullptr, nullptr); + if (ier != 0) + throw std::runtime_error("finufft_setpts failed"); + + ier = finufft_execute(plan, c_complex.data(), fhat_complex.data()); + if (ier != 0) + throw std::runtime_error("finufft_execute failed"); + + finufft_destroy(plan); + + // Extract positive frequencies and compute PSD + std::vector psd(n_pos_freqs); + double win_ss = 0.0; + for (size_t i = 0; i < hann.size(); ++i) + win_ss += hann[i] * hann[i]; + + if (win_ss > 0.0) + { + // Match computeSpectrogram: magnitude / sqrt(fs * Σw²) + double scale_factor = 1.0 / std::sqrt(median_fs * win_ss); + // Positive half starts at index nfft_padded/2 + for (int k = 0; k < n_pos_freqs; ++k) + { + int idx = nfft_padded / 2 + k; + double real = fhat_complex[idx].real(); + double imag = fhat_complex[idx].imag(); + psd[k] = std::sqrt(real * real + imag * imag) * scale_factor; + } + } + + spectra.push_back(psd); + window_centres.push_back(win_start + win_dur / 2.0); + + win_start += hop_dur; + } + + // Build result + SpectrogramResult result; + result.freqs = freqs; + + if (spectra.empty()) + { + result.times.clear(); + result.Sxx.clear(); + } + else + { + result.times.resize(window_centres.size()); + for (size_t i = 0; i < window_centres.size(); ++i) + result.times[i] = window_centres[i] - t_start; + result.Sxx = spectra; + } + + // If target_fs is specified, resample to a common frequency grid via cubic spline + if (target_fs > 0.0 && !result.Sxx.empty()) + { + double fmax = target_fs / 2.0; + double freq_spacing = 1.0 / secperseg; + + // Build target frequency grid + std::vector target_freqs; + for (double f = 0; f <= fmax + freq_spacing / 2; f += freq_spacing) + target_freqs.push_back(f); + + // Resample each time window using cubic spline + std::vector> resampled_Sxx(result.Sxx.size()); + + for (size_t t = 0; t < result.Sxx.size(); ++t) + { + const auto &psd_original = result.Sxx[t]; + + // Compute cubic spline coefficients + std::vector freqs_dbl(freqs.begin(), freqs.end()); + auto coeffs = computeNaturalCubicSplineCoeffs(freqs_dbl, psd_original); + + // Evaluate at target frequencies + auto resampled = evaluateCubicSpline(freqs_dbl, psd_original, coeffs, target_freqs); + resampled_Sxx[t] = resampled; + } + + result.freqs = target_freqs; + result.Sxx = resampled_Sxx; + } + + return result; + } + // Compute Short-Time Fourier Transform returning complex values // Returns shape: (n_times, n_frequencies, 2) where last dimension is [real, imag] static std::vector>> computeShortTimeFT( @@ -888,6 +1092,156 @@ class SensorProcessor return {resampledTime, std::make_tuple(rx, ry, rz)}; } + // ── Cubic-spline helpers ────────────────────────────────────────── + // Solve for natural cubic spline coefficients given knot values y + // at (not-necessarily-uniform) knot positions t. + // Returns vectors b, c, d such that on interval [t_i, t_i+1]: + // S(x) = y_i + b_i*(x-t_i) + c_i*(x-t_i)^2 + d_i*(x-t_i)^3 + struct CubicSplineCoeffs + { + std::vector b, c, d; + }; + + static CubicSplineCoeffs computeNaturalCubicSplineCoeffs( + const std::vector &t, + const std::vector &y) + { + size_t n = t.size(); + // n >= 2 guaranteed by caller + + std::vector h(n - 1); + for (size_t i = 0; i < n - 1; ++i) + h[i] = t[i + 1] - t[i]; + + // Solve tridiagonal system for c (natural boundary: c[0]=c[n-1]=0) + std::vector alpha(n, 0.0); + for (size_t i = 1; i < n - 1; ++i) + { + alpha[i] = 3.0 * ((y[i + 1] - y[i]) / h[i] - (y[i] - y[i - 1]) / h[i - 1]); + } + + std::vector c(n, 0.0); + std::vector l(n, 1.0), mu(n, 0.0), z(n, 0.0); + + for (size_t i = 1; i < n - 1; ++i) + { + l[i] = 2.0 * (t[i + 1] - t[i - 1]) - h[i - 1] * mu[i - 1]; + mu[i] = h[i] / l[i]; + z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i]; + } + + // Back-substitution + for (int i = static_cast(n) - 2; i >= 0; --i) + { + c[i] = z[i] - mu[i] * c[i + 1]; + } + + std::vector b(n - 1), d(n - 1); + for (size_t i = 0; i < n - 1; ++i) + { + d[i] = (c[i + 1] - c[i]) / (3.0 * h[i]); + b[i] = (y[i + 1] - y[i]) / h[i] - h[i] * (c[i + 1] + 2.0 * c[i]) / 3.0; + } + + // Trim c to n-1 entries to match b, d + c.resize(n - 1); + + return {b, c, d}; + } + + // Evaluate cubic spline at a set of sorted query points + static std::vector evaluateCubicSpline( + const std::vector &t, // knot positions (sorted) + const std::vector &y, // knot values + const CubicSplineCoeffs &coeffs, + const std::vector &queries) // query positions (sorted) + { + size_t n = t.size(); + size_t m = queries.size(); + const auto &b = coeffs.b; + const auto &c = coeffs.c; + const auto &d = coeffs.d; + + std::vector result(m); + size_t seg = 0; + + for (size_t j = 0; j < m; ++j) + { + double q = queries[j]; + + // Clamp to data range + if (q <= t[0]) + { + result[j] = y[0]; + continue; + } + if (q >= t[n - 1]) + { + result[j] = y[n - 1]; + continue; + } + + // Advance segment pointer (queries are sorted) + while (seg < n - 2 && t[seg + 1] < q) + seg++; + + double dx = q - t[seg]; + result[j] = y[seg] + dx * (b[seg] + dx * (c[seg] + dx * d[seg])); + } + + return result; + } + + // ── Cubic-spline resampling ─────────────────────────────────────── + static std::pair, std::tuple, std::vector, std::vector>> + resampleAccelerometerCubic( + const std::vector ×tamps, + const std::vector &x, + const std::vector &y, + const std::vector &z, + double targetFs) + { + size_t n = timestamps.size(); + if (n < 2) + { + return {std::vector(), std::make_tuple(std::vector(), std::vector(), std::vector())}; + } + + long startTimeUs = timestamps.front(); + long endTimeUs = timestamps.back(); + double durationSec = (endTimeUs - startTimeUs) / 1000000.0; + int numSamples = static_cast(durationSec * targetFs); + if (numSamples < 1) + numSamples = 1; + + // Build uniform output time grid (microseconds) + double intervalUs = 1000000.0 / targetFs; + std::vector resampledTime(numSamples); + std::vector queryT(numSamples); + for (int i = 0; i < numSamples; ++i) + { + resampledTime[i] = startTimeUs + static_cast(i * intervalUs); + queryT[i] = static_cast(resampledTime[i]); + } + + // Convert timestamps to double for spline computation + std::vector tDouble(n); + for (size_t i = 0; i < n; ++i) + tDouble[i] = static_cast(timestamps[i]); + + // Compute spline coefficients for each axis + auto coeffsX = computeNaturalCubicSplineCoeffs(tDouble, x); + auto coeffsY = computeNaturalCubicSplineCoeffs(tDouble, y); + auto coeffsZ = computeNaturalCubicSplineCoeffs(tDouble, z); + + // Evaluate + auto rx = evaluateCubicSpline(tDouble, x, coeffsX, queryT); + auto ry = evaluateCubicSpline(tDouble, y, coeffsY, queryT); + auto rz = evaluateCubicSpline(tDouble, z, coeffsZ, queryT); + + return {resampledTime, std::make_tuple(rx, ry, rz)}; + } + static std::pair, std::vector> computeJerk( const std::vector ×tamps, const std::tuple, std::vector, std::vector> &accelerometerData) @@ -1208,6 +1562,56 @@ py::dict resampleAccelerometer_wrapper( return result_dict; } +py::dict resampleAccelerometerCubic_wrapper( + py::array_t timestamps, + py::array_t x, + py::array_t y, + py::array_t z, + double targetFs) +{ + auto ts_buf = timestamps.request(); + auto x_buf = x.request(); + auto y_buf = y.request(); + auto z_buf = z.request(); + + if (ts_buf.size != x_buf.size || ts_buf.size != y_buf.size || ts_buf.size != z_buf.size) + { + throw std::runtime_error("Input arrays must have the same length"); + } + + std::vector ts_vec(static_cast(ts_buf.ptr), + static_cast(ts_buf.ptr) + ts_buf.size); + std::vector x_vec(static_cast(x_buf.ptr), + static_cast(x_buf.ptr) + x_buf.size); + std::vector y_vec(static_cast(y_buf.ptr), + static_cast(y_buf.ptr) + y_buf.size); + std::vector z_vec(static_cast(z_buf.ptr), + static_cast(z_buf.ptr) + z_buf.size); + + auto result = SensorProcessor::resampleAccelerometerCubic(ts_vec, x_vec, y_vec, z_vec, targetFs); + auto &resampled_time = result.first; + auto &[rx, ry, rz] = result.second; + + py::array_t out_timestamps(resampled_time.size()); + py::array_t out_x(rx.size()); + py::array_t out_y(ry.size()); + py::array_t out_z(rz.size()); + + std::copy(resampled_time.begin(), resampled_time.end(), + static_cast(out_timestamps.request().ptr)); + std::copy(rx.begin(), rx.end(), static_cast(out_x.request().ptr)); + std::copy(ry.begin(), ry.end(), static_cast(out_y.request().ptr)); + std::copy(rz.begin(), rz.end(), static_cast(out_z.request().ptr)); + + py::dict result_dict; + result_dict["timestamps"] = out_timestamps; + result_dict["x"] = out_x; + result_dict["y"] = out_y; + result_dict["z"] = out_z; + + return result_dict; +} + py::dict computeJerk_wrapper( py::array_t timestamps, py::array_t x, @@ -1361,6 +1765,58 @@ py::array_t computeShortTimeFT_wrapper( return output; } +py::dict computeSpectrogramNUFFT_wrapper( + py::array_t timestamps, + py::array_t signal, + double secperseg, + double secoverlap, + double target_fs = 0.0) +{ + auto ts_buf = timestamps.request(); + auto sig_buf = signal.request(); + std::vector timestamps_vec(static_cast(ts_buf.ptr), + static_cast(ts_buf.ptr) + ts_buf.size); + std::vector signal_vec(static_cast(sig_buf.ptr), + static_cast(sig_buf.ptr) + sig_buf.size); + + auto result = SensorProcessor::computeSpectrogramNUFFT(timestamps_vec, signal_vec, secperseg, secoverlap, target_fs); + + // Convert frequencies to NumPy array + py::array_t freqs(result.freqs.size()); + std::copy(result.freqs.begin(), result.freqs.end(), + static_cast(freqs.request().ptr)); + + // Convert times to NumPy array + py::array_t times(result.times.size()); + std::copy(result.times.begin(), result.times.end(), + static_cast(times.request().ptr)); + + // Convert Sxx to 2D NumPy array (times x freqs) + size_t n_times = result.Sxx.size(); + size_t n_freqs = (n_times > 0) ? result.Sxx[0].size() : 0; + py::array_t Sxx({n_times, n_freqs}); + + if (n_times > 0 && n_freqs > 0) + { + auto Sxx_buf = Sxx.request(); + double *Sxx_ptr = static_cast(Sxx_buf.ptr); + for (size_t t = 0; t < n_times; ++t) + { + for (size_t f = 0; f < n_freqs; ++f) + { + Sxx_ptr[t * n_freqs + f] = result.Sxx[t][f]; + } + } + } + + py::dict result_dict; + result_dict["freqs"] = freqs; + result_dict["times"] = times; + result_dict["Sxx"] = Sxx; + + return result_dict; +} + py::dict computeMotionFeatures_wrapper( py::array_t jerkSignal, double fs, @@ -1581,6 +2037,14 @@ PYBIND11_MODULE(_core, m) py::arg("z"), py::arg("targetFs")); + m.def("resample_accelerometer_cubic", &resampleAccelerometerCubic_wrapper, + "Resample accelerometer data using natural cubic spline interpolation", + py::arg("timestamps"), + py::arg("x"), + py::arg("y"), + py::arg("z"), + py::arg("targetFs")); + m.def("compute_jerk", &computeJerk_wrapper, "Compute jerk (derivative of acceleration) from accelerometer data", py::arg("timestamps"), @@ -1601,6 +2065,14 @@ PYBIND11_MODULE(_core, m) py::arg("nperseg"), py::arg("noverlap")); + m.def("compute_spectrogram_nufft", &computeSpectrogramNUFFT_wrapper, + "Compute spectrogram from non-uniformly sampled data using finufft", + py::arg("timestamps"), + py::arg("signal"), + py::arg("secperseg"), + py::arg("secoverlap"), + py::arg("target_fs") = 0.0); + m.def("compute_short_time_ft", &computeShortTimeFT_wrapper, "Compute Short-Time Fourier Transform returning complex values (n_times, n_frequencies, 2)", py::arg("signal"),