Clarify velocity input requirement for smoothness module and improve jerk computation

Author: sanketsabharwal

examples/smoothness.py

@@ -22,12 +22,23 @@ def extract_wrist_xy(results, keypoint_idx, width, height):
     return None
 
 def main():
+    """
+    Example demonstrating smoothness analysis of hand movement.
+
+    This script:
+    1. Tracks hand position using MediaPipe
+    2. Computes velocity from position changes (dx, dy)
+    3. Feeds velocity to the Smoothness analyzer
+    4. Outputs SPARC and Jerk RMS metrics in real-time
+
+    Note: The Smoothness module expects velocity signal.
+    """
     mp_pose = mp.solutions.pose
     pose = mp_pose.Pose(min_detection_confidence=0.5, min_tracking_confidence=0.5)
 
     # Enable filter to stabilize signal
     smoother = Smoothness(rate_hz=30, use_filter=True)
-    sliding_window = SlidingWindow(50, 1)
+    sliding_window = SlidingWindow(50, 1)  # Store velocity values
     cap = cv2.VideoCapture(0)
 
     LEFT_WRIST_IDX = 15
@@ -53,16 +64,23 @@ def main():
                 if xy and prev_xy and prev_time:
                     dt = now - prev_time
                     if dt > 0:
+                        # Compute position changes
                         dx = xy[0] - prev_xy[0]
                         dy = xy[1] - prev_xy[1]
-                        velocity = np.sqrt(dx**2 + dy**2) / dt
+
+                        # IMPORTANT: Smoothness module expects VELOCITY as input
+                        # We compute velocity magnitude from position changes
+                        velocity = np.sqrt(dx**2 + dy**2) / dt  # pixels/second
 
                         # Clamp unrealistic velocity spikes (in pixels/sec)
                         if velocity < 1000:
+                            # Feed velocity (not position!) to the smoothness analyzer
                             sliding_window.append([velocity])
-                            sparc, jerk = smoother(sliding_window)
-                            if sparc is not None:
-                                print(f"SPARC: {sparc:.3f}, Jerk RMS: {jerk:.1f}")
+
+                            # Get smoothness metrics (SPARC and Jerk RMS)
+                            result = smoother(sliding_window)
+                            if not np.isnan(result['sparc']):
+                                print(f"SPARC: {result['sparc']:.3f}, Jerk RMS: {result['jerk_rms']:.1f}")
 
                 prev_xy = xy
                 prev_time = now

pyeyesweb/low_level/smoothness.py

@@ -1,8 +1,11 @@
 """Movement smoothness analysis module.
 
-This module provides tools for quantifying the smoothness of movement signals
+This module provides tools for quantifying the smoothness of movement velocity signals
 using multiple metrics including SPARC (Spectral Arc Length) and Jerk RMS.
-Designed for real time analysis of motion capture or sensor data.
+Designed for real-time analysis of motion capture or sensor data.
+
+IMPORTANT: This module expects velocity signals as input, not position signals.
+Velocity should be computed from position data before analysis (for example velocity = sqrt(dx^2 + dy^2) / dt).
 
 Smoothness metrics are important indicators of movement quality in:
 1. Motor control assessment
@@ -20,10 +23,11 @@
 
 
 class Smoothness:
-    """Compute movement smoothness metrics from signal data.
+    """Compute movement smoothness metrics from velocity signal data.
 
     This class analyzes movement smoothness using SPARC (Spectral Arc Length)
-    and Jerk RMS metrics. It can optionally apply Savitzky-Golay filtering
+    and Jerk RMS metrics. This class expects velocity signals as input,
+    not position signals. It can optionally apply Savitzky-Golay filtering
     to reduce noise before analysis.
 
     SPARC implementation is based on Balasubramanian et al. (2015) "On the analysis
@@ -54,16 +58,17 @@ class Smoothness:
     >>> from pyeyesweb.data_models.sliding_window import SlidingWindow
     >>> import numpy as np
     >>>
-    >>> # Generate sample movement data (simulated velocity profile)
+    >>> # Generate sample velocity data
+    >>> # For example, from tracking hand movement: velocity = sqrt(dx^2 + dy^2) / dt
     >>> t = np.linspace(0, 2, 200)
-    >>> movement_data = np.sin(2 * np.pi * t) + 0.1 * np.random.randn(200)
+    >>> velocity_data = np.sin(2 * np.pi * t) + 0.1 * np.random.randn(200)
     >>>
     >>> smooth = Smoothness(rate_hz=100.0, use_filter=True)
     >>> window = SlidingWindow(max_length=200, n_columns=1)
     >>>
-    >>> # Add movement data
-    >>> for value in movement_data:
-    ...     window.append([value])
+    >>> # Add velocity data to the window
+    >>> for velocity in velocity_data:
+    ...     window.append([velocity])
     >>>
     >>> result = smooth(window)
     >>> print(f"SPARC: {result['sparc']:.3f}, Jerk RMS: {result['jerk_rms']:.3f}")
@@ -103,20 +108,21 @@ def _filter_signal(self, signal):
         return apply_savgol_filter(signal, self.rate_hz)
 
     def __call__(self, sliding_window: SlidingWindow):
-        """Compute smoothness metrics from windowed signal data.
+        """Compute smoothness metrics from windowed velocity signal data.
 
         Parameters
         ----------
         sliding_window : SlidingWindow
-            Buffer containing signal data to analyze.
+            Buffer containing velocity signal data to analyze.
+            Input should be velocity values, not position values.
 
         Returns
         -------
         dict
             Dictionary containing:
             - 'sparc': Spectral Arc Length (closer to 0 = smoother).
                       Returns NaN if insufficient data.
-            - 'jerk_rms': RMS of jerk (third derivative).
+            - 'jerk_rms': RMS of jerk (computed from velocity input).
                          Returns NaN if insufficient data.
         """
         if len(sliding_window) < 5:
@@ -132,6 +138,7 @@ def __call__(self, sliding_window: SlidingWindow):
         normalized = normalize_signal(filtered)
 
         sparc = compute_sparc(normalized, self.rate_hz)
-        jerk = compute_jerk_rms(filtered, self.rate_hz)
+        # Note: Since the smoothness module expects velocity input (as shown in examples), we specify signal_type='velocity' to compute proper jerk
+        jerk = compute_jerk_rms(filtered, self.rate_hz, signal_type='velocity')
 
         return {"sparc": sparc, "jerk_rms": jerk}

pyeyesweb/utils/math_utils.py

@@ -133,36 +133,56 @@ def compute_sparc(signal, rate_hz=50.0):
     return -arc
 
 
-def compute_jerk_rms(signal, rate_hz=50.0):
-    """Compute RMS of jerk (third derivative) from a signal.
+def compute_jerk_rms(signal, rate_hz=50.0, signal_type='velocity'):
+    """Compute RMS of jerk (rate of change of acceleration) from a signal.
 
-    Jerk is the rate of change of acceleration. Lower RMS jerk values
-    indicate smoother movement.
+    Jerk is the third derivative of position or the first derivative of acceleration. Lower RMS jerk values indicate smoother movement.
 
     Parameters
     ----------
     signal : ndarray
-        1D movement signal (position or velocity).
+        1D movement signal.
     rate_hz : float, optional
         Sampling rate in Hz (default: 50.0).
+    signal_type : str, optional
+        Type of input signal: 'position' or 'velocity' (default: 'velocity').
+        - 'position': Computes third derivative to get jerk
+        - 'velocity': Computes second derivative to get jerk
 
     Returns
     -------
     float
         Root mean square of jerk.
-        Returns NaN if signal has less than 2 samples.
+        Returns NaN if signal has insufficient samples for the required derivatives.
 
     Notes
     -----
-    This implementation uses first-order finite differences to approximate
-    the derivative. For position signals, this computes jerk directly.
+    Uses numpy.gradient for smooth derivative approximation with central differences
+    where possible, providing better accuracy than forward differences.
     """
-    if len(signal) < 2:
-        return np.nan
     rate_hz = validate_numeric(rate_hz, 'rate_hz', min_val=0.0001)
-    dt = 1.0 / rate_hz
-    jerk = np.diff(signal) / dt
-    return np.sqrt(np.mean(jerk ** 2))
+
+    # Define derivative orders needed for each signal type
+    derivative_orders = {
+        'velocity': 2,  # if signal type is velocity we can get velocity -> acceleration -> jerk
+        'position': 3   # if signal type is position we can get position -> velocity -> acceleration -> jerk
+    }
+
+    if signal_type not in derivative_orders:
+        raise ValueError(f"signal_type must be 'position' or 'velocity', got '{signal_type}'")
+
+    n_derivatives = derivative_orders[signal_type]
+    min_samples = n_derivatives + 1
+
+    if len(signal) < min_samples:
+        return np.nan
+
+    # Apply derivatives using numpy.gradient for better accuracy
+    result = np.asarray(signal)
+    for _ in range(n_derivatives):
+        result = np.gradient(result, 1.0/rate_hz)
+
+    return np.sqrt(np.mean(result ** 2))
 
 
 def normalize_signal(signal):