# Debiasing Welch’s Method for Spectral Density Estimation

Lachlan C. Astfalck1,2, Adam M. Sykulski3, and Edward J. Cripps1

###### Abstract

Welch’s method provides an estimator of the power spectral density that is
statistically consistent. This is achieved by averaging over periodograms
calculated from overlapping segments of a time series. For a finite length
time series, while the variance of the estimator decreases as the number of
segments increase, the magnitude of the estimator’s bias increases: a bias-
variance trade-off ensues when setting the segment number. We address this
issue by providing a a novel method for debiasing Welch’s method which
maintains the computational complexity and asymptotic consistency, and leads
to improved finite-sample performance. Theoretical results are given for
fourth-order stationary processes with finite fourth-order moments and
absolutely continuous fourth-order cumulant spectrum. The significant bias
reduction is demonstrated with numerical simulation and an application to
real-world data, where several empirical metrics indicate our debiased
estimator compares favourably to Welch’s. Our estimator also permits irregular
spacing over frequency and we demonstrate how this may be employed for signal
compression and further variance reduction. Code accompanying this work is
available in the R and python languages.

Keywords: Welch’s method; bias reduction; tapering; non-parametric estimation;
periodogram. ††1 School of Physics, Mathematics & Computing, The University of
Western Australia, 35 Stirling Hwy, Crawley, Australia ††2 Oceans’ Graduate
School, The University of Western Australia, 35 Stirling Hwy, Crawley,
Australia ††3 Department of Mathematics, Imperial College London, South
Kensington, London SW7 2BX, U.K.

##  1 Introduction

The periodogram is the fundamental non-parametric estimator of the power
spectral density; however, it suffers from two deficiencies: the periodogram
is statistically inconsistent and is biased with finite samples. By averaging
multiple periodograms calculated from partitioned or overlapping segments of a
time series, Welch (1967) overcomes the issue of inconsistency to provide a
statistically consistent estimator of the power spectral density, eponymously
known as Welch’s estimator. Welch’s estimator is simple to implement and easy
to understand: all reasons that has led to its ubiquity in the applied
sciences. However, for finite length time series, the issue of bias remains.
Each segment’s periodogram is biased and Welch’s estimator inherits this bias.
Further, as the magnitude of the bias is inversely related to the length of
the segments, there is the classical inverse relationship between bias and
variance of the estimator; this leads to the undesirable property that we
become increasingly more certain in an estimator that is increasingly more
biased. In this paper, we introduce a debiased analogue to Welch’s estimator
that we entitle the debiased Welch estimator.

Historically, a number of methodologies have been developed to mitigate the
effects of bias. One such methodology is tapering (Thomson, 1982), where a
taper is applied to the data; although, typically this reduces correlation
between distant frequencies at the cost of increasing correlation between
nearby frequencies, and so bias remains. More recently, Sykulski et al. (2019)
provided methodology to debias spectral density estimates, but this
methodology relies on the availability of the parametric form of the power
spectral density. In this paper, we extend the results in Sykulski et al.
(2019) to non-parameteric estimation by constructing a natural basis
representation to the power spectral density. From this basis representation
we intentionally bias the bases, fit the biased bases to a biased estimate of
the power spectral density, and recover debiased coefficients that are used to
form our debiased non-parametric estimator. Tapering is ubiquitous in the
applied sciences and augments both parametric and non-parametric estimation of
the power spectral density, hence our method is developed to be applicable
with any choice of taper.

Novel insight is provided across two theorems. The first presents theoretical
results that make explicit how the properties of Welch’s estimator are a
function of the number of observations, number of segments, and the segment
lengths. This theorem makes the finite sample bias explicit and motivates the
construction of a new debiased alternative which we propose via a basis
representation. To obtain a good representation of the power spectral density,
however, we must specify a large number of basis functions and so estimation
of the coefficients via direct optimisation of the Whittle likelihood, as in
Sykulski et al. (2019), becomes computationally infeasible. The major
contribution of this work is therefore housed in the second theorem, where we
introduce methodology that leverages results established in the first theorem,
to provide asymptotic results that enable fast and direct computation of the
debiased Welch estimator. We build theory for time series requiring fourth-
order stationarity with finite fourth-order moments and absolutely continuous
fourth-order cumulant spectrum. The theory is established for all tapers,
subject to the standard conditions of squared summability of the tapering
sequence, thereby extending many classical results that use the raw
periodogram, for example, those in Anderson (1971), Brillinger (1975) and
Katznelson (1968).

The theoretical developments are also demonstrated by simulation. First, we
use a canonical auto-regressive model of order four studied throughout the
spectral estimation literature, see for example Percival et al. (1993). We use
this model to contrast the performance of the standard and debiased Welch
estimators as a function of the number of segments and the segment length.
Next, we demonstrate the debiased Welch estimator in a real-world case study
of coastal wave measurements. Here, we provide an example of where parametric
methods struggle to capture the complexity of the process, thus necessitating
the use of non-parametric methods. Finally, unlike Welch’s estimator, the
debiased Welch estimator does not require even spacing over frequency; we use
this to demonstrate how we may simultaneously debias and compress the signal.
We show this for a simulation of a stochastic process with log-linear decay in
the power spectral density over frequency, a pertinent case-study for many
applications. Code accompanies our results on multiple software platforms;
both the R and python implementations are available at github.com/TIDE-ITRH.

##  2 Spectral analysis

###  2.1 Fundamentals and definitions

Define by $\\{X_{t}\\}$ a discrete real-valued stochastic process, sampled at
interval $\Delta$ and indexed by $t\in\mathcal{Z}$, where $\mathcal{Z}$ is the
set of integers. The process $\\{X_{t}\\}$ may be either thought of as a true
discrete-time stochastic process, such as an auto-regressive model, or as
uniformly sampled observations of a continuous-time stochastic process
$\\{X(t)\\}$ in which case $X_{t}=X(\Delta t)$. We assume fourth-order
stationarity, with finite fourth-order moments and absolutely continuous
fourth-order cumulant spectrum, $\kappa(\alpha,\beta,\zeta)$. Define the auto-
covariance sequence of $\\{X_{t}\\}$ as
$\gamma(\tau)=\mathrm{E}[X_{t}X_{t-\tau}]$ for $\tau\in\mathcal{Z}$. Following
Wiener–Khinchin’s theorem, a necessary and sufficient condition for
$\gamma(\tau)$ to be the auto-covariance function of $\\{X_{t}\\}$ is that
there exists a power spectral distribution function $F(\omega)$ such that

$\gamma(\tau)=\frac{1}{2\pi}\int_{-\pi/\Delta}^{\pi/\Delta}e^{i\omega\tau\Delta}\;\mathrm{d}F(\omega),$
(1)

for $\omega\in[-\pi/\Delta,\pi/\Delta]$, where (1) is a Riemann-Stieltjes
integral and $\omega$ is defined in radians. Absolute continuity of
$\kappa(\alpha,\beta,\zeta)$ implies absolute summability in $\gamma(\tau)$.
Consequently, there exists a uniformly continuous power spectral density (see,
for instance, Section 4.2 Shumway and Stoffer, 2000), $f(\omega)$, such that
$f(\omega)\mathrm{d}\omega=\mathrm{d}F(\omega)$, and $f(\omega)$ forms a
Fourier pair with $\gamma(\tau)$ so that

$\gamma(\tau)=\frac{1}{2\pi}\int_{-\pi/\Delta}^{\pi/\Delta}f(\omega)e^{i\omega\tau\Delta}\;\mathrm{d}\omega,\hskip
8.53581ptf(\omega)=\Delta\sum_{\tau=-\infty}^{\infty}\gamma(\tau)e^{-i\omega\tau\Delta}.$
(2)

As $f(\omega)$ is uniformly continuous over this interval, and according to
the Vitali convergence theorem, $f(\omega)$ is Riemann integrable in $\omega$.

Say we observe from $\\{X_{t}\\}$ an $n$-dimensional observation
$\mathrm{X}_{n}=(x_{0},\dots,x_{n-1})$. The most common estimator of
$f(\omega)$ from $\mathrm{X}_{n}$ is the periodogram, $I_{n}(\omega)$, defined
as the squared modulus of the discrete Fourier transform of $\mathrm{X}_{n}$,

$I_{n}(\omega)=|J_{n}(\omega)|^{2},\hskip
8.53581ptJ_{n}(\omega)=\Delta^{1/2}n^{-1/2}\sum_{t=0}^{n-1}x_{t}e^{-i\omega
t\Delta}.$ (3)

Alternatively, we may calculate the periodogram via the discrete Fourier
transform of $\check{\gamma}(\tau)$, a biased estimator to $\gamma(\tau)$,

$I_{n}(\omega)=\Delta\sum_{\tau=-n+1}^{n-1}\check{\gamma}(\tau)e^{-i\omega\tau\Delta},\hskip
8.53581pt\check{\gamma}(\tau)=\frac{1}{n}\sum_{t=0}^{n-|\tau|-1}x_{t+|\tau|}x_{t}.$

The periodogram is a biased and inconsistent estimator of $f(\omega)$. It has
expectation $\mathrm{E}[I_{n}(\omega)]=f(\omega)*\mathcal{F}_{n}(\omega)$
where

$\mathcal{F}_{n}(\omega)=\frac{\Delta}{2\pi
n}\frac{\sin^{2}(n\omega\Delta/2)}{\sin^{2}(\omega\Delta/2)}$ (4)

is the Fejér kernel, and variance $\mathrm{var}[I_{n}(\omega)]\approx
f^{2}(\omega)$ which does not reduce with $n$, thus causing the inconsistency.
The periodogram is asymptotically unbiased as $\mathcal{F}_{n}(\omega)$ tends
to a Dirac delta function as $n\rightarrow\infty$. Often, practitioners will
appeal to this to argue that the periodogram is a suitable estimator; however,
the effects due to bias are related not only to $n$, but also to the dynamic
range of $f(\omega)$. Even with massive $n$, for common $f(\omega)$
encountered in the sciences, the periodogram can be too biased to be useful
(e.g. see the discussion in Thomson, 1982).

A standard technique to alleviate the bias in the periodogram is to apply a
taper, or window, to the data. Tapering pre-multiplies the data with a taper
to form the product $h_{t}X_{t}$ for $t=0,\dots,n-1$ where
$\sum_{t=0}^{n-1}h_{t}^{2}=1$. The sequence $\\{h_{t}\\}$ is called the data
taper. The tapered periodogram, otherwise known as the modified periodogram,
is defined as

$I_{n}(\omega;h)=|J_{n}(\omega;h)|^{2},\hskip
8.53581ptJ_{n}(\omega;h)=\Delta^{1/2}\sum_{t=0}^{n-1}h_{t}x_{t}e^{-i\omega
t\Delta}.$

The expectation of the tapered periodogram is
$\mathrm{E}[I_{n}(\omega;h)]=f(\omega)*\mathcal{H}_{n}(\omega)$ where
$\mathcal{H}_{n}(\omega)$ is the spectral window for $\\{h_{t}\\}$ defined by

$\mathcal{H}_{n}(\omega)=\Delta\left|\sum_{t=0}^{n-1}h_{t}e^{-i\omega
t\Delta}\right|^{2}.$

Tapering can be effective in removing bias; however, it is not a panacea. The
data taper is generally chosen so as to reduce bias in the form of broadband
blurring in $I_{n}(\omega;h)$; this comes with the trade-off that bias from
narrowband blurring is increased. Further, tapering, in effect, reduces the
amount of data and so the tapered periodogram has a higher variance as
compared to the standard periodogram. This variance can be reduced by applying
multiple tapers (multi-tapering), but this then reintroduces bias such that
the level of the bias-variance trade-off must be in practice controlled by the
practitioner. Setting $h_{t}=n^{-1/2}$ recovers the periodogram as
$\mathcal{H}_{n}(\omega)=\mathcal{F}_{n}(\omega)$. For this reason all
periodograms may be considered to be tapered, and so we will use the notation
$I_{n}(\omega;h)$ to refer to both tapered and standard periodograms
throughout.

###  2.2 Debiased parametric estimation

Rather than using $I_{n}(\omega;h)$ as a direct estimate to $f(\omega)$, it is
common to instead use $I_{n}(\omega;h)$ to fit some parametric form
$f(\omega;\vartheta)$, with Fourier pair $\gamma(\omega;\vartheta)$. To avoid
the unwieldy computation often associated with time-domain methods requiring
matrix inversion, it is standard to approximate the exact maximum likelihood
with a discretized form of the Whittle likelihood

$l_{\mathrm{W}}(\vartheta)=-\sum_{\omega\in\Omega_{n}}\left(\log
f(\omega;\vartheta)+\frac{I_{n}(\omega;h)}{f(\omega;\vartheta)}\right),$ (5)

where $\Omega_{n}=2\pi(\Delta n)^{-1}(-\lfloor
n/2\rfloor,\dots,-1,0,1,\dots,\lceil n/2\rceil-1)$ is the set of Fourier
frequencies for an observed stochastic process of length $n$. Due to the
inherent bias in $I_{n}(\omega;h)$, parameter estimation from using (5) can be
biased (for example, see Dahlhaus, 1988). Sykulski et al. (2019) proposed a
solution to this problem, whereby the pseudolikelihood

$l_{\mathrm{D}}(\vartheta)=-\sum_{\omega\in\Omega_{n}}\left(\log\check{f}_{n}(\omega;\vartheta)+\frac{I_{n}(\omega;h)}{\check{f}_{n}(\omega;\vartheta)}\right)$
(6)

is defined. Here, we substitute
$\check{f}_{n}(\omega;\vartheta)=E[I_{n}(\omega;h)]=f(\omega;\vartheta)*\mathcal{H}_{n}(\omega)$
for $f(\omega;\vartheta)$ in (5). We use the caron diacritic throughout this
manuscript to indicate intentionally biased quantities. Sykulski et al. (2019)
further showed that the convolution
$f(\omega;\vartheta)*\mathcal{H}_{n}(\omega)$ can be efficiently calculated as

$\check{f}_{n}(\omega;\vartheta)=2\Delta\times\mathrm{Re}\left\\{\sum_{\tau=0}^{n-1}\left(\sum_{t=0}^{n-\tau-1}h_{t}h_{t+\tau}\right)\gamma(\tau;\vartheta)e^{-i\omega\tau\Delta}\right\\}-\Delta\times\gamma(0;\vartheta)$
(7)

using fast Fourier transforms. When $h_{t}=n^{-1/2}$, and so $I_{n}(\omega;h)$
is the periodogram, $\sum_{t=0}^{n-\tau-1}h_{t}h_{t+\tau}=1-\tau/n$, and (7)
is the standard form for the debiased Whittle likelihood as presented in
Equation (9) of Sykulski et al. (2019). This technique debiases estimates of
$\vartheta$, and hence $f(\omega;\vartheta)$ and $\gamma(\omega;\vartheta)$,
in parametric estimation.

###  2.3 Welch’s method for spectral density estimation

Welch’s method segments a time series $x_{0},\dots,x_{n-1}$ into $M$
partitioned or overlapping blocks of length $L$, and applies a length $L$ data
taper to each block. We define the periodogram of the $m$th block, for
$m=0,\dots,M-1$, as

$I_{L}^{m}(\omega;h)=\Delta\left|\sum_{t=0}^{L-1}h_{t}x_{(t+mS)}e^{-i\omega
t\Delta}\right|^{2},$

where $S$ is an integer-valued shift factor satisfying $0<S\leq L$, so that
$p=(L-S)/L$ is the percentage overlap of the blocks, and the relationship
between $M$, $L$ and $p$ is

$ML(1-p)+Lp=n,$ (8)

such that setting two of $\\{M,L,p\\}$ fixes the third for a given $n$.
Welch’s estimator is defined as

$\bar{I}_{L}(\omega;h)=\frac{1}{M}\sum_{m=0}^{M-1}I_{L}^{m}(\omega;h).$ (9)

The degree of overlap, $p$, and the choice of data taper, $\\{h_{t}\\}$, are
modelling choices often informed by either desired properties of the analysis
or preferences within scientific disciplines. For instance, it is standard in
the engineering sciences to assume a 50% overlap ($p=0.5$) with a Hamming
taper: see the MATLAB documentation of the function pwelch. Welch’s estimator
has expectation
$\mathrm{E}[\bar{I}_{L}(\omega;h)]=f(\omega)*\mathcal{H}_{L}(\omega)$ and
variance

$\displaystyle\mathrm{var}\left[\bar{I}_{L}(\omega;h)\right]=$
$\displaystyle\frac{1}{M^{2}}\sum_{m=0}^{M-1}\mathrm{var}\left[I_{L}^{m}(\omega;h)\right]+\frac{2}{M^{2}}\sum_{m<m^{\prime}}\mathrm{cov}\left[I_{L}^{m}(\omega;h),I_{L}^{m^{\prime}}(\omega;h)\right]$
$\displaystyle=$
$\displaystyle\frac{\mathrm{var}[I_{L}(\omega;h)]}{M}\left(1+\frac{2}{M}\sum_{m<m^{\prime}}\rho(\omega;|m-m^{\prime}|,h)\right),$
(10)

where
$\rho(\omega;|m-m^{\prime}|,h)=\mathrm{cor}\left[I_{L}^{m}(\omega;h),I_{L}^{m^{\prime}}(\omega;h)\right]$
and $\mathrm{var}[I_{L}(\omega;h)]=\mathrm{var}[I_{L}^{m}(\omega;h)]$ for all
$m$. We establish some important properties of Welch’s estimator in Theorem 1.

######  Theorem 1.

Assume a fourth-order stationary stochastic process $\\{X_{t}\\}$ with finite
fourth-order moments and absolutely continuous fourth-order cumulant spectrum.
Welch’s estimator, defined in (9), of the observed process
$\mathrm{X}_{n}=(x_{0},\dots,x_{n-1})$, with data taper $\\{h_{t}\\}$ that
satisfies $\sum h_{t}^{2}=1$, has the following properties.

  1. a)

The estimator, $\bar{I}_{L}(\omega;h)$, satisfies
$\bar{I}_{L}(\omega;h)=f(\omega)*\mathcal{H}_{L}(\omega)+\mathcal{O}_{p}(M^{-1/2})$.

  2. b)

The variance $\mathrm{var}[\bar{I}_{L}(\omega;h)]$ converges asymptotically,
in $L$, to

$\lim_{L\rightarrow\infty}\mathrm{var}[\bar{I}_{L}(\omega;h)]=c\;\mathrm{var}[I_{L}(\omega;h)]$

for all $\omega\neq 0,\omega_{\mathcal{N}}$. Here, $c$ is a constant of
proportionality that is constant over $\omega$, and $\omega_{\mathcal{N}}$
denotes the Nyquist frequency.

  3. c)

The quantity $M^{-1/2}\bar{I}_{L}(\omega;h)$ tends to a Gaussian distribution
as $M\rightarrow\infty$.

Proof. Proofs are provided in the Supplementary Material.

######  Remark 1.

As $n$ increases, then it is not possible to increase both $L$, the block
length, and $M$, the number of blocks, linearly with $n$ for a fixed overlap
$p$, see (8); a trade-off is required. For fixed $L$ and linearly increasing
$M$, Theorem 1a establishes that Welch’s estimate has the optimal convergence
rate of $\mathcal{O}_{p}(M^{-1/2})=\mathcal{O}_{p}(n^{-1/2})$, but the
estimator is inconsistent for all processes except white noise due to the bias
inherent in $f(\omega)*\mathcal{H}_{L}(\omega)$ for fixed $L$ as discussed. On
the other extreme, for fixed $M$ and linearly increasing $L$, the bias
asymptotically vanishes but the variance is $\mathcal{O}_{p}(1)$ (Theorem 1a)
and therefore the Welch estimate suffers the same inconsistency as the
periodogram (Theorem 1b). Consistency is only achieved by increasing both $L$
and $M$ with $n$ at reduced rates.

##  3 Debiasing Welch’s Method

###  3.1 Constructing biased bases

To construct a debiased and consistent estimate to $f(\omega)$ based on
Welch’s method, we use concepts proposed by Sykulski et al. (2019) but extend
them to a non-parametric setting. Specifically, here we construct from a
family of bases $B_{k}(\omega)$, a family of biased bases
$\check{B}_{k}(\omega)$. We assume the functional form
$f(\omega;\vartheta)=\sum_{k=1}^{K}a_{k}B_{k}(\omega)$ and set
$\check{f}(\omega;\vartheta)=\sum_{k=1}^{K}a_{k}\check{B}_{k}(\omega)$ in (6)
to provide debiased inference for the parameters
$\vartheta=(a_{1},\dots,a_{K})$. Here, the $B_{k}(\omega)$ are basis functions
specified from some family of bases chosen to represent some desirable
property of $f(\omega)$. We set
$\int_{-\infty}^{\infty}B_{k}(\omega)\mathrm{d}\omega=1$. Due to the linearity
of the integral in (2), assuming
$f(\omega;\vartheta)=\sum_{k=1}^{K}a_{k}B_{k}(\omega)$ implies
$\gamma(\tau;\vartheta)=\sum_{k=1}^{K}a_{k}\rho_{k}(\tau)$ where each
$\rho_{k}(\tau)$ are the auto-correlation functions calculated from
$B_{k}(\omega)$ as

$\rho_{k}(\tau)=\frac{1}{2\pi}\int_{-\pi/\Delta}^{\pi/\Delta}B_{k}(\omega)e^{i\omega\tau\Delta}\mathrm{d}\omega.$
(11)

Due to the distributivity of the convolution operator the expected periodogram
is
$\check{f}(\omega;\vartheta,h)=\mathrm{E}[I_{n}(\omega;h)]=\sum_{k=1}^{K}a_{k}\check{B}_{k}(\omega;h)$
where each $\check{B}_{k}(\omega;h)$ is calculated exactly using a fast
Fourier transform by substituting $\rho_{k}(\tau)$ for
$\gamma(\tau;\vartheta)$ in (7).

###  3.2 A Riemann approximation to the spectral density

Section 3.1 describes a general approach to bias any basis family. Here, we
provide specific details with respect to a Riemannian approximation of
$f(\omega)$. As $f(\omega)$ is Riemann integrable then we may approximate
$f(\omega)$ by summating contiguous rectangular functions; this approximation
becomes exact as the width of the rectangular functions tends to zero. Define
the rectangular function by

$\mathrm{rect}(\omega)=\begin{cases}1&|\omega|<1/2\\\ 1/2&|\omega|=1/2\\\
0&\text{elsewhere}\end{cases}$

and the symmetric rectangular function with centre $\omega^{\mathrm{c}}$ and
width $\delta$ as

$\mathrm{symrect}(\omega;\omega^{\mathrm{c}},\delta)=\mathrm{rect}\left(\frac{\omega-\omega^{\mathrm{c}}}{\delta}\right)+\mathrm{rect}\left(\frac{\omega+\omega^{\mathrm{c}}}{\delta}\right),$
(12)

for $\omega^{\mathrm{c}}\geq\delta/2$. As we are concerned with real-valued
time series, $f(\omega)=f(-\omega)$. Define
$(\omega_{0},\omega_{1},\dots,\omega_{K})$ as a partition of
$[0,\omega_{\mathcal{N}}]$. The mid-point Riemann approximation to $f(\omega)$
is

$f(\omega)\approx\sum_{k=1}^{K}f\left(\omega_{k}^{\mathrm{c}}\right)\mathrm{symrect}\left(\omega;\omega_{k}^{\mathrm{c}},\delta_{k}\right),\quad\omega\in[-\pi/\Delta,\pi/\Delta],$
(13)

with centres $\omega_{k}^{\mathrm{c}}=\omega_{k-1}+\delta_{k}/2$ and widths
$\delta_{k}=\omega_{k}-\omega_{k-1}$. The approximation in (13) converges to
equality as $K\rightarrow\infty$ and all $\delta_{k}\rightarrow 0$. This
approach does not require us to space the bases regularly; any partition can
be defined, allowing for irregularly spaced bases. For now, so as to mimic the
behaviour of Welch’s estimate, we assume the bases to be regularly spaced in
$\omega$, and so $\delta_{k}=\delta$ for all $k$. We define the centres and
widths of the symmetric rectangular functions so that for a maximum value of
$K=\lceil(L-1)/2\rceil$, $\omega_{k}^{\mathrm{c}}=k\delta-\delta/2$ and
$\delta=\pi/(\Delta K)$. Setting
$B_{k}(\omega)=\mathrm{symrect}\left(\omega;\omega_{k}^{\mathrm{c}},\delta\right)$,
and following the results of Tobar (2019), we calculate via (11)

$\rho_{k}\left(\tau;\omega_{k}^{\mathrm{c}},\delta\right)=\delta^{-1}\mathrm{sinc}\left(\tau\delta\right)\cos\left(\omega_{k}^{\mathrm{c}}\tau\right),$
(14)

where $\mathrm{sinc}(\tau)=(\pi\tau)^{-1}\sin(\pi\tau)$ and
$\mathrm{sinc}(0)=1$. To calculate $\check{B}_{k}(\omega;h)$ we substitute
$\rho_{k}\left(\tau;\omega_{k}^{\mathrm{c}},\delta\right)$ in (7).

###  3.3 Model fitting

For a given choice of bases $B_{k}(\omega)$ the final step is to estimate
$\vartheta=(a_{1},...,a_{k})$ to recover the spectral estimate
$f(\omega;\vartheta)=\sum_{k=1}^{K}a_{k}B_{k}(\omega)$. The preferred debiased
approach would be to specify
$\check{f}(\omega;\vartheta)=\sum_{k=1}^{K}a_{k}\check{B}_{k}(\omega;h)$ and
solve for the $a_{k}$ by maximising the debiased Whittle likelihood defined in
(6). For the approximation to $f(\omega)$ in (13) to be reasonable, however,
we are required to define a potentially large number of bases, $K$, and so it
is likely that maximising (6) is computationally infeasible using naive
optimisation routines.

To proceed we observe from Theorem 1c that $\bar{I}_{L}(\omega;h)$ converges,
with $M$, to be Gaussian, and so we model
$\bar{I}_{L}(\omega;h)=\sum_{k=1}^{K}a_{k}\check{B}_{k}(\omega;h)+\epsilon(\omega)$,
for heteroskedastic and Gaussian $\epsilon(\omega)$ with
$\mathrm{var}[\epsilon(\omega)]=\mathrm{var}[\bar{I}_{L}(\omega;h)]$. Define
the vector
$\check{\mathrm{B}}(\omega;h)=(\check{B}_{1}(\omega;h),\dots,\check{B}_{K}(\omega;h))$
as the column vector of all the biased bases at $\omega$. Weighted least
squares can be used to estimate a solution of the form

$\underset{\vartheta}{\arg\min}\;\sum_{\omega\in\Omega_{L}}\left\\{\mathrm{var}[\bar{I}_{L}(\omega;h)]^{-1}\left(\bar{I}_{L}(\omega;h)-\vartheta^{\mathrm{\scriptscriptstyle
T}}\check{\mathrm{B}}(\omega;h)\right)^{2}\right\\}.$ (15)

As $\mathrm{var}[\bar{I}_{L}(\omega;h)]$ is dependent on $f(\omega)$, the true
power spectral density we are trying to estimate, the solution to (15) is not
analytical. As an approximation to (15), we substitute
$\bar{I}_{L}(\omega;h)^{2}$ for $\mathrm{var}[\bar{I}_{L}(\omega;h)]$ and
instead solve

$\hat{\vartheta}=\underset{\vartheta}{\arg\min}\;\sum_{\omega\in\Omega_{L}}\left\\{\bar{I}_{L}(\omega;h)^{-2}\left(\bar{I}_{L}(\omega;h)-\vartheta^{\mathrm{\scriptscriptstyle
T}}\check{\mathrm{B}}(\omega;h)\right)^{2}\right\\}.$ (16)

The motivation for (16) is as follows. From Theorem 1b, as
$\lim_{L\rightarrow\infty}\mathrm{var}\left[\bar{I}_{L}(\omega;h)\right]=c\;\mathrm{var}\left[I_{L}(\omega;h)\right]$,
over all $\omega$, we substitute $\mathrm{var}\left[I_{L}(\omega;h)\right]$
for $\mathrm{var}\left[\bar{I}_{L}(\omega;h)\right]$ in (15). Further, as we
establish in Theorem 2 below, $\bar{I}_{L}(\omega;h)^{2}$ converges in $M$ and
$L$ to $\mathrm{var}[I_{L}(\omega;h)]$, and so we approximate
$\mathrm{var}[I_{L}(\omega;h)]$ by $\bar{I}_{L}(\omega;h)^{2}$, resulting in
(16). From $\hat{\vartheta}=(\hat{a}_{1},\dots,\hat{a}_{K})$, we define a
debiased functional estimate of $f(\omega)$ as
$f(\omega;\hat{\vartheta})=\sum_{k=1}^{K}\hat{a}_{k}B_{k}(\omega)$. For the
symmetric rectangular basis function of Section 3.2, we define the debiased
Welch estimator at a discrete set of frequencies as
$\hat{I}(\omega_{k})=\hat{a}_{k}B_{k}(\omega_{k})=\hat{a}_{k}/\delta_{k}$ for
$\omega_{k}\in\\{\omega_{1}^{\mathrm{c}},\dots,\omega_{K}^{\mathrm{c}}\\}$;
this is analogous to the frequencies for which Welch’s estimator is generally
defined.

######  Theorem 2.

Assume a fourth-order stationary process, $\\{X_{t}\\}$, with finite fourth-
order moments and absolutely continuous fourth-order cumulant spectrum,
$\kappa(\alpha,\beta,\zeta)$. The variance of the periodogram of any observed
length-$L$ segment of $\mathrm{X}_{n}=(x_{0},\dots,x_{n-1})$, for all
$\omega\neq 0,\omega_{\mathcal{N}}$, and with data taper $\\{h_{t}\\}$ that
satisfies $\sum h_{t}^{2}=1$, is

  1. a)

$\mathrm{var}[I_{L}(\omega;h)]=\bar{I}_{L}(\omega;h)^{2}+\mathcal{O}_{p}\left(\frac{1}{\sqrt{M}}+\frac{\log^{2}L}{L^{2}}\right)$,
when $\kappa(\alpha,\beta,\zeta)=0$ for all $\alpha,\beta,\zeta$; and

  2. b)

$\mathrm{var}[I_{L}(\omega;h)]=\bar{I}_{L}(\omega;h)^{2}+\mathcal{O}_{p}\left(\frac{1}{\sqrt{M}}+\frac{1}{L}\right)$,
otherwise.

Proof. The proof is provided in the Supplementary Material.

######  Remark 2.

From Theorem 1 of Sykulski et al. (2019), we find that the debiased Welch
estimator, $\hat{I}(\omega_{k})$, when computed through (6), is unbiased so
that $\mathrm{E}[\hat{I}(\omega_{k})]=f(\omega_{k};\vartheta)$. The weighted
least squares routine in (16) follows asymptotic arguments in Theorems 1 and 2
that result in an estimator that is asymptotically unbiased and
$\mathcal{O}_{p}(M^{-1/2})$ as long as $L$ and $M$ both increase with $n$, and
significantly debiased for finite samples compared to the standard Welch
estimator, as we will demonstrate in the following sections through simulation
and applications.

###  3.4 Computation

We summarize the computation of the debiased Welch estimator in Algorithm 1 in
the Supplementary Material, with two variants provided. The first variant,
shown in lines 3–8 of Algorithm 1, provides the option to space the $K$ bases
evenly or unevenly over frequency. Even spacing of the bases provides a closer
analogue to Welch’s estimator; here, the choice of integer valued
$K\leq\lceil(L-1)/2\rceil$ is the only additional user decision required by
the debiased Welch estimator. When there is large bias present in Welch’s
estimator, choosing $K$ close to $\lceil(L-1)/2\rceil$ can lead to an ill-
defined solution. In general, we recommend $K=\lceil(L-1)/4\rceil$ as a good
initial estimate which, in effect down-samples the Fourier frequencies by an
order of two. For uneven spacing, the centres and widths are all user-
specified. We discuss uneven spacing further in Section 6 with an example
provided for high-frequency slope estimation. The second variant of the model,
shown in lines 12–15 of Algorithm 1, provides an option to impose non-
negativity on $\vartheta$ in (16), constraining the solution space to
$\hat{\vartheta}\in\mathcal{R}^{+}$. The debiased Welch estimator is not
guaranteed to yield a non-negative solution; due to finite sampling error,
small negative values may result. This is particularly so when the choice of
$K$ is too high and so the solution to the linear system in (16) is unstable.
Efficiently computing non-negative least-squares estimates has been long
understood; in the accompanying code we implement results from Lawson and
Hanson (1974).

The computational order of the debiased Welch estimator, for both even and
uneven spacing, is $\mathcal{O}(\max\\{n\log L,K^{3}\\})$. Here, the first
term is the computational order of Welch’s estimator, and the second term is
the computational order to do the debiasing in (16). The number of bases $K$
scales with $L$ and so $\mathcal{O}(K^{3})\propto\mathcal{O}(L^{3})$ for even
spacing. From Theorem 2, we observe that the optimal rate with which to scale
$L$ with $n$ is at maximum $n^{1/3}$, and so the debiasing in (16) scales at
most $\mathcal{O}(n)$. Thus, the debiased Welch estimator retains the same
computational order of the standard Welch estimator, $\mathcal{O}(n\log L)$.

##  4 Simulations of autoregressive processes

To test empirical performance of the debiased Welch estimator, we run a
simulation study using a canonical model for biased spectral density
estimation: the ar(4) model of Percival et al. (1993). The model is specified
as $X_{t}=\sum_{j=1}^{4}\phi_{j}X_{t-j}+\epsilon_{t}$ where
$\epsilon_{t}\sim\mathcal{N}(0,\sigma^{2})$ and
$\\{\phi_{1},\phi_{2},\phi_{3},\phi_{4},\sigma\\}=\\{2.7607,-3.8106,2.6535,-0.9238,1\\}$.
To show units in subsequent plots, we measure $\\{X_{t}\\}$ in metres, $t$ in
seconds and assume $\Delta=1$. Throughout, we calculate standard periodograms,
corresponding to $h_{t}=L^{-1/2}$. Similar results with other tapers were
found, with the bias in Welch’s estimator being reduced but not eliminated.
First, we plot in Figure 1 the standard (left) and debiased (right) Welch
estimates with black solid lines, calculated from a single realisation from
the ar(4) model. We generate a time series of length $32,768$ ($2^{15}$) and
set $M=64$ ($2^{6}$) and $L=512$ ($2^{9}$) with zero overlap for both
estimates. In both plots, the true spectral density is shown by the dashed
grey line, and the expected value of Welch’s estimator, calculated from (7)
using the auto-covariance function of the ar(4) model, is shown by the solid
grey line. Visually, the debiased Welch estimate provides a closer estimate of
the true spectral density, especially in regions of low power where we know
the bias to be significant.

Figure 1: Welch estimates of an ar(4) random time series. The ar(4) model is
defined formally in the text; the true spectral density is shown by the dashed
grey line and the expected Welch estimate is shown by the solid grey line. The
Welch estimates are shown by solid black lines with the standard Welch
estimate shown left, and the debiased Welch estimate shown right.

We now explore how the bias, standard deviation and root-mean-square error of
the Welch estimators behave as a function of $M$ for fixed $L$. For each
$M\in\\{8,16,32,64,128,256\\}$ and $p\in\\{0,0.5\\}$ we generate 1000 random
time series from the ar(4) model with fixed segment length $L=1024$
($2^{10}$). For each $M$ and $p$, we calculate, over the ensemble, the
empirical mean absolute bias, standard deviation, and root-mean-square error,
at each frequency. We aggregate the scores by computing the mean log value
across frequency and for each metric; this makes the error process of the
metrics independent of the expected power at each frequency and hence
aggregates the frequencies with equal importance. Results are shown in Figure
2. The dashed lines show the results from the Welch estimator, the solid lines
from the debiased Welch estimator, the grey lines with an overlap of $p=0.5$
and the black lines with an overlap of $p=0$. The results of this simulation
study support a number of our theoretical claims. The black dot-dashed line
shows the $M^{-1/2}$ convergence rate which is consistent with the theory
established in Theorem 1a for the Welch estimator, and is maintained for the
debiased Welch estimator. The absolute bias tends to zero for the debiased
Welch estimator as the approximation in Theorem 2 improves with increasing
$M$, justifying the computation of the debiased Welch estimates using (16).
The absolute bias does not converge to zero for the standard Welch estimator
due to the finite-sample bias of periodograms of fixed length $L$ as already
established in Theorem 1a.

Figure 2: Empirical bias, standard deviation and root mean squared error of
the standard and debiased Welch estimates over $M$, for fixed $L$. Each data
point represents the metrics calculated over an 1000 member ensemble for each
frequency and the aggregated by the mean log values. The dashed lines
represent the Welch estimates, the solid lines the debiased Welch estimates,
the grey lines with an overlap of 0.5 and the black lines with an overlap of
0. Table 1:  Metrics calculated over an ensemble of 1000 members for fixed $M$
and varying $L$.

| $L=512$ | $L=1024$ | $L=2048$
---|---|---|---
| Welch | D.Welch | Welch | D.Welch | Welch | D.Welch
RMSE | 17 (13) | 19(15) | 12 (9.8) | 13 (10) | 8.8 (6.7) | 9.2 (6.9)
RMSLE | 17 (1.6) | 3.1 (2.5) | 9.9 (1.2) | 1.8 (1.4) | 5.6 (0.76) | 1.0 (0.74)
MAPE | 64 (8.6) | 13 (4.3) | 34 (4.5) | 9.5 (1.8) | 19 (2.3) | 8.1 (0.67)

  * •

RMSE, Root Mean Squared Error; RMSLE, Root Mean Squared Log Error; MAPE, Mean
Absolute Percentage Error.

Finally, we examine the performance of the Welch estimators as a function of
$L$, for fixed $M$. Similarly to before, we generate 1000 random time series,
but with fixed $M=64$ and $L\in\\{512,1024,2048\\}$, all with zero overlap. We
calculate three metrics over frequency, for each time series: the root-mean-
squared-error; the root-mean-squared-log-error, that is, the root-mean-
squared-error calculated over the log values; and the mean-absolute-
percentage-error. The results are summarized in Table 1 by the mean values and
standard deviation, in brackets, calculated over the ensemble. Root-mean-
squared-error is largely determined by the higher power regions in the
spectral estimates. In these regions we do not expect much deviation between
the estimators, as bias due to blurring predominately occurs in the low power
regions. This is reflected in the results where both Welch estimators have
similar root-mean-squared-error values across all $L$. Conversely, the root-
mean-squared-log-error is power normalising, that is, the performance of the
spectral estimates is calculated evenly across all frequencies, regardless of
the spectral density. Here, we see a clear increase in performance by the
debiased Welch estimator across all simulations with highly significant
separation between the scores. This is also reflected in the mean-absolute-
percentage-error, a metric commonly used in evaluating spectral density
estimation (see Percival et al., 1993).

##  5 Application to coastal wave monitoring

Figure 3: Contemporaneous measurements of sea-surface height from acoustic
(grey) and pressure (black-dashed) sensors, measured in Western Australia at
Torbay (top) and Oean Beach (bottom). The right hand plots show the
attenuation of the pressure signal over frequency and are dependent on the
mean water depths.

We now demonstrate the use of the debiased Welch estimator applied to a
measured data-set of coastal sea-surface heights. The data are measured from
two locations with differing mean water depths on the Southern coast of
Western Australia: Ocean Beach and Torbay. Contemporaneous measurements of
sea-surface height are taken from acoustic and pressure sensors. Acoustic
sensors provide measurements of the true process via acoustic ranging and
pressure sensors measure pressure variation on the sea-bed and translate this
to a measurement of sea surface height. Shallow water wave theory states that
pressure is attenuated through depth as

$\Lambda_{p}(z;k)=\frac{\cosh(k(d+z))}{\cosh(kd)},$ (17)

where $d$ is the undisturbed depth, $-z$ is depth from the sea-surface to the
sensor and, here, $k$ is wavenumber. We can also define (17) as a function of
frequency using the dispersion relation (Dean and Dalrymple, 1991). We show
two 5-minute chunks of the measured time series from both sites on the left of
Figure 3; the solid grey and dashed black lines are measurements from the
acoustic and pressure sensors, respectively. On the right of Figure 3 we show
the attenuation of the pressure signal over frequency. The high frequency
signal in the Torbay measurements see larger attenuation than in the Ocean
Beach measurements due to the deeper mean water level. Due to the attenuation
of the high frequency signal, in both records, we expect a high dynamic range
in the pressure signal, and hence non-ignorable spectral bias due to blurring.

Figure 4: Welch estimates of sea-surface height spectra from acoustic and
pressure sensors located at Torbay (top) and Ocean Beach (bottom). The left
plots show the standard Welch estimates of the acoustic (grey) and pressure
(black) measurements. The middle plots show the debiased Welch estimate of the
pressure data (black) and, for comparison, the Welch estimate (grey). The
right plots show the standard (grey) and debiased (black) Welch estimates of
the pressure measurements inverted through the attenuation function. The thin
grey line is the acoustic signal used as a notion of ground truth.

We calculate power spectral density estimates of a three-hour record of sea-
surface heights at each location; the results are shown in Figure 4. First, we
calculate the Welch estimates, shown in the left hand panels. Here, the grey
and black lines are from the acoustic measurements, respectively. As expected,
there is good agreement between the estimators at low frequencies, and the
pressure signal is attenuated at higher frequencies. The spectral shape of the
Ocean Beach record is interesting: there are two spectral peaks around
$0.06\;$Hz and $0.12\;$Hz, implying the presence of both a well formed swell
(low-frequency) and wind-sea (high-frequency), likely due to a local storm.
This behaviour is hard to capture parametrically and is a good example to
motivate the use of non-parametric methods. The black lines in the middle
plots are the debiased Welch estimates calculated from the pressure signal,
and the grey lines are the Welch estimates, for reference. We do not display
the debiased estimate for the acoustic signal as the bias correction is very
small here due to the small dynamic range. Finally, the right hand plot shows
the Welch estimates of the pressure measurements inverted through the
attenuation function to give an estimate of the true wave power spectra; the
thick grey and black lines show the inverted Welch and debiased Welch
estimates, respectively. These estimates are not shown above $0.18\;$Hz and
$0.28\;$Hz for Torbay and Ocean Beach, respectively, as this is where the
noise floor of the pressure sensor dominates the signal. The thin grey line is
the acoustic Welch estimate, also shown in the left hand plots, and is used
for a notion of the ground truth. Due to in-situ pre-processing of the
pressure signal, whereby the signal is demeaned of tidal effects, there is
some bias between the two signals at low-frequencies below $0.05\;$Hz across
the entire data-set. Overall, the debiased Welch estimator gives a better
estimate to the true process over broader range of frequencies. As compared to
pressure sensors, acoustic sensors are more expensive, and due to larger power
requirements, do not allow for long field deployments. Increasing the signal
extraction from pressure sensors motivates an increase in use which will, in
turn, enable cheaper and longer measurement campaigns.

##  6 Unevenly spaced bases for signal compression

To be analogous to Welch’s estimator, in Section 3 we define the debiased
Welch estimator with even basis spacing over frequency. As noted therein, we
do not require this, and can alternatively define the bases by an irregular
partition over frequency. For spectra whose value varies slowly over certain
neighbourhoods of frequency, we argue that a wider basis spacing would lead to
little reduction of accuracy in estimating the true process. Allowing for
wider, and irregularly spaced bases serves two purposes: first we achieve a
degree of signal compression, as compared to the Welch estimator; and second,
we reduce the variance of the estimator as there are more data informing the
coefficient of each basis. The code provided alongside this manuscript accepts
uneven bases by providing the $\omega^{\mathrm{c}}_{k}$ and $\delta_{k}$, as
defined in (13), as user defined inputs for all $k$.

To demonstrate this concept, we calculate Welch and debiased Welch estimates
of a random sample from a Matérn process parameterized as

$\gamma(\tau)=\frac{2\sigma^{2}}{\Gamma(\alpha-1/2)2^{\alpha-1/2}}|\tau|^{\alpha-1/2}\kappa_{\alpha-1/2}(|\tau|),\hskip
14.22636ptf(\omega)=\frac{\sigma^{2}}{(\omega^{2}+\lambda^{2})^{\alpha}},$

with $\sigma=1$, $\lambda=0.1$ and $\alpha=11/6$, where $\Gamma(\cdot)$
denotes the gamma function and $\kappa_{\alpha-1/2}(\cdot)$ the order
$\alpha-1/2$ modified Bessel function of the second kind. The true power
spectral density and the expected value of Welch’s estimator for the Matérn
process are shown in the plots in Figure 5 by the dashed and solid grey lines,
respectively. The Matérn power spectral density is asymptotically, in
$\omega$, log-linear and so we space the bases of the debiased Welch estimator
evenly in log-frequency. The location of the bin-centres,
$\omega_{k}^{\mathrm{c}}$, are shown by the rug-plot on the frequency axes in
Figure 5. Asymptotically log-linear processes are common in many sciences, for
example, background ocean spectra and turbulent dissipation in oceanography
(Garrett and Munk, 1972; Zulberti et al., 2022), stellar rotation periods and
red-shift in astronomy (Foreman-Mackey et al., 2017) and long-memory processes
(Palma, 2007), and so this example provides a useful case-study to many
applied fields. We sample a length $32,768$ $(2^{15})$ time series with
$\Delta=1$ and calculate the Welch estimate, with parameters $M=32$ $(2^{5})$,
$L=1024$ $(2^{10})$ and no overlap, shown in the left plots of Figure 5 in
black. The debiased Welch estimate, with log-spaced bases, is calculated and
shown in the right plots of Figure 5 in black. The bottom plots of Figure 5
are identical to the top, except with frequency shown on a log scale. Not only
do we reduce empirical bias in the estimate, we also compress the signal by a
factor of $\sim 50$ and reduce the standard deviation of the estimate from the
true spectral density.

Figure 5: Welch estimates of a Matérn power spectrum with uneven basis
spacing. The top plots show frequency on a linear scale and the bottom plots
on a log scale. The standard Welch estimate is shown left and debiased Welch
estimate, with uneven basis spacing, is shown right with black solid lines.
Locations of the Fourier frequencies and basis centres for the Welch and
debiased Welch estimates are shown by the rug plots on the frequency axes. In
all plots, the true power spectral distribution and expected Welch estimate
are shown by the dashed and solid grey lines, respectively.

## Acknowledgement

All authors are supported by the ARC ITRH for Transforming energy
Infrastructure through Digital Engineering (TIDE), Grant No. IH200100009. We
further thank Dr. Michael Cuttler and Dr. Jeff Hanson for providing us with
the wave data examined in Section 5.

## References

  * Anderson (1971) Anderson, T.W., 1971.  The statistical analysis of time series.  John Wiley & Sons.
  * Brillinger (1975) Brillinger, D.R., 1975.  Time series: data analysis and theory.  Holt, Rinehart and Winston Inc.
  * Dahlhaus (1988) Dahlhaus, R., 1988.  Small sample effects in time series analysis: a new asymptotic theory and a new estimate.  The Annals of Statistics , 808–841.
  * Dean and Dalrymple (1991) Dean, R.G., Dalrymple, R.A., 1991\.  Water wave mechanics for engineers and scientists. volume 2 of Advanced Series on Ocean Engineering.  World Scientific Publishing Co Pte Ltd.
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., Angus, R., 2017\.  Fast and scalable Gaussian process modeling with applications to astronomical time series.  The Astronomical Journal 154, 220\.
  * Garrett and Munk (1972) Garrett, C., Munk, W., 1972\.  Space-time scales of internal waves.  Geophysical Fluid Dynamics 3, 225–264.
  * Katznelson (1968) Katznelson, Y., 1968.  An introduction to harmonic analysis.  John Wiley & Sons, Inc.
  * Lawson and Hanson (1974) Lawson, C.L., Hanson, R.J., 1974\.  Solving least squares problems.  Prentice Hall.
  * Palma (2007) Palma, W., 2007.  Long-memory time series: theory and methods.  John Wiley & Sons.
  * Percival et al. (1993) Percival, D.B., Walden, A.T., et al., 1993\.  Spectral analysis for physical applications.  Cambridge University Press.
  * Shumway and Stoffer (2000) Shumway, R.H., Stoffer, D.S., 2000\.  Time series analysis and its applications. volume 3.  Springer.
  * Sykulski et al. (2019) Sykulski, A.M., Olhede, S.C., Guillaumin, A.P., Lilly, J.M., Early, J.J., 2019.  The debiased Whittle likelihood.  Biometrika 106, 251–266.
  * Thomson (1982) Thomson, D.J., 1982.  Spectrum estimation and harmonic analysis.  Proceedings of the IEEE 70, 1055–1096.
  * Tobar (2019) Tobar, F., 2019.  Band-limited Gaussian processes: The sinc kernel.  Advances in Neural Information Processing Systems 32\.
  * Welch (1967) Welch, P., 1967.  The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms.  IEEE Transactions on audio and electroacoustics 15, 70–73.
  * Zulberti et al. (2022) Zulberti, A., Jones, N., Rayson, M., Ivey, G., 2022\.  Mean and turbulent characteristics of a bottom mixing-layer forced by a strong surface tide and large amplitude internal waves.  Journal of Geophysical Research: Oceans 127, e2020JC017055.


# Exoskeleton-Mediated Physical Human-Human Interaction for a Sit-to-Stand
Rehabilitation Task

Lorenzo Vianello1, Emek Barış Küçüktabak1,2, Matthew Short1,3, Clément
Lhoste1, Lorenzo Amato1,4,
Kevin Lynch2, and Jose Pons1,2,3 1 Legs and Walking Lab of Shirley Ryan
AbilityLab, Chicago, IL, USA2 Center for Robotics and Biosystems of
Northwestern University, Evanston, IL, USA3 Department of Biomedical
Engineering, Northwestern University, Evanston, IL, USA4 The Biorobotics
Institute, Scuola Superiore Sant’Anna, 56025 Pontedera, Italy and Department
of Excellence in Robotics & AI, Scuola Superiore Sant’Anna, 56127 Pisa, Italy

###### Abstract

Sit-to-Stand (StS) is a fundamental daily activity that can be challenging for
stroke survivors due to strength, motor control, and proprioception deficits
in their lower limbs. Existing therapies involve repetitive StS exercises, but
these can be physically demanding for therapists while assistive devices may
limit patient participation and hinder motor learning. To address these
challenges, this work proposes the use of two lower-limb exoskeletons to
mediate physical interaction between therapists and patients during a StS
rehabilitative task. This approach offers several advantages, including
improved therapist-patient interaction, safety enforcement, and performance
quantification. The whole body control of the two exoskeletons transmits
online feedback between the two users, but at the same time assists in
movement and ensures balance, and thus helping subjects with greater
difficulty. In this study we present the architecture of the framework,
presenting and discussing some technical choices made in the design.

##  I Introduction

Transitioning from sitting to standing is fundamental to living independently
and is one of the most frequently executed functional tasks. This task,
commonly referred to as Sit-to-Stand (StS), can be difficult to perform for
post-stroke individuals due to deficits in strength, motor control and
proprioception of the lower limb. Compared to healthy individuals, post-stroke
patients present many problems related to their StS ability: asymmetric force
and weight distribution, reduced peak vertical reaction force, increased time
to complete the movement, and a larger mediolateral center of pressure
displacement [1]. In addition, there is a correlation between the ability to
perform StS tasks and clinical balance scales in post-stroke patients [2].

Typically, the therapist’s role in StS tasks is to instruct and guide the
patient toward more correct movements via physical assistance and helps those
with greater disabilities by providing partial to full weight support. These
approaches require a large amount of physical effort from the therapist which
can, over time, lead to work-related injuries [3].

For this reason, assistive devices are usually offered to patients who lack
the strength to support their own weight [4]. Nevertheless, many of these
devices require the patient to play mostly a passive role, resulting in a
smaller activation of motor learning mechanisms [5]. Moreover, they usually do
not allow the therapist to provide external physical guidance or assistance.
Eliminating the physical patient-therapist interaction significantly limits
the usage of therapist’s knowledge and ability, and therefore the efficacy of
the rehabilitation.

Figure 1: Sit-to-Stand rehabilitative framework: two exoskeletons are used to
mediate the interaction between two users performing a tracking task in which
the position of the two Center of Masses is used as a reference. The
controller of the two exoskeletons produces a threefold function: 1) they
virtually connect the movements of the two users; 2) they produce vertical
assistance on the two by limiting fatigue; and 3) they preserve balancing.

An infrastructure with two exoskeletons can be used to mediate physical
interaction between a patient and therapist similar to the previously proposed
human-robot-human systems [6]. Such an infrastructure has several advantages.
First, it allows a therapist to feel and guide their patient on multiple
joints and contact points simultaneously. Second, the exoskeleton can be
programmed to further assist and enforce defined safety criteria which are
particularly relevant during StS tasks. Third, interfacing with exoskeletons
facilitates the objective quantification of motor performance over many
repetitions, increasing reproducibility. Last but not least, this
infrastructure could enhance tele-rehabilitation practices, allowing
therapists and patients to physically interact in different remote settings
(e.g., at home, in the clinic).

To the best of our knowledge, there hasn’t been any prior presentation of a
rehabilitation framework using overground lower-limb exoskeletons to mediate
physical interaction between therapists and patients for StS task. In fact,
while overground systems offer the benefit of promoting natural movements
while engaging the balance, they can be challenging to control [7, 8, 9]. The
whole body controller of the exoskeleton must not only be designed to consider
the degree of impairment of the patient, but at the same time to meet safety
prerequisites such as maintaining balance [10].

In our prior research, we introduced a control algorithm that controls the
forces between a floating-base lower-limb exoskeleton and its wearer while
taking into account actuation constraints and safety requirements such as
velocity and torque limits [11]. Building on this formalization, we introduced
the practical application of exoskeleton-mediated physical interaction between
two users [12]. This approach combines the benefits of a whole-body
controller, including constraint and force regulation, with spring-damper
rendering to couple the joints of two users. Exoskeleton-mediated physical
interaction aims to explore how task performance and motor learning effects,
initially assessed on simpler robots with just one joint, can be scaled to
address more complex, bilateral tasks [6]. Indeed, robot-mediated physical
interaction is acknowledged for its potential to facilitate knowledge transfer
between individuals while also enabling the recording of movements and the
repetition of exercises [13, 14].

In this context, we propose a StS framework for physical rehabilitation that
virtually connects two lower-limb exoskeletons operated by human users. To
promote repetitive training, we designed a task in which users control the
vertical position of their Center of Mass (CoM) with the help of continuous
visual feedback. In this work, we compared two haptic interaction strategies:
(1) virtual connections applied directly between the joints of the two
participants and (2) virtual connections linking the vertical position of the
two CoMs. The haptic interaction is represented as a spring-damper system, and
several compliance profiles are tested. The proposed framework can scale the
assistance provided to patients according to their functional level and uses a
balancing recovery strategy We present the structure of our framework and the
results of this preliminary testing.

##  II METHODS

###  II-A Exoskeleton Whole Body Model

For the StS task, we considered the exoskeleton in a state of double stance.
The equation of motion of the exoskeleton during this state is described as:

$\bm{M}_{\text{ds}}(\bm{q})\ddot{\bm{q}}+\bm{b}_{\text{ds}}(\bm{q},\dot{\bm{q}})+\bm{g}_{\text{ds}}(\bm{q},\bm{\alpha})=S^{T}(\bm{\tau}_{\text{motor}}+\bm{\tau}_{\text{friction}})+\bm{\tau}_{\text{int}}.$
(1)

where $\bm{q}=[\theta_{0},\theta_{1},\theta_{2},\theta_{3},\theta_{4}]$ are
the generalized coordinates corresponding to the backpack, hip, and knee
angles. The matrix and vectors $\bm{M}_{ds}\in\mathbb{R}^{5\times
5},\bm{b}_{ds}\in\mathbb{R}^{5},\bm{g}_{ds}\in\mathbb{R}^{5}$, are the mass
matrix, Coriolis and centrifugal torques, and gravitational torques
respectively during double stance. The variable
$\bm{\tau}_{\text{motor}}\in\mathbb{R}^{4}$ corresponds to the motor torques
sent to the driver, $\bm{\tau}_{\text{friction}}\in\mathbb{R}^{4}$ is the
friction component of the torque, $S=[0_{4\times 1},\mathbb{I}_{4\times 4}]$
is a selection matrix of actuated joints, and
$\bm{\tau}_{\text{int}}\in\mathbb{R}^{5}$ is a vector of interaction torques
applied to the exoskeleton by the user.

The parameter $\bm{\alpha}$ is used to quantify the amount of vertical
assistance the exoskeleton provides to the user expressed in $kg$. This
parameter is used to delay or avoid the onset of fatigue in StS exercise. In
fact, fatigue may impair motor learning during rehabilitation. This parameter
can be easily scaled during rehabilitation by increasing or decreasing
assistance depending on the patient’s level of impairment.

###  II-B Human-Exoskeleton Interaction

The model presented in Eq. (1) is used to control the joint acceleration such
that it follows a desired interaction profile. With this purpose, we
implemented the virtual mass controller presented in [11] to simulate a
desired virtual mass or inertia on the joints of the exoskeleton:

$\ddot{\bm{q}}^{*}=\bm{M}_{\text{virt}}^{-1}(\bm{\tau}_{\text{int}}-\bm{\tau}_{\text{int}}^{*}),$
(2)

where $\bm{\tau}_{\text{int}}^{*}$ is the desired interaction torque. For
instance, whenever we want to minimize the interaction between the exoskeleton
and the user, we apply zero desired interaction torque:

$\bm{\tau}_{\text{int}}^{*}=0.$ (3)

This condition is typically referred to as transparent control; in this paper,
it is used as a baseline for motor performance.

In many rehabilitation exercises, the exoskeleton is commanded to interact
with the user in a certain way (e.g., providing assistance along a trajectory
[15]). In our proposed framework of virtual interaction between two
exoskeletons ($A$, $B$), the desired interaction with the user is defined as a
virtual spring-damper between the configuration of the two exoskeletons. This
behavior can act directly at joint level by defining:

$\begin{split}\bm{\tau}^{\text{int}}_{A}&=K_{\text{virt,q}}(\theta_{B}-\theta_{A})+C_{\text{virt,q}}(\dot{\theta}_{B}-\dot{\theta}_{A})\\\
\bm{\tau}^{\text{int}}_{B}&=-\bm{\tau}^{\text{int}}_{A},\end{split}$ (4)

where $\theta_{A,B},\dot{\theta}_{A,B}$ are respectively the joints angle and
velocities of the robots $A$ and $B$. The variables $K_{\text{virt,q}}$ and
$C_{\text{virt,q}}$ are the constants of the rendered spring and damper
between the users, respectively.

In addition to the joint space interaction, it is also possible to design the
interaction in the task space. In other words, virtual elements can be
rendered between the two points of the exoskeletons in the sagittal plane. In
this work, we rendered spring and damper between the vertical positions ($z$)
of their CoM.

$\begin{split}\bm{F}^{\text{int}}_{A}&=K_{\text{virt,z}}(z_{B}-z_{A})+C_{\text{virt,z}}(\dot{z}_{B}-\dot{z}_{A})\\\
\bm{F}^{\text{int}}_{B}&=-\bm{F}^{\text{int}}_{A}\\\
\bm{\tau}^{\text{int}}_{A}&=\mathbb{J}_{\text{CoM}}^{z}(\bm{\theta}_{A})^{T}\bm{F}^{\text{int}}_{A}\\\
\bm{\tau}^{\text{int}}_{B}&=\mathbb{J}_{\text{CoM}}^{z}(\bm{\theta}_{B})^{T}\bm{F}^{\text{int}}_{B},\end{split}$
(5)

where $z_{A,B},\dot{z}_{A,B}$ are respectively the vertical positions and
velocities of the robots $A$ and $B$ and
$\mathbb{J}_{\text{CoM}}^{z}(.)\in\mathbb{R}^{1\times 5}$ is the vertical
component of CoM Jacobian. This type of interaction allows freedom of movement
in the null space of the vertical component of motion and a larger selection
of the postures that can meet the given position of the CoM.

###  II-C Quadratic Optimal Controller

To track the desired generalized acceleration presented in Eq. (2) under
physical and safety constraints, we use real-time constrained optimization
with the OSQP library [16] to solve for motor torques and joint accelerations,
$\textbf{x}=(\ddot{\bm{q}},\bm{\tau}_{\text{motor}})\in\mathbb{R}^{9}$:

$\min_{\textbf{x}}\lVert f(\textbf{x})\rVert\qquad\text{ s.t.
}\textbf{A}\textbf{x}=\textbf{b}\textit{ and
}\textbf{D}\textbf{x}\leq\textbf{f}.$ (6)

The objective function $\lVert f(\textbf{x})\rVert$ is the 2-norm of the
difference between the actual and desired generalized accelerations,

$f(\textbf{x})=[\mathbb{I}_{5\times 5},\;0_{5\times
4}]\textbf{x}-\ddot{\bm{q}}^{*}.$ (7)

The equality constraint ensures that optimized variables satisfy the equation
of motion (Eq. (1)),

$\begin{split}\textbf{A}=\textbf{A}_{\text{EOM}}&=[\bm{M}_{\text{ds}},\>-S^{\top}],\\\
\textbf{b}=\textbf{b}_{\text{EOM}}&=[-\bm{b}_{\text{ds}}-\bm{g}_{\text{ds}}+\bm{\tau}_{\text{int}}+S^{\top}\bm{\tau}_{\text{friction}}].\end{split}$
(8)

Inequality constraints arise from physical and safety constraints on motor
torque, power and velocity as presented in our previous work [11]. In this
paper we extend our previous formulation including a constraint on the joints
positions that avoids reaching the robot’s mechanical joint limits and
creating unexpected interactions. This is achieved by constraining the joint
accelerations such that the joint position limits will not be violated in
$\beta\Delta t$ amount of time, where $\Delta t$ is the control period.

$\begin{split}\textbf{D}_{{\bm{q}}^{+}}=[\mathbb{I}_{5\times 5},\>0_{5\times
4}],\;\textbf{f}_{{\bm{q}}^{+}}&=2\frac{{\bm{q}}_{\text{max}}-{\bm{q}}-\dot{\bm{q}}\beta\Delta
t}{\beta^{2}\Delta t^{2}}\\\ \textbf{D}_{{\bm{q}}^{-}}=[-\mathbb{I}_{5\times
5},\>0_{5\times
4}],\;\textbf{f}_{{\bm{q}}^{-}}&=-2\frac{{\bm{q}}_{\text{min}}-{\bm{q}}-\dot{\bm{q}}\beta\Delta
t}{\beta^{2}\Delta t^{2}},\end{split}$ (9)

We also added a safety constraint on the Cartesian position of the CoM that
ensures balancing in the sagittal plane which is presented in the following
subsection.

Figure 2: Balancing constraint for the StS task: Fig. (A) shows how the pose
(position and orientation) of the CoM is constrained in both the horizontal
and vertical planes; Fig. (B) shows how the position is constrained in the
forward movement to avoid falling; Fig. (C) shows how the backward movement is
constrained but still leaves the possibility of performing the StS task.

###  II-D Balancing Constraint

Human-exoskeleton balancing is a critical aspect of exoskeleton control to
provide seamless coordination between the wearer and the exoskeleton for
enhanced stability and mobility [17]. The main solutions proposed are based on
fall recognition and balance recovery and satisfaction [18]. Many of these
algorithms in the literature are based on strategies already used to control
humanoid robots [19, 20, 21]. In fact, balancing is usually satisfied by
constraining the position of the CoM within a support region defined by the
position of the feet and the activity being performed [22, 23, 24, 25]. Many
of the algorithms presented for applications in rehabilitation give little
freedom of motion to the patient and tend to prefer positional control
strategies [26]. This, on the one hand, allows greater safety for the wearer,
but on the other hand limits their contribution by enforcing passive
repetitions of a movement. In this scenario, Huo et al. [10] presented a
system that prioritizes the wearer’s control and delivers adequate assistance
for power and balance.

In this work, we employ a similar approach in order to satisfy the balancing
of the exoskeleton-human couple. We constrain the motion of the Divergent
Component of Motion (DCM) to remain inside the support region of the sagittal
plane of the exoskeleton. Since all the actuation of the exoskeleton is on the
sagittal plane, we implemented constraints only on that plane. Similarly to
[27] we define the DCM as

$\bm{p}_{\text{DCM}}=\bm{p}_{\text{CoM}}+(\sqrt{h/g})\dot{\bm{p}}_{\text{CoM}}$
(10)

where $\bm{p}_{\text{CoM}}=(x_{\text{CoM}},z_{\text{CoM}})$ is the position of
the CoM in the sagittal plane (Fig. 2A) while $h$ is the height of the linear
inverted pendulum which approximates the exoskeleton-human couple and is
assumed as the height of the CoM of the exoskeleton as the human is at an
upright posture. This aspect characterizes the pendulum’s susceptibility to
losing balance, with a larger $h$ indicating greater instability. In the case
of patients, this parameter may be increased, increasing the demand for
assistance, especially for those facing greater difficulties.

To constraint the position of the DCM we implemented the following inequality
constraint ($Dx\leq f$):

$\begin{split}\textbf{D}_{\text{CoM}^{+}}&=[\mathbb{J}_{\text{CoM}},\>0_{2\times
4}],\qquad\textbf{D}_{\text{CoM}^{-}}=-\textbf{D}_{\text{CoM}^{+}},\\\
\textbf{f}_{\text{DCM}^{+}}&=2\frac{\bm{p}^{+}-\bm{p}_{\text{DCM}}-\dot{\bm{p}}_{\text{DCM}}\beta\Delta
t}{\beta^{2}\Delta t^{2}}-\dot{\mathbb{J}}_{\text{CoM}}\dot{\bm{q}},\\\
\textbf{f}_{\text{DCM}^{-}}&=-2\frac{\bm{p}^{-}-\bm{p}_{\text{DCM}}-\dot{\bm{p}}_{\text{DCM}}\beta\Delta
t}{\beta^{2}\Delta t^{2}}+\dot{\mathbb{J}}_{\text{CoM}}\dot{\bm{q}},\\\
\end{split}$ (11)

where $\bm{p}^{+},\bm{p}^{-}$ are the maximum and minimum limits for the DCM
position, respectively. This formalization is similar to the joint position
limit expressed in Eq. (9). The joint accelerations are constrained such that
the CoM pose does not exceed the defined limits in $\beta\Delta t$ amount of
time.

Horizontal DCM position is defined as constant along the horizontal ($x$) axes
while the vertical component ($z$) is defined as a function of the horizontal
position. Similarly, to limit the over extension of the backpack angle we made
the upper limit of the backpack angle ($\bm{q}_{bp}^{+}$) variable from the
horizontal position. Specifically, the two constraints are defined as linear
functions that with a maximum at the maximum equilibrium position of the user
($x_{eq}$) and two minima in the limits of the horizontal position of the
sagittal plane. This is formalized as follows and visualized in Fig. 2A.

$\begin{split}\bm{p}_{z}^{+}&=z_{max}-||a_{z}(x-x_{eq})||,\\\
\bm{q}_{bp}^{+}&=\theta_{max}-||a_{\theta}(x-x_{eq})||,\\\ \end{split}$ (12)

We defined the constraint in that way to compensate for the robot’s
uncontrollable movements (non-actuated ankle joint). A more in-depth
discussion of the implications of such a constraint is covered in Sec. III.

Figure 3: Composition of trials with different visual feedback profiles: each trial consists of an equal time (30 s) of performing the motion alone, resting, and performing the motion while subjected to virtual interaction. Three types of reference were tested: (A) simple sine, (B) composite sine, and (C) discrete tasks in the range of motion. TABLE I: Values of the gains used in the experimentation. Parameter | Value | Parameter | Value
---|---|---|---
$\bm{K}_{\text{virt,q}}^{\text{soft}}$ | $\mathbb{I}_{4\times 4}\times 30[\frac{\text{Nm}}{\text{rad}}]$ | $\bm{C}_{\text{virt,q}}^{\text{soft}}$ | $\mathbb{I}_{4\times 4}\times 4[\frac{\text{Nms}}{\text{rad}}]$
$\bm{K}_{\text{virt,q}}^{\text{stiff}}$ | $\mathbb{I}_{4\times 4}\times 70[\frac{\text{Nm}}{\text{rad}}]$ | $\bm{C}_{\text{virt,q}}^{\text{stiff}}$ | $\mathbb{I}_{4\times 4}\times 10[\frac{\text{Nms}}{\text{rad}}]$
${K}_{\text{virt,z}}$ | 1000 $[\frac{\text{N}}{\text{m}}]$ | ${C}_{\text{virt,z}}$ | 50 $[\frac{\text{Ns}}{\text{m}}]$
$p^{+}$ | 0.2 [m] | $\bm{p}^{-}$ | 0.3 [m]

###  II-E Experimental Protocol

In this paper, we performed exploratory work to define a StS exercise using
exoskeletons to mediate patient-therapist interaction. Two users wore lower-
limb exoskeletons (X2, Fourier Intelligence, Singapore). The presented
controller was implemented on a ROS and C++ based software stack named CANOpen
Robot Controller (CORC) [28]. The gains used in the controller are presented
in Table I.

The distance between the chair and the ankle was adjusted according to the
participant’s preference, while the height of the chair was chosen according
to knee height. Dyads were seated side by side, facing two monitors which
displayed visual feedback for the tracking task (Fig. 1). To validate the
usability of the exoskeleton by therapists, users wore the exoskeleton without
external help. At this point, users familiarized themselves for two minutes
with the StS task wearing the exoskeleton. During this time, the range of
motion was recorded for visual reference and to define balancing constraints.

Figure 4: Tracking performance improved when two pilot users were coupled via
joint space connection compared to no connection. The left panel shows the
total error across repetitions of the StS exercise (mean $\pm$ SD) for two
users during solo and dyad trials under two connection stiffnesses. Right
panels show the error separated into bias and random across conditions.

Each trial consisted of the alternating succession of 3 periods lasting 30
seconds (Fig. 3): one period in which users performed the exercise alone
(transparent exoskeleton, $\bm{\tau}^{\text{int}}=0$), one period in which the
exoskeleton is used to virtually connect the two users (Eq. (4), (5)), and a
rest period. Each experimental condition is composed by 8 trials with 30
seconds of rest between each trial in which the two users performed the task
in parallel. During the tracking task both users were provided with visual
feedback of the target vertical CoM position (red point) and their own
instantaneous vertical CoM position (green dot) as shown in Fig. 1.
Participants were instructed to match their CoM position to the target
position as accurately as possible while the target moved along a vertical
bar. The upper and lower bounds of the target movement were determined by the
range of motion (ROM) of the two participants. Specifically, we used the
distance between the maximum of the minima and the minimum of maxima to define
a common ROM. We chose to use this convention to model the therapist-patient
case, in which the therapist is more likely to adapt to the patient’s joint
limits. Three types of references were tested with the common ROM as the
amplitude: (1) a simple sine function, (2) multi-sine function with different
frequencies, and (3) a step function with discrete values in the ROM,

$\begin{split}z_{\text{des}}=\text{ROM}\times[&\text{sin}(0.5(2\pi)(t+t_{r})+\phi_{1})],\\\
z_{\text{des}}=\text{ROM}\times[&0.33\text{sin}(0.50(2\pi)t)\\\
+&0.33\text{sin}(0.20(2\pi)t)\\\ +&0.33\text{sin}(0.16(2\pi)t)],\\\
\end{split}$ (13)

where $z_{\text{des}}$ is the target vertical CoM position at a given time.

In this preliminary study, we compared different configurations to define the
usability of the StS infrastructure. With this intent, we tested different
interaction profiles (i.e., stiffness and damping) as shown in Tab. I. We
present the results and discuss measures of performance which are relevant for
a StS exercise. In particular, we are interested in a user’s tracking errors
while they attempt to accurately follow the visually-displayed target during
the exercise. We use the root-mean-square of the difference between each
individual’s actual trajectory and the target trajectory to represent the
total error during each tracking trial. In addition, we distinguish between
the type of error, bias (i.e., systematic) and random, as they each contribute
to the total tracking error. Bias error gives an indication of repeated
mistakes by the user (e.g., undershooting the peak target position) across
tracking trials. Random error indicates the consistency of the movement across
trials, as low random error would suggest the user smoothly and repeatably
pursues the target from trial to trial. We present these three categories of
error (i.e., total, bias and random) for each of the stiffness configurations
during the joint space interaction to show how haptic coupling influences
tracking performance for an example dyad.

##  III RESULTS

Comparing performances with and without haptic interaction during the tracking
task, our preliminary results are consistent with previous work in dyadic
tracking at the ankle [13]. In fact, not only does the interaction seem to
improve one’s ability to follow a reference, but also the degree of
improvement appears related to the magnitude of virtual stiffness of the
connection. It is not the intent of this paper to validate these findings for
lower-limb exoskeletons but we report here a change in user performance for
one dyad (Fig. 4).

In these initial observations, we did not observe any particular differences
between applying a virtual connection in the task space or joint space. This
is likely due to the small variability with which StS movements are performed
at knee and hip levels. Fig. 5C shows the limited amount of movement in the
null space for a single joint (right knee) both during the solo and dyad
conditions. We expect that greater differences may be observed in users with
differences in anatomy, strength and flexibility (e.g., therapist and
patient). For instance, it is possible that interaction in the task space
could reduce the quality of information the therapist can receive from the
patient (i.e., lack of feedback in each joint) [29]. In future experiments, we
intend to test higher stiffness values and compare users with different
kinematic and dynamic properties.

In Sec.II-D, we presented the balancing constraint we implemented to
compensate for the robot’s movement limitations (i.e., non-actuated ankle
joint). In fact, ankle movements can lead to solutions that violate the
constraint and thus to instability and unnatural movements for the user such
as hyper-extension of the hip. The constraints on the vertical position and
backpack angle are designed specifically to avoid these occurrences. For
example, a forward movement of the ankle that violates the horizontal
constraint generates solutions that lead the user to a safe position by
lowering the CoM position and backpack angle. Backward movements, on the other
hand, are less constrained to allow execution of the StS movement. Constraints
on joint velocities and the backpack remain active, which limits dangerous
movements. The results on the horizontal position of the CoM can be seen in
Fig. 2B,C.

The influence of vertical assistance from the robot was also evaluated in this
study. Its contribution in guiding vertical movement and aiding descent was
assessed through subjective feedback and comparison of CoM trajectories. Users
reported difficulty in distinguishing the contribution of the assistance from
that of the virtual connection. This is evident when comparing Fig. 5A and
Fig. 5B: with assistance, in both solo and dyadic conditions, the users
experienced a vertical force exerted by the exoskeleton, prompting them to
adapt their approach. They responded by quickening their upward movements and
postponing their downward movements. Based on these results, we believe that
assistance from the robot is critical to limit the onset of fatigue, but it is
better suited for discrete tracking tasks (Fig.3C) in which users have more
degree of freedom on the choice of the trajectory to execute.

Figure 5: Distribution of Sit-to-Stand motion: (A) task space ($z_{CoM}$)
movement distribution under different interaction conditions; (B) task space
($z_{CoM}$) movement distribution while both exoskeletons produce an upward
assistance component; (C) joint space ($\bm{q}$) movement distribution of the
right knee angle.

In this study we present three different types of tracking tasks for a StS
exercise (Fig. 3). Each has advantages: sinusoidal functions (Fig. 3A) were
reported to be less-fatiguing because they allow users to pause in comfort
positions (i.e., standing or sitting). However, the predictive nature of a
single frequency sine function is likely too simple to analyze motor learning.
Complex sinusoidal functions (Fig. 3B), on the other hand, promote variability
of movement, which is fundamental to rehabilitative exercises. The drawback of
multi-sine functions is that they are particularly fatiguing, requiring users
to slow down and hold flexed postures to follow the speed of the reference.
The discrete sinusoidal trajectory (Fig. 3C and Fig. 6) reduces these
limitations, but is more difficult to evaluate from a performance standpoint
for the variability with which users can choose to change their approach
trajectory. Based on these results, we intend to perform additional pilot
testing with more users to construct a rehabilitation task that appropriately
challenges users while limiting fatigue over many repetitions.

Figure 6: Example of a StS tracking task with discrete targets.

##  IV Conclusion

In this paper, we presented a framework for rehabilitation exercises specific
to Sit-to-Stand movements based on the concept of exoskeleton mediated
interaction. To share the task between the two participants, we proposed three
types of tracking tasks. These tasks are aimed at promoting repeatability and
at the same time movement variability, which are essential to activate motor
learning in patients. We then extended the controller for exoskeleton mediated
haptic interaction by including scalable and task-specific Sit-to-Stand
assistance and support to ensure balance during movement execution. In the
future, we intend to test this approach first on a healthy population to
validate the learning ability of the Sit-to-Stand task and then on a
population of therapists and patients to test its usability and contribution
in the context of physical rehabilitation.

## Acknowledgments

This work was supported by the National Science Foundation / National Robotics
Initiative (Grant No: 2024488). We would like to thank Tim Haswell for his
technical support on the hardware improvements of the ExoMotus-X2 exoskeleton.

## References

  * [1] A. Boukadida, F. Piotte, P. Dehail, and S. Nadeau, “Determinants of sit-to-stand tasks in individuals with hemiparesis post stroke: a review,” _Annals of physical and rehabilitation medicine_ , vol. 58, no. 3, pp. 167–172, 2015.
  * [2] P.-T. Cheng, M.-Y. Liaw, M.-K. Wong, F.-T. Tang, M.-Y. Lee, and P.-S. Lin, “The sit-to-stand movement in stroke patients and its correlation with falling,” _Archives of physical medicine and rehabilitation_ , vol. 79, no. 9, pp. 1043–1046, 1998.
  * [3] B. McCrory, J. M. Burnfield, A. R. Darragh, J. L. Meza, S. L. Irons, P. Chernyavskiy, A. M. Link, and G. Brusola, “Work injuries among therapists in physical rehabilitation,” in _Proceedings of the Human Factors and Ergonomics Society Annual Meeting_ , vol. 58, no. 1. SAGE Publications Sage CA: Los Angeles, CA, 2014, pp. 1072–1076.
  * [4] J. M. Burnfield, B. McCrory, Y. Shu, T. W. Buster, A. P. Taylor, and A. J. Goldman, “Comparative kinematic and electromyographic assessment of clinician-and device-assisted sit-to-stand transfers in patients with stroke,” _Physical Therapy_ , vol. 93, no. 10, pp. 1331–1341, 2013.
  * [5] M. Fray, S. Hignett, A. Reece, S. Ali, and L. Ingram, “An evaluation of sit to stand devices for use in rehabilitation,” in _Proceedings of the 20th Congress of the International Ergonomics Association (IEA 2018) Volume VII: Ergonomics in Design, Design for All, Activity Theories for Work Analysis and Design, Affective Design 20_. Springer, 2019, pp. 774–783.
  * [6] E. B. Küçüktabak, S. J. Kim, Y. Wen, K. Lynch, and J. L. Pons, “Human-machine-human interaction in motor control and rehabilitation: a review,” _Journal of NeuroEngineering and Rehabilitation_ , 2021. [Online]. Available: https://doi.org/10.1186/s12984-021-00974-5
  * [7] C. Camardella, F. Porcini, A. Filippeschi, S. Marcheschi, M. Solazzi, and A. Frisoli, “Gait phases blended control for enhancing transparency on lower-limb exoskeletons,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 3, pp. 5453–5460, 2021.
  * [8] M. Sharifi, J. K. Mehr, V. K. Mushahwar, and M. Tavakoli, “Adaptive cpg-based gait planning with learning-based torque estimation and control for exoskeletons,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 4, pp. 8261–8268, 2021.
  * [9] R. M. Andrade, S. Sapienza, E. E. Fabara, and P. Bonato, “Trajectory tracking impedance controller in 6-dof lower-limb exoskeleton for over-ground walking training: Preliminary results,” in _2021 international symposium on medical robotics (ISMR)_. IEEE, 2021, pp. 1–6.
  * [10] W. Huo, H. Moon, M. A. Alouane, V. Bonnet, J. Huang, Y. Amirat, R. Vaidyanathan, and S. Mohammed, “Impedance modulation control of a lower-limb exoskeleton to assist sit-to-stand movements,” _IEEE Transactions on Robotics_ , vol. 38, no. 2, pp. 1230–1249, 2021.
  * [11] E. B. Küçüktabak, Y. Wen, S. J. Kim, M. Short, D. Ludvig, L. Hargrove, E. Perreault, K. Lynch, and J. Pons, “Haptic transparency and interaction force control for a lower-limb exoskeleton,” _arXiv preprint arXiv:2301.06244_ , 2023.
  * [12] E. B. Küçüktabak, Y. Wen, M. Short, E. Demirbaş, K. Lynch, and J. Pons, “Virtual physical coupling of two lower-limb exoskeletons,” _arXiv preprint arXiv:2307.06479_ , 2023.
  * [13] M. R. Short, D. Ludvig, E. B. Küçüktabak, Y. Wen, L. Vianello, E. J. Perreault, L. Hargrove, K. Lynch, and J. Pons, “Haptic Human-Human Interaction During an Ankle Tracking Task: Effects of Virtual Connection Stiffness,” 5 2023. [Online]. Available: https://www.techrxiv.org/articles/preprint/Haptic_Human-Human_Interaction_During_an_Ankle_Tracking_Task_Effects_of_Virtual_Connection_Stiffness/23234069
  * [14] S. J. Kim, Y. Wen, D. Ludvig, E. B. Küçüktabak, M. R. Short, K. Lynch, L. Hargrove, E. J. Perreault, and J. L. Pons, “Effect of dyadic haptic collaboration on ankle motor learning and task performance,” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , 2022\.
  * [15] R. L. Medrano, G. C. Thomas, C. G. Keais, E. J. Rouse, and R. D. Gregg, “Real-time gait phase and task estimation for controlling a powered ankle exoskeleton on extremely uneven terrain,” _IEEE Transactions on Robotics_ , 2023.
  * [16] B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd, “OSQP: an operator splitting solver for quadratic programs,” _Mathematical Programming Computation_ , vol. 12, no. 4, pp. 637–672, 2020. [Online]. Available: https://doi.org/10.1007/s12532-020-00179-2
  * [17] M. F. Hamza, R. A. R. Ghazilla, B. B. Muhammad, and H. J. Yap, “Balance and stability issues in lower extremity exoskeletons: A systematic review,” _Biocybernetics and Biomedical Engineering_ , vol. 40, no. 4, pp. 1666–1679, 2020.
  * [18] I. Farkhatdinov, J. Ebert, G. Van Oort, M. Vlutters, E. Van Asseldonk, and E. Burdet, “Assisting human balance in standing with a robotic exoskeleton,” _IEEE Robotics and automation letters_ , vol. 4, no. 2, pp. 414–421, 2019.
  * [19] J. Ramos and S. Kim, “Dynamic locomotion synchronization of bipedal robot and human operator via bilateral feedback teleoperation,” _Science Robotics_ , vol. 4, no. 35, p. eaav4282, 2019.
  * [20] L. Penco, B. Clément, V. Modugno, E. M. Hoffman, G. Nava, D. Pucci, N. G. Tsagarakis, J.-B. Mouret, and S. Ivaldi, “Robust real-time whole-body motion retargeting from human to humanoid,” in _2018 IEEE-RAS 18th International Conference on Humanoid Robots (Humanoids)_. IEEE, 2018, pp. 425–432.
  * [21] L. Vianello, L. Penco, W. Gomes, Y. You, S. M. Anzalone, P. Maurice, V. Thomas, and S. Ivaldi, “Human-humanoid interaction and cooperation: a review,” _Current Robotics Reports_ , pp. 1–14, 2021.
  * [22] O. Harib, A. Hereid, A. Agrawal, T. Gurriet, S. Finet, G. Boeris, A. Duburcq, M. E. Mungai, M. Masselin, A. D. Ames _et al._ , “Feedback control of an exoskeleton for paraplegics: Toward robustly stable, hands-free dynamic walking,” _IEEE Control Systems Magazine_ , vol. 38, no. 6, pp. 61–87, 2018\.
  * [23] F. Xu, X. Lin, H. Cheng, R. Huang, and Q. Chen, “Adaptive stair-ascending and stair-descending strategies for powered lower limb exoskeleton,” in _2017 IEEE International Conference on Mechatronics and Automation (ICMA)_. IEEE, 2017, pp. 1579–1584.
  * [24] C. Deng and Z. Li, “Human-guided robotic exoskeleton cooperative walking for climbing stairs,” in _2018 3rd International Conference on Advanced Robotics and Mechatronics (ICARM)_. IEEE, 2018, pp. 60–65.
  * [25] G. Menga and M. Ghirardi, “Lower limb exoskeleton for rehabilitation with improved postural equilibrium,” _Robotics_ , vol. 7, no. 2, p. 28, 2018.
  * [26] G. Rodriguez Tapia, I. Doumas, T. Lejeune, and J.-G. Previnaire, “Wearable powered exoskeletons for gait training in tetraplegia: a systematic review on feasibility, safety and potential health benefits,” _Acta Neurologica Belgica_ , vol. 122, no. 5, pp. 1149–1162, 2022.
  * [27] Y. Ishiguro, T. Makabe, Y. Nagamatsu, Y. Kojio, K. Kojima, F. Sugai, Y. Kakiuchi, K. Okada, and M. Inaba, “Bilateral humanoid teleoperation system using whole-body exoskeleton cockpit tablis,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 6419–6426, 2020.
  * [28] J. Fong, E. B. Küçüktabak, V. Crocher, Y. Tan, K. Lynch, J. Pons, and D. Oetomo, “CANopen Robot Controller (CORC): An open software stack for human robot interaction development,” _The International Symposium on Wearable Robotics (WeRob2020)_ , vol. 2020, 2020.
  * [29] E. B. Brokaw, R. J. Holley, and P. S. Lum, “Comparison of joint space and end point space robotic training modalities for rehabilitation of interjoint coordination in individuals with moderate to severe impairment from chronic stroke,” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , vol. 21, no. 5, pp. 787–795, 2013.