Version 0.16.8

meson.build

@@ -1,7 +1,7 @@
 project(
     'scikit-digital-health',
     'c',
-    version: '0.16.7',
+    version: '0.16.8',
     license: 'MIT',
     meson_version: '>=1.1',
 )

src/skdh/activity/core.py

@@ -287,20 +287,22 @@ def _setup_plotting(self, save_name):  # pragma: no cover
         self.plot_fname = save_name
 
         self._t60 = arange(0, 24.1, 1 / 60)
-    
+
     def _update_date_results(
         self, results, time, day_n, day_start_idx, day_stop_idx, day_start_hour
     ):
         # add 15 seconds to make sure any rounding effects for the hour don't adversely
         # effect the result of the comparison
         start_dt = self.convert_timestamps(time[day_start_idx])
 
-        window_start_dt = start_dt + Timedelta(15, unit='s')
+        window_start_dt = start_dt + Timedelta(15, unit="s")
         if start_dt.hour < day_start_hour:
-            window_start_dt -= Timedelta(1, unit='day')
+            window_start_dt -= Timedelta(1, unit="day")
 
         results["Date"][day_n] = window_start_dt.strftime("%Y-%m-%d")
-        results["Day Start Timestamp"][day_n] = start_dt.strftime("%Y-%m-%d %H:%M:%S.%f")
+        results["Day Start Timestamp"][day_n] = start_dt.strftime(
+            "%Y-%m-%d %H:%M:%S.%f"
+        )
         results["Day End Timestamp"][day_n] = self.convert_timestamps(
             time[day_stop_idx]
         ).strftime("%Y-%m-%d %H:%M:%S.%f")

src/skdh/base.py

@@ -281,7 +281,7 @@ def convert_timestamps(self, t):
             Datetime object
         """
         # set if we want to return aware time
-        kw = {'unit': 's', 'utc': self.tz_name is not None}
+        kw = {"unit": "s", "utc": self.tz_name is not None}
         dt = to_datetime(t, **kw)
         if self.tz_name is not None:
             dt = dt.tz_convert(self.tz_name)

src/skdh/context/ambulation_classification.py

@@ -111,7 +111,7 @@ def predict(self, time, accel, fs=None, tz_name=None, **kwargs):
             (N, 3) array of acceleration, in units of 'g', collected at 20hz.
         fs : float, optional
             Sampling rate. Default None. If not provided, will be inferred.
-        
+
         Other Parameters
         ----------------
         tz_name : {None, str}, optional
@@ -127,7 +127,15 @@ def predict(self, time, accel, fs=None, tz_name=None, **kwargs):
             Results dictionary for downstream pipeline use including start and stop times for ambulation bouts.
 
         """
-        super().predict(expect_days=False, expect_wear=False, time=time, accel=accel, fs=fs, tz_name=tz_name, **kwargs)
+        super().predict(
+            expect_days=False,
+            expect_wear=False,
+            time=time,
+            accel=accel,
+            fs=fs,
+            tz_name=tz_name,
+            **kwargs,
+        )
         # check that input matches expectations, downsample to 20hz if necessary
         time_ds, accel_ds, fs = self._check_input(time, accel, fs)
 

src/skdh/context/gait_classification.py

@@ -74,7 +74,7 @@ def predict(self, *, time, accel, fs=None, tz_name=None, **kwargs):
         fs : float, optional
             Sampling frequency in Hz of the accelerometer data. If not provided,
             will be computed form the timestamps.
-        
+
         Other Parameters
         ----------------
         tz_name : {None, str}, optional

src/skdh/context/motion_detect.py

@@ -8,7 +8,18 @@
 from warnings import warn
 from importlib import resources
 
-from numpy import round, all, linalg, mean, diff, absolute, max, asarray, concatenate, clip
+from numpy import (
+    round,
+    all,
+    linalg,
+    mean,
+    diff,
+    absolute,
+    max,
+    asarray,
+    concatenate,
+    clip,
+)
 from scipy.signal import butter, sosfiltfilt
 
 from skdh.base import BaseProcess, handle_process_returns
@@ -81,7 +92,7 @@ def predict(self, time, accel, fs=None, tz_name=None, **kwargs):
             (N, 3) array of acceleration, in units of 'g'.
         fs : float, optional
             Sampling rate. Default None. If not provided, will be inferred.
-        
+
         Other Parameters
         ----------------
         tz_name : {None, str}, optional
@@ -97,7 +108,13 @@ def predict(self, time, accel, fs=None, tz_name=None, **kwargs):
 
         """
         super().predict(
-            expect_days=False, expect_wear=False, time=time, accel=accel, fs=fs, tz_name=tz_name, **kwargs
+            expect_days=False,
+            expect_wear=False,
+            time=time,
+            accel=accel,
+            fs=fs,
+            tz_name=tz_name,
+            **kwargs,
         )
         # check input requirements are met
         time, accel, fs = self._check_input(time, accel, fs)
@@ -135,7 +152,11 @@ def predict(self, time, accel, fs=None, tz_name=None, **kwargs):
         time_1s = time[0 :: int(fs)][0 : len(predictions)]
 
         # Group results
-        results = {"movement_detected": predictions, "movement_time": time_1s, "tz_name": tz_name}
+        results = {
+            "movement_detected": predictions,
+            "movement_time": time_1s,
+            "tz_name": tz_name,
+        }
 
         # Get starts and stops of movement bouts
         # Run length encoding of predictions @ 1s

src/skdh/features/__init__.py

@@ -34,6 +34,7 @@
     DominantFrequency
     DominantFrequencyValue
     PowerSpectralSum
+    RangePowerSum
     SpectralFlatness
     SpectralEntropy
     ComplexityInvariantDistance

src/skdh/features/lib/extensions/__init__.py

@@ -13,6 +13,7 @@
     dominant_frequency,
     dominant_frequency_value,
     power_spectral_sum,
+    range_power_sum,
     spectral_entropy,
     spectral_flatness,
 )
@@ -42,6 +43,7 @@
     "dominant_frequency",
     "dominant_frequency_value",
     "power_spectral_sum",
+    "range_power_sum",
     "spectral_entropy",
     "spectral_flatness",
     "signal_entropy",

src/skdh/features/lib/extensions/ffeatures.f95

@@ -416,6 +416,62 @@ subroutine power_spectral_sum_1d(n, x, fs, nfft, low_cut, hi_cut, pss) bind(C, n
 end subroutine
 
 
+! --------------------------------------------------------------------
+! SUBROUTINE  range_power_sum_1d
+!     Compute the sum of the spectral power in the range specified.
+! 
+!     Input
+!     n         : integer(long)
+!     x(n)      : real(double), array to compute signal entropy for
+!     nfft      : integer(long), number of points to use in the FFT computation
+!     fs        : real(double), sampling frequency in Hz
+!     low_cut   : real(double), low frequency cutoff for the range to use
+!     hi_cut    : real(double), high frequency cutoff for the range to use
+!     normalize : integer(int), normalize to the power sum across the whole range
+! 
+!     Output
+!     pss : real(double)
+! --------------------------------------------------------------------
+subroutine range_power_sum_1d(n, x, fs, nfft, low_cut, hi_cut, normalize, rps) bind(C, name="range_power_sum_1d")
+    use, intrinsic :: iso_c_binding
+    use real_fft, only : execute_real_forward
+    implicit none
+    integer(c_long), intent(in) :: n, nfft
+    real(c_double), intent(in) :: x(n), fs, low_cut, hi_cut
+    integer(c_int), intent(in) :: normalize
+    real(c_double), intent(out) :: rps
+    ! local
+    real(c_double), parameter :: log2 = log(2._c_double)
+    integer(c_long) :: i, ihcut, ilcut, ier
+    real(c_double) :: sp_norm(nfft + 1)
+    real(c_double) :: sp_hat(2 * nfft + 2), y(2 * nfft)
+
+    ! find the cutoff indices for the high and low cutoffs
+    ihcut = min(floor(hi_cut / (fs / 2) * (nfft - 1) + 1, c_long), nfft + 1)
+    ilcut = max(ceiling(low_cut / (fs / 2) * (nfft - 1) + 1, c_long), 1_c_long)
+
+    if (ihcut > nfft) then
+        ihcut = nfft
+    end if
+
+    rps = 0._c_double  ! ensure starts at 0
+
+    y = 0._c_double
+    y(:n) = x
+    sp_hat = 0._c_double
+    call execute_real_forward(2 * nfft, y, 1.0_c_double, sp_hat, ier)
+
+    sp_norm = sp_hat(1:2*nfft+2:2)**2 + sp_hat(2:2*nfft+2:2)**2
+    ! sp_norm = sp_norm / sum(sp_norm(ilcut:ihcut)) + 1.d-10
+
+    if (normalize == 1_c_int) then
+        rps = sum(sp_norm(ilcut:ihcut)) / sum(sp_norm(1:nfft))
+    else
+        rps = sum(sp_norm(ilcut:ihcut))
+    end if
+end subroutine
+
+
 ! --------------------------------------------------------------------
 ! SUBROUTINE  spectral_entropy_1d
 !     Compute the spectral entropy of the specified frequency range

src/skdh/features/lib/extensions/frequency.c

@@ -10,6 +10,7 @@
 extern void dominant_freq_1d(long *, double *, double *, long *, double *, double *, double *);
 extern void dominant_freq_value_1d(long *, double *, double *, long *, double *, double *, double *);
 extern void power_spectral_sum_1d(long *, double *, double *, long *, double *, double *, double *);
+extern void range_power_sum_1d(long *, double *, double *, long *, double *, double *, int *, double *);
 extern void spectral_entropy_1d(long *, double *, double *, long *, double *, double *, double *);
 extern void spectral_flatness_1d(long *, double *, double *, long *, double *, double *, double *);
 extern void destroy_plan(void);
@@ -239,6 +240,92 @@ PyObject * power_spectral_sum(PyObject *NPY_UNUSED(self), PyObject *args){
 }
 
 
+PyObject * range_power_sum(PyObject *NPY_UNUSED(self), PyObject *args)
+{
+    PyObject *x_;
+    long padlevel;
+    double fs=0., low_cut=0., hi_cut=12.;
+    int norm;
+    int fail = 0;
+
+    if (!PyArg_ParseTuple(args, "Odlddi:range_power_sum", &x_, &fs, &padlevel, &low_cut, &hi_cut, &norm)) return NULL;
+
+    if (fs <= 0.)
+    {
+        PyErr_SetString(PyExc_ValueError, "Sampling frequency cannot be negative");
+        return NULL;
+    }
+    if (hi_cut <= low_cut)
+    {
+        PyErr_SetString(PyExc_ValueError, "High frequency cutoff must be greater than low frequency cutoff");
+        return NULL;
+    }
+
+    PyArrayObject *data = (PyArrayObject *)PyArray_FromAny(
+        x_, PyArray_DescrFromType(NPY_DOUBLE), 1, 0,
+        NPY_ARRAY_ENSUREARRAY | NPY_ARRAY_CARRAY_RO, NULL
+    );
+    if (!data) return NULL;
+    // catch size 0 inputs
+    if (PyArray_SIZE(data) == 0)
+    {
+        PyErr_SetString(PyExc_ValueError, "Input data size must be larger than 0.");
+        Py_XDECREF(data);
+        return NULL;
+    }
+
+    int ndim = PyArray_NDIM(data);
+
+    npy_intp *ddims = PyArray_DIMS(data);
+    npy_intp *rdims = (npy_intp *)malloc((ndim - 1) * sizeof(npy_intp));
+    if (!rdims)
+    {
+        Py_XDECREF(data); return NULL;
+    }
+    for (int i = 0; i < (ndim - 1); ++i)
+    {
+        rdims[i] = ddims[i];
+    }
+
+    PyArrayObject *res = (PyArrayObject *)PyArray_Empty(ndim-1, rdims, PyArray_DescrFromType(NPY_DOUBLE), 0);
+    free(rdims);
+
+    long nfft = (long)pow(2, ceil(log((double)ddims[ndim-1]) / log(2.)) - 1 + padlevel);
+
+    if (!res) fail = 1;
+    if (!fail)
+    {
+        double *dptr = (double *)PyArray_DATA(data);
+        double *rptr = (double *)PyArray_DATA(res);
+
+        long stride = ddims[ndim-1];
+        int nrepeats = PyArray_SIZE(data) / stride;
+
+        for (int i = 0; i < nrepeats; ++i)
+        {
+            range_power_sum_1d(&stride, dptr, &fs, &nfft, &low_cut, &hi_cut, &norm, rptr);
+            dptr += stride;
+            rptr ++;
+        }
+    }
+
+    if (fail)
+    {
+        Py_XDECREF(data);
+        Py_XDECREF(res);
+        // destroy the FFT plan created in the fortran module
+        destroy_plan();
+        return NULL;
+    }
+
+    Py_XDECREF(data);
+    // destroy the FFT plan created in the fortran module
+    destroy_plan();
+
+    return (PyObject *)res;
+}
+
+
 PyObject * spectral_entropy(PyObject *NPY_UNUSED(self), PyObject *args){
     PyObject *x_;
     long padlevel;
@@ -393,6 +480,7 @@ static struct PyMethodDef methods[] = {
     {"dominant_frequency",   dominant_frequency,   1, NULL},  // last is test__doc__
     {"dominant_frequency_value",   dominant_frequency_value,   1, NULL},
     {"power_spectral_sum",   power_spectral_sum,   1, NULL},
+    {"range_power_sum", range_power_sum, 1, NULL},
     {"spectral_entropy",   spectral_entropy,   1, NULL},
     {"spectral_flatness",   spectral_flatness,   1, NULL},
     {NULL, NULL, 0, NULL}          /* sentinel */