# Matchings in the hypercube with specified edges

Joshua Erde Institute of Discrete Mathematics, Graz University of Technology,
Steyergasse 30, 8010 Graz, Austria. Email<EMAIL_ADDRESS>
###### Abstract

Given a matching $M$ in the hypercube $Q^{n}$, the _profile_ of $M$ is the
vector $\bm{x}=(x_{1},\ldots,x_{n})\in\mathbb{N}^{n}$ such that $M$ contains
$x_{i}$ edges whose endpoints differ in the $i$th coordinate. If $M$ is a
perfect matching, then it is clear that $||\bm{x}||_{1}=2^{n-1}$ and it is
easy to show that each $x_{i}$ must be even. Verifying a special case of a
conjecture of Balister, Győri, and Schelp, we show that these conditions are
also sufficient.

##  1 Introduction

Balister, Győri, and Schelp [1] considered the following partitioning problem:
we are given a sequence $\bm{d_{1}},\ldots,\bm{d_{2^{n-1}}}$ of elements of
$\mathbb{F}_{2}^{n}$ and our aim is to partition $\mathbb{F}_{2}^{n}$ into
pairs $\\{\bm{a_{i}},\bm{b_{i}}\\}$ such that
$\bm{a_{i}}-\bm{b_{i}}=\bm{d_{i}}$ for each $i\in\left[2^{n-1}\right]$. They
noted that, since we are working over $\mathbb{F}_{2}^{n}$, for such a
partition

$\sum_{i=1}^{2^{n-1}}\bm{d_{i}}=\sum_{i=1}^{2^{n-1}}\bm{a_{i}}-\bm{b_{i}}=\sum_{i=1}^{2^{n-1}}\bm{a_{i}}+\bm{b_{i}}=\sum_{\bm{v}\in\mathbb{F}_{2}^{n}}\bm{v}=\bm{0}.$
(1)

They conjectured that this necessary condition was also sufficient.

######  Conjecture 1.1 (Balister, Győri, and Schelp).

Let $n\geq 2$ and let $\bm{d_{1}},\ldots,\bm{d_{2^{n-1}}}$ be elements of
$\mathbb{F}_{2}^{n}$ such that $\sum_{i=1}^{2^{n-1}}\bm{d_{i}}=\bm{0}$. Then
there exists a partition of $\mathbb{F}_{2}^{n}$ into pairs
$\\{\bm{a_{i}},\bm{b_{i}}\\}$ such that $\bm{a_{i}}-\bm{b_{i}}=\bm{d_{i}}$ for
each $i\in\left[2^{n-1}\right]$.

Similar problems of partitioning groups into pairs with prescribed differences
have also been considered in $\mathbb{F}_{p}$ for $p>2$ prime [4, 7],
$\mathbb{F}_{p}^{n}$ [3] and cyclic groups [5].

Balister, Győri, and Schelp [1] showed that Conjecture 1.1 holds under certain
assumptions on the differences $\bm{d_{i}}$ \- namely if the first half of the
differences are identical, and the second half come in identical pairs. More
recently, Kovács [6] gave further evidence towards Conjecture 1.1.

######  Theorem 1.2 ([6]).

If either

  1. (a)

the number of distinct $\bm{d}_{i}$ is at most $n-2\log_{2}n-1$; or

  2. (b)

$n$ is sufficiently large and at least a $\frac{28}{29}$ proportion of the
$\bm{d}_{i}$ are identical,

then Conjecture 1.1 holds

In this note we consider another natural restriction of Conjecture 1.1, where
we insist that each difference $\bm{d_{i}}$ must be some unit vector, i.e., of
the form $\bm{e_{i}}$ for some $i$. We note that, since there are precisely
$n$ unit vectors, this is not quite covered by Theorem 1.2 (a). In this case,
Conjecture 1.1 has a simple reformulation in terms of matchings in the
hypercube $Q^{n}$.

Given a matching $M\subseteq Q^{n}$ we say the _profile_ of $M$ is the vector
$\bm{x}=(x_{1},\ldots,x_{n})\in\mathbb{N}^{n}$ such that $M$ contains $x_{i}$
edges whose endpoints differ in the $i$th coordinate, and we are interested in
which tuples $\bm{x}\in\mathbb{N}^{n}$ are achievable as the profile of some
matching. We call such tuples _admissable_. In this case, the obvious
necessary condition for a tuple to be the profile of a perfect matching is
that each unit vector appears an even number of times as a matching edge,
i.e., each $x_{i}$ is even, in which case we say $\bm{x}$ is _even_ , and
Conjecture 1.1 would imply that this condition is also sufficient. The aim of
this short note is to verify the conjecture in this case. Let us call a tuple
$\bm{x}\in\mathbb{N}^{n}$ _perfect_ if
$||\bm{x}||_{1}:=\sum_{i}x_{i}=2^{n-1}$.

######  Theorem 1.3.

Let $n\geq 2$. A perfect tuple $\bm{x}\in\mathbb{N}^{n}$ is admissable if and
only if it is even.

After distributing an initial version of the manuscript it was bought to the
author’s attention, that Conjecture 1.1 and related questions have also been
considered independently in the context of _batch codes_ in information
theory, and that Theorem 1.3 also follows from work of Hollmann, Khathuria,
Riet, and Skachek [2], which implies Conjecture 3.1 in the case where all the
vectors $\bm{d_{i}}$ lie in some fixed hyperplane $H$. However, since we feel
the reformulation in terms of matchings in $Q^{n}$ leads to interesting
questions (see Section 3), and our proof is short and simple, we hope the work
in the note may still be of interest.

The structure of the paper is as follows. In Section 2 we prove Theorem 1.3
and in Section 3 we discuss some related questions.

##  2 Main result

We first note that by symmetry, the ordering of the coordinates of $\bm{x}$ is
irrelevant, and so we may assume from this point forwards that the coordinates
appear in non-decreasing order, i.e., $x_{1}\leq x_{2}\leq\ldots\leq x_{n}$.
Given $\bm{x},\bm{z}\in\mathbb{N}^{n}$ we say that _$\bm{z}$ precedes
$\bm{x}$_, which we write as $\bm{z}\preceq\bm{x}$ if $z_{i}\leq x_{i}$ for
all $i\in[n]$. Clearly the set of admissable tuples is down-closed in the
partial order $\preceq$. In particular, Theorem 1.3 is equivalent to the
statement that every even tuple is admissable.

The following simple observation underpins our proof.

######  Lemma 2.1.

If $\bm{x}\in\mathbb{N}^{n}$ is admissable and
$x_{n+1}=2^{n}-2||\bm{x}||_{1}$, then $(2\bm{x},x_{n+1})$ is a perfect
admissable tuple.

###### Proof.

By assumption, there is a matching $M_{1}\subseteq Q^{n}$ with profile
$\bm{x}$. We form a matching $M_{2}\subseteq Q^{n+1}$ by taking disjoint two
copies of $M_{1}$, that is

$M_{2}:=\\{\big{(}(\bm{u},0),(\bm{v},0)\big{)},\big{(}(\bm{u},1),(\bm{v},1)\big{)}\colon(\bm{u},\bm{v})\in
M_{1}\\}$

and note that $M_{2}$ has profile $2\bm{x}$.

Let $X\subseteq V\left(Q^{n}\right)$ be the set of vertices which are not
covered by $M_{1}$, where we note that $x_{n+1}=:|X|=2^{n}-2||\bm{x}||_{1}$.
We form $M_{3}$ by adding to $M_{2}$ the edges in direction $\bm{e_{n+1}}$
joining the two copies of each vertex in $X$ in $Q^{n+1}$. That is,

$M_{3}:=M_{2}\cup\\{\big{(}(\bm{v},0),(\bm{v},1)\big{)}\colon\bm{v}\in X\\}.$

It is clear that $M_{3}\subseteq Q^{n}$ is a matching with profile
$(2\bm{x},x_{n+1})$ and since $2||\bm{x}||_{1}+|X|=2^{n}$ it follows that
$(2\bm{x},x_{n+1})$ is a perfect admissable tuple. ∎

$\xrightarrow{\hskip 28.45274pt}$

Figure 1: Example showing how to build matching with profile $(2,2,2,2)$ from
a matching with profile $(1,1,1)$

Lemma 2.1 lends itself well to an inductive argument - given a perfect even
tuple $\bm{x}\in\mathbb{N}^{n+1}$, if $x_{i}\equiv 0\mod 4$ for all $i\in[n]$,
then
$\bm{y}:=\left(\frac{x_{1}}{2},\frac{x_{2}}{2},\ldots,\frac{x_{n}}{2}\right)$
is an even tuple in $\mathbb{N}^{n}$ and we can deduce the admissability of
$\bm{x}$ from that of $\bm{y}$. However, it is not clear what to do if some of
the coordinates of $\bm{y}$ are odd.

The key idea here is that, if $x_{n+1}$ is sufficiently large, then $\bm{y}$
will not be perfect, and so we might hope to find an even tuple
$\bm{y^{\prime}}$ which $\preceq$-dominates $\bm{y}$ with
$||\bm{y^{\prime}}||_{1}\leq 2^{n-1}$. The admissability of $\bm{y^{\prime}}$
would then imply the admissability of $\bm{y}$ and thus the admissability of
$\bm{x}$.

So, for example, to show that $(0,2,4,4,6)$ is admissable, by Lemma 2.2 it is
sufficient to show that $(0,1,2,2)$ is admissable, however
$(0,1,2,2)\preceq(2,2,2,2)$ and we know (cf. Figure 1) that $(2,2,2,2)$ is
admissable.

Let us make the preceding idea formal with the following lemma. Given
$\bm{x}\in\mathbb{N}^{n}$, let us write $o(\bm{x})$ for the number of
coordinates of $\bm{x}$ except the last which are congruent to $2\mod 4$.

######  Lemma 2.2.

Let $n\geq 2$, and let $\bm{x}\in\mathbb{N}^{n+1}$ be even with
$||\bm{x}||_{1}\leq 2^{n}$ and $x_{n+1}\geq 2o(\bm{x})$. If
$\bm{y}:=\left(\frac{x_{1}}{2},\frac{x_{2}}{2},\ldots,\frac{x_{n}}{2}\right)$,
then there exists a perfect even $\bm{y^{\prime}}\in\mathbb{N}^{n}$ such that
$\bm{y}\preceq\bm{y^{\prime}}$

###### Proof.

For each $i\in[n]$ let $o_{i}=1$ if $x_{i}\equiv 2\mod 4$ and let $o_{i}=0$
otherwise. Let $\bm{z}=\bm{y}+\bm{o}$. Note in particular that
$\bm{y}\preceq\bm{z}$ and that by construction $\bm{z}$ is even.

Since $||\bm{o}||_{1}=o(\bm{x})\leq\frac{x_{n+1}}{2}$, it follows that

$||\bm{z}||\leq\sum_{i=1}^{n}\frac{x_{i}}{2}+o_{i}\leq\frac{||\bm{x}||_{1}}{2}\leq
2^{n-1}.$

Hence, since $2^{n-1}$ is even, there is some perfect even
$\bm{y^{\prime}}\succeq\bm{z}\succeq\bm{y}$ as claimed. ∎

###### Proof of Theorem 1.3.

We first note that by (1) every perfect admissable tuple is even.

For the other direction, let us suppose towards a contradiction that the
theorem does not hold, and let $\bm{x}\in\mathbb{N}^{n+1}$ be a perfect even
tuple which is not admissable, with $n$ as small as possible. Let
$\bm{y}=\left(\frac{x_{1}}{2},\frac{x_{2}}{2},\ldots,\frac{x_{n}}{2}\right)$.
Note that, $\bm{x}=(2\bm{y},x_{n+1})$ and since $\bm{x}$ is perfect,
$x_{n+1}=2^{n}-2||\bm{y}||_{1}$. Furthremore, since $\bm{x}$ is not
admissable, by Lemma 2.1 neither is $\bm{y}$.

If $\bm{x}$ is such that $x_{n+1}\geq 2n$, then since $o(\bm{x})\leq n$, it
follows from Lemma 2.2 that there exists a perfect even
$\bm{y^{\prime}}\succeq\bm{y}$. However, since $\bm{y}$ is not admissable,
neither is $\bm{y^{\prime}}\in\mathbb{N}^{n}$, contradicting our choice of
$\bm{x}$.

Hence we way assume that $x_{n+1}<2n$. However, if $n+1\geq 8$, then

$x_{n+1}\geq\frac{2^{n}}{n+1}\geq 2n,$ (2)

and hence we may assume that $n+1\leq 7$.

This reduces the proof to a finite case check, and in fact more precise
applications of Lemma 2.2 will deal with all but three of the remaining cases.
In the table below we list all perfect even tuples $\bm{x}\in\mathbb{N}^{n+1}$
with $3\leq n+1\leq 7$ such that $x_{n+1}<2n$. For tuples listed in blue,
$x_{n+1}\geq 2o(\bm{x})$, and so by Lemma 2.2 these tuples cannot be minimal
counterexamples. The remaining tuples are listed in red.

$n+1=3$ | $(0,2,2)$
---|---
$n+1=4$ |  $(0,0,4,4)$, $(0,2,2,4)$, $(2,2,2,2)$
$n+1=5$ |  $(0,0,4,6,6)$, $(0,2,2,6,6)$, $(0,2,4,4,6)$,
|  $(2,2,2,4,6)$, $(0,4,4,4,4)$, $(2,2,4,4,4)$
$n+1=6$ |  $(0,0,8,8,8,8)$,$(0,2,6,8,8,8)$,$(0,4,4,8,8,8)$,$(2,2,4,8,8,8)$,
|  $(0,4,6,6,8,8)$, $(2,2,6,6,8,8)$,$(2,4,4,6,8,8)$,
|  $(4,4,4,4,8,8)$, $(0,6,6,6,6,8)$, $(2,4,6,6,6,8)$,
|  $(4,4,4,6,6,8)$, $(2,6,6,6,6,6)$, $(4,4,6,6,6,6)$
$n+1=7$ |  $(4,10,10,10,10,10,10)\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$,$(6,8,10,10,10,10,10)$,$(8,8,8,10,10,10,10)$

So, the only three cases we have to construct by hand are $(0,2)$, $(2,2,2,2)$
and $(6,6,6,6,6,2)$. At this point it is slightly easier to apply Lemma 2.1 to
see that is remains to show that $(1)$, $(1,1,1)$ and $(3,3,3,3,1)$ are
admissable, where in fact we show the stronger statement that $(3,3,3,3,3)$ is
admissable.

The first is trivial. The second is only slightly less trivial and is given
already in Figure 1. The third is a fun exercise and one possible solution is
given in Figure 2.

Figure 2: A matching in $Q^{5}$ with profile $(3,3,3,3,3)$.

∎

##  3 Discussion

As we noted in Section 2, if $\bm{x}$ is an admissable tuple and
$\bm{y}\preceq\bm{x}$, then $\bm{y}$ is also admissable. In particular,
Theorem 1.3 also tells us a lot about the admissable profiles of non-perfect
matchings.

However, whilst Theorem 1.3 fully characterises the _perfect_ admissable
tuples, there are admissable tuples which are not $\preceq$-dominated by any
perfect admissable tuples, for example two of the base cases in the proof of
Theorem 1.3. Apart from the parity restrictions on perfect admissable tuples,
there doesn’t seem to be any restriction, and we have checked by hand that all
tuples $\bm{x}\in\mathbb{N}^{n}$ with $||\bm{x}||_{1}<2^{n-1}$ are admissable
for $n\leq 4$.

######  Conjecture 3.1.

Let $n\in\mathbb{N}$. Every tuple $\bm{x}\in\mathbb{N}^{n}$ with
$||\bm{x}||_{1}<2^{n-1}$ is admissable.

It would also be interesting to characterise which tuples are the profiles of
Hamilton cycles in $Q^{d}$. Since each Hamilton cycle can be split into two
matchings, it is clear that any such tuple must be the sum of two perfect
admissable tuples, and so in particular must also be even. Furthermore,
clearly no direction can appear more than $2^{n-1}$ times and each direction
much appear at least twice. We conjecture that these are the only
restrictions.

######  Conjecture 3.2.

Let $n\in\mathbb{N}$. Every even tuple $\bm{x}\in\mathbb{N}^{n}$ with
$||\bm{x}||_{1}=2^{n}$, $\max_{i}\\{x_{i}\\}\leq 2^{n-1}$ and
$\min_{i}\\{x_{i}\\}\geq 2$ is the profile of a Hamilton cycle.

A different generalisation would be to consider decomposing the hypercube into
$2$-dimensional faces with specified orientations. More concretely, given a
decomposition $D$ of $E\left(Q^{n}\right)$ into $4$-cycles, the _profile_ of
this decomposition is the weighting $w$ of $K_{n}$ where the weight $w(ij)$ of
an edge is the number of $4$-cycles in $D$ where the antipodal points differ
in coordinates $i$ and $j$. It is easy to see that (1) implies that for any
$i$ the total number of faces with an edge in direction $i$ must be even,
i.e., $\sum_{j\neq i}w(ij)=0\mod 2$. Indeed, if we consider the intersection
of such a decomposition with the subcube whose $i$th coordinate is $0$ then we
can easily construct a perfect matching with $\sum_{j\neq i}w(ij)$ edges in
direction $i$. However, it is a simple exericse to check that the weighting
$w$ of $K_{4}$ given by $w(12)=w(34)=2$ and $w(ij)=0$ otherwise is not the
profile of any such decomposition, and so at least in small dimensions there
are other obstructions.

######  Question 1.

What are the admissable profiles for decompositions of $Q^{n}$ into
$4$-cycles?

Another interesting question would be to characterise the admissable profiles
for matchings of the _middle layer graph_ , the induced subgraph $M_{n}$ of
$Q^{2n+1}$ on the middle two layers. A similar argument as in (1) shows that
if $\bm{x}$ is the profile of a perfect matching on $M_{n}$ then
$x_{i}\equiv\binom{2n+1}{n+1}\mod 2$ for each $i\in[2n+1]$. An application of
Lucas’ theorem shows that these binomial coefficients are odd if and only if
$n+1$ is a power of $2$.

However, there is another necessary condition beyond this parity condition and
the trivial condition that $||\bm{x}||_{1}=\frac{|V(M_{n})|}{2}$. Since only
$\binom{2n}{n}$ vertices in each layer are incident to an edge in a fixed
direction, we require that $x_{i}\leq\binom{2n}{n}$ for each $i\in[2n+1]$. In
fact, since the $\binom{2n}{n-1}$ vertices in the lower layer not adjacent to
an edge in direction $\bm{e_{i}}$ are only adjacent to the vertices in the
upper layer which are adjacent to an edge in direction $\bm{e_{i}}$, and
therefore must be matched to a subset of these vertices in a direction not
equal to $\bm{e_{i}}$, we can have at most
$\binom{2n}{n}-\binom{2n}{n-1}=\binom{2n-1}{n-1}$ edges in direction
$\bm{e_{i}}$.

For $n=1$, it is easy to verify that these two conditions, together with the
trivial condition that
$||\bm{x}||_{1}=\frac{|V(M_{n})|}{2}=\binom{2n+1}{n+1}$, are also sufficient,
since $(1,1,1)$ is the only tuple satisfying these conditions.

######  Question 2.

What are the admissable profiles of perfect matchings of $M_{n}$?

Finally, another intriguing generalisation would be to look at matchings on
other polytopes. For example, each of the edges in the $n$-dimensional
permutahedron Perm$(n)$, embedded in $\mathbb{R}^{n+1}$ as the convex hull of
all _permutation vectors_ , is parallel to some vector of the form
$\bm{e_{i}}-\bm{e_{j}}$ with $i,j\in[n+1]$. Let us define the _profile_ of a
matching of Perm$(n)$ to be the weighting $w$ of $K_{n+1}$ where $w(ij)$ is
the number of edges parallel to the vector $\bm{e_{i}}-\bm{e_{j}}$.

######  Question 3.

What are the admissable profiles of perfect matchings of Perm$(n)$?

In the case $n=2$ there are only $2$ matchings, each of which use an edge in
each direction once. For $n=3,4$ computational evidence suggests that each
direction must be used an even number of times, and that if view the profile
$w$ as a vector in $\mathbb{R}^{\binom{n+1}{2}}$ then the set of admissable
vectors are the even interior points of some convex polytope. For $n=3$, this
polytope is determined by the inequalities

$\displaystyle w(12)+w(34)=w(13)+w(24)=w(14)+w(23)=4;$ $\displaystyle
w(ij)\geq 0\text{ for all }i,j,$

which essentially say that each pair of orthogonal directions has to be used
exactly four times in the matching. For $n=4$ we were able to calculate a list
of $1080$ inequalities which determine the polytope, but it is not clear if
there is a more enlightening description.

### Acknowledgements

The author would like to thank Maurício Collares for pointing out the
relevance of Lucas’ theorem, as well as for providing computational evidence
to motivate the questions in Section 3. The author was supported in part by
the Austrian Science Fund (FWF) [10.55776/P36131]. For the purpose of open
access, the author has applied a CC-BY public copyright licence to any Author
Accepted Manuscript version arising from this submission.

## References

  * [1] P. N. Balister, E. Győri, and R. H. Schelp.  Coloring vertices and edges of a graph by nonempty subsets of a set.  Eur. J. Comb., 32(4):533–537, 2011.
  * [2] H. D. L. Hollmann, K. Khathuria, A. Riet, and V. Skachek.  On some batch code properties of the simplex code.  Des. Codes Cryptography, 91(5):1595–1605, 2023.
  * [3] R. N. Karasev and F. V. Petrov.  Partitions of nonzero elements of a finite field into pairs.  Isr. J. Math., 192:143–156, 2012.
  * [4] D. Kohen and I. Sadofschi.  A new approach on the seating couples problem.  arXiv preprint arXiv:1006.2571, 2010.
  * [5] D. Kohen and I. Sadofschi Costa.  On a generalization of the seating couples problem.  Discrete Math., 339(12):3017–3019, 2016.
  * [6] B Kovács.  Finding pairwise disjoint vector pairs in $\mathbb{F}_{2}^{n}$ with a prescribed sequence of differences.  In European Conference on Combinatorics, Graph Theory and Applications, number 12, 2023.
  * [7] E. Preissmann and M. Mischler.  Seating couples around the king’s table and a new characterization of prime numbers.  Am. Math. Mon., 116(3):268–272, 2009.


# DV-QKD Coexistence With 1.6 Tbps Classical Channels Over Hollow Core Fibre

Obada Alia, Rodrigo S Tessinari, Sima Bahrani, Thomas D Bradley, Hesham Sakr,
Kerrianne Harrington, John Hayes, Yong Chen, Periklis Petropoulos, David
Richardson, Francesco Poletti,
George T Kanellos, Reja Nejabati, Dimitra Simeonidou O. Alia, S. Bahrani, G.
Kanellos, R. Nejabati and D. Simeonidou are with the High Performance Network
Group School of Computer Science, Electrical & Electronic Engineering and
EngineeringMaths (SCEEM), University of Bristol BS8 1TH, Bristol, UK, e-mail:
(obada.alia@bristol.ac.uk).R. S. Tessinari is with Toshiba Europe Limited,
Cambridge Research Laboratory, 208 Science Park, Milton Road, Cambridge CB4
0GZ, UK. Previously with the High Performance Network Group at the University
of BristolT. D. Bradley is with the High-Capacity Optical Transmission
Laboratory, Eindhoven University of Technology, 5600 MB, Eindhoven, The
Netherlands.H. Sakr, K. Harrington, J. Hayes, Y. Chen, P. Petropoulos, D.
Richardson and F. Poletti are with the Optoelectronics Research Centre,
University of Southampton, Southampton SO17 1BJ, UK.

###### Abstract

The feasibility of coexisting a quantum channel with carrier-grade classical
optical channels over Hollow Core Nested Antiresonant Nodeless Fibre (HC-NANF)
is experimentally explored for the first time in terms of achievable quantum
bit error rate (QBER), secret key rate (SKR) as well as classical signal bit
error rates (BER). A coexistence transmission of 1.6 Tbps is achieved for the
classical channels simultaneously with a quantum channel over a 2 km-long HC-
NANF with a total coexistence power of 0 dBm. To find the best and worst
wavelength position for the classical channels, we simulated different
classical channels bands with different spacing between the quantum and
classical channels considering the crosstalk generated from both Raman
scattering and four-wave-mixing (FWM) on the quantum channel. Following our
simulation, we numerically estimate the best (Raman spectrum dip) and worst
locations (Raman spectrum peak) of the classical channel with respect to its
impact on the performance on the quantum channel in terms of SKR and QBER. We
further implemented a testbed to experimentally test both single mode fibre
(SMF) and HC-NANF in the best and worst-case scenarios. In the best-case
scenario, the spacing between quantum and classical is 200 GHz (1.6 nm) with
50 GHz (0.4 nm) spacing between each classical channel. The SKR was preserved
without any noticeable changes when coexisting the quantum channel with eight
classical channels at 0 dBm total coexistence power in HC-NANF compared to a
significant drop of 73% when using SMF at $-24$ dBm total coexistence power
which is 250 times lower than the power used in HC-NANF. In the worst-case
scenario using the same powers, and with 1 THz (8 nm) spacing between quantum
and classical channels, the SKR dropped 10% using the HC-NANF, whereas in the
SMF the SKR plummeted to zero.

###### Index Terms:

Hollow Core Nested Antiresonant Nodeless Fibre, Quantum Key Distribution,
Bandwidth variable transceivers, Quantum and Classical Channel Coexistence.

##  I Introduction

Quantum communication networks present a revolutionary step in the field of
quantum communication [1]. In the last two decades, quantum networks evolved
from relayed point-to-point quantum communication links - with main focus in
quantum security applications such as QKD - to multiuser quantum networks with
partial or full connectivity beyond point-to-point links functioning in real
life environment [2]. However, one of the main goals of quantum communication
is to provide a worldwide connectivity via complex networks – similar to the
current internet – with an ultimate security based on the laws of physics
rather than computational complexity [3]. The Quantum Internet vision is
moving towards interconnecting seamlessly multiple nodes and enabling
applications beyond QKD such as blind and distributed quantum computing [4].
For such applications to reach their maximum capabilities, a seamless medium
allowing the coexistence of high power classical channels with quantum
channels would be beneficial for the integration of quantum technologies with
the current optical infrastructure.

QKD is an implementation of quantum cryptography technologies that utilises
the basic principles of quantum mechanics such as the uncertainty principles
to generate information secure symmetric keys that could be used to encrypt
the decrypt key information [1]. QKD has been considered as the ultimate
physical layer security since the quantum-grade encryption mechanism proves to
be information theoretical secure [5]. In discrete variable QKD (DV-QKD) [5],
the process of generating the symmetric key is undertaken by transmitting
encoded single photons between the sender and receiver, Alice and Bob
respectively over an authenticated quantum channel. This quantum channel is
also supported by classical communication channels for basis reconciliation,
error correction and privacy amplification protocols [6].

The integration of quantum technologies with the current classical optical
networking infrastructure would help in the wider deployment of QKD schemes.
However, since the commercially deployed DV-QKD schemes rely on the exchange
of single or few photons between Alice and Bob and are limited by the power
budget of the optical link used for quantum channel, deployment of quantum
channels coexisting with classical optical communication channels is therefore
a very challenging task. To further elaborate on the challenges of such
schemes, the optical power of classical communication channels is orders of
magnitude higher than the power of the quantum channel for quantum
communication and the additional noise generated from such channels due to
optical non-linear effects like Raman scattering and four wave mixing (FWM)
[7] degrades the quantum channel performance extensively and was studied in
[8, 9, 10, 11, 12]. Therefore, coping with optical nonlinearity represents a
major challenge for QKD systems. A high isolation (>110 dB) [13] between the
quantum and classical channel is required and usually is implemented by sharp
filters or cascading multiple filters to partly isolate the noise from the
classical channels and distinguish the quantum signal from adjacent channel
crosstalk; hence improving the quantum signal. Nevertheless, non-linear
processes in glass fibres also generate excessive photons at the same
wavelength as the quantum channel (in-band noise) that cannot be mitigated or
filtered, dramatically affecting its performance.

The coexistence of quantum and classical channels where they are both in the
same optical wavelength bands (C-band) or in different optical wavelength
bands (O-band for the quantum and C-band for the classical) has been
extensively studied and demonstrated in [9, 14, 8, 15, 16]. Coexisting the
quantum and classical channels in the C-band has many advantages, such as
lower fibre loss in the C-band comparing to the O-band and better
compatibility with the current fibre infrastructure. However, the nonlinear
effects such as Raman scattaring becomes more dominant, hence a lower total
coexistence power is achievable as recorded from previous implementations,
namely $-18.5$ dBm, $-22$ dBm and $-25.5$ dBm in [9, 8, 15] respectively. In
contrast, when coexisting the quantum and classical channels in different
bands, the total coexistence power could be as high as $-2$ dBm as shown in
[14].

Here, we consider a radical solution to this challenge by adopting a
transmission fibre technology, in which nonlinearities are simply not present.
Antiresonant Hollow Core Fibres (AR-HCF) are optical fibres where light is
guided in a hollow core via anti-resonance effects within the glass membranes
surrounding the core [17]. AR-HCF surpass the performance of glass core fibres
as they have ultra-low optical mode overlap with the glass (and hence reduced
optical nonlinearity and Rayleigh scattering), a lower latency, a very low
chromatic dispersion, an unprecedented polarisation purity and ultra-low back
scattering [18, 19, 20]. Significant improvement in optical performance of
HCFs were made possible by inducing new features such as negative curvature
core surround [21], then by fibre topologies with a single layer of non-
touching tubes [22], and finally by the addition of small nested tubes - a
design known as hollow core Nested Antiresonant Nodeless Fibre (HC-NANFs)
[23]. HC-NANFs provide a solution for coexistence scheme as they do not only
provide several attractive advantages compared to glass core fibres, such as
reduced optical nonlinearity and ultra-lower latency, but also have the
ability to reach a total loss value lower than that of conventional solid
single-mode fibres [23]. These desirable qualities enable the transmission of
classical channels at high optical powers while coexisting with quantum
channels over the same medium [24]. The losses of HC-NANF have improved
immensely from 1.3 dB/km [25] (used in this experiment), to 0.28 dB/km [26]
and recently to 0.22 dB/km [27].

In this paper, we demonstrate the coexistence of a 16 dB power budget
commercial DV-QKD system (Clavis3 [28]) and 8 x 200 Gbps 16-QAM carrier-grade
classical optical channels at an extremely high total coexistence power of 0
dBm over a 2 km HC-NANF, revealing minimal effects on the quantum channel
performance. We also numerically estimate the best and worst wavelength
location for coexisting an 8-channels classical band with a fixed quantum
channel and compare the QKD performance in terms of SKR and QBER for both
scenarios in HC-NANF and SMF. The paper is organized as follows: in Section II
we discuss the experimental system setup followed by Section III where we
discuss the numerical analysis for the classical channels position. Section IV
present the results and is divided into three parts, optical spectrum for
quantum/classical coexistence IV-A, quantum system characterisation IV-B and
coexistence analysis based on the quantum channel performance IV-C. Finally,
in Section V we conclude the paper.

##  II Experimental System Setup

Figure 1: Experimental Testbed for the Coexistence of 1.6 Tbps classical
channels and DV-QKD channel over 2km HC-NANF and SMF. Inset: scanning electron
micrograph (SEM) of the HC-NANF cross section.

Fig. 1 shows the experimental system setup used to demonstrate the coexistence
of eight classical channels with a DV-QKD channel over a 2 km-long HC-NANF and
SMF (as inset figure, Fig. 1 also illustrates the cross-sectional diagram of
the HC-NANF used [25]). The testbed facilitates the experimental evaluation of
the nonlinear effects on the performance of the quantum channel created by the
presence of classical channels spectrally close to the quantum channel using
both HC-NANF and SMF. For the classical channels, two optical packet DWDM
platforms are used with bandwidth-variable transponders (BVTs) [29]. Each of
these units include four BVT ports reconfigurable to coherent 100 Gbps (PM-
QPSK), 150 Gbps (8-QAM) or 200 Gbps (16-QAM) and each port can be tuned to any
of the 100 wavelengths in the C-band included in the ITU-T grid with 50 GHz
offset. Adaptable soft-decision forward error correction (SD-FEC) is also
available in the transponders to enable maximum transmission capacity with
minimum errors. Since the HC-NANF has ultra-low optical nonlinearity, it
allows the transmission of classical channels at an extremely high power (over
0 dBm) with minimal effect on the quantum channel. Therefore for the HC-NANF
since we are using high total coexistence power, all the ports were configured
with the 16-QAM modulation for a maximum capacity of 200 Gbps per channel
resulting in 1.6 Tbps of transmission overall and the SD-FEC considered was
15% with a detector sensitivity of -26 dBm. As for the SMF, due to the
nonlinear effects in this type of fibre the coexistence power was reduced
compared to the HC-NANF case, and a substantially higher detector sensitivity
is required to receive an error-free classical signal. Therefore, for the SMF
a PM-QPSK modulation was used for a maximum capacity of 100 Gbps per channel
resulting in 800 Gbs of transmission overall and the SD-FEC considered was 25%
with a detector sensitivity of -35 dBm.
For the quantum channel, IDQuantique DV-QKD systems are used (Clavis3 QKD
Platform [28]). These systems are implemented to run with the COW (Coherent
One-Way) protocol [30]. The aim of the COW protocol is to make the
implementation of a QKD system as simple as possible using only two single
photon detectors and to allow a significant increase of the final secret key
rate. Moreover in the Clavis3 system, Alice and Bob encode bits with quantum
state carried by single photons (qubits) using time-bin encoding method by
creating a pair of coherent pulses propagating in the same spatial mode and
separated by a given time. The first pulse is called the early pulse and the
second one is called the late pulse. Hence, the measurement method to analyse
this basis is simply to measure the time of detection of the optical pulse. If
one detection occurs in the early time-bin, the qubit value is a $|0\rangle$
state, whereas if it occurs in the late time-bin, the qubit is a $|1\rangle$
state.

TABLE I: Parameters for HC-NANF Coexistence Testbed Parameters | Value
---|---
Classical Channels
Number of Channels | 8
|  Classical Channel
---
Frequencies scenario 1
| 193.50 THz, 193.45 THz,
---
193.40 THz, 193.35 THz,
193.30 THz, 193.25 THz,
193.20 THz, 193.15 THz
|  Classical Channel
---
Frequencies scenario 2
| 192.70 THz, 192.65 THz,
---
192.60 THz, 192.55 THz,
192.50 THz, 192.45 THz,
192.40 THz, 192.35 THz
Grid Spacing | 50 GHz
Modulation Format | 16-QAM
| Optical Signal-to-Noise
---
Ratio (OSNR)
20 dB
Transmission rate per channel | 200 Gbps
Total transmission rate | 1.6 Tbps
Pre-FEC Level | 15%
Detector sensitivity* | -26 dBm
Quantum Channel
DV-QKD Wavelength | 1547.72 nm
DV-QKD Frequency | 193.70 THz
QKD Protocol | COW
Maximum Loss | 16 dB
Optical Band Pass/Rejection Filter (OBRF)
| Insertion loss band pass port
---
0.5 dB
| Center wavelength
---
band pass port
1547.72 nm
Bandwidth band pass port | 100 GHz
*Corresponding to 16-QAM Modulation @200 Gbps and back-to-back.

As shown in Fig. 1, the eight coherent output ports of the two BVTs are
multiplexed using a wavelength selective switch (WSS) shown as MUX1 with 5 dB
of insertion loss for a total throughput of 1.6 Tbps. The WSS is used as a
multiplexer and a band pass filter to couple the classical channels into a
single fibre and provide a 30 dB isolation. The WSS combined output is
connected to the input port of a tunable band pass filter (TBPF1) with 60 dB
of isolation and extremely sharp filter edges followed by a Gaussian shape
tunable band pass filter (TBPF2) with 30 dB of isolation to further suppress
the noise generated by the classical channels. The classical channels are then
coupled to the quantum channel through a WDM multiplexer (DWDM add drop
filter) MUX2 in a co-propagation coexisting configuration and travels through
a 2 km of HC-NANF or SMF. After the coexistence within the selected fibre, the
output is connected into a WDM demultiplexer tuned to 1547.72 nm (DEMUX1) that
passes the quantum channel while rejecting all other channels (classical
channels). After the WDM demultiplexer, the quantum signal is passed through a
double stage filtering using two fixed WDM band-pass filters centred at the
quantum channel wavelength (1547.72 nm) to eliminate any noise generated by
the classical channels. The output of the second filter is connected to the
DV-QKD Bob unit, which undertakes the single-photon detection and processing
of the encoded photons allowing the completion of the COW protocol of the DV-
QKD systems. The output of the reflection port of the WDM demultiplexer is
directed to an optical isolator with insertion loss of 1 dB to prevent the
tunable laser used by the BVT Rx as a local oscillator from returning to the
Bob-QKD unit and interfering with the QKD measurements. It also prevents the
Amplified Spontaneous Emission (ASE) noise generated by the Erbium-Doped Fibre
Amplifier (EDFA) which is used to set the classical signals to the suitable
detectable power levels before demultiplexing via a 1x8 splitter. Furthermore,
an optical isolator is used after the Alice-QKD unit to prevent the noise
generated from the classical channels returning to the Alice-QKD unit and
interfering with the quantum signal. After the optical isolator at the Alice
side, a variable optical attenuator (VOA) is used to adjust the losses of the
quantum channel to the minimum operational level of 10 dB to optimise the
functionality of the Clavis3 QKD system and to prevent over saturating the
single photon detector in the Bob-QKD unit. It also allows us to obtain a
direct comparison of the quantum/classical coexistence between the SMF and HC-
NANF in terms of losses. In this testbed, the total end-to-end (Alice to Bob)
quantum channel loss is 10.5 dB. A summary of main parameters of the testbed
is shown in Table I.

##  III Channel Spacing/Position for Quantum/Classical Coexistence

For the classical channels, we consider a classical band of 8 channels and 50
GHz spacing between each channel, coexisting with a fixed quantum channel
(193.7 THz) in SMF. In [31, 32, 33], optimal wavelength assignment in hybrid
quantum-classical systems have been investigated. Our calculation is based on
the fact that the classical channels are located in the C-band and the quantum
channel is at ITU channel 37 (193.7 THz). Furthermore, it is assumed that the
quantum channel is in the anti-stokes region of the Raman spectrum of the
classical channels. We denote the coexistence power by $I_{\rm c}$.
Considering both Raman scattering and four-wave-mixing (FWM), we numerically
calculate the best and worst locations of a classical band with respect to its
impact on the performance of the quantum channel in terms of SKR and QBER.
Fig. 2 shows the best and worst cases for the spacing between the classical
band and the quantum channel for different values of $I_{\rm c}$. For low
values of $I_{\rm c}$, Raman scattering is the dominant source of noise.
Therefore, the optimum location of the classical band is at the dip of the
Raman spectrum which is near the quantum channel (100-200 GHz spacing),
whereas the worst case is at the peak of the Raman spectrum (1.0-1.4 THz
spacing). As $I_{\rm c}$ increases, both Raman noise and FWM crosstalk
increase as well, however, FWM crosstalk grows much faster and becomes more
dominant at higher coexistence power $\approx$ -7 dBm. Since FWM products
close to the classical band are stronger, for higher values of $I_{\rm c}$ the
space between the quantum channel and the classical band should be chosen such
that the FWM products with high power $f_{ijk}=f_{i}\pm f_{j}\pm f_{k}$ do not
lie in the quantum channel filter bandwidth [34]. In Fig. 2, the photon counts
for FWM noise and Raman noise in the special case of 200 GHz (1.6 nm) spacing
are also depicted. Here, the time duration of the detector and the detection
efficiency are assumed to be $100~{}{\rm ps}$ and $0.7$, respectively.

Figure 2: Best and worst case scenarios for the classical band spacing,
considering different values for the coexistence power in the SMF. The photon
counts for Raman noise and FWM noise in the special case of classical band
spacing equal to 200 GHz (1.6 nm) are also shown in the figure.

##  IV Experimental Results

###  IV-A Optical Spectrum of Quantum and Classical Channels

Fig. 3 shows the spectrum of $8~{}\times~{}200$ Gbps classical channels
coexisting with a quantum channel obtained by using a tap coupler as shown in
Fig. 1 MON1 & MON2 for the best and worst case scenarios. Fig. 3 also
demonstrates the optical filters profile (pink) of the tunable band pass
filters (TBPF1 & TBPF2) used to filter out the noise generated by the
classical channel. It also shows the optical filter profile of the Band-pass
filter (BPF) used before the Bob-QKD unit (red) to suppress the nonlinear
effects from passing to the QKD detectors. The filter profiles were generated
by using the output of an EDFA as a noise source. Furthermore, the spectrum of
the eight classical channels after the filtering (black) MON1 and the
amplification stage (blue) MON2 is shown before demultiplexing and coherent
detection at the BVT receiver.

Figure 3: Spectrum of the combined transmission of quantum and classical
channels with optical filter profiles of both scenarios.

|
---|---
(a) | (b)
|
(c) | (d)

(a)

(b)

(c)

(d)

Figure 5:  a) Average SKR versus total coexistence power using HC-NANF. b)
Average SKR versus total coexistence power using SMF. c) Average QBER versus
total coexistence power using HC-NANF. d) Average QBER versus total
coexistence power using SMF. 200 GHz spacing between quantum and classical
channels. The solid lines in the figures are trend-line fit.

###  IV-B Co-propagation of Quantum and Classical Channel Coexistence without
Fibre - System Characterisation

To investigate the performance of the filtering stages in our system, we
implemented the best-case scenario using the same coexistence scheme in Fig. 1
without using fibres; hence removing most of the nonlinearites in the system.
Since no fibres are used, the degradation in the quantum channel performance
is mainly due to the photon leakage to the Bob-QKD unit from filtering
deficiencies. As shown in Fig. 4, the average SKR stays within the expected
range between 2100-2300 bps for a total coexistence power of -12 dBm which is
higher than the maximum total coexistence power of -24 dBm used in SMF
verifying that the degradation in the SMF fibre case is mostly due to
nonlinear effects i.e. Raman scattering. This comparison is to demonstrate the
effectiveness of the filtering stages in preventing photon leakage to the Bob
QKD-unit for the power levels which were used in SMF. Moreover, the 1,2,4 and
8 channels symbols shown in Fig. 4, corresponds to the number of channels and
their aggregated power. For example, 2 channels at -9 dBm corresponds to a
power of -12 dBm per channel, 4 channels at -9 dBm corresponds to a power of
-15 dBm per channel, and finally 8 channels at -9 dBm corresponds to a power
of -18 dBm per channel.

Figure 4: Average SKR versus total coexistence power without fibre - system
characterisation at 200 GHz spacing.

###  IV-C Co-propagation of Quantum and Classical Channel Coexistence over HC-
NANF and SMF with different Coexistence Power and Channels Spacing

|
---|---
(a) | (b)
|
(c) | (d)

(a)

(b)

(c)

(d)

Figure 6: a) Average SKR versus total coexistence power using HC-NANF. b)
Average SKR versus total coexistence power using SMF. c) Average QBER versus
total coexistence power using HC-NANF. d) Average QBER versus total
coexistence power using SMF. 1 THz spacing between quantum and classical
channels. The solid lines in the figures are trend-line fit.

Fig. 5 and Fig. 6 show the experimental evaluation of SKR and QBER of the
quantum channel with different spectral spacing between the classical and
quantum channels based on the calculation in Section III. For a spacing of 1.6
nm (200 GHz) in SMF (the best case scenarios for SMF considering both Raman
scaterring and FWM), Fig. 5b and Fig. 5d show that the SKR value drops from
$\approx$ 2300 bps to $\approx$ 600 bps while the QBER values increase from
$\approx$ 3% to $\approx$ 4.2% when coexisting 8 classical channels at a total
coexistence power of -24 dBm. This significant 73% drop in the SKR and 40%
rise in the QBER is due to high optical nonlinearity in SMF causing a high
noise leakage to the Bob-QKD unit. Using the same spacing of 1.6 nm in HC-
NANF, Fig. 5a and Fig. 5c show that, the SKR and QBER values without the
presence of any classical channel (no coexistence) are similar to the values
when coexisting 8 classical channels at a total coexistence power of 0 dBm.
This is due to the ultra-low optical nonlinearity in HC-NANF. For a spacing of
8 nm (1 THz) in SMF (the worst case scenarios for SMF considering Raman
scaterring) and using the same coexistence power of -24 dBm, the SKR drops to
zero as the QBER exceeds the operational threshold of 5.2% of the QKD as shown
in Fig. 6b and Fig. 6d. This is because the 8 nm spacing between the quantum
and classical channels resulted in the peak of the Raman scattering spectrum
to be located within the bandwidth of the BPF filter used to filter the
quantum channel; hence the system stops generating keys due to excessive Raman
noise leakage over the quantum channel. With the same spacing of 8 nm and
coexistence power of 0 dBm in HC-NANF, Fig. 6a and Fig. 6c show that, the SKR
drops from $\approx$ 2300 bps to $\approx$ 2050 bps while the QBER values
increase from $\approx$ 2.2% to $\approx$ 3.1% using HC-NANF.

In both scenarios, the back-to-back losses of the quantum channel is similar
when using both SMF and HC-NANF (10.5 dB) and the highest coexistence power (0
dBm) in the HC-NANF is 250 times higher than the highest coexistence power
(-24 dBm) in SMF. Even though the coexistence power is much higher is HC-NANF,
the ultra-low nonlinear effects of the due to its hollow core preserve or have
a slight effect on the SKR and QBER of the quantum channel. This proves the
suitability of HC-NANF as an excellent transmission medium for high power
coexistence of quantum and classical channels.

##  V Conclusion

The coexistence of a DV-QKD channel and 8 × 200 Gbps 16-QAM classical channels
was successfully demonstrated over a 2-km long Hollow Core Nested Antiresonant
Nodeless Fibre for a record-high coexistence transmission of 1.6 Tbps. The
best and worst spectral positions of classical channels were determined via
numerical calculation of the impact of the nonlinear effects such as Raman
scattering and four-wave mixing on the quantum channel performance. In the
best case scenario with 1.6 nm spacing between the quantum and classical
channels (200 GHz) at -24 dBm total coexistence power in SMF, the SKR dropped
by 73%, whereas, at 0 dBm total coexistence power in HC-NANF (250 times higher
than the power used in SMF), the SKR was preserved. In the worst case scenario
8 nm spacing between the quantum and classical channels (1 THz - Raman
spectrum peak at the quantum channel wavelength) and using the same powers,
the SKR dropped to 0% in SMF meaning no keys were generated, whereas in HC-
NANF the SKR dropped only by 10%. This significant difference in the QKD
performance proves the advantage of using HC-NANF to provide a seamless
coexistence of quantum and classical channels for future quantum networks with
a minimal effect on the quantum channel performance regardless of the number
and power of the classical channels.

## Acknowledgement

This work was funded by EU funded project UNIQORN (820474) and the EPSRC
Airguide Photonics Collaboration Fund Award (ref: 517129) (EP/P030181/1). Part
of the research leading to this work has been supported by the Quantum
Communication Hub funded by the EPSRC grant ref. EP/T001011/1 and the ERC
LightPipe project (grant n∘ 682724).

## References

  * [1] S. Pirandola, _et al._ , “Advances in quantum cryptography,” _Advances in Optics and Photonics_ , vol. 12, no. 4, pp. 1012–1236, 2020\.
  * [2] M. Mehic, _et al._ , “Quantum key distribution: a networking perspective,” _ACM Computing Surveys (CSUR)_ , vol. 53, no. 5, pp. 1–41, 2020.
  * [3] S. Valerio, _et al._ , “The security of practical quantum key distribution,” _Rev. Mod. Phys._ , Sep 2009.
  * [4] R. V. Meter, _et al._ , “The path to scalable distributed quantum computing,” _Computer_ , vol. 49, no. 9, pp. 31–42, 2016.
  * [5] N. Gisin, _et al._ , “Quantum cryptography,” _Reviews of modern physics_ , vol. 74, no. 1, p. 145, 2002.
  * [6] D. Elkouss, _et al._ , “Information reconciliation for quantum key distribution,” _arXiv preprint arXiv:1007.1616_ , 2010.
  * [7] A. R. Chraplyvy, “Limitations on lightwave communications imposed by optical-fiber nonlinearities,” _Journal of Lightwave Technology_ , vol. 8, no. 10, pp. 1548–1557, 1990.
  * [8] P. Eraerds, _et al._ , “Quantum key distribution and 1 gbps data encryption over a single fibre,” _New Journal of Physics_ , vol. 12, no. 6, p. 063027, 2010.
  * [9] K. Patel, _et al._ , “Coexistence of high-bit-rate quantum key distribution and data on optical fiber,” _Physical Review X_ , vol. 2, no. 4, p. 041010, 2012.
  * [10] D. Zavitsanos, _et al._ , “Coexistence of discrete-variable qkd with wdm classical signals in the c-band for fiber access environments,” in _21st ICTON_. IEEE, 2019, pp. 1–5.
  * [11] T. F. da Silva, _et al._ , “Impact of raman scattered noise from multiple telecom channels on fiber-optic quantum key distribution systems,” _Journal of lightwave technology_ , vol. 32, no. 13, pp. 2332–2339, 2014\.
  * [12] D. Subacius, _et al._ , “Backscattering limitation for fiber-optic quantum key distribution systems,” _Applied Physics Letters_ , vol. 86, no. 1, p. 011103, 2005.
  * [13] N. Peters, _et al._ , “Dense wavelength multiplexing of 1550 nm qkd with strong classical channels in reconfigurable networking environments,” _New Journal of physics_ , vol. 11, no. 4, p. 045012, 2009.
  * [14] T. Chapuran, _et al._ , “Optical networking for quantum key distribution and quantum communications,” _New Journal of Physics_ , vol. 11, no. 10, p. 105001, 2009.
  * [15] J. F. Dynes, _et al._ , “Ultra-high bandwidth quantum secured data transmission,” _Scientific reports_ , vol. 6, no. 1, pp. 1–6, 2016.
  * [16] E. Hugues-Salas, _et al._ , “11.2 tb/s classical channel coexistence with dv-qkd over a 7-core multicore fiber,” _Journal of Lightwave Technology_ , vol. 38, no. 18, pp. 5064–5070, 2020.
  * [17] S.-f. Gao, _et al._ , “Hollow-core conjoined-tube negative-curvature fibre with ultralow loss,” _Nature communications_ , vol. 9, no. 1, pp. 1–6, 2018.
  * [18] H. Sakr, _et al._ , “Hollow core optical fibres with comparable attenuation to silica fibres between 600 and 1100 nm,” _Nature communications_ , vol. 11, no. 1, pp. 1–10, 2020.
  * [19] A. Taranta, _et al._ , “Exceptional polarization purity in antiresonant hollow-core optical fibres,” _Nature Photonics_ , vol. 14, no. 8, pp. 504–510, 2020.
  * [20] V. Michaud-Belleau, _et al._ , “Backscattering in antiresonant hollow-core fibers: over 40 db lower than in standard optical fibers,” _Optica_ , vol. 8, no. 2, pp. 216–219, 2021.
  * [21] Y. Wang, _et al._ , “Low loss broadband transmission in hypocycloid-core kagome hollow-core photonic crystal fiber,” _Optics letters_ , vol. 36, no. 5, pp. 669–671, 2011.
  * [22] A. N. Kolyadin, _et al._ , “Light transmission in negative curvature hollow core fiber in extremely high material loss region,” _Optics express_ , vol. 21, no. 8, pp. 9514–9519, 2013.
  * [23] F. Poletti, “Nested antiresonant nodeless hollow core fiber,” _Optics express_ , vol. 22, no. 20, pp. 23 807–23 828, 2014.
  * [24] O. Alia, _et al._ , “1.6 tbps classical channel coexistence with dv-qkd over hollow core nested antiresonant nodeless fibre (hc-nanf),” in _2021 European Conference on Optical Communication (ECOC)_. IEEE, 2021, pp. 1–4.
  * [25] T. D. Bradley, _et al._ , “Record low-loss 1.3 db/km data transmitting antiresonant hollow core fibre,” in _European Conference on Optical Communication (ECOC)_. IEEE, 2018, pp. 1–3.
  * [26] G. T. Jasion, _et al._ , “Hollow core nanf with 0.28 db/km attenuation in the c and l bands,” in _Optical Fiber Communication Conference_. Optical Society of America, 2020, pp. Th4B–4.
  * [27] H. Sakr, _et al._ , “Hollow core nanfs with five nested tubes and record low loss at 850, 1060, 1300 and 1625nm,” in _2021 Optical Fiber Communications Conference and Exhibition (OFC)_. IEEE, 2021, pp. 1–3.
  * [28] IDQuantiqueSA, “Clavis3 qkd platform,” https://www.idquantique.com/quantum-safe-security/products/clavis3-qkd-platform-rd/, May 2021.
  * [29] F. Engineering, “An open approach for switching, routing, and transport,” https://engineering.fb.com/2016/11/01/connectivity/an-open-approach-for-switching-routing-and-transport/, Nov. 2016.
  * [30] D. Stucki, _et al._ , “Continuous high speed coherent one-way quantum key distribution,” _Optics express_ , vol. 17, no. 16, pp. 13 326–13 334, 2009\.
  * [31] S. Bahrani, _et al._ , “Wavelength assignment in hybrid quantum-classical networks,” _Scientific reports_ , vol. 8, no. 1, pp. 1–13, 2018.
  * [32] Y. Ou, _et al._ , “Field-trial of machine learning-assisted quantum key distribution (qkd) networking with sdn,” in _2018 European Conference on Optical Communication (ECOC)_. IEEE, 2018, pp. 1–3.
  * [33] S. Bahrani, _et al._ , “Resource optimization in quantum access networks,” in _ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2019, pp. 7988–7992.
  * [34] L. He, _et al._ , “The four wave mixing effects in quantum key distribution based on conventional wdm network,” in _2014 12th International Conference on Optical Internet 2014 (COIN)_. IEEE, 2014, pp. 1–2.


# A Piecewise Monotonic Gait Phase Estimation Model for Controlling a Powered
Transfemoral Prosthesis in Various Locomotion Modes

Xinxing Chen, Chuheng Chen, Yuxuan Wang, Bowen Yang, Teng Ma, Yuquan Leng, and
Chenglong Fu Manuscript received: February 24, 2022; Revised May 17, 2022;
Accepted July 12, 2022. This paper was recommended for publication by Editor
Pietro Valdastri upon evaluation of the Associate Editor and Reviewers’
comments. This work was supported by the National Natural Science Foundation
of China [Grant U1913205, 62103180, and 52175272]; Guangdong Innovative and
Entrepreneurial Research Team Program [Grant 2016ZT06G587]; the China
Postdoctoral Science Foundation (2021M701577) ; the Science, Technology and
Innovation Commission of Shenzhen Municipality [ZDSYS20200811143601004 and
KYTDPT20181011104007]; the Stable Support Plan Program of Shenzhen Natural
Science Fund [Grant 20200925174640002]; and Centers for Mechanical Engineering
Research and Education at MIT and SUSTech.The authors are all with Shenzhen
Key Laboratory of Biomimetic Robotics and Intelligent Systems, Shenzhen,
518055, China, and also with Guangdong Provincial Key Laboratory of Human-
Augmentation and Rehabilitation Robotics in Universities, Southern University
of Science and Technology, Shenzhen, 518055, China. Corresponding author:
Chenglong Fu. Email<EMAIL_ADDRESS>is also with Department of
Biomedical Engineering, National University of Singapore.Digital Object
Identifier (DOI): see top of this page.

###### Abstract

Gait phase-based control is a trending research topic for walking-aid robots,
especially robotic lower-limb prostheses. Gait phase estimation is a challenge
for gait phase-based control. Previous researches used the integration or the
differential of the human’s thigh angle to estimate the gait phase, but
accumulative measurement errors and noises can affect the estimation results.
In this paper, a more robust gait phase estimation method is proposed using a
unified form of piecewise monotonic gait phase-thigh angle models for various
locomotion modes. The gait phase is estimated from only the thigh angle, which
is a stable variable and avoids phase drifting. A Kalman filter-based smoother
is designed to further suppress the mutations of the estimated gait phase.
Based on the proposed gait phase estimation method, a gait phase-based joint
angle tracking controller is designed for a transfemoral prosthesis. The
proposed gait estimation method, the gait phase smoother, and the controller
are evaluated through offline analysis on walking data in various locomotion
modes. And the real-time performance of the gait phase-based controller is
validated in an experiment on the transfemoral prosthesis.

###### Index Terms:

Prosthetics and exoskeletons, human-robot collaboration, gait pattern
analysis, gait phase estimation, virtual constraint

##  I Introduction

Gait phase is a variable that describes the progress of a gait cycle. Compared
with conventional finite-state control methods for powered lower-limb
prosthesis [1, 2, 3, 4], which have to tune several joint impedance
controllers to mimic the damping and stiffness behaviors in separate gait
stages, gait phase-based controllers determine the desired joint angles with
unified models related to the estimated gait phase and require less manually-
tuned hyper-parameters. A more significant superiority of the gait phase-based
controllers is that they allow for semi-volitional control of the prosthesis
because the users have the autonomy in adjusting their own gait phase.
Therefore, gait phase-based control has been a trending research topic in the
field of walking-aid robots [5, 6, 7, 8, 9, 10].

Gait phase-based control was initially proposed by Gregg et al. [6]. A human-
inspired representation of the gait phase was utilized to parameterize the
virtual constraint models of the joint angles of a powered prosthesis. The
virtual constraints encoded the desired joint angles as functions of the
monotonically increasing gait phase in each gait cycle.

A common challenge for gait phase-based control of walking-aid robots is gait
phase estimation. Since the gait phase is highly correlated with a human’s
real-time pose during walking, accurate gait phase estimation can help
determine the desired state of the robots through the virtual constraints in
assistive walking. The gait phase is usually estimated from the human’s lower-
limb motion state measured by mechanical sensors, such as inertial measurement
units (IMUs) [11, 12, 13] and pressure sensors [14, 15]. The thigh angle was
one of the most widely-used variables for gait estimation, because it can
robustly represent ipsilateral lower-limb motion patterns during unsteady
walking [16].

Gregg et al. [6, 17] estimated the gait phase by constructing a scaled and
shifted phase plane of the thigh angle versus its integral. The phase angle
calculated from the phase plane was then normalized and served as the gait
phase. However, the computation of the phase angle presents a challenge for
real-time control of the prosthesis. Besides, the integral of the thigh angle
may introduce accumulative measurement errors in gait phase estimation. Hone
et al. [18] also reported that the thigh angle profile and its integral were
phase-shifted sinusoids and an additional phase correction method was required
for accurate gait phase estimation. The gait phase can also be estimated from
the pair of thigh angle and its angular velocity [19, 20, 21], but the
measured angular velocity is prone to be noisy and make estimation inaccurate.
The thigh angle can also be directly utilized to estimate the gait phase and
it is reported to be more stable than the integral and the angular velocity[5,
8, 22], but there was an obvious saturation of the estimated gait phase in the
late swing phase in the previous works, which was caused by the retraction of
the thigh before the heel-strike event and affected the estimation accuracy.

Due to the above-mentioned challenges, most existing gait phase estimation
methods were dedicatedly designed for a certain locomotion mode (i.e., terrain
type), adopting various smoothing methods, phase plane shifting methods, and
different definitions of the gait cycle to compensate for the measurement
noise and the thigh retraction problem, yet a unified gait phase estimation
method for various locomotion modes is rarely seen[23] and has not been
validated through online robotic experiments.

To provide a solution for estimating gait phase in various locomotion modes,
this paper proposes a unified piecewise monotonic model to describe the
relationship between the gait phase and the human’s thigh angle. Different
from some previous gait phase estimation researches, the retraction stage is
considered in the proposed model, because it is a common feature when
traversing rough terrains[24], and can improve the stability of biped
walk[25]. The inverse function of the model can then be utilized to estimate
the gait phase from only the thigh angle. A Kalman filter-based smoother
inspired by [26] was designed to mitigate the mutations of the estimated gait
phase. A powered transfemoral prosthesis is controlled to track the desired
knee and ankle angles, which are determined by their virtual constraints and
the estimated gait phase.

The main contributions of this paper are twofold:

(1) Towards gait phase estimation in various locomotion modes, a unified
piecewise monotonic function of the thigh angle was proposed to estimate the
gait phase, which can reduce the parameter tuning workload and make it easier
to adapt the controller of the prosthesis to various terrain types.

(2) The thigh retraction stage, which was commonly observed when humans
traverse through rough terrains, was taken into consideration to model a
reverse motion of the swing leg before heel-strike. The added retraction stage
helps improve the accuracy and adaptation of the proposed method in various
locomotion modes and mitigates the potential saturation problem when using
only the thigh angle.

The proposed method is evaluated on a public dataset [27] and our collected
lower-limb joint angle data of five able-bodied subjects. The results show
that the proposed gait phase estimation method can achieve low estimation
errors in various locomotion modes, including transition modes, which have
rarely been analyzed in previous researches, especially for robotic prosthesis
[28]. Compared with other two gait phase estimation methods using the pair of
$\phi$ and $\int\phi$ and using the pair of $\phi$ and $\dot{\phi}$,
respectively, the RMSE of the smoothed gait phase estimated by the proposed
method can reach $37.6\%$ lower than the other two methods.

The real-time performance of the proposed phase estimation method is validated
through the experiments on a transfemoral prosthesis. The results demonstrate
that the proposed gait phase estimation method and the controller can achieve
good real-time gait phase estimation results and joint angle tracking
performance.

The rest of the paper is organized as follows: Sec. II describes the proposed
method. The proposed models and phase smoother are evaluated offline with
datasets in Sec. III. Sec. III also presents the walking experiment conducted
on a powered transfemoral prosthesis and provides discussions and insights on
the obtained results. Sec. IV concludes the paper and presents outlooks on the
gait phase-based control method.

Figure 1: The proposed gait phase estimation method, the Kalman filter-based
smoother, and the joint angle tracking controller for the powered transfemoral
prosthesis. In (A), the gait cycle is divided into three stages, which are
presented with different background colors. In each gait stage, the thigh
rotates to the direction shown as the blue arrow. The thigh angle is
illustrated with a colored solid curve. The gait phase (illustrated with the
red dashed curve) is estimated from the thigh angle $\phi$ with the inverse
function of the piecewise monotonic thigh angle-gait phase model. In (B), the
black dashed-dotted line is a straight line increasing from 0 to 1, which is
regarded as the actual gait phase in a gait cycle from time 0 to $t_{e}$. The
red dashed curve is the estimated gait phase. The purple dashed curve is the
gait phase smoothed by the Kalman filter-based smoother. In (C), the smoothed
gait phase is used to calculate the desired joint angles. The prosthesis
tracks the desired joint angles with a PD-like impedance controller.

##  II Methods

In this section, the proposed gait phase estimation method based on a
piecewise monotonic model was presented. A Kalman filter was designed to
smooth the estimated gait phase. The desired joint angles of the powered
prosthesis are determined by the estimated gait phase and the virtual
constraints fitted from walking patterns of able-bodied subjects. Then a PD-
like impedance controller is applied to the powered prosthesis to track the
desired knee and ankle angles. Fig. 1 presents an overview of the proposed
gait phase estimation method and the controller for the prosthesis. Sec. II-A
presents the mathematical derivation of the proposed phase variable. Sec. II-B
shows how the Kalman filter-based smoother is designed. Sec. II-C introduces
the virtual constraints between the estimated gait phase and the joint angles
on various terrain types. The PD-like joint angle tracking controller of the
powered prosthesis is also presented in Sec. II-C.

###  II-A Gait Phase Estimation

As illustrated in Fig. 1, the thigh angle is not a monotonic variable during a
gait cycle, which makes it hard to estimate the gait phase only from the thigh
angle. Previous researches figured out a solution which introduced another
variable (e.g., the integral or the differential of the thigh angle) to map
the pair of two variables to the gait phase and solve the one-to-many problem.
However, we have discussed in the Introduction that these supplementary
variables are not stable in non-rhythmic walking and can bring in gait phase
estimation drifting. In spite of this, the thigh angle itself is stable
because it is directly measured from the human’s pose with much fewer drifting
and noise issues.

In this paper, in order to estimate the gait phase from only the non-monotonic
thigh angle, an intuitive solution is to divide the gait cycle into a few
stages to make the thigh angle monotonic in each stage. As shown in Fig. 1, a
gait cycle can be divided into three stages: stance, swing and retraction. In
the stance stage, the thigh rotates backwards around the hip joint, thus the
thigh angle $\phi$ decreases monotonically from $\phi_{0}$ to $\phi_{1}$; once
the swing stage starts, the thigh angle increases monotonically from
$\phi_{1}$ to $\phi_{2}$; in the retraction stage, the thigh rotates backwards
again until the heel-strike event, and the thigh angle decreases from
$\phi_{2}$ to $\phi_{0}$. Since the thigh angle is piecewise monotonic, as
long as the gait stage can be divided accurately, the gait phase $s$ can be
mapped to a certain thigh angle $\phi$.

In this paper, a set of biased sigmoid functions serve as the mapping
functions from $s$ to $\phi$:

$\phi=f_{n}(s)=\frac{h_{n}}{1+e^{-k_{n}\left(s-s_{0n}\right)}}+b_{n},n\in\\{\text{st,
sw, re}\\},$ (1)

where st, sw, re refer to the stance stage, the swing stage, and the
retraction stage, respectively. $h_{n},k_{n},s_{0n}$ and $b_{n}$ are
hyperparameters related to the shape of the functions and the optimal
hyperparameters $\left(h_{n},k_{n},s_{0n},b_{n}\right)_{*}$ can be determined
by minimizing the gap between the fitted and the measured thigh angle curve:

$\left(h_{n},k_{n},s_{0n},b_{n}\right)_{*}=\mathop{\arg\min}\limits_{\left(h_{n},k_{n},s_{0n},b_{n}\right)}\sum_{t}\left(f(s_{t})-\phi_{t}\right)^{2}.$
(2)

This set of functions are selected because they fulfill the following
requirements for the mapping functions:

(1) The mapping functions should be monotonic in each gait stage: For any two
gait phases $s_{i}$ and $s_{j}$ in the same gait stage $n$,

$\displaystyle\dot{f}_{n}(s_{i})\cdot\dot{f}_{n}(s_{j})$
$\displaystyle=\frac{k_{n}^{2}h_{n}^{2}e^{-k_{n}\left(s_{i}-s_{0n}\right)}e^{-k_{n}\left(s_{j}-s_{0n}\right)}}{\left(1+e^{-k_{n}\left(s_{i}-s_{0n}\right)}\right)^{2}\left(1+e^{-k_{n}\left(s_{j}-s_{0n}\right)}\right)^{2}}$
$\displaystyle\geqslant 0,n\in\\{\text{st, sw, re}\\},$ (3)

therefore, the selected functions are piecewise monotonic.

(2) The inverse function of the mapping functions should exist, so that the
gait phase can be estimated from the thigh angle: The inverse function of the
selected functions can be formulated as:

$\displaystyle\hat{s}=f_{n}^{-1}(\phi)=$ $\displaystyle
s_{0n}-\frac{1}{k_{n}}\ln\left(\frac{h_{n}}{\phi-b_{n}}-1\right),$
$\displaystyle\phi\in(b_{n},b_{n}+h_{n}),n\in\\{\text{st, sw, re}\\},$ (4)

where $\hat{s}$ is the estimated gait phase.

(3) The inverse function in each gait stage should be continuous to avoid
mutations of the estimated gait phase: The inverse functions are derivable in
each gait stage, and the derivative functions can be expressed as Eq. (5),
therefore they are continuous.

$\displaystyle\dot{f}_{n}^{-1}(\phi)=$
$\displaystyle\frac{h_{n}}{k_{n}(b_{n}+h_{n}-\phi)(\phi-b_{n})},$
$\displaystyle\phi\in(b_{n},b_{n}+h_{n}),n\in\\{\text{st, sw, re}\\}.$ (5)

With the above piecewise monotonic gait phase-thigh angle model, the gait
phase can be estimated from Eq. (4), as illustrated by the red dashed curve in
Fig. 1(a).

###  II-B Kalman Filter-based Smoother

The gait phase estimated by Eq. (4) can be noisy due to thigh angle
measurement noises and mutations may occur when gait stages switch (as
illustrated in Fig. 1(b)) because the thigh angle curve in each gait cycle is
not exactly the same as the fitted gait phase-thigh angle model in Eq. (1) and
(2).

To mitigate the above issues, a Kalman filter-based smoother is designed to
filter the noisy estimated phase by assuming that gait phase changing velocity
does not vary largely in each gait cycle. The smoothed gait phase
$\tilde{s}_{t}$ at time step $t$ can be calculated from estimated gait phase
$\hat{s}_{t-1}$ and the phase changing velocity $\dot{\hat{s}}_{t-1}$:

$\left\\{\begin{array}[]{l}\boldsymbol{x}_{t}=\left[\begin{array}[]{l}\tilde{s}_{t}\\\
\dot{\tilde{s}}_{t}\end{array}\right]=\boldsymbol{A}\boldsymbol{x}_{t-1}+\boldsymbol{B}_{t-1}\\\
\boldsymbol{y}_{t-1}=\left[\begin{array}[]{l}\hat{s}_{t-1}\\\
\dot{\hat{s}}_{t-1}\end{array}\right]=\boldsymbol{C}\boldsymbol{x}_{t-1}+\boldsymbol{D}_{t-1}\end{array}\right.,$
(6)

where $\boldsymbol{x}_{t}$ is the smoothed state vector at time step $t$,
$\boldsymbol{y}_{t-1}$ is the state vector estimated by Eq. (4) at time step
$t-1$. $\boldsymbol{A}=\left[\begin{array}[]{cc}1&\delta t\\\
0&1\end{array}\right]$ is the state transition matrix,
$\boldsymbol{C}=\left[\begin{array}[]{cc}1&0\\\ 0&1\end{array}\right]$ is the
observation matrix, $\boldsymbol{B}\sim N(0,\boldsymbol{Q})$ and
$\boldsymbol{D}\sim N(0,\boldsymbol{R})$ are the process noises and
measurement errors, respectively, where
$\boldsymbol{Q}=\left[\begin{array}[]{cc}0&0\\\
0&{\sigma_{B}}^{2}\end{array}\right]$ and
$\boldsymbol{R}=\left[\begin{array}[]{cc}{\sigma_{D1}}^{2}&0\\\
0&{\sigma_{D2}}^{2}\end{array}\right]$ are the covariance matrices. Since the
gait phase changing velocity $\dot{\tilde{s}}$ is expected to be stable,
$\sigma_{B}$ is set to be a small positive value in this paper, while
$\sigma_{D1}$ and $\sigma_{D2}$ are determined based on experience.

A Kalman filter is utilized to update the above matrices through the following
equations:

$\displaystyle\boldsymbol{x}_{t}$
$\displaystyle=\boldsymbol{A}\boldsymbol{x}_{t-1},$ (7)
$\displaystyle\boldsymbol{P}_{t}$
$\displaystyle=\boldsymbol{A}\boldsymbol{P}_{t-1}\boldsymbol{A}^{\mathsf{T}}+\boldsymbol{Q},$
(8) $\displaystyle\boldsymbol{K}_{t}$
$\displaystyle=\boldsymbol{P}_{t}\boldsymbol{C}^{\mathsf{T}}\left(\boldsymbol{C}\boldsymbol{P}_{t}\boldsymbol{C}^{\mathsf{T}}+\boldsymbol{R}\right)^{-1},$
(9) $\displaystyle\boldsymbol{x}_{t}$
$\displaystyle=\boldsymbol{x}_{t}+\boldsymbol{K}_{t}\left(\boldsymbol{y}_{t}-\boldsymbol{C}\boldsymbol{x}_{t}\right),$
(10) $\displaystyle\boldsymbol{P}_{t}$
$\displaystyle=\left(\boldsymbol{I}-\boldsymbol{K}_{t}\boldsymbol{C}\right)\boldsymbol{P}_{t},$
(11)

where $\boldsymbol{P_{t}}$ is the covariance matrix of the state, which
indicates the uncertainty of the state $\boldsymbol{x}_{t}$, and
$\boldsymbol{K}$ is the Kalman Gain.

###  II-C Joint Angle Tracking Control

The models between the gait phase and the desired joint angles (i.e., the
virtual constraints of the joint angles) can be fitted with collected walking
data. In this paper, the virtual constraints
$\alpha_{\text{d}i}=g_{i}(\tilde{s}),i\in\\{1,2\\}$ are built by fitting the
actual joint angle values at all the knot points in a gait cycle to a cubic
spline interpolation of the gait phase, where $i=1$ is for the knee and $2$ is
for the ankle, and $\alpha_{\text{d}i}$ is the desired joint angle.

Once the smoothed gait phase $\tilde{s}$ was obtained through the Kalman
filter-based smoother, the desired knee angle $\alpha_{\text{d}1}$ and the
ankle angle $\alpha_{\text{d}2}$ for the powered prosthesis can be determined.
The desired joint angle and angular velocity then serve as the inputs of a PD-
like impedance controller:

$\tau_{i}=K_{\text{p}i}\left(\alpha_{\mathrm{d}i}-\alpha_{i}\right)+K_{\mathrm{di}}\left(\dot{\alpha}_{\mathrm{d}i}-\dot{\alpha}_{i}\right),i\in\\{1,2\\},$
(12)

where $K_{\text{p}i}$ and $K_{\text{d}i}$ are adjustable parameters,
$\alpha_{i}$ and $\dot{\alpha}_{i}$ are the measured joint angle and angular
velocity of the powered prosthesis. $\tau_{i}$ represents the current sent to
the $i$th joint motor. By tuning $K_{\text{p}i}$ and $K_{\text{d}i}$, the
powered prosthesis can track the desired joint angles accurately while
performing human leg-like stiffness and damping behaviors.

##  III Performance Evaluation

The proposed method was evaluated on both a public dataset [27] and our
collected lower-limb kinematic data of five able-bodied subjects. Data
collection was approved and performed under the supervision of the Sustech
Medical Ethics Committee (approval number: 20210009, date: 2021/3/2). The
errors between the estimated gait phase and the actual gait phase were
calculated, and the results were compared with two previously proposed gait
phase estimation methods: (1) using the pair of $\phi$ and $\int\phi$ and (2)
using the pair of $\phi$ and $\dot{\phi}$. The virtual constraints of the knee
and ankle angles were evaluated by calculating the errors between the desired
joint angles and the real joint angles of the subjects. The proposed joint
angle tracking controller was validated on a powered transfemoral prosthesis.

TABLE I: Details of eight subjects No. | Gender | Height (cm) | Weight (kg) | Age | Source
---|---|---|---|---|---
1 | Female | 164 | 53 | 20 | dataset
2 | Female | 158 | 55 | 23 | dataset
3 | Male | 173 | 66 | 25 | dataset
4 | Female | 168 | 54 | 24 | campus
5 | Male | 170 | 66 | 24 | campus
6 | Male | 180 | 90 | 25 | campus
7 | Female | 164 | 54 | 28 | campus
8 | Male | 183 | 71 | 24 | campus

###  III-A The lower-limb motion datasets

In the public dataset [27], subjects were instructed to walk with five
locomotion modes, including : level walking (LW), ramp ascent (RA) and descent
(RD), stairs ascent (SA) and descent (SD). The whole-body joint angles of the
subjects in the above locomotion modes, including eight transition modes (LW-
RA, RA-LW, LW-RD, RD-LW, LW-SA, SA-LW, LW-SD and SD-LW) were recorded in the
dataset. The walking data of three subjects from this dataset (refer to
subject 1 $\sim$ 3 in Table. I) were utilized to evaluate the proposed gait
phase estimation and smoothing method.

We also recruited five able-bodied subjects (subject 4 $\sim$ 8 in Table. I)
with various heights, weights and genders on campus to walk with the above
five locomotion modes and recorded their lower-limb kinematic data. The data
collection devices and environments are shown in Fig. 2. Three IMUs are
attached to the thigh, shank, and heel to record the Euler angles and
accelerations every 0.04 s. The subjects were instructed to walk on a flat
treadmill with 0.8 m/s and 1.2 m/s (LW_slow and LW_fast), an inclined
treadmill with 0.8 m/s (RA and RD), and stairs (SA and SD).

(a) The IMUs attached to the leg of the subjects

(b) The flat treadmill for LW_fast and LW_slow modes

(c) The treadmill with $8^{\circ}$ inclination for RA and RD modes

(d) The stairs for SA and SD modes

Figure 2: Lower-limb motion data collection devices and environments. (a)
shows the IMUs attached to the legs of the subjects. (b) $\sim$ (d) show the
treadmill and the stairs where the data were collected.

###  III-B Offline Evaluation on the Datasets

The lower-limb motion data was processed to divide time-series data into
separate gait cycles. The beginning of the gait cycles is characterized by a
large vertical acceleration measured by the foot IMU. In this paper, the gait
phase $s$ maps the progress of the current gait cycle linearly to $[0,100\%]$,
and its actual value at time $t$ is calculated by:

$s_{t}=(t-t_{h})/T*100\%,$ (13)

where $t_{h}$ is the time of the last heel-strike event, $T$ is the duration
of the current gait cycle. The hyperparameters of the proposed gait phase
model in Eq. (4) were fitted from the actual gait phase $s$ and the mean thigh
angle profile $\Phi$ in all the gait cycles . The virtual constraint models
for the knee and ankle angles in each locomotion mode were also built from the
actual gait phase and the mean joint angles.

The gait phase estimated from the resulting models and the Kalman filtered-
based smoother was evaluated by the mean errors (ME), the root mean squared
errors (RMSE), and the Pearson correlation coefficients [29] ($r$) between the
estimated gait phase $\hat{s}$, the smoothed gait phase $\tilde{s}$, and the
actual gait phase $s$ at all the time steps in all the gait cycles. The
virtual constraints for the knee and ankle angles are evaluated by ME, RMSE,
and $r$ between the desired angles $\alpha_{\text{d}i}$ determined by the
virtual constraints and the actual joint angles $\alpha_{i}$ ($i\in{1,2}$).

The evaluation results on the public datasets are presented in Table. II. The
estimation accuracy was evaluated for all the five rhythmic locomotion modes
and eight transition modes. Despite that the errors of the estimated $\hat{s}$
are low enough for every locomotion mode, the smoothed gait phase $\tilde{s}$
are even closer to the actual gait phase. This demonstrates the advantages of
the Kalman filter-based smoother.

TABLE II: Performance evaluation results on the public dataset Mode |  | $\hat{s}(\%)$ |  | $\tilde{s}(\%)$ |  | $\alpha_{1}\left({}^{\circ}\right)$ |  | $\alpha_{2}\left({}^{\circ}\right)$
---|---|---|---|---|---|---|---|---
| Mean | RMSE | $r$ |  | Mean | RMSE | $r$ |  | Mean | RMSE | $r$ |  | Mean | RMSE | $r$
LW |  | -0.16 | 1.77 | 1.00 |  | -0.22 | 1.25 | 1.00 |  | 0.09 | 3.05 | 0.99 |  | -0.11 | 2.59 | 0.94
RA |  | -0.01 | 1.32 | 1.00 |  | 0.00 | 0.34 | 1.00 |  | 0.02 | 2.38 | 0.99 |  | -0.09 | 1.93 | 0.98
RD |  | -0.68 | 6.16 | 0.98 |  | -0.67 | 2.97 | 1.00 |  | 1.20 | 5.53 | 0.96 |  | -0.69 | 3.64 | 0.94
SA |  | -0.00 | 1.17 | 1.00 |  | 0.01 | 0.43 | 1.00 |  | 0.13 | 2.39 | 1.00 |  | -0.05 | 2.28 | 0.96
SD |  | 0.03 | 2.21 | 1.00 |  | 0.01 | 0.77 | 1.00 |  | 0.04 | 3.27 | 0.99 |  | 0.06 | 2.54 | 0.99
LW-RA |  | 0.02 | 1.45 | 1.00 |  | -0.01 | 0.56 | 1.00 |  | 0.04 | 5.13 | 0.95 |  | -0.13 | 3.73 | 0.94
RA-LW |  | -0.02 | 2.64 | 1.00 |  | 0.06 | 0.74 | 1.00 |  | 0.09 | 3.60 | 0.98 |  | -0.08 | 2.20 | 0.97
LW-RD |  | -0.47 | 5.83 | 0.98 |  | -0.48 | 2.57 | 1.00 |  | 1.06 | 6.92 | 0.94 |  | -0.59 | 2.94 | 0.96
RD-LW |  | -0.72 | 3.48 | 0.99 |  | -1.04 | 2.93 | 1.00 |  | 1.02 | 5.81 | 0.96 |  | -0.58 | 4.29 | 0.95
LW-SA |  | 0.00 | 1.13 | 1.00 |  | 0.03 | 0.58 | 1.00 |  | 0.04 | 7.30 | 0.96 |  | -0.04 | 4.17 | 0.88
SA-LW |  | -0.04 | 0.85 | 1.00 |  | -0.03 | 0.27 | 1.00 |  | 0.02 | 2.56 | 0.99 |  | 0.01 | 1.70 | 0.99
LW-SD |  | 0.06 | 1.65 | 1.00 |  | 0.07 | 0.59 | 1.00 |  | 0.02 | 4.62 | 0.98 |  | -0.09 | 3.44 | 0.98
SD-LW |  | -0.02 | 3.74 | 0.99 |  | 0.03 | 2.54 | 1.00 |  | -0.40 | 8.69 | 0.94 |  | -0.20 | 4.56 | 0.91

On our collected dataset, the proposed method was compared with two previous
gait phase estimation methods, denoted by $\phi$ & $\int\phi$ and $\phi$ &
$\dot{\phi}$, respectively. These two methods were chosen for comparison
because they were typically and widely used in recent gait phase estimation
works[18, 6, 17, 16, 21] and have been validated through robotic experiments
on prostheses, thus can represent the state of the art. Some other references
in this paper[23, 26] were not considered for comparison because the former
one has not been tested on a prosthesis and the latter one was designed only
for stance phase.

The method $\phi$ & $\int\phi$ estimated the gait phase as:

$\hat{s}=\operatorname{atan2}((\int\phi+\Gamma)z,(\phi+\gamma))/2\pi,$ (14)

where the parameters $(z,\gamma,\Gamma)$ were calculated for each locomotion
modes respectively, using the mean thigh angle profile $\Phi$ derived from the
collected data:
$z=\frac{\left|\Phi_{\max}-\Phi_{\min}\right|}{\left|\int\Phi_{\max}-\int\Phi_{\min}\right|}$,
$\gamma=-\left(\frac{\Phi_{\max}+\Phi_{\min}}{2}\right)$, and
$\Gamma=-\left(\frac{\int\Phi_{\max}+\int\Phi_{\min}}{2}\right)$.

The method $\phi$ & $\dot{\phi}$ estimated the gait phase as:

$\hat{s}=\operatorname{atan2}((\dot{\phi}+\Lambda)k,(\phi+\lambda))/2\pi,$
(15)

where the parameters were also calculated for each locomotion modes
respectively:
$k=\frac{\left|\Phi_{\max}-\Phi_{\min}\right|}{\left|\dot{\Phi}_{\max}-\dot{\Phi}_{\min}\right|}$,
$\lambda=-\left(\frac{\Phi_{\max}+\Phi_{\min}}{2}\right)$, and
$\Lambda=-\left(\frac{\dot{\Phi}_{\max}+\dot{\Phi}_{\min}}{2}\right)$. The
above equations calculating the parameters of the two methods were all from
the cited references.

The comparison results are presented in Table. III. It can be seen that our
proposed method can achieve lower estimation errors for the gait phase and
also the desired joint angles. To better illustrate the gait phase and joint
angle estimation results with our proposed method and two other methods, the
results of subject 8 were presented in Fig. 3. It can be seen that the gait
phase curves and the joint angle curves estimated from the proposed method are
closer to the actual values compared with the other two methods. The
estimation performance of $\phi$ & $\dot{\phi}$ is the worst, because the
thigh angle velocity can be very noisy and less stable than the thigh angle
itself.

TABLE III: Performance evaluation results on our collected dataset Mode |  | Method |  | $\tilde{s}(\%)$ |  | $\alpha_{1}\left({}^{\circ}\right)$ |  | $\alpha_{2}\left({}^{\circ}\right)$
---|---|---|---|---|---|---|---|---
|  | Mean | RMSE | $r$ |  | Mean | RMSE | $r$ |  | Mean | RMSE | $r$
LW_slow |  |  $\phi$ & $\int\phi$ |  | 2.55 | 7.12 | 0.97 |  | -2.61 | 10.52 | 0.87 |  | 0.14 | 4.96 | 0.99
|  $\phi$ & $\dot{\phi}$ |  | 9.18 | 20.32 | 0.85 |  | -0.00 | 25.58 | 0.25 |  | -0.98 | 7.26 | 0.98
| Ours: $\phi$ |  | -0.46 | 4.44 | 0.99 |  | -0.17 | 5.61 | 0.96 |  | 0.10 | 3.95 | 0.99
LW_fast |  |  $\phi$ & $\int\phi$ |  | 1.62 | 5.87 | 0.98 |  | -2.77 | 11.05 | 0.87 |  | 0.28 | 5.62 | 0.99
|  $\phi$ & $\dot{\phi}$ |  | 5.77 | 15.89 | 0.90 |  | -2.92 | 24.46 | 0.37 |  | -0.61 | 9.09 | 0.97
| Ours: $\phi$ |  | -0.41 | 4.78 | 0.99 |  | -0.26 | 6.2 | 0.96 |  | -0.05 | 4.58 | 0.99
RA |  |  $\phi$ & $\int\phi$ |  | -1.52 | 4.42 | 0.99 |  | -1.35 | 6.99 | 0.95 |  | 0.49 | 5.21 | 0.98
|  $\phi$ & $\dot{\phi}$ |  | 3.1 | 13.36 | 0.91 |  | 0.58 | 12.58 | 0.85 |  | -0.44 | 6.86 | 0.97
| Ours: $\phi$ |  | -0.23 | 2.85 | 1.0 |  | -0.22 | 4.38 | 0.98 |  | 0.02 | 3.85 | 0.99
RD |  |  $\phi$ & $\int\phi$ |  | 4.45 | 11.84 | 0.93 |  | -3.58 | 16.15 | 0.76 |  | 0.16 | 6.3 | 0.98
|  $\phi$ & $\dot{\phi}$ |  | 8.11 | 23.19 | 0.8 |  | -6.61 | 25.21 | 0.48 |  | -1.01 | 7.62 | 0.97
| Ours: $\phi$ |  | 3.46 | 18.27 | 0.83 |  | -2.8 | 13.59 | 0.84 |  | 0.08 | 5.9 | 0.98
SA |  |  $\phi$ & $\int\phi$ |  | -3.63 | 9.97 | 0.95 |  | -2.87 | 12.6 | 0.88 |  | 0.23 | 6.89 | 0.76
|  $\phi$ & $\dot{\phi}$ |  | 2.37 | 17.34 | 0.84 |  | 1.29 | 17.84 | 0.75 |  | -0.68 | 8.44 | 0.65
| Ours: $\phi$ |  | 0.75 | 10.04 | 0.94 |  | -0.21 | 7.64 | 0.96 |  | -0.14 | 5.44 | 0.86
SD |  |  $\phi$ & $\int\phi$ |  | 19.49 | 24.25 | 0.88 |  | 5.53 | 33.55 | 0.26 |  | -0.99 | 17.88 | 0.46
|  $\phi$ & $\dot{\phi}$ |  | -16.66 | 30.56 | 0.58 |  | -8.51 | 32.44 | 0.19 |  | -1.84 | 22.04 | 0.05
| Ours: $\phi$ |  | 5.14 | 17.73 | 0.82 |  | 9.59 | 22.19 | 0.71 |  | 1.99 | 11.2 | 0.76
Figure 3: The gait phase and joint angle of subject 8 estimated by our
proposed method and two other methods. The solid blue curves and the shaded
region are the mean and standard deviation of actual values. The orange and
the green curves and the shaded region are the estimation results of the two
benchmark methods. The red curves and the shaded region are the results of our
proposed method.

###  III-C Real-time Evaluation on the Prosthesis

In order to validate the proposed gait phase estimation method and the gait
phase-based controller, an experiment was conducted on a powered transfemoral
prosthesis shown in Fig. 4. Subject 8 was instructed to wear the prosthesis
with an L-shaped holder and walk on the treadmill with 0.8 m/s speed. The
process of the experiment was recorded and presented in the supplementary
video of this paper. An IMU was attached to the L-shaped holder to measure the
thigh angle of the subject. The proposed gait phase estimation and smoothing
method was applied to estimate the gait phase in real time. The prosthesis was
driven by the proposed controller to track the joint angles determined by the
virtual constraints. The current commands $\tau$ calculated with Eq. (12) were
sent to the motor driver every 0.03 s. The knee angle and the ankle angle of
the prosthesis were measured by the encoders of the motors.

Figure 4: The powered transfemoral prosthesis. Two brushless DC electric
motors actuated by the Elmo drivers serve as the knee and the revolute pair of
the ankle. The able-bodied subject wears the prosthesis with an L-shaped
holder. A height-increase insole was used in the left shoe of the subject to
compensate for the difference between the lengths of his leg and the
prosthesis.

An example series of thigh angle and the real-time estimated gait phase and
the joint angles are presented in Fig. 5. The ground truth of $s$ in the
experiment is calculated by Eq. (13). It can be seen that the smoothed
estimated gait phase $\tilde{s}$ is highly aligned to the actual gait phase.
The RMSE of $\tilde{s}$ is $7.50\%$, and the correlation coefficient between
$\tilde{s}$ and $s$ is 0.99. The knee angle and the ankle angle of the
prosthesis also show good rhythmicity and follow the subject’s motion pattern
in the LW_slow locomotion mode.

Figure 5: An example series of recorded data and the smoothed estimated gait
phase in the experiment. The desired joint angles are determined by the
virtual constraint models. The actual thigh angle is measured by an IMU. The
actual knee and ankle angle of the prosthesis is measured by the encoders.

###  III-D Discussions

From the above results, the superiority of the proposed gait phase estimation
method can be proved. The feasibility of the proposed gait phase-based
controller has also been validated.

In the offline evaluation of the proposed method, the proposed method can
achieve high estimation accuracy in various locomotion modes, including
transition modes, which have never been analyzed in other gait phase
estimation papers. The gait phase estimation errors can reach 37.6% lower than
other gait phase estimation methods. This result can inspire the researchers
to investigate how the gait phase-based controllers can drive the walking-aid
robots in complex environments, e.g., outdoor environments where the
locomotion modes often change.

It should be noted that in the evaluation results on our collected dataset,
the gait phase estimation accuracy in the RD and SD modes was worse than in
other modes. It can be explained that human is in an unstable state going down
the ramp and the stairs, and the joint angles can vary a lot during the single
support stage [30, 31]. This is also a reason why most accidental falls happen
during ramp descent and stair descent.

From the experiment on the thigh prosthesis, it can be seen that the estimated
gait phase matched the actual gait phase well, however, gait phase estimation
results may be further improved, especially at the early stance stage in each
gait cycle. The Kalman filter-based smoother designed in this paper is
initialized at every heel-strike event, thus leading to inaccurate gait phase
changing velocity estimation at the early stance stage. The smoother can be
improved by estimating a moving average of the gait phase changing velocity in
previous steps, and using the estimated gait phase changing velocity to update
the Kalman filter.

##  IV Conclusions and Outlooks

In this paper, a unified form of piecewise monotonic models was proposed to
reveal the relationship between the gait phase and the thigh angle in various
locomotion modes. The inverse function of the model was used to estimate the
gait phase with only the thigh angle. The estimated gait phase was further
smoothed by a Kalman filter-based smoother to mitigate estimation errors and
mutations. A joint angle tracking control method was proposed for a
transfemoral prosthesis.

The proposed gait phase estimation method was evaluated on a public dataset
and our collected dataset and showed high estimation accuracy. Because the
measured thigh angle is more stable and less noisy compared with its integral
and differential, the proposed method using only the thigh angle can achieve
better performance than two existing gait phase estimation methods using the
pair of $\phi$ and $\int\phi$ and using the pair of $\phi$ and $\dot{\phi}$,
respectively.

The proposed gait phase-based impedance controller showed good real-time
performance in the experiment. But it was only validated in a continuous and
rhythmic locomotion mode, further evaluations should be conducted for various
locomotion modes, especially the RD mode, the SD mode, and the transition
modes. To improve the performance of the powered prosthesis in complex
environments, vision-based environment perception can contribute to updating
the parameters of the gait phase estimation model to adapt to various types of
terrains.

##  V ACKNOWLEDGMENT

The authors would like to thank Dr. Kuangen Zhang from the University of
British Columbia for his contribution to the mathematical reasoning of the
piecewise monotonic gait phase model during his visit at Southern University
of Science and Technology from 2020 to 2021.

## References

  * [1] S. K. Au and H. M. Herr, “Powered ankle-foot prosthesis,” _IEEE Robot. Autom. Mag._ , vol. 15, no. 3, pp. 52–59, 2008.
  * [2] B. E. Lawson, J. Mitchell, D. Truex, A. Shultz, E. Ledoux, and M. Goldfarb, “A robotic leg prosthesis: Design, control, and implementation,” _IEEE Robot. Autom. Mag._ , vol. 21, no. 4, pp. 70–81, 2014.
  * [3] B. E. Lawson, H. A. Varol, A. Huff, E. Erdemir, and M. Goldfarb, “Control of stair ascent and descent with a powered transfemoral prosthesis,” _IEEE Trans. Neural Syst. Rehabil. Eng._ , vol. 21, no. 3, pp. 466–473, 2012.
  * [4] K. Zhang, J. Luo, W. Xiao, W. Zhang, H. Liu, J. Zhu, Z. Lu, Y. Rong, C. W. de Silva, and C. Fu, “A subvision system for enhancing the environmental adaptability of the powered transfemoral prosthesis,” _IEEE Trans. Cybern._ , vol. 51, no. 6, pp. 3285–3297, 2020.
  * [5] S. Rezazadeh, D. Quintero, N. Divekar, E. Reznick, L. Gray, and R. D. Gregg, “A phase variable approach for improved rhythmic and non-rhythmic control of a powered knee-ankle prosthesis,” _IEEE Access_ , vol. 7, pp. 109 840–109 855, 2019.
  * [6] D. Quintero, D. J. Villarreal, D. J. Lambert, S. Kapp, and R. D. Gregg, “Continuous-phase control of a powered knee–ankle prosthesis: Amputee experiments across speeds and inclines,” _IEEE Trans. Robot._ , vol. 34, no. 3, pp. 686–701, 2018.
  * [7] S. Rezazadeh, D. Quintero, N. Divekar, and R. D. Gregg, “A phase variable approach to volitional control of powered knee-ankle prostheses,” _2018 IEEE/RSJ Int. Conf. Intell. Robot. Syst. (IROS)_. IEEE, 2018, pp. 2292–2298.
  * [8] T. K. Best, K. R. Embry, E. J. Rouse, and R. D. Gregg, “Phase-variable control of a powered knee-ankle prosthesis over continuously varying speeds and inclines,” _2021 IEEE/RSJ Int. Conf. Intell. Robot. Syst. (IROS)_. IEEE, 2021, pp. 6182–6189.
  * [9] I. Kang, P. Kunapuli, and A. J. Young, “Real-time neural network-based gait phase estimation using a robotic hip exoskeleton,” _IEEE Trans. Med. Robot._ , vol. 2, no. 1, pp. 28–37, 2019.
  * [10] I. Kang, D. D. Molinaro, S. Duggal, Y. Chen, P. Kunapuli, and A. J. Young, “Real-time gait phase estimation for robotic hip exoskeleton control during multimodal locomotion,” _IEEE Robot. Autom. Lett._ , vol. 6, no. 2, pp. 3491–3497, 2021.
  * [11] C. A. Cifuentes, A. Frizera, R. Carelli, and T. Bastos, “Human–robot interaction based on wearable imu sensor and laser range finder,” _Robot. Auton. Syst._ , vol. 62, no. 10, pp. 1425–1439, 2014.
  * [12] X. Chen, K. Zhang, H. Liu, Y. Leng, and C. Fu, “A probability distribution model-based approach for foot placement prediction in the early swing phase with a wearable IMU sensor,” _IEEE Trans. Neural Syst. Rehabil. Eng._ , vol. 29, pp. 2595–2604, 2021.
  * [13] X. Chen, K. Zhang, Y. Leng, and C. Fu, “A multi-task learning method for human motion classification and person identification,” _2021 IEEE Int. Conf. Adv. Robot. Mechatron. (ICARM)_. IEEE, 2021, pp. 132–137.
  * [14] K. Kong and M. Tomizuka, “Smooth and continuous human gait phase detection based on foot pressure patterns,” _2008 IEEE Int. Conf. Robot. Automat. (ICRA))_. IEEE, 2008, pp. 3678–3683.
  * [15] E. Martini, T. Fiumalbi, F. Dell’Agnello, Z. Ivanić, M. Munih, N. Vitiello, and S. Crea, “Pressure-sensitive insoles for real-time gait-related applications,” _Sensors_ , vol. 20, no. 5, p. 1448, 2020.
  * [16] D. J. Villarreal, H. A. Poonawala, and R. D. Gregg, “A robust parameterization of human gait patterns across phase-shifting perturbations,” _IEEE Trans. Neural Syst. Rehabil. Eng._ , vol. 25, no. 3, pp. 265–278, 2017.
  * [17] D. J. Villarreal, D. Quintero, and R. D. Gregg, “Piecewise and unified phase variables in the control of a powered prosthetic leg,” _2017 Int. Conf. Rehabil. Robot. (ICORR)_ , 2017, pp. 1425–1430.
  * [18] W. Hong, N. A. Kumar, and P. Hur, “A phase-shifting based human gait phase estimation for powered transfemoral prostheses,” _IEEE Robot. Autom. Lett._ , vol. 6, no. 3, pp. 5113–5120, 2021.
  * [19] M. A. Holgate, T. G. Sugar, and A. W. Bohler, “A novel control algorithm for wearable robotics using phase plane invariants,” _2009 IEEE Int. Conf. Robot. Automat. (ICRA)_ , 2009, pp. 3845–3850.
  * [20] J. Lee, W. Hong, and P. Hur, “Continuous gait phase estimation using LSTM for robotic transfemoral prosthesis across walking speeds,” _IEEE Trans. Neural Syst. Rehabil. Eng._ , vol. 29, pp. 1470–1477, 2021.
  * [21] R. Macaluso, K. Embry, D. J. Villarreal, and R. D. Gregg, “Parameterizing human locomotion across quasi-random treadmill perturbations and inclines,” _IEEE Trans. Neural Syst. Rehabil. Eng._ , vol. 29, pp. 508–516, 2021.
  * [22] T. Ma, Y. Wang, X. Chen, C. Chen, Z. Hou, H. Yu, and C. Fu, “A piecewise monotonic smooth phase variable for speed-adaptation control of powered knee-ankle prostheses,” _IEEE Robot. Autom. Lett._ , vol. 7, no. 3, pp. 8526–8533, 2022.
  * [23] K. Zhang, “Sensor fusion for predictive control of human-prosthesis-environment dynamics in assistive walking,” Ph.D. dissertation, University of British Columbia, 2022.
  * [24] D. G. E. Hobbelen and M. Wisse, “Swing-leg retraction for limit cycle walkers improves disturbance rejection,” _IEEE Trans. Robot._ , vol. 24, no. 2, pp. 377–389, 2008.
  * [25] M. Wisse, C. G. Atkeson, and D. K. Kloimwieder, “Swing leg retraction helps biped walking stability,” _5th IEEE-RAS Int. Conf. Hum. Robot., 2005._ IEEE, 2005, pp. 295–300.
  * [26] N. Thatte, T. Shah, and H. Geyer, “Robust and adaptive lower limb prosthesis stance control via extended kalman filter-based gait phase estimation,” _IEEE Robot. Autom. Lett._ , vol. 4, no. 4, pp. 3129–3136, 2019.
  * [27] J. Camargo, A. Ramanathan, W. Flanagan, and A. Young, “A comprehensive, open-source dataset of lower limb biomechanics in multiple conditions of stairs, ramps, and level-ground ambulation and transitions,” _J. Biomech._ , vol. 119, p. 110320, 2021.
  * [28] I. Kang, D. D. Molinaro, S. Duggal, Y. Chen, P. Kunapuli, and A. J. Young, “Real-time gait phase estimation for robotic hip exoskeleton control during multimodal locomotion,” _IEEE Robot. Autom. Lett._ , vol. 6, no. 2, pp. 3491–3497, 2021.
  * [29] P. Schober, C. Boer, and L. A. Schwarte, “Correlation coefficients: appropriate use and interpretation,” _Anesth. Analg._ , vol. 126, no. 5, pp. 1763–1768, 2018.
  * [30] J. E. Zachazewski, P. O. Riley, D. E. Krebs, _et al._ , “Biomechanical analysis of body mass transfer during stair ascent and descent of healthy subjects,” _J. Rehabil. Res. Dev._ , vol. 30, pp. 412–412, 1993.
  * [31] I. V. Boyaninska, “Balance control when ascending and descending stairs of different architectural designs,” Master’s thesis, University of Toronto (Canada), 2018.