1 Introduction↩︎

The electromyography (EMG)-based robotic rehabilitation outperforms conventional methods, such as constraint-induced movement and physical therapy. It improves motor recovery, reduces spasticity, and enhances patient engagement by maximizing their voluntary action [1]. Estimating and predicting grip force from real-time EMG signals enables accessing the control variable—force—allowing robotic assistance to supplement patients’ voluntary effort adaptively.

The non-invasive sEMG technique measures muscle activity by detecting electrical signals generated by motor units (MU), which are activated by motor neurons. The EMG signal is obtained as an interference of firing signals from each MU, known as motor unit activation potential (MUAP). Factors such as electrode configuration (placement, size, inter-distance), tissue properties (fat layers, skin conductivity), temporal/spectral firing patterns, and cross-talk from neighboring muscles affect signal variability [2].

In [3], pinching force was predicted using a 6-channel sEMG sleeve with RMS features and the gene expression programming algorithm, achieving RMSE errors of 7.5–8.5% and cross-correlation coefficients up to 95%. The study examined four MVC levels (20%, 40%, 60%, and 80%) but excluded predictions near the transient states between these levels. In [4], sEMG with four bipolar electrodes and finger force signals were used to predict grip force with 90% accuracy during the transient phase, utilizing five EMG and eleven force features. Optimal sensing required 2–4 features across three positions. Another study [5] identified the brachioradialis muscle as optimal for grip force sensing based on EMG signal strength. In [6], EMG signals were acquired using an eight-channel Myo’s armband, and an LSTM network predicted normalized pinching force 1, 3, and 5 seconds ahead directly from the EMG data. Based on [7], which used eight sEMG sensors and 24 healthy subjects, extrinsic muscle coordination reliably predicted grip and pinch force levels due to its greater sensitivity to force changes compared to intrinsic muscles.

Some studies, including [8] and [9], have explored predicting gripping force from the transient phase of EMG signals, which captures the initial burst of muscle activity. In [8], a 192-channel sEMG setup was used with 12 participants to predict grasp force, achieving absolute errors as low as 2.5% of MVC using ten features and regularized linear regression. Meanwhile, [9], utilized 8 sEMG sensors placed on the forearms of 16 participants. The model used ten EMG features to predict gripping force 330 ms ahead with the elastic net regression, resulting in errors of 2% of MVC.

The analysis above, along with findings in [10], show that while the majority of the research has achieved favorable results in terms of accuracy, it relies on a high—and constantly increasing—number of necessary sEMG sensors for the accurate prediction of grasping force. Moreover, forecasts during the transient state are rarely considered. Additionally, existing methods often overlook advanced signal processing for extracting meaningful information from noisy EMG signals, which is a crucial step in causal grip force modeling. While Koopman operator theory (KOT) [11] has demonstrated success in rehabilitation (e.g., Functional Electrical Stimulation [12]) and multi-modal physiological dynamics reconstruction via deep Koopman auto-encoders [13], we present the first online framework for grip force estimation and short-term forecasting using a single sEMG sensor pair. The main contributions of this work are:

• We devised and optimized a composition of EMG signal processing methods, achieving high peak cross-correlations between EMG and grip force signals using a single sensing position on the forearm.

• We devised a novel Koopman-based data-driven approach with problem-specific observables for estimating and short-term predicting grip force in real-time during both transient and plateau phases.

• The methodology enables fast execution, thus facilitating real-time implementation.

2 Materials and methods↩︎

The first part of the study involves collecting time-series data on human⁶ muscle activity through sEMG and measuring hand-grip force with a hand dynamometer. The dynamometer’s shape necessitated a medium wrap type of grasp, following the taxonomy developed by [14] and refined by [15]. To minimize prediction errors, we considered sEMG positioning on the forearm to maximize signal measurement. Advanced signal processing was applied to extract relevant features from the EMG signals and establish a robust correlation with grip force. Offline optimization of decision variables was also conducted.

We use KOT to represent nonlinear dynamics with linear operators, distinguishing between "dynamic" and "static" types [16]. Both act on observables—functions mapping the system’s state to scalar or vector values—rather than directly on the state. The dynamic operator advances observables in time, making it suitable for short-term grip force forecasting, while the static operator maps observables between different spaces, estimating current grip force from EMG signals.

An important step is the signal filtering, essential for proper inferences using EMG signals, which are often affected by noise and cross-talk between muscles. The signals are typically processed using notch filtering at 50 Hz to eliminate ground noise [4], bandpass (BP) filtering between 10 and 500 Hz [5], [7]–[10], and BP filtering that captures only power spectrum peaks between 20 and 60 Hz [6]. In contrast, we focused on extracting meaningful features (comparable to observables in KOT) from the sEMG signal—those highly correlated with grip force—that can be applied within the KOT framework. This strategy aligns with the structural approach to sEMG modeling in [2], which states that recruiting MUs is fundamental for generating muscle force. Greater MU recruitment and higher firing rates amplify force output, with firing rates increasing almost linearly with muscle force. These firing patterns and the interference between the active MUs influence the characteristics of the measured sEMG. We aim to isolate the difficult-to-identify signal components—specifically the recruited MUs, their interference patterns, and firing rates—that predominantly contribute to measured grip force during a medium wrap, to separate meaningful signals from unwanted noise. To achieve this, we utilized Fast Fourier Transformation (FFT) and sensitivity analysis (SA) to develop a spectral mask that selectively targets specific spectral components.

The ultimate goal is to develop a fully integrated module that accurately estimates and predicts exerted grip force from EMG signals using KOT, encompassing advanced signal processing and personalized calibration protocols.

2.1 Experimental Design↩︎

Two complementary wireless sensing devices were utilized: an EMG sensing unit and a grip dynamometer. The EMG sensing unit records electrical activity on the skin surface associated with muscle contractions, while the dynamometer measures gripping force. The EMG signals were sampled at 1 kHz, while the dynamometer data was acquired at 200 Hz. Detailed device specifications and acquisition parameters are provided in Appendix 5.1.1. The reference dynamometer was calibrated as described in Appendix 5.1.2.

Various electrode positions were tested, and two forearm positions (see Fig. 1) were selected for further analysis based on signal strength observed in a simple screening experiment. A two-factor randomized block design (RBD) [17] was chosen, as described in Table ¿tbl:tab:rbd-design?. Electrode positioning followed configurations from prior studies [5], [10]. The experiment involved 13 male participants (ages 22–24), with subject variability controlled through blocking to isolate measurement position effects on estimation and prediction errors. Randomization was implemented within blocks. Skin preparation included hair removal and alcohol cleansing at measurement positions. Participants were given time to familiarize themselves with the required force levels. For all trials, the initial 5 s were used to zero the dynamometer signal at each position.

2.2 Procedural Parameter Optimization and sEMG Procesing Steps↩︎

The optimization procedure maximizes the mean peak cross-correlation between processed EMG signal \(\boldsymbol{e}_\textrm{proc}\) and the synchronously collected grip force \(\boldsymbol{g}\). This metric is selected for its effectiveness in capturing signal similarity across time lags, which vary between subjects and measurements. Algorithm \(\mathop{\mathrm{\mathrm{\small PeakCrossCorr}}}\) (see Algorithm 10 in Appendix 5.2.1) computes the per-run peak cross-correlation. Multi-step SA and optimization identify key spectral components and narrow decision vector bounds near the optimum using scatterplot projections with smoothed trends and means, generated via Latin hypercube (LH) sampling. This integrated approach ensures robust cross-correlations, consistent performance, and stability by avoiding unstable maxima despite small variations.

First, EMG signals are processed in mini-batches of approximately 0.5 seconds (496 data points) with a sampling rate of 992.97 Hz. This batch size enables quasi-stationary analysis while maintaining sensitivity to physiologically relevant variations in grip force. This aligns with [18], who demonstrated that while initial stimulus processing activates the motor system early (100 ms–130 ms), complex decision-making processes are reflected in grip force patterns only after 300 ms–350 ms. The selected batch size optimizes the tradeoff between frequency and temporal resolution. It is sufficient for reliable frequency analysis via Short-Time Fourier Transform (STFT) and capturing meaningful changes in signal dynamics. We transform each batch using FFT (resolution: 2.002 Hz), followed by the application of an optimized spectral mask \(\mathbf{w}_\text{mask}\) to modify frequency components selectively. After inverse FFT, the signal is rectified by taking absolute value and smoothed using an exponential moving average (MA) with optimized window size \(w_\text{EMA}\) and decay factor \(\alpha_\text{EMA}\). The exponential MA window \(\mathbf{wnd}_{\text{EMA}}\) includes previous batch data for the first \(w_\text{EMA}-1\) points to ensure smooth transitions between batches. The batch processing pipeline for optimal EMG signal extraction is illustrated in Fig. 2, with the detailed Algorithm 11: \(\mathop{\mathrm{\mathrm{\small BatchProcEMG}}}\) provided in Appendix 5.2.2.

The high-dimensional optimization problem was formulated to maximize the mean peak cross-correlation \(\boldsymbol{r}_\text{peak,all}\) across all 52 experiment runs using a decision vector comprised of 250 variables. The entries and initial bounds of the decision vector:

• Spectral mask \(\mathbf{w}_\text{mask}\) with 248 entries corresponding to each frequency bin (without DC offset) in range 0–5,

• Exponential MA window size \(w_\text{EMA}\) in range 2–495,

• Exponential MA decay factor \(\alpha_\text{EMA}\) in range 0–0.05.

The optimal decision vector is obtained by maximizing the objective function: \[\begin{align} \underset{\boldsymbol{w}_\textrm{mask}, w_\textrm{EMA}, \alpha_\textrm{EMA}}{\arg\max} \mathop{\mathrm{\mathrm{\small OptimProb}}}\big( \boldsymbol{w}_\textrm{mask}, w_\text{EMA}, \alpha_\text{EMA} \big), \end{align}\] where the \(\mathop{\mathrm{\mathrm{\small OptimProb}}}\) is described in Algorithm 3. The previously described algorithms for EMG batch processing and computing peak cross-correlation are applied to each of the 52 experimental runs, extracting paired EMG and grip force data \((\mathbf{e}_\textrm{raw}, \mathbf{g})\) from the complete dataset \((\mathbf{e}_\textrm{raw,all}, \mathbf{g}_{\text{all}})\).

The DC offset was set to zero and excluded from the optimization process. Initial mask bounds allowed either complete removal of a spectral component’s amplitude or amplification by up to five times. In contrast, the decay factor bounds were designed to support simple MA (if set to 0) and exponential MA (if greater than 0). Broad bounds on the window size also enabled fine-tuning of the applied smoothing.

Due to differences in sampling rates—where the EMG signal is sampled \sim 5 times faster than the dynamometer—the grip force signal \(\boldsymbol{g}\) was resampled using EMG timestamps. Intermediate points were estimated using linear regression to ensure proper signal alignment for cross-correlation computation. We hypothesize that not all 248 FFT frequency bins significantly impact the cross-correlation. To address this, SA is performed in conjunction with optimization to narrow the problem’s scope and identify the most influential spectral components on the optimization problem in Algorithm 3. A preliminary Sobol SA is conducted separately for positions P1 and P2, utilizing grouped spectral mask variables to provide a high-level understanding of how cross-correlation sensitivity is affected by smoothing and filtering. Details regarding the preliminary SA procedure, parameters, and sample size are outlined in Appendix 5.2.3.

Fig. 4 (a) presents the obtained first-order and total-order sensitivity indices (SIs) and their sums. SIs sums are \sim 1, indicating that they can reliably approximate the percentage of output variance attributable to each variable, and interactions between variables can be safely ignored. When sampling the problem within the initial bounds, the contributions of smoothing parameters and the spectral mask to the variance in mean peak cross-correlation depend on the measurement position. Smoothing window size \(w_\text{EMA}\) is the most influential factor, explaining the majority of output variance across both positions. The decay factor contributes \sim 20 %, while the spectral mask contribution varies most: 14 % on P1 to 37 % on P2. Additional LH sampling with 10000 decision vector samples was conducted within the same bounds to investigate the partial contributions of smoothing parameters further. Fig. 4 (b) displays the resulting projections onto the decay factor \(\alpha_\text{EMA}\) and window size \(w_\text{EMA}\). Increasing \(w_\text{EMA}\) improves the peak cross-correlation until it reaches a plateau, supporting the decision to narrow the window size bounds to the 200–495 range for subsequent steps. Conversely, increasing \(\alpha_\text{EMA}\) shows a decreasing trend in peak cross-correlation, so we narrowed the range to 0–0.01. The grouped SA was repeated with the narrowed parameter bounds, where the spectral mask contributed to 88 % of the variance in P1 and 95 % in P2. This highlights the necessity of ungrouping the mask for more detailed analysis in the subsequent steps.

For iterative, multi-step, ungrouped SA, a more efficient RBD-FAST SA method [19] was used. Details regarding the simultaneous SA and optimization procedure steps are outlined in Appendix 5.2.4. The top spectral components with narrowed bounds after each SA and optimization step are presented in Fig. 5. In a step-by-step analysis, a reverse funnel-shaped pattern can be observed (indicated by the two red arrows), with one starting point at low-frequency components and the other around 50 Hz. This pattern highlights the transition from the most influential spectral components to the less influential ones. After the \(19^\text{th}\) step of ungrouped SA, the cumulative contribution of the SIs corresponding to spectral mask for frequencies >= 204 Hz amounted to less than 2 %. Therefore, we reduced the spectral mask \(\boldsymbol{w}_\textrm{mask}\) to 101 components (2 Hz–202 Hz, step \sim 2 Hz), consistent with [20], which shows most sEMG power is below 250 Hz. Last two SA and optimization steps (\(20^\text{th}, 21^\text{st}\)) reveal interesting trends:

• Increasing decay factor \(\alpha_\text{EMA}\) decreases cross-correlation, and optimal values are close to zero (0–5 × 10⁻⁴), indicating a need for simple MA.

• Window size \(w_\text{EMA}\) optimum is in the 275–330 range.

The description of the optimal spectral mask is detailed in Section 3. With the optimal set of signal processing procedures and parameters, next section outlines the Koopman-driven methodology for estimating current and predicting future values of the grip force from the processed EMG signal.

2.3 Koopman-driven Estimation and Prediction of Hand Grip Force↩︎

The framework begins with a one-time 20 s–30 s calibration experiment using sEMG sensors and a dynamometer, with grip force levels as in Table ¿tbl:tab:rbd-design?. The processed EMG signal \(\boldsymbol{e}_{\text{proc}}\) is used to train the estimation model. Data is processed batch-by-batch for real-time application to estimate current grip force \(\boldsymbol{g}\) under the window. This estimation is used to train a prediction model that forecasts grip force over a 0.5 s horizon.

2.3.1 Koopman Operator Theory↩︎

The state-space representation of a dynamical system involves an \(n\)-dimensional manifold \(\mathcal{M}\), where states \(\mathbf{x}\) evolve discretely over time: \[\begin{align} \mathbf{x}[i+1] = \mathbf{F}\big(\mathbf{x}[i]\big), \end{align} \label{eq:dynamicalMap}\tag{1}\] where \(\mathbf{F}\mathbin{:} \mathcal{M} \to \mathcal{M}\) is the potentially nonlinear state transition function, and \(\mathbf{x}[i+1]\) represents the time-shifted state. To address the challenges of modeling complex nonlinear dynamics, we adopt an operator-theoretic perspective on the dynamics of observables [11]. The nonlinear system (1 ) is mapped onto observable dynamics \(\phi(\mathbf{x})\), typically complex-valued, but here we focus on real-valued \(\phi\mathbin{:} \mathcal{M} \to \mathbb{R}\). Collecting all possible observables constitutes a vector space that is generally infinite-dimensional. The Koopman operator \(\mathcal{K}\), describing observable evolution over \(\Delta t\), is defined as [21]: \[\begin{align} \phi\big(\mathbf{x}[i+1]\big) = \mathcal{K}\phi\big(\mathbf{x}[i]\big) = \phi\Big(\mathbf{F}\big(\mathbf{x}[i]\big)\Big). \end{align} \label{eq:dynamical95evolution}\tag{2}\] \(\mathcal{K}\) remains linear even for nonlinear underlying system [11].

2.3.2 Estimating Grip Force↩︎

The Koopman operator can also describe “static” nonlinear maps between different spaces \(\mathcal{M} \to \mathcal{N}\) [16], a property we leverage to estimate current-batch grip force from processed EMG. Through lifting and proper choice of observables, we can describe static nonlinear maps using spaces of observables and a linear mapping operator \(\mathcal{K}_\text{e}\mathbin{:} \mathcal{O}_\mathcal{M} \to \mathcal{O}_\mathcal{N}\). The processed EMG signal \(\boldsymbol{e}_\text{proc}\) and measured grip force \(\boldsymbol{g}\) are lifted using yet-to-be-determined vectors of functions \(\mathbf{\phi}\) and \(\psi\), respectively. The input and output matrices of the lifted variables, where \(\big(\boldsymbol{e}_\text{proc}[i], \boldsymbol{g}[i]\big)\) represents a single data pair realization, can be expressed as: \[\begin{align} E = \phi\big(\boldsymbol{e}_\text{proc}\big), \quad G = \psi\big(\boldsymbol{g}\big). \end{align}\] The approximation of the static Koopman operator is obtained by minimizing the Frobenius norm [16]: \[\begin{align} \min_{\mathcal{K}_\text{e}} \| G - \mathcal{K}_\text{e}E \|_F \quad \to \quad \overline{\mathcal{K}}_\text{e} = GE^{\dagger}, \end{align} \label{eqn:static-koopman}\tag{3}\] where \(^{\dagger}\) is the SVD-based Moore-Penrose pseudoinverse. The estimation Koopman operator \(\overline{\mathcal{K}}_\text{e}\), trained on the complete personalized calibration experiment, is applied sequentially to each \sim 0.5 s batch window to generate real-time grip force estimations. Empirical analysis showed 8-fold downsampling (\(\require{physics} 993 \to \qty{124}{\Hz}\)) of processed EMG (Fig. 9 (b)) and normalization to \([0, 1]\) significantly reduced grip force estimation error. We employed Hankel lifting with time-delay embedding, a technique well-suited for capturing the system’s temporal dynamics [22]. The Hankel data matrix \(E\), constructed from \(\boldsymbol{e}_\text{proc}\) with \(N\) data points, incorporates state \(\mathbf{e}_0\) in the first row and \(d\) time-delayed observables \(\mathbf{e}_{\textrm{td}(i)}\) in subsequent rows: \[\begin{align} {!}{\( E = \begin{bmatrix} \mathbf{e}_0 \\ \mathbf{e}_{\textrm{td}(1)} \\ \vdots \\ \mathbf{e}_{\textrm{td}(d-1)} \\ \mathbf{e}_{\textrm{td}(d)} \end{bmatrix} = \begin{bmatrix} \boldsymbol{e}_\text{proc}[1] & \boldsymbol{e}_\text{proc}[2] & \cdots & \boldsymbol{e}_\text{proc}[N-d] \\ \boldsymbol{e}_\text{proc}[2] & \boldsymbol{e}_\text{proc}[3] & \cdots & \boldsymbol{e}_\text{proc}[N-d+1] \\ \vdots & \vdots & \ddots & \vdots \\ \boldsymbol{e}_\text{proc}[d] & \boldsymbol{e}_\text{proc}[d+1] & \cdots & \boldsymbol{e}_\text{proc}[N-1] \\ \boldsymbol{e}_\text{proc}[d+1] & \boldsymbol{e}_\text{proc}[d+2] & \cdots & \boldsymbol{e}_\text{proc}[N] \end{bmatrix}. \)} \end{align} \label{eqn:hankel-lifting}\tag{4}\] The lifting in (4 ) was applied to the downsampled EMG matrix \(E\) and the grip force matrix \(G\), using \(d=60\) time delays determined empirically. A complementary nonlinear lifting was used to map EMG plateaus and transitions to grip force patterns, producing gridded indicator observables. These observables, derived from time-delay embedded data, discretize the embedding space into a 3D grid to capture temporal dynamics. The grid is defined on a Cartesian plane with three time delays—\(\mathbf{e}_{\textrm{td}(1)}, \mathbf{e}_{\textrm{td}(1+\tau_1)},\, \mathbf{e}_{\textrm{td}(1+\tau_2)}\)—with each subregion treated as a binary observable (Fig. 6). For any subregion \(S_{ijk}\), an indicator function assigns a value of one to points within the subregion and zero otherwise. Let \(S_{ijk}\) be a subregion in \([0, 1]^3\), defined by the lower (LO) and upper (UP) grid bounds \(\boldsymbol{b}_{ijk} = \big[b_{i,\text{LO}}, b_{i,\text{UP}}, b_{j,\text{LO}}, b_{j,\text{UP}}, b_{k,\text{LO}}, b_{k,\text{UP}}\big]\):

• For \(\mathbf{e}_{\textrm{td}(1)}\): lower limit \(b_{i,\text{LO}}\), upper limit \(b_{i,\text{UP}}\).

• For \(\mathbf{e}_{\textrm{td}(1+\tau_1)}\): lower limit \(b_{j,\text{LO}}\), upper limit \(b_{j,\text{UP}}\).

• For \(\mathbf{e}_{\textrm{td}(1+\tau_2)}\): lower limit \(b_{k,\text{LO}}\), upper limit \(b_{k,\text{UP}}\).

We conducted a targeted experimental investigation to optimize the grid for observables and minimize the model’s estimation error. Since finer grids yield too many indicator observables that hinder real-time training, we limited the grid to 22 divisions \(\mathbf{b}_{ijk}\), yielding 21 subregions \(S_{ijk}\) per time delay. A power function (exponent 1.8) was applied to adjust grid spacing, aligning processed EMG plateaus (red, Fig. 9 (b)) with 5 grip force levels (blue, Fig. 9 (b) and Table ¿tbl:tab:rbd-design?). The exponent was determined via regression analysis. The gridded indicator observable \(\mathbf{e}_{\mathrm{I},S_{ijk},\tau_1,\tau_2} \in \{0,1\}\) is defined as: \[\begin{align} \mathbf{e}_{\mathrm{I},S_{ijk},\tau_1,\tau_2} = & \begin{cases} & \big(b_{i,\text{LO}} \leq \mathbf{e}_{\textrm{td}(1)} \leq b_{i,\text{UP}}\big) \, \text{ and }\\ & \big(b_{j,\text{LO}} \leq \mathbf{e}_{\textrm{td}(1+\tau_1)} \leq b_{j,\text{UP}}\big) \,\text{ and } \\ & \big(b_{k,\text{LO}} \leq \mathbf{e}_{\textrm{td}(1+\tau_2)} \leq b_{k,\text{UP}}\big), \end{cases} \\ & \text{for} \;i, j, k = 0, \dots, 20. \end{align} \label{eqn:gridded-indicator-obs}\tag{5}\] Lifting in this manner can produce many empty (all-zero) or sparse (mostly zero) observables, potentially leading to overfitting of the estimation model. We applied a constraint to retain only observables with at least 0.1 % density during algorithm testing to mitigate this. Gridding was performed simultaneously on time delays \(\mathbf{e}_{\textrm{td}(1)}\), \(\mathbf{e}_{\textrm{td}(30)}\), and \(\mathbf{e}_{\textrm{td}(60)}\) (\(\tau_1=29,\;\tau_2=59 \to \mathbf{e}_{\mathrm{I},S_{ijk},29,59}\)), initially producing \(21^3 = \num{9261}\) observables (Fig. 6), most of which were discarded due to high sparsity. The final lifted data matrices were assembled by stacking time-delayed observables augmented with gridded indicators for processed EMG and zero rows for grip force. \[{!}{\( \begin{align} E & = \phi\big(\boldsymbol{e}_\text{proc}\big) = \big( \mathbf{e}_{0}, \mathbf{e}_{\text{td}(1)}, \ldots, \mathbf{e}_{\text{td}(60)}, \ldots, \mathbf{e}_{\mathrm{I},S_{ijk},29,59}, \ldots \big)^T \\ G & = \psi\big(\boldsymbol{g}\big) = \big( \mathbf{g}_{0}, \mathbf{g}_{\text{td}(1)}, \ldots, \mathbf{g}_{\text{td}(60)}, \ldots, \mathbf{0}, \ldots \big)^T. \end{align} \)} \label{eqn:final-lifting-estimation}\tag{6}\] Using (3 ) and (6 ), the Koopman estimation model \(\overline{\mathcal{K}}_\text{e}\) was trained on the full 20 s–30 s calibration experiment in approximately 1.5 s. The resulting grip force approximations \(\boldsymbol{g}_\text{e}\) were thresholded to a minimum value of −1 N.

2.3.3 Forecasting Grip Force↩︎

We propose a novel methodology for short-term grip force forecasting over a future horizon of 0.5 s, leveraging previously estimated grip force \(\boldsymbol{g}_\text{e}\) and the Koopman operator for dynamic systems in 2 . The Koopman operator is well-suited for this task due to its proven adaptability in dynamically changing environments [16], such as object grasping. Our goal is to determine the Koopman operator \(\mathcal{K}_\text{p}\) that predicts the system’s evolution by mapping observables within the same space (\(\mathcal{K}_\text{p}: \mathcal{O}_\mathcal{N} \to \mathcal{O}_\mathcal{N}\)) for forecasting. This approach enables fast training on data mini-batches, ensuring adaptation to the system’s latest state. To approximate the potentially infinite-dimensional \(\mathcal{K}_\text{p}\), we employ Dynamic Mode Decomposition (DMD) [21] to compute its spectral properties as Ritz pairs (\(\lambda_j\) - Ritz values and \(\mathbf{z}_j\) - Ritz vectors). \[\begin{align} \mathcal{K}_\text{p}Z & = Z\Lambda, \\ Z & = \left( \mathbf{z}_1, \ldots, \mathbf{z}_j, \ldots, \mathbf{z}_r \right), \\ \Lambda & = \mathop{\mathrm{diag}}{\left( \lambda_1, \ldots, \lambda_j, \ldots, \lambda_r \right)}, \end{align} \label{eqn:dmd-ritz-pairs}\tag{7}\] where \(r\) represents the total number of Koopman modes after DMD. To solve 7 and obtain the spectral properties of the Koopman operator, we used a lifted input matrix of grip force estimations \(\boldsymbol{g}_\text{e}\) and employed the pyKMD suite⁷, which implements the Refined Rayleigh-Ritz Data-Driven Modal Decomposition (DDMD_RRR) method alongside QR compression [23]. Enhanced DDMD_RRR refines Ritz vectors \(\mathbf{z}_j\) and improves spectral accuracy by minimizing residuals in a data-driven setting, while QR compression ensures efficient real-time execution. We applied Hankel-DMD time-delay embedding [24] to \(\boldsymbol{g}_\text{e}\), as done for the estimation Koopman operator in 4 , effectively capturing temporal dynamics for grip force prediction. The number of time delays was kept within 4–10, fine-tuned during optimization. Before lifting, \(\boldsymbol{g}_\text{e}\) was smoothed using Locally Weighted Scatterplot Smoothing (LOWESS) [25] with a single iteration and window size calculated as: (smoothing coefficient \(\times\) batch size) = (1.1–1.9 \(\times 496\)). This smoothing reduced spikes and prediction errors by incorporating data from current and previous batches, with the smoothing coefficient fine-tuned as a hyperparameter. Additional lifting produced an interaction (combination) matrix \(G_{\textrm{e,int}}\) from the \(\ln\)-transformed time-delay observables of \(\boldsymbol{g}_\text{e}\) (see (4 )) to reduce prediction error. A vertical shift of \(\require{physics} +\qty{10}{\N}\) was applied to ensure all values were positive for the \(\ln\) transformation. The interaction component of the forecast input matrix is: \[{!}{\( \begin{align} G_{\textrm{e,int}} = \begin{bmatrix} \ln{\boldsymbol{g}_{\textrm{e}}[1]} \ln{\boldsymbol{g}_{\textrm{e}}[2]} & \ln{\boldsymbol{g}_{\textrm{e}}[2]} \ln{\boldsymbol{g}_{\textrm{e}}[3]} & \cdots & \ln{\boldsymbol{g}_{\textrm{e}}[N-d]} \ln{\boldsymbol{g}_{\textrm{e}}[N-d+1]} \\ \ln{\boldsymbol{g}_{\textrm{e}}[1]} \ln{\boldsymbol{g}_{\textrm{e}}[3]} & \ln{\boldsymbol{g}_{\textrm{e}}[2]} \ln{\boldsymbol{g}_{\textrm{e}}[4]} & \cdots & \ln{\boldsymbol{g}_{\textrm{e}}[N-d]} \ln{\boldsymbol{g}_{\textrm{e}}[N-d+2]} \\ \vdots & \vdots & \ddots & \vdots \\ \ln{\boldsymbol{g}_{\textrm{e}}[d-1]} \ln{\boldsymbol{g}_{\textrm{e}}[d+1]} & \ln{\boldsymbol{g}_{\textrm{e}}[d]} \ln{\boldsymbol{g}_{\textrm{e}}[d+2]} & \cdots & \ln{\boldsymbol{g}_{\textrm{e}}[N-2]} \ln{\boldsymbol{g}_{\textrm{e}}[N]} \\ \ln{\boldsymbol{g}_{\textrm{e}}[d]} \ln{\boldsymbol{g}_{\textrm{e}}[d+1]} & \ln{\boldsymbol{g}_{\textrm{e}}[d+1]} \ln{\boldsymbol{g}_{\textrm{e}}[d+2]} & \cdots & \ln{\boldsymbol{g}_{\textrm{e}}[N-1]} \ln{\boldsymbol{g}_{\textrm{e}}[N]} \end{bmatrix} \end{align} \)} \label{eqn:log-interaction-lifting}\tag{8}\] The complete lifted input matrix is obtained by stacking 8 with the time-delay embedding of \(\mathbf{g}_\text{e}\), yielding a sequence of snapshots \(\mathbf{g}_{\textrm{e,lift}(i)}\): \[\begin{align} G_{\textrm{e,lift}} = \begin{pmatrix} G_{\textrm{e,int}} \\ G_{\textrm{e,td}} \end{pmatrix} = \big( \mathbf{g}_{\textrm{e,lift}(1)}, \ldots, \mathbf{g}_{\textrm{e,lift}(N-d)} \big). \end{align} \label{eqn:lift-in-data}\tag{9}\] Before applying DMD to 9 , a snapshot thinning step was performed by removing nearby columns in \(G_{\textrm{e,lift}}\) while preserving the time-delay structure, effectively reducing computational burden [22]. Thinning demonstrated no loss of accuracy during subsequent hyperparameter tuning within the range 3–8. Additionally, it enabled forecasting at 16 Hz–41 Hz, sufficient for capturing decision-making patterns during gripping [18]. After obtaining the spectral decomposition of the Koopman operator \(\mathcal{K}_\text{p}\), we compute the Koopman amplitudes \(\boldsymbol{\alpha} = [\alpha_1, \alpha_2, \dots, \alpha_\ell]^\top \in \mathbb{C}^\ell\), to reconstruct the input data matrix or predict the next state. This step is formulated as a least-squares minimization problem: \[\begin{align} \underset{\boldsymbol{\alpha} \in \mathbb{C}^\ell}{\min} \sum_{i=1}^{N-d} \left\| \mathbf{g}_{\textrm{e,lift}(i)} - \sum_{j=1}^{\ell} \mathbf{z}_j \alpha_j \lambda_j^{i-1} \right\|_2^2, \end{align} \label{eqn:coeffs-ls-min}\tag{10}\] where \(\ell\) denotes the number of modes kept after mode reduction \(r \to \ell\). Implemented solver from [26] can efficiently solve (10 ) using normal equations or, if the problem is ill-conditioned, a QR factorization-based solver, both of which are contained within the pyKMD framework. After obtaining predictions, as a final step, too low or too high values were thresholded to the minimum and maximum grip force from the calibration experiment. Finally, forecasting future snapshots in horizon \(h\) can be approximated using: \[\begin{align} \mathbf{g}_{\textrm{e,lift}(N-d+h)} \approx \sum_{j=1}^{\ell} \mathbf{z}_j \alpha_j \lambda_j^{N-d+h-1}, \quad h = 1, \dots. \end{align} \label{eqn:predict-snapshot}\tag{11}\] The reported error metric for estimating and predicting grip force is the Weighted Mean Absolute Percentage Error (wMAPE): \[\text{wMAPE} = \frac{\sum_{i=1}^{N} |\hat{g_i} - g_i|}{\sum_{i=1}^{N} |g_i|} \label{eqn:wmape}\tag{12}\] This relative metric was chosen because it effectively handles close-to-zero values by normalizing absolute errors with the sum of actual values, making it suitable for comparison across measurements with varying absolute grip magnitudes.

Final hyperparameter tuning for the prediction model was conducted using a grid search across five parameters: number of Koopman modes after reduction \(\ell\), number of time delays \(d\), batch window modifier coefficient for smoothing, batch window modifier for prediction, and thinning step. Based on the hyperparameter tuning results shown in Fig. 7, with runs featuring extremely high errors excluded, it can be concluded that the minimum wMAPE error is achieved with four Koopman modes, 7–10 time delays, batch window modifier coefficient for smoothing in the range of 1.1–1.2, thinning step in the range of 7–8, and a batch window size modifier coefficient in the range of 1.2–1.4. To select the final hyperparameters for the prediction algorithm, the median wMAPE value across all measurements was also computed to enhance metric robustness. The hyperparameters resulting in the minimal sum of mean and median wMAPE are presented in Table ¿tbl:tab:hyperparameter-tuning?. The final error metric using these optimal hyperparameters is discussed in Section 3.

3 Results & Experimental Validation↩︎

The analysis of the most sensitive frequencies from the multi-step SA (Fig. 5) shows the need to first act on reducing the lower frequency spectral components (\(\leq\)​14 Hz) and those in the 46 Hz–50 Hz range to improve the peak cross-correlations. Once the bounds of these initial sensitive spectral components are narrowed to near-optimal values, the range of sensitive frequencies expands, filling the gap between 14 Hz–46 Hz and extending towards higher-frequency components (arrows in Fig. 5). Typically, parameter optimization would be performed within the narrowed decision vector bounds. However, it is unnecessary in this case due to the descriptive statistics of the resulting mean peak cross-correlation computed from all LH samples after the final SA step. Descriptive statistics for different positions are:

• P1: mean 0.956, SD 4.13 × 10⁻⁴, min. 0.954, max. 0.958,

• P2: mean 0.960, SD 5.35 × 10⁻⁴, min. 0.958, max. 0.962.

Because the standard deviation of the mean peak cross-correlations indicates extremely low variability, the average values between each decision variable’s upper and lower bounds will be used as the optimal mask values, as shown in Fig. 8. Sections of the optimal spectral mask align with findings from [20], where lower-frequency components associated with wire movements are either attenuated or removed. The attenuation is nearly linear, with DC component and 2 Hz fully filtered out, 10 Hz reduced by 50 %, and 18 Hz left unaffected. Frequencies from 20 Hz–48 Hz require amplification following an inverted U-shaped curve, starting and ending at 25 %, with a peak amplification of 150 % between 32 Hz–42 Hz.

The inverted U is followed by a sharp dip at the 50 Hz component, indicating that electrical ground noise impacts the recorded signal. The 50 Hz spectral component proves particularly challenging to process, as neither retaining it at nominal amplitudes (with a BP filter) nor completely removing it (using a notch filter) is sufficient. Our optimizations show that maintaining 50 Hz at 37.5 % of its amplitude preserves enough of the signal’s power for grip force modeling. This finding supports the conclusions of [20], which state that a significant portion of the EMG signal’s power spectrum is concentrated at 50 Hz and should not be entirely filtered out.

The next section of the mask targets amplifying the mid-frequency spectra, where most of the signal’s power is concentrated. As the amplitudes of the spectral components decrease, the mask progressively amplifies them, increasing approximately linearly from 50 % at 52 Hz to 450 % at 110 Hz. This amplification trend plateaus at 425 %–450 % around 110 Hz and extends up to 202 Hz. A possible explanation for the need to amplify mid-frequency spectral components may lie in three physiological phenomena related to the spatial low-pass filtering of recorded EMG signals: MU structure, volume conduction, and electrode positioning [2]. First, the MU structure introduces spatial smoothing at the signal source due to the scattered arrangement of muscle fibers within a MU. Volume conduction refers to the signal’s attenuation and distortion as it propagates through biological tissues from the source to the skin surface, where it is measured. Finally, electrode positioning contributes to low-pass filtering, as the electrodes capture an averaged signal over the area they cover, further reducing higher-frequency content.

We conducted a stepwise ablation study on the spectral mask plateau (110 Hz–202 Hz) to validate our conclusions. We systematically excluded frequency bands in 4 Hz increments, starting from 110 Hz–202 Hz and progressing to 198 Hz–202 Hz. Each step resulted in a lower-than-optimal mean peak cross-correlation between processed EMG and measured grip force.

Our SA revealed the minimal influence of higher-frequency components (204 Hz–498 Hz). We performed an additional ablation study by varying the mask from 0–5 in this range. As no significant changes in mean peak cross-correlation were observed, we set the mask to 0 for the 204 Hz–498 Hz range, effectively shielding the system from potential noise in this band. This approach to processing sEMG signals is innovative, and some of the resulting observations warrant further investigation for a clearer understanding.

By applying the signal processing methods outlined in Section 2.2, the processed EMG signals, shown in red in Fig. 9 (b), demonstrate a clear high cross-correlation with the measured grip force. Table ¿tbl:tab:correlations-summary? summarizes peak cross-correlations for all measurements at both sensing positions, showing very strong correlations. The analysis reveals that EMG signals lag behind grip force measurements, with delays ranging from 0 ms–156 ms.

The RBD experimental design allows for investigating the effects of subject and sensing location on the wMAPE estimation metric. Using the procedure described in Section 2.3.2, wMAPE was calculated for 52 experiment runs (see Table ¿tbl:tab:rbd-estimate-wmape?). The mean estimation wMAPE of 5.5 % reflects fit error, as the model was trained on each subject’s full dataset and evaluated on 0.5 s batches without a test set, given Koopman operators’ deterministic linearity. These results demonstrate effective grip force estimation from a single sensing position. In addition to single-measurement error metrics, we computed means and effects across blocks containing all measurements for each position and subject. The effects were computed as the difference between the block mean and the overall mean wMAPE. Although the means and effects varied considerably across subjects, the mean for P1 was only 0.2 % worse than the overall mean, while the mean for P2 was 0.2 % better.

An analysis of variance (ANOVA) was conducted on the blocked RBD with estimation wMAPE to assess the significance of the position effect using a 5 % significance level. The results, presented in Table ¿tbl:tab:rbd-aov?, show that, while the effect of the subject on wMAPE was significant (\(p\)-value = 0.015), the effect of position (when the subject effect is removed) was not (\(p\)-value = 0.422). It can be concluded that the placement of EMG electrodes along the flexor carpi ulnaris muscle, whether in position 1 or 2, does not significantly influence the estimation error. Three representative examples of estimated grip force, with errors corresponding approximately to the first, second, and third quartiles, are highlighted in yellow in Fig. 9 (b). The complete set of graphs is available in the first author’s GitHub repository⁸.

Following a similar approach to the estimations, the methodology from Section 2.3.3 was applied to short-term grip force forecasting in 0.5 s mini-batches over a 0.5 s horizon. The adaptive forecasting model was trained on the previous \(\require{physics} 1.3 \times \qty{0.5}{s}\) batches and tested on the subsequent \(\require{physics} 1 \times \qty{0.5}{s}\) batch. wMAPE values were computed for all measurements and reported as RBD results in Table ¿tbl:tab:rbd-predict-wmape?. Using the optimal hyperparameters from Table ¿tbl:tab:hyperparameter-tuning?, we achieved approximately 17.9 % overall forecast wMAPE. While the effect of the subject varied, the effect of position on wMAPE remained within the  ± 0.1 % range. The ANOVA analysis in Table ¿tbl:tab:rbd-aov? confirmed the non-significant effect of position (\(p = \num{0.853}\)) and subject (\(p = \num{0.06}\)) on forecast error, further demonstrating robustness to electrode placement along the flexor carpi ulnaris muscle. Examples of forecasts with error metric values approximately corresponding to the first, second, and third quartiles are shown as red dots in Fig. 9 (c), with the smoothed grip force estimation, which serves as the input signal for forecasting, displayed in yellow.

Fig. 9 contrasts all described sEMG and grip signals, from raw and processed EMG to estimated and measured grip force, resulting in smoothed estimates and short-term batch predictions.

4 Conclusion↩︎

This work focuses on developing real-time procedures for sensing hand muscle activity and estimating grip force, with potential applications in controlling rehabilitation devices. The methodology involves sensing muscle activity (physiological signals) using non-invasive surface electromyography (sEMG) sensors. The signals are cross-correlated with synchronously collected hand grip force, determined using a force-calibrated hand dynamometer as part of a properly designed experiment. We devised the composition of the sEMG signal processing steps to achieve a strong mean peak cross-correlation (\sim 96 %) between the EMG signal and grip force. We optimized decision variables (spectral mask, smoothing window size, and decay rate) to maximize cross-correlation. Finally, using processed EMG from a single sensing location, the static Koopman operator enables grip force estimation, while the dynamic operator facilitates short-term prediction. The algorithm follows these steps: processing raw EMG, transforming and lifting the data for estimation, applying the estimation model, smoothing the estimates and lifting them for prediction, continuously retraining the prediction model online, and generating forecasts.

While creating a single large model that would generalize across different individuals and sessions would be challenging due to the inherent variability in EMG signals [2], the proposed Koopman methodology requires less than 30 s to perform a one-time, patient- and session-specific calibration experiment using both sEMG sensors and a hand dynamometer, followed by 1.5 s for estimation model training. Experimental assessment demonstrated the forecasting model’s rapid training and execution, requiring only 30 ms to generate 0.5 s predictions after receiving the latest 0.5 s data mini-batch. This confirms the framework’s readiness for real-time implementation. Our method demonstrated robustness across sensing locations, wMAPE of 5.5 % for estimations and 17.9 % for predictions. To contextualize, we compared our results to a baseline [27] which uses frequency-domain BP filtering (STFT with 0.5 s window and 50 % overlap), RMS smoothing (0.5 s window), and ordinary least squares regression. This baseline showed significantly higher wMAPE of 24.4 % for estimations and requires an extra FFT-to-IFFT step, complicating real-time implementation.

Future work will integrate developed Robot Operating System (ROS) modules for sensor interaction, data acquisition, and signal processing with the Koopman real-time force estimation and prediction module. Advancements in simplifying hand kinematics from [28] will be incorporated to develop a comprehensive system, including a human agent model with various grasp modalities. This approach aims to accelerate the development of robotic rehabilitation devices, enhancing their effectiveness and adaptability to individual patient needs.

5 Appendix↩︎

5.1 Experiment Setup↩︎

5.1.1 Hardware Components↩︎

Shimmer3 EMG Unit⁹ was used for muscle activity sensing. A maximum internal gain of 12 was set, and an internal calibration procedure was followed. Two electrodes at a single sensing location, along with a third reference electrode, were used to address the challenge of the small EMG signal relative to noise. This setup employs Common Mode Rejection to eliminate shared noise through signal subtraction, effectively preserving and amplifying the localized EMG signal. For grip force measuring, a Vernier Go Direct® Hand Dynamometer¹⁰ was used, with a force range 0 N–550 N, resolution 0.05 N, and uncertainty 1.96 SD.

5.1.2 Dynamometer Calibration↩︎

Dynamometer calibration was performed in-house according to the ASTM E74 norm. Calibration weights of OIML classes F1, M1, and M3 were used, and the laboratory temperature was maintained at \sim 23 °C. Calibration forces spanned: 5 N, 20 N, 50 N, 100 N, 150 N, 200 N, 250 N, 300 N, 350 N, 400 N, 450 N, 500 N, & 550 N. The instrument was first preloaded from min to max force to establish hysteresis. Each measurement was preceded by reducing the instrument to zero between successive loadings, and the entire process was repeated three times. For measured forces greater than 50 N, the uncertainty as a percentage of the measured force is less than 4 %, while for forces below 50 N, the percentage can be as high as 10 %. The obtained calibration equation that relates the raw dynamometer signal \(g_\textrm{r}\) to the reported grip force \(g\) is:

\[g = 1.063 g_\textrm{r} - 2.588 \times 10^{-4} g_\textrm{r}^2 - 9.003 \times 10^{-8} g_\textrm{r}^3 + 7.615 \times 10^{-10} g_\textrm{r}^4 \label{eqn:calibration}\tag{13}\] The most prominent term in (13 ) is linear, and \sim 1, due to demands of the ASTM E74 norm, a fourth-degree polynomial was deemed necessary. The ROS and Python were utilized to integrate devices and enable the acquisition of time-series raw EMG and calibrated dynamometer signals. The implemented modules are available in the first author’s GitHub repositories, namely godirect_ros¹¹ and shimmer_ros¹².

5.2 Sensitivity Analysis and Parameter Optimization↩︎

5.2.1 Peak Cross-correlation↩︎

Algorithm for computing peak cross-correlation between processed EMG signal \(\boldsymbol{e}_\textrm{proc}\) and grip force \(\boldsymbol{g}\) from single experiment run:

5.2.2 Batch Processing EMG Signal↩︎

Algorithm for computing processed EMG signal \(\mathbf{e}_{\text{proc}}\) from raw EMG data \(\mathbf{e}_\text{raw}\) by sequentially processing mini-batches \(\mathbf{e}_\textrm{raw,batch}\) (496 points each) using optimized spectral mask \(\mathbf{w}_\textrm{mask}\), exponential moving average window size \(w_\text{EMA}\), and decay factor \(\alpha_\text{EMA}\).

5.2.3 Preliminary SA↩︎

Grouped Sobol variance-based SA was conducted independently for P1 and P2 EMG measurements. For each position, \(2^{16}\) samples of 250-dimensional decision vectors were generated using Saltelli’s sampling scheme. Bootstrapping with \(2^{16}\) resamples was applied to account for uncertainty, creating additional datasets by sampling the original one with replacement. SIs were calculated for each bootstrapped dataset, and empirical distributions were generated. From these distributions, 95 % confidence intervals were derived.

5.2.4 Iterative Multi-step Ungrouped SA↩︎

RBD-FAST SA method returned only first-order SIs, with \(2^{16}\) decision vector samples and resampling with bootstrapping using 8192 (\(2^{13}\)) samples for CI computation. After each step, the projections on the most sensitive spectral components, or smoothing parameters, were visualized using LH sampling with 10000 samples, and manual bounds narrowing was performed. The first step identified the 2 Hz spectral component as the most influential on cross-correlation variance, with an SI in the 78 %–83 % range. The LH sample projection onto the 2 Hz component showed that increasing its mask variable significantly reduces the mean peak cross-correlation, prompting the narrowing of its variable bounds to the 0–0.5 range. Additionally, LH sampling was projected onto the next three spectral components with the highest SI from the first SA step. Similar trends to those observed for the 2 Hz component were noted, prompting the narrowing of bounds to 0–1 for 4 Hz, 0–2 for 6 Hz, and 0–3 for 50 Hz. Another 20 steps (\(2^\text{nd}-21^\text{th}\)) of iterative SA and optimization were performed similarly (see Fig. 5). Figures displaying LH sampling projections, along with smooth and mean lines that support the reasoning behind the narrowing of bounds, are in the first author’s GitHub repository¹³.

Acknowledgment↩︎

The authors thank students for conducting initial experiments and providing the experimental data. This work is supported by the Air Force Office of Scientific Research under award number FA9550-22-1-0531 and by the University of Rijeka under Grant uniri-iskusni-tehnic-23-47.

References↩︎

---

1. This work was supported by the U.S. Air Force Office of Scientific Research (AFOSR) under Award FA9550-22-1-0531 and by the University of Rijeka under Grant uniri-iskusni-tehnic-23-47. (Corresponding author: Ervin Kamenar).↩︎

2. Tomislav Bazina is with the University of Rijeka, Faculty of Engineering, 51000 Rijeka, Vukovarska 58, Croatia (e-mail: tomislav.bazina@uniri.hr).↩︎

3. Ervin Kamenar is with the University of Rijeka, Faculty of Engineering, 51000 Rijeka, Vukovarska 58, Croatia (e-mail: ekamenar@uniri.hr).↩︎

4. Igor Mezić is with the University of California, Santa Barbara, CA, 93101 USA (e-mail: mezic@ucsb.edu).↩︎

5. The authors are with AIMdyn, Inc., Santa Barbara, CA, 93101, USA (e-mail: tomislav@aimdyn.com; mezici@aimdyn.com; mfonoberova@aimdyn.com; ervin@aimdyn.com).↩︎

6. Informed consent was obtained from all human subjects involved. The research was performed under the oversight of the University of Rijeka Faculty of Engineering Ethics Committee, reference number 2409.17340.↩︎

7. https://apps.aimdyn.com/↩︎

8. https://github.com/tbazina/grip-emg-optimize↩︎

9. https://shimmersensing.com/product/shimmer3-emg-unit/↩︎

10. https://vernier.com/product/go-direct-hand-dynamometer/↩︎

11. https://github.com/tbazina/godirect_ros↩︎

12. https://github.com/tbazina/shimmer_ros↩︎

13. https://github.com/tbazina/grip-emg-optimize↩︎