Advancements in magnetomicrometry, state-estimation and neuromuscular control

Abstract

Systems, devices, or methods described herein may involve a human augmentation device configured to estimate the state of at least one implanted magnetic target. In some example embodiments, the device may comprise an array of sensors configured to measure one or more target field characteristics. In some example embodiments, the device may comprise a processing element configured to acquire measurements of the one or more target field characteristics and apply information-theoretic processing of the measurements to estimate in real-time the state and error covariance of at least one target state vector. In some example embodiments, the device may comprise a processing element configured to control an actuator of the human augmentation device.

Description

Advancements in Magnetomicrometry, State-estimation and Neuromuscular Control RELATED APPLICATION

[0001] This application claims the benefit of U.S. Provisional Application No.63 / 516,760, filed on July 31, 2023.The entire teachings of the above application are incorporated herein by reference. BACKGROUND

[0002] The design and use of magnetometer arrays and global optimization algorithms that enable high-bandwidth tracking of implanted magnets in muscle-tendon structures and other tissues is taught in US 2020 / 0305765 and WO 2022 / 087297, the entire teachings of each being incorporated herein by reference.

[0003] The platforms presented in US 2020 / 0305765 and WO 2022 / 087297 suffer from performance gaps relating to: limited precision and accuracy when magnet depths are greater than 2 cm to 3 cm and / or more than two magnet targets are presented; and excessive bulk, mass and power consumption—the latter adding to the mass when the additional battery capacity needed to power the electronics is considered. SUMMARY

[0004] A human augmentation device may be configured to estimate the state of at least one implanted magnetic target. The device may comprise an array of sensors, each being configured to measure one or more target field characteristics. The device may comprise a first processing element coupled to the sensors, the processing element being configured to: acquire measurements of the one or more target field characteristics; apply information-theoretic processing of the measurements to estimate in real-time at least one of the state and error covariance of at least one target state vector, the target state vector comprising at least one of: target position, velocity, orientation, angular rate, or target strength; and update the at least onetarget state employing a state update innovation method that incorporates at least one of: (1) a posteriori covariance; (2) an observation matrix relating the change in a sensor characteristic to the change in the state of the at least one target; (3) a measurement variance; or (4) the innovation representing a difference between the measurement and its expected value. The device may comprise a second processing element configured to control an actuator of the human augmentation device.

[0005] In some example embodiments, the first processing element reduces the error covariance by employing at least one stochastic dynamic model of the at least one target state, the model taking the form of at least one of: a Markov model driven by white noise; a dynamic system described by differential or difference equations driven by noise; an oscillator comprising time varying quadrature amplitudes and frequency with the amplitude and frequency driven by a Markov process itself driven by white noise; or a disturbance field that is varying in an inertial frame that is driven by a Markov model, itself driven by white noise or by a dynamic system described by differential or difference equations driven by noise.

[0006] In some example embodiments, the information-theoretic processing of the measurements is accomplished by an Information Filter implemented in which the time-varying observation matrix used in the information gain calculation is derived through real-time computation of the Jacobian that relates the change of a signal characteristic to a change in target state.

[0007] In some example embodiments, the filter employs knowledge of the measurement noise characteristics, the characteristics comprising: measurement noise variance unique to at least one of the signal characteristics and to the at least one sensor component in the array; and measurement noise correlation in time.

[0008] In some example embodiments, the a posteriori state covariance is assumed to have substantially achieved steady-state.

[0009] In some example embodiments, the Jacobian is computed at a nominal location of the target.

[0010] In some example embodiments, the innovation contribution is a sum of Information Vectors computed by one of a plurality of multiple processors within the processing elementwhere each processor is connected to a sub-population of sensors and an Information Vector unique to a sensor is computed.

[0011] In some example embodiments, the observation matrix is computed via interpolation in relation to its value at a nominal location of the target, the type of interpolation being selected from: a Linear interpolation, a Bilinear interpolation, a Trilinear interpolation, or an interpolation based on computation of the Hessian matrix evaluated at the nominal location.

[0012] In some example embodiments, the processing element normalizes and applies known corrections to each in the array of sensor outputs to correct for at least one of the following error sources: translation of the sensor with respect to its nominal position within the sensor coordinate frame comprising at least one of a component of a Cartesian translation; orientation of the sensor with respect to the nominal orientation within the sensor coordinate frame comprising at least one of a rotation about x, a rotation about y and a rotation about z; orientation offset of the unit vectors corresponding to the principal coordinate axes along which signal characteristics are acquired by the sensor, where the offsets are two-degree of freedom offsets represented by a rotation about two of the coordinate-frame axes; sensitivity and / or non- linearity of the measurement of the signal characteristics on at least one of principal coordinate axes; or bias offset in the measurement of at least one of the signal characteristics.

[0013] In some example embodiments, the sensor array and the processing element are packaged onto at least one of one side or both sides of at least one of: a rigid-flex PCBA; a rigid PCBA; or interconnected rigid or rigid-flex PCBAs.

[0014] In some example embodiments, the first and second processing elements are packaged in the same assembly.

[0015] In some example embodiments, the first and second processing elements are packaged into one or more FPGAs.

[0016] In some example embodiments, the communications between the second processing element and the human augmentation device is carried out on a network.

[0017] In some example embodiments, the network synchronizes time using at least one of a Network Time Protocol and Precision Time Protocol.

[0018] In some example embodiments, the network is an Ethernet, EtherCAT, or of a Time- Sensitive Network (TSN) type.

[0019] In some example embodiments, the communications is run on a multi-drop synchronous serial interface network.

[0020] In some example embodiments, the device further comprising an inertial measurement unit (IMU) arranged to report the angular rate of the device. In some example embodiments, the first or second processing element is configured to: predict temporal changes in the magnetic field owing to the instantaneous angular rate; and cancels the predicted component of the magnetic field measurement changes so as to reduce the covariance of the target state estimates.

[0021] In some example embodiments, the initial vector of the at least one target state is informed by a global optimization algorithm that minimizes the prediction residuals.

[0022] In some example embodiments, an information-theoretic smoothing algorithm is applied to the state estimates retrospectively and those smoothed state estimates determine a target state with substantially reduced covariance.

[0023] In some example embodiments, the rigid-flex PCBA can be used to create a tile structure where each tile comprises at least one magnetometer.

[0024] In some example embodiments, the oscillatory response is detected through use of an adaptive information-theoretic tracker that applies synchronous demodulation and phase detection to track the evolution of amplitude and frequency of the oscillation.

[0025] In some example embodiments, the innovation update is computed as the product of an innovation gain matrix and the innovation.

[0026] In some example embodiments, the actuator comprises at least one of an electric motor, a hydraulic motor, or a pneumatic motor activated by the processing element.

[0027] An apparatus employed as part of a neuromuscular controller of a device aimed at augmenting human biomechanical function may comprise an array of sensors each designed to measure one or more target field characteristics. The apparatus may comprise a processing element coupled to the sensors, the processing element configured to: acquire measurements of the field characteristics; and employ a neural processing unit to determine at least one of a set of response parameters comprising a muscle-tendon force, a time derivative of a muscle-tendon force, an elongation or a time derivative of an elongation; the at least one of a set of response parameters being determined based upon a trained model of the relationship between the targetsand the response parameters. The apparatus may comprise an actuator that communicates with the processing element. In some example embodiments, the processing element commands the actuator to deliver the determined dynamic response.

[0028] In some example embodiments, the neural processing unit is trained by a series of synthesized Monte Carlo runs representing a stochastic model of the muscle-tendon unit dynamics and the target dynamics therein, wherein the stochastic model captures the range of variation in at least one of signal characteristics, magnet strength, range of magnet motion and geometric relation between the magnets and the sensing array; or disturbance state in the target state that will present to the apparatus when in use.

[0029] In some example embodiments, the neural processing unit employs a Convolutional Neural Network (CNN) to process at least one of sensed target field characteristic on at least one target.

[0030] In some example embodiments, the actuator is muscle tissue that is activated by artificial muscle stimulation delivered to the tissue by the processing element.

[0031] In some example embodiments, the actuator comprises at least one of an electric motor, a hydraulic motor, or a pneumatic motor activated by the processing element.

[0032] A self-calibrating apparatus for use in tracking the state of magnetic targets may comprise at least one addressable coil that can project a reference magnetic field with known characteristics. The apparatus may comprise an array of sensors each designed to measure at least one component of the magnetic field characteristics arising from the presence of the magnetic targets and the coil. The apparatus may comprise a processing element connected to the array and the coil, the processing element programmed to compensate for at least one of a sensor bias, sensitivity and / or non-linearity, or geometric alignment as a means of improving at least one of precision and accuracy of the target state tracking, the target state tracking derived from the sensor measurements, where the target state comprises at least one of a target position, velocity, or orientation within an apparatus coordinate frame.

[0033] In some example embodiments, the optimization is a global optimization.

[0034] In some example embodiments, finding the coordinate transformation that minimizes the innovation employs an information-theoretic tracking algorithm.

[0035] In some example embodiments, the global optimization is used to create the a priori information that serves to initialize the information-theoretic tracking algorithm.

[0036] In some example embodiments, the tiled-sensor array and the modulated coil are integrated into the same rigid printed circuit assembly or into the same rigid printed circuit subassembly.

[0037] In some example embodiments, the tiled-sensor array and the modulated coil are integrated into the same rigid-flex circuit assembly or into the same rigid-flex printed circuit subassembly.

[0038] In some example embodiments, the global optimization is accomplished by gradient descent and not a Jacobian. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the Office upon request and payment of the necessary fee.

[0040] The foregoing will be apparent from the following more particular description of example embodiments, as illustrated in the accompanying drawings in which like reference characters refer to the same parts throughout the different views. The drawings are not necessarily to scale, emphasis instead being placed upon illustrating embodiments.

[0041] FIG.1 shows, for the case of limb amputation, an example embodiment of a magnetomicrometric-based Neuromuscular Controller modelled as a bi-directional communications channel in an Information-Theoretic Context (diagram not to scale).

[0042] FIG.2 shows Equation 1.2.

[0043] FIG.3 shows Equation 1.3.

[0044] FIG.4 shows Equation 1.4.

[0045] FIG.5 shows Equation 1.5.

[0046] FIG.6 shows Equation 1.6.

[0047] FIG.7 shows Equation 1.7, which relates to capturing the two-state Markov Model Dynamics.

[0048] FIG.8 shows a table with answer to key questions impacting the Magnetomicrometric-based Neuromuscular Controller Design provided by the Information- Theoretic (Information Filter) Analyses of a Typical First-Generation Magnetomicrometric Configuration.

[0049] FIG.9A shows an overview of an example embodiment of a Magnetometer Array Information (MMAI) Processor.

[0050] FIG.9B shows details of an example embodiment of a Magnetic Field Information Vector Processor (MFIVP) Architecture.

[0051] FIG.9C shows details of an example embodiments of an Information Aggregator Processor (IAP) Architecture. sEMG refers to Surface EMG.

[0052] FIG.10A shows numerical calculations associated with adjusting Hij* to account for large magnet excursions from the nominal operating point that invalidates use of a Jacobian excursion-invariant value.

[0053] FIG.10B shows a schematic of Simple Instrumented Muscle Tendon

[0054] FIG.10C shows an example derivation of a Tendon Force, ftendon, and Muscle Length, ^lmuscle. Here the force of the muscle-tendon unit is computed from the measured elongation of a tendon about its nominal length using the positions of magnets 1 and 2, as well as information on tendon modulus. Muscle belly length is measured from the positions of magnets 2 and 3 implanted within or on the origin and insertion tendons, respectively, at the muscle belly output locations near the musculotendinous junctions.

[0055] FIG.11 shows Equations 10a, 10b, and 10c. Equation 10a: Calculation of the output vector, E{ ^yk} and the tendon and muscle displacements, E{ ^p21k} and E{ ^p32k}. The expectation operator E{} is replaced in the following with the ^ designation. Equation 10b shows a calculation of the tendon and muscle unit vectors. Equation 10c shows a calculation of the muscle-tendon Force and Muscle Elongation, including 1- ^ confidence interval for these.

[0056] FIG.12 shows an example implementation of a Phase-Lock Loop to detect lateral oscillations of unknown frequency and amplitude.

[0057] FIG.13 shows Equations 1.13, which relate to an example Information-theoretic Model of a one-magnet system with a disturbance field that is generally rotating in accordance with rotations measured by an IMU.

[0058] FIG.14 shows an example embodiment of Algorithm Deployment in a Network of Magnetometer Arrays.

[0059] FIG.15 shows an example embodiment of a Deep Learning Array Information Processor (DLAIP) Model-based architecture.

[0060] FIGS.16A-16B show an example embodiment of an x-axis magnet tracking at a nominal magnet tracking depth of 2.5 cm on a 10 second DLAIP validation run—the simulation driven by a two-state Markov model with ^=100 msec.

[0061] FIGS.16C-16D show an example embodiment of a y-axis magnet tracking at a nominal tracking depth of 2.5 cm on a 10 second DLAIP validation run—the simulation driven by a two-state Markov model with ^=100 msec.

[0062] FIGS.16E-16F show an example embodiment of a z-axis magnet tracking at a nominal tracking depth of 2.5 cm on a 10 second validation run—the simulation driven by a two- state Markov model with ^=100 msec.

[0063] FIGS.17A-17B show an example embodiment of a x-axis magnet tracking at a nominal tracking depth of 5 cm on a 10 second DLAIP validation run—the simulation driven by a two-state Markov model with ^=100 msec.

[0064] FIGS.17C-17D show an example embodiment of a y-axis magnet tracking at a nominal tracking depth of 5 cm on a 10 second DLAIP validation run—the simulation driven by a two-state Markov model with ^=100 msec.

[0065] FIGS.17E-17F show an example embodiment of a z-axis magnet tracking at a nominal tracking depth of 5 cm on a 10 second DLAIP validation run—the simulation driven by a two-state Markov model with ^=100 msec.

[0066] FIG.18A shows an example Calibration Test Fixture employing the Modulated-Field Calibration Framework. SOM refers to system-on-module. NIST refers to National Institute fof Standards and Technology. NVRAM refers to non-volatile random-access memory.

[0067] FIG.18B illustrates representative Calibration Test Fixture Implementation details showing row and column selection circuitry, current paths into and out of the coil and the equivalent circuit showing the topology of the current drive, row / column drive and coaxial stacking of coil layers. PCBA refers to Printed-Circuit Board Assembly.

[0068] FIG.18C is a top view showing an example arrangement of coils and magnetometer ICs in a region of the PCBA.

[0069] FIG.19A shows an example Coordinate Frame Calibration Arrangement Using Active Field Generation on a separate assembly or integrated into the MMAI.

[0070] FIG.19B shows example Coordinate-frame Mathematics to Align Two Sensor (MMAI) Frames with respect to an Auxiliary Active-Magnetic Field Frame (Note: The auxiliary frame may be integrated onto an MMAI or separately in proximity to the network of MMAI).

[0071] FIG.20 shows an example embodiment of Rigid-flex MMAI Packaging.

[0072] FIG.21 shows an example muscle-tendon magnet deployment for estimation of muscle force, elongation and velocity of agonist muscle and antagonist muscle by magnetometer array for thumb actuation of a limb.

[0073] FIG.22 shows an example embodiment in which a magnet is implanted into or onto bone within an arm, and serves as a reference for measurement of muscle fascicle length changes using magnets embedded within muscles.

[0074] FIG.23 is a cutaway view of a limb showing an example embodiment in which reference bone magnets are implanted into bones. These bone attached magnets serve as a reference for muscle fascicle length change measurements.

[0075] FIG.24 shows an amputated limb that comprises a reference bone magnet implanted into bone. This bone attached magnet serves as a reference for muscle fascicle length change measurements.

[0076] FIG.25 is a cutaway view of the amputated residuum and shows reference bone magnets implanted into bones. These bone-attached magnets serve as a reference for muscle fascicle length change measurements.

[0077] FIG.26 shows an example of a reference bone magnet implanted into tibial bone. This bone-attached magnet serves as a reference for muscle fascicle length change measurements.

[0078] FIG.27 shows an example embodiment of a self-calibrating, tiled sensor array integrated into socket that measures magnet displacements in muscle fascicles and relative to bone.

[0079] FIG.28 shows an example embodiment of a reference bone magnet implanted into tibial bone. This bone-attached magnet serves as a reference for muscle fascicle length change measurements.

[0080] FIG.29 shows an example embodiment of Magnetometer arrays packaged onto rigid fiber-composite arms that are affixed to an osseointegrated implant.

[0081] FIG.30 shows an example embodiment of magnetometer arrays packaged onto rigid fiber-composite arms that are affixed to an osseointegrated implant.

[0082] FIG.31 shows an example embodiment of magnetometer arrays attached onto rigid fiber-composite arms, that are affixed to an osseointegrated implant.

[0083] FIG.32 shows an example embodiment for Free-Space Control of a Robotic Prosthesis via Muscle Magnetomicrometry.

[0084] FIG.33 shows Equations 11a-11r of Appendix A.

[0085] FIG.34 shows an example generic pipeline of estimating error distribution of the estimated pose. DETAILED DESCRIPTION

[0086] A description of example embodiments follows.

[0087] Systems, devices, and methods described herein address previously mentioned deficiencies in the prior art. Some embodiments involve real-time model-based, information- theoretic and AI-based magnet and muscle-dynamics state estimation methods in distributed, scalable, and embedded hardware that augment or replace the global optimization algorithms. Some embodiments involve calibration infrastructure employing programmable, modulated magnetic field references deployed in board-level calibration and tested in manufacturing and embedded in the real-time magnet tracking and neuromuscular control hardware implementations. Some embodiments involve advanced systems integration and packaging employing high-precision, self-calibrating, deformable magnetometer arrays designed for use in dense, multi-muscle-tendon control applications that place a premium on precision and accuracy.

[0088] The systems, devices, and methods taught in this disclosure advance dynamic tissue tracking for applications related to the rehabilitation of persons with neuromechanical pathology, including but not limited to amputation, muscle weakness, and / or muscle paralysis. Thedisclosure also relates to the application of augmentation technology for persons without neuromechanical pathology. Documentation of this work and these advances is contained in Appendix A: Further Magnetomicrometric Advances in Robotic Control.

[0089] In some example embodiments, systems, devices, and / or methods herein may involve a system in which targets, which may be magnetic targets, are implanted in a muscle-tendon complex of a limb. The targets are sensed by a sensing array mounted to the limb. An electronic processor processes the signals from the sensing array to compute muscle-tendon signals used to command an actuator.

[0090] In some example embodiments, systems, devices, and / or methods herein may involve a system in which a sensing array is mounted to an osseointegrated implant to an amputated residuum.

[0091] In some example embodiments, systems, devices, and / or methods herein may involve implanting a target on bone of the limb as a reference for computations relative to targets mounted to a muscle-tendon complex of the limb.

[0092] In some example embodiments, systems, devices, and / or methods herein may involve the application of information-theoretic processing to the signals received from a sensing array.

[0093] In some example embodiments, systems, devices, and / or methods herein may involve calibrating a magnetometer array that senses a magnetic field created by an array of addressable coil structures.

[0094] In some example embodiments, systems, devices, and / or methods herein may involve a distributed network of measurement devices that applies information-theoretic-fusion of measurements to determine target state.

[0095] In some example embodiments, systems, devices, and / or methods herein may involve a neural processing unit that determines state of targets.

[0096] In some example embodiments, systems, devices, and / or methods herein may involve a neural processing unit that determines a set of response parameters based on a trained model of the relationship between targets and the response parameters.

[0097] In some example embodiments, systems, devices, and / or methods herein may involve a tiled-array of sensors that measure a magnetic field from an addressable coil. The current inthe coil can be modulated and calculations performed to compensate for deformation of the tiled array while tracking target state.

[0098] In some example embodiments, systems, devices, and / or methods herein may involve a tiled-array in which modulated fields are used to determine the coordinate frame of each node in relation to any other node.

[0099] A neuromuscular control apparatus to augment or replace the function of a biological limb of a wearer for augmenting muscle function may comprise: (a) targets implanted within at least one biological limb of a wearer with at least one of a muscle-tendon complex of the biological limb having two tendons and one muscle; (b) at least three (3) targets arranged to measure muscle force and elongation of the at least one muscle-tendon complex, wherein at least one target is implanted in a first tendon, and at least two targets are implanted in a second tendon of the at least one muscle-tendon complex; (c) an actuator attached to or within to augment the function of the limb; (d) a sensing array mounted to the limb of one or more elements designed to measure one or more components of a signal emitted by the targets; (e) a first processing element connected to the sensors to process the signals from the sensing arrays and to compute muscle-tendon signals comprising at least one of the force, elongation and / or time-derivatives associated with the at least one muscle-tendon complex; and (f) a second processing element to command an actuating signal to drive the actuator in accordance with the muscle tendon signals in or to augment or replace the function of the at least one muscle-tendon complex.

[0100] The term actuator shall be taken to mean a mechatronic component that produces force, torque, or displacement in a controlled way in response to an input. The mechatronic component can comprise at least one of an electric, hydraulic, or pneumatic motor with inputs supplied by one or more processing elements in the control apparatus. In other embodiments, the actuator’s mechatronic component can be a muscle tissue wherein the actuator is activated by open-loop or closed-loop artificial muscle stimulation delivered to the tissue by the processing element. In some embodiments, the muscle stimulation can be applied using Functional Electrical Stimulation (FES). In some embodiments, the muscle stimulation can be applied using optogenetics.

[0101] In some example embodiments, at least two targets are implanted in the second tendon, at least one target is proximal and at least one target is distal to its bone attachment.

[0102] In some embodiments, the first processing element and the second processing element are integrated.

[0103] In some example embodiments, the force of the muscle-tendon complex is computed as the relative displacement between the two targets in the second tendon multiplied by a tendon stiffness; and the elongation of the muscle is computed as the relative displacement between the magnet implanted in the first tendon and the nearest magnet within the second tendon.

[0104] In some example embodiments, the targets are magnets with properties that are estimated by the first processor; and the sensor is a magnetometer array, each magnetometer measuring at least one component of the magnetic field projected by the three or more targets.

[0105] A neuroprosthetic limb to replace the function of a missing biological limb of a wearer may comprise: (a) an amputated limb residuum with the abutment of an osseointegrated implant passing from the residual bone through the skin envelope to an external mechatronic limb; (b) a target implanted within at least one muscle-tendon complex of the residual limb; (c) an actuator designed to actuate the mechatronic limb; (d) a sensing array mounted to the osseointegrated implant distal to the amputated residuum comprising one or more elements designed to measure one or more components of a signal emitted by the target; (e) a first processing element connected to the sensors to process the signals from the sensing arrays and to compute muscle-tendon signals comprising at least one of force, elongation and time-derivatives associated with the at least one muscle-tendon complex; and (f) a second processing element to command an actuating signal to drive the actuator in accordance with the muscle tendon signals to replace the function of a missing biological limb of a wearer.

[0106] A neuromuscular control apparatus to augment or replace the function of a biological limb of a wearer may comprise: (a) a target implanted within at least one muscle-tendon complex of the biological limb; (b) a target implanted on a bone of the limb; (c) an actuator attached to or within the limb to augment the function of the limb; (d) a sensing array mounted to the limb of one or more elements designed to measure one or more components of a signal emitted by the targets; (e) a first processing element connected to the sensors to process the signals from the sensing array and to compute muscle-tendon signals comprising at least one of the distance between the bone target and the muscle-tendon target and time-derivatives associated with the distance between the bone target and the muscle-tendon target; (f) a second processing elementto command an actuating signal to drive the actuator in accordance with the muscle tendon signals to augment or replace the function of the biological limb.

[0107] An apparatus designed to estimate the state of at least one implanted target within the body (the at least one implanted target comprising magnetic material positioned near the apparatus) may comprise: (a) an array of sensors each designed to measure one or more target field characteristics; (b) a processing element connected to the sensors that: (i) acquires measurements of the field characteristics; (ii) applies information-theoretic processing of the measurements to estimate in real-time the state and error covariance of the at least one target state vector, the state vector comprising at least one of (1) target position, (2) velocity, (3) orientation, (4) angular rate, and / or (5) target strength; and (iii) updates the at least one of target states employing a state update innovation method that incorporates at least one of: (1) a posteriori covariance; (2) an observation matrix relating the change in a sensor characteristic to the change in the state of the at least one target; (3) a measurement variance; and / or (4) the innovation method representing a difference between the measurement and its expected value; and (c) a processing element that communicates to a human augmentation device at least one of a: (i) component of a target state; or (ii) a function of the target states that can include at least one of a (1) force, a (2) torque, or (3) an elongation, or their derivatives.

[0108] In some example embodiments, the processing element reduces the error covariance by employing at least one stochastic dynamic model of the at least one target state, the model taking the form of at least one of: a Markov model or a dynamic system represented by a set of differential or difference equations driven by one or more noise signals, an oscillator comprising time varying quadrature amplitudes and frequency with the amplitude and frequency driven by a Markov process itself driven by white noise, or a disturbance field that is varying in an inertial frame that is driven by a Markov model driven by white noise, or driven by a dynamic system represented by a set of differential or difference equations driven by one or more noise signals.

[0109] In some example embodiments, the information-theoretic processing of the measurements is accomplished by an information filter implemented where the time-varying observation matrix used in the information gain calculation is derived through real-time computation of the Jacobian that relates the change of a signal characteristic to a change in targetstate. In some cases, the information filter is, through linearization about the estimated state or other operating point, approximated by an extended Kalman filter.

[0110] In some example embodiments, the filter employs knowledge of the measurement noise characteristics, the characteristics comprising: measurement noise variance unique to at least one of the signal characteristic and to the at least one sensor component in the array; and measurement noise correlation in time.

[0111] In some example embodiments, the a posteriori state covariance is assumed to have achieved steady-state.

[0112] In some example embodiments, the Jacobian is computed at a nominal location of the target.

[0113] In some example embodiments, the innovation contribution is a sum of information vectors computed by one of a plurality of multiple processors within the processing element where each processor is connected to a sub-population of sensors and an information vector unique to a sensor is computed.

[0114] In some example embodiments, the Jacobian is computed via interpolation in relation to its value at a nominal location of the target, the type of interpolation being selected from: Linear interpolation; Bilinear interpolation; Trilinear interpolation; or an interpolation based on computation of the Hessian matrix evaluated at the nominal location.

[0115] In some example embodiments, the processing element normalizes and applies known corrections to each in the array of sensor outputs to correct for at least one of the following error sources: (a) Translation of the sensor with respect to its nominal position within the sensor coordinate frame comprising at least one of a component of a Cartesian translation; (b) Orientation of the sensor with respect to the nominal orientation within the sensor coordinate frame comprising at least one of a rotation about x, a rotation about y and a rotation about z; (c) Orientation offset of the unit vectors corresponding to the principal coordinate axes along which signal characteristics are acquired by the sensor, where the offsets are two-degree of freedom offsets represented by a rotation about two of the coordinate frame axes; (d) Sensitivity of the measurement of the signal characteristics on at least one of principal coordinate axes; or (e) Bias offset in the measurement of at least one of the signal characteristics.

[0116] In some example embodiments, the sensor array and the processing element are packaged onto at least one of one side or both sides of at least one of: (a) a rigid-flex PCBA; (b) a rigid PCBA; or (c) an interconnected rigid or rigid-flex PCBAs.

[0117] In some example embodiments, the first and second processing elements are packaged in the same assembly.

[0118] In some example embodiments, the first and second processing elements are packaged into one or more FPGAs.

[0119] In some example embodiments, the communications between the second processing element and the human augmentation device is carried out on a network.

[0120] In some example embodiments, the network synchronizes time using at least one of a network time protocol and precision time protocol.

[0121] In some example embodiments, the network is an Ethernet or EtherCAT or other time-sensitive network (TSN) with time synchronization capability.

[0122] In some example embodiments, the communications is accomplished in a multi-drop synchronous serial interface network.

[0123] An apparatus for calibrating a magnetometer array may comprise: (a) a PCBA printed with an array of addressable multi-layer coil structures that can carry current so as to create a magnetic field that can be sensed by the magnetometer array, the PCBA arranged with row and column switches to address the coil that intersects the row and column; (b) a power electronic assembly designed to deliver a calibrated current to the selected coil; (c) a processing element designed to: (i) sequence and drive row and column switch closure, (ii) modulate the coil current command, (iii) sequence the sampling of the magnetometer array measurements, (iv) process the magnetometer array measurement samples so as to compute the calibration parameters, these signature parameters comprising at least one of the following for each sensor in the array: (1) Cartesian translation with respect to the PCBA coordinate system, (2) 3x3 normalization matrix that corrects for scale and orientation variation with respect to the sensor coordinate frame on the PCBA, (3) Bias in each component of the magnetometer measurement, (4) Noise variance of each component of the magnetometer measurement, 5) Signal sensitivity and non-linearity (v) Store the resulting signature in a non-volatile memory connected to the magnetometer array.

[0124] In some example embodiments, the magnetometer array under-test is housed in an enclosure that substantially eliminates stray magnetic fields, these that may be created by at least one of the following disturbances: (a) the earth’s geomagnetic field, (b) field generated by other equipment, or (c) RF and EMI generated by other equipment

[0125] In some example embodiments, the apparatus employs: (a) an inertial measurement unit (IMU) function arranged to report the angular rate of the measurement apparatus; (b) signal processing functions within the processing element that: (i) predict temporal changes in the magnetic field owing to the instantaneous angular rate; (ii) cancels the predicted component of the magnetic field measurement changes so as to reduce the covariance of the target state estimates.

[0126] An apparatus may comprise a distributed, heterogeneous network of measurement devices each designed to estimate the state of at least one target in an array of targets within proximity to the network in its coordinate frame, the apparatus may be designed to apply information-theoretic-fusion of the measurements to determine the target state in a single unified (global) coordinate frame tied to at least one coordinate frame in the network.

[0127] In some example embodiments, the network of measurement devices comprises magnetometer arrays and the targets comprise at least one of: (a) magnetic dipoles; or (b) magnetic material with known field characteristics in relation to the geometry of the magnetic material.

[0128] In some example embodiments, the target state estimates are shared across the network to correct for far-field influences of targets so as to reduce the target state covariance, the target state comprising at least one of: (a) Target Location; (b) Target orientation; (c) Target strength; or (d) Time-derivatives of a-c.

[0129] In some example embodiments, the initial target states are determined by a global optimization algorithm that minimizes the prediction residual—the residual computed as the difference between the magnetic field measurements at each node in the network and the expected value of those measurements based upon magnet state estimates.

[0130] In some example embodiments, the initial target states are informed by simulated relationships between nominal location of magnets and their magnetic characteristics and those of the network of magnetometers.

[0131] In some example embodiments, the coordinate transformation between each node- referenced coordinate system employed in the network and a limb-referenced coordinate system is calibrated through use of a common magnetic field stimulus that can be measured at each node in the network.

[0132] In some example embodiments, the common magnetic field stimulus is applied by an apparatus with an addressable array of coils where the: (a) coil location and current amplitude is known within the coordinate system of the addressable array of coils; (b) amplitude of the magnetic stimulus is sufficient to be at any time detectable in at least two nodes in the network and where (c) the coil apparatus is moved within the proximity of the network so that each element of the network can detect its influence during the calibration sequence; and (d) the applied stimulus location and amplitude is shared across the network so that each node can determine the coordinate transformation between it and the local coordinate system of the addressable coil apparatus.

[0133] In some example embodiments, the common magnetic field is applied passively in the form of an array of at least one magnetic dipole of sufficient strength to influence the measurements of the magnetometer array in each node in the network.

[0134] In some embodiments, each node in the magnetometer array network has at least one addressable printed-circuit coil which can project a modulated magnetic field that at least one other node can detect so as to be able to determine the coordinate transformation between the two nodes.

[0135] In some example embodiments, the apparatus employs processing to calibrate a tiled array of magnetometers deployed in a fixed but non-planar arrangement, where the processing computes the coordinate transformation between each magnetometer in the tiled array to a common apparatus reference coordinate frame by measuring the magnetic field in response to a plurality of coil excitations at each member of a population of the addressable coils.

[0136] In some example embodiments, the initial vector of the at least one target state is informed by a global optimization algorithm that minimizes the prediction residuals.

[0137] In some embodiments, each node in the network can pair with at least one other node such that the coordinate transformation between any node and any other node in the network can be determined.

[0138] In some example embodiments, the on-board coil stimulus can be used to confirm calibration of the magnetometer array packaged on the board.

[0139] In some example embodiments, the on-board coil stimulus can be used as a risk control measure to confirm function during use.

[0140] In some embodiments, an information-theoretic smoothing algorithm is applied to the state estimates retrospectively and those smoothed state estimates determine a target state with substantially reduced covariance.

[0141] In some example embodiments, the rigid-flex PCBA can be used to create a tile structure where each tile comprises at least one magnetometer.

[0142] An apparatus designed to estimate the state of at least one implanted target within the body positioned near the apparatus may comprise: (a) An array of sensors each designed to measure one or more target field characteristics; and (b) A processing element connected to the sensors that: (i) acquires measurements of the field characteristics; (ii) employs a neural processing unit to determine at least one of the at least one of the target states comprising: (1) target locations; (2) target orientation; (3) target strength; or (4) time-derivatives of 1-3; and (iii) communicates the target states to an external processing element.

[0143] An apparatus employed as part of a neuromuscular controller of a device aimed at augmenting muscle-tendon function may comprise: (a) An array of sensors each designed to measure one or more target field characteristics; (b) A processing element connected to the sensors that: (i) Acquires measurements of the field characteristics; (ii) Employs a neural processing unit to determine at least one of a set of response parameters comprising muscle- tendon force, elongation or time derivatives of these based upon a trained model of the relationship between the targets and the response parameters; and (c) An actuator that communicates with the processing element where the processing unit commands the actuator to deliver the determined dynamic response.

[0144] In some example embodiments, the neural processing unit is trained by a series of synthesized Monte Carlo runs representing a stochastic model of the muscle-tendon unit dynamics and the target dynamics therein, where the stochastic model captures the range of variation in at least one of signal characteristics; magnet strength; range of magnet motion andgeometric relation between the magnets and the sensing array; disturbance state and in target state that will present to the apparatus when in use.

[0145] In some example embodiments, at least one target is implanted in the bone and the remaining two targets are mounted in the tendons distal from the bone attachment to enable at least one of the processing elements to: (a) Reference target location and velocity to a bone- centered coordinate system; (b) Compute muscle force and elongation; or (c) Communicate with the actuator to command it to deliver the dynamic response from the computation.

[0146] In some embodiments, the biological limb has been adapted for osseointegration and the sensing array is attached to the bone to minimize the relative motion between the sensing array and the biological limb.

[0147] In some example embodiments, the neural processing unit employs a convolutional neural network to process at least one of sensed target field characteristic on at least one target.

[0148] A self-calibrating apparatus for use in tracking the state of magnetic targets within an apparatus-referenced (“World”) coordinate frame may comprise: (a) At least one of an addressable coil that can be driven by a programmable current source to project a magnetic field with known geometric characteristics; (b) A tiled-array of sensors each designed to measure at least one component of the magnetic field characteristics arising from the presence of the magnetic targets and the coil, the tiled array assembled in at least one of a rigid or deformable substrate; and (c) A processing element connected to the tiled array and the coil and programmed to: (i) modulate the current in the coil so as to vary the polarity and strength of the magnetic field; (ii) demodulate the sensor signals in accordance with the coil modulation; (iii) calculate for each sensor the modulated field component while also retaining a measure of the modulation- free magnetic field signals—these with the modulated component removed; (iv) predict the magnetic field arising from the modulated field from coil using a priori information about the coordinate transformation between the coil on a sensor in the array; (v) compute the difference between the demodulated magnetic field measurement and the magnetic field prediction on the sensor—the difference serving as the innovation; (vi) compute the model-based, analytic Jacobian or other sensitivity matrix that determines the marginal impact of each component of the coordinate transformation parameters on the sensor measurement; (vii) find the coordinate transformation between the coil and the sensor that minimizes the innovation; and (viii) Use thecoordinate transformations on each element of the sensor array as the basis to estimate the state of the magnet targets within an apparatus coordinate frame—the state comprising the location, orientation and magnet strength and time-derivatives therein—and the state of the disturbance field and time-derivatives using an optimization algorithm that minimizes the difference between the modulation-free magnetic field measurements and the predicted magnetic field from the target and disturbance field estimates the above processing aimed at compensating in real-time for deformation of the tiled array while tracking magnet state.

[0149] In some example embodiments, the optimization is a global optimization.

[0150] In some example embodiments, finding the coordinate transformation that minimizes the innovation employs an information-theoretic tracking algorithm.

[0151] In some example embodiments, the global optimization is used to create the a priori information that serves to initialize the information-theoretic tracking algorithm.

[0152] In some example embodiments, the optimization algorithm that minimizes the modulation-free magnetic field from the target and disturbance field estimates is an information- theoretic tracking algorithm.

[0153] In some example embodiments, the information-theoretic tracking algorithm is initialized by a global optimization algorithm.

[0154] In some example embodiments, the tiled-sensor array and the modulated coil are integrated into the same rigid printed circuit assembly or into the same rigid printed circuit subassembly.

[0155] In some example embodiments, the tiled-sensor array and the modulated coil are integrated into the same rigid-flex circuit assembly or into the same rigid-flex printed circuit subassembly.

[0156] In an aspect, systems and / or methods may comprise a network of rigid, semi-rigid or deformable tiled-array apparatus for tracking the state of magnets in a common framework of coordinate transformations, at least one node in the network comprising an active, modulated field generator, sensing and computation that compensates in real-time for coordinate frame distortion in at least one of the tiled-array apparatus, the one or more modulated fields being used as a means to determine in real-time the coordinate frame of each node in relation to any othernode, for the purpose of tracking magnet state within proximity of the network in a common coordinate system.

[0157] In some example embodiments, the modulated fields projected by one or more of the nodes are at different frequencies, the frequencies projected by the nodes are known to the other nodes, enabling demodulating in real-time the modulated signals.

[0158] In some example embodiments, the influences of all the magnet targets in proximity to the network are integrated by using an information-theoretic fusion of measurements from all or a subset of the tiled-arrays on the network.

[0159] In some example embodiments, the oscillatory response is detected through use of an adaptive information-theoretic tracker that applies synchronous demodulation and phase detection to track the evolution of amplitude and frequency of the oscillation.

[0160] In some example embodiments, the implanted magnetic dipole orientation in pitch and roll is substantially aligned so as to minimize one or more components of the target state estimation covariance, where the target state, as referenced to a coordinate frame, is defined by at least one of: (a) Target location along Cartesian x, y and z; (b) Target velocity along Cartesian x, y and z; (c) Angular rotation in pitch and roll; or (d) Time-derivatives of a) and b) above.

[0161] In some embodiments, the global optimization is accomplished by gradient descent and not a Jacobian.

[0162] In some example embodiments, the first processing element: (a) detects the lateral vibration of at least one of the tendon-mounted targets and where the first processing element employs an information filter processing function comprising an adaptive phase-lock loop that estimates the instantaneous vibration amplitude and frequency associated with a muscle activation signal; and (b) incorporates the muscle activation signal into the muscle-tendon kinetics calculation.

[0163] In some example embodiments, the adaptive phase-lock loop employs is implemented in an information-theoretic context using an information filter that incorporates a white-noise driven stochastic model of the vibration dynamics.

[0164] In some example embodiments, the stochastic model is a Markov model.

[0165] In some example embodiments, the innovation update is computed as the product of an innovation gain matrix and the innovation.

[0166] In some example embodiments, the sensing array employs a PCBA employing buried-capacitance onto which sensors are mounted on one or more sides.

[0167] In some example embodiments, the PCBA is at least one of a rigid, rigid-flex or stacked, interconnected PCBAs.

[0168] In some example embodiments, the apparatus includes a reference magnetic dipole to verify and adjust the calibration.

[0169] In some example embodiments, the addressable coil is modulated at fixed, reference frequency and the sensor measurements are demodulated using the reference frequency to compute signature parameters not associated with sensor bias.

[0170] In some example embodiments, the second processing element senses when the wearer is inactive and reduces battery power consumption in at least one of the actuator, sensing array, or the first processing element.

[0171] In some example embodiments, the reduced battery power consumption is achieved by at least one of: (a) Disconnecting the battery voltage from select electronics subsystems; (b) Setting FPGA and processor electronics in at least one of an Idle or Standby state; or (c) at least one of the first or second processing elements monitors lack of substantive changes in at least one of a magnetometer or IMU measurement that signals wearer inactivity.

[0172] Each example embodiment is a specific description of one or more broader concepts as will be evident to those skilled in the art. Further, the embodiments set forth in the appendix and throughout the application may be utilized individually or be combined in any combination.

[0173] 1. Information-theoretic Magnet and Muscle-Dynamics State Estimation Methods

[0174] In some embodiments, the efferent and afferent channels of the magnetomicrometry- based neuromuscular controller—installed on a wearer of a human augmentation device (Wearer)—are formally treated as a bi-directional communication channel shown in FIG.1. Information-theoretic methods are invoked to perform signal processing operations involved in the efferent-afferent communications.

[0175] FIG.1 shows, for the case of limb amputation, a magnetomicrometric-based Neuromuscular Controller modelled as a bi-directional communications channel in an Information-Theoretic Context (Diagram not to scale). The magnetic fields arising from magnets implanted within muscle-tendon units adjacent to the magnetometer arrays are viewed as theinput to an efferent communication channel in a neuromuscular controller for a wearable robot or for the artificial stimulation control of living muscle tissue. The first magnetometer signal processor in the efferent channel corrects for scale, bias, and coordinate-frame irregularities in each sensor integrated circuit (IC) and calculates the magnet and disturbance state. The second processor uses magnet state to compute muscle-dynamics information that captures the wearer’s volitional intent, or for the acquisition of muscle-tendon signals for use in a feedback control loop for muscle tissue control. For the application of augmentative wearable robotics, the muscle-dynamics information commands the robotic dynamics through a transformation from the linear muscle-tendon space to the rotary space of the wearable robotic joint. The outputs of this transformation are control targets such as joint torque and / or displacement including time- derivatives of these—to the multi-degree of freedom (DOF) robotic actuator. Cutaneous and proprioceptive afferent feedback signals from the wearable robot into the nervous system of the user comprising artificial muscle stimulation are delivered by ancillary hardware. The neuromuscular controller function is critical to achieving volitional integrity and afferent, reflexive feedback in this bi-directional communication channel.

[0176] Some embodiments, such as the one shown in FIG.1 include three magnetic beads. Alternative embodiments may include a different number of magnetic beads.

[0177] Some embodiments described herein use information theory to guide the design— including number of magnetometer ICs, their placement, number of magnets, the orientation of the magnets, the signal processing and control algorithms—to maximize the precision of the magnet and muscle-dynamics state estimates.

[0178] From a mathematics perspective we equate precision with the information matrix, Sk|k, where: Sk|k= ^k|k-1Eq 1.1

[0179] Where Xk|kmeans the value of X at t=tkgiven all information acquired up to and including t=tkand ^k|kis the state covariance evaluated as:

[0180] Large values in the diagonal elements in the covariance matrix imply a lack of precision in those state estimate components. Conversely, large values in the diagonal elementsin the information matrix, Sk|kimply a high degree of precision in the state estimates of magnet location and velocity.

[0181] Note: The expected value notation—E{x}—shall be taken to mean the expected value of a variable evaluated over the probability, p(x). The “hat” symbol over a variable shall be taken to be equivalent to the expected value notation.

[0182] As shown in Eq 1.2 (see FIG.2), the measurement, zik, comprising the magnetic field measurements at time k from sensor, i is a non-linear function, h(xk) of the vector of composite magnet states, xjk, and is corrupted by additive stationary white noise, vkat time k with known measurement variance, Ri. See Appendix A, Eq.8-15.

[0183] When evaluated at an operating point, x0k, the measurement variation, ^zkrelates to the variation of state, ^xk, as defined in Equation 1.3 (see FIG.3).

[0184] Note: The E{vikvip}=0 ^ constraint ensures that the measurement noise is “white”— meaning that it is not correlated in time. That’s a fundamental assumption in the optimality of the information filter. Also note that the measurement variance has a sensor dependance. Tests have shown the measurement variance across measurement components (x,y,z) and sensors is quite large. So, discounting of measurements of high variance is an important implementation aspect to maintain optimality.

[0185] ^hi(x) is the Jacobian which can be evaluated analytically per Appendix A. (See Appendix A Equations 8-15) The optimal least-squares a posteriori estimate of ^xkhas an information matrix, Sk|k, defined by the recursive relationship in Equation 1.4 (see FIG.4) and a recursive relation for the information vector defined as in Eq 1.5 (see FIG.5).

[0186] Equation 1.6 (see FIG.6) shows that in the information-theoretic framework, the measurement update, at any time, k, from the N sensors in the magnetometer array (MMA), is simply the sum of the information matrix and information vectors respectively from each sensor.

[0187] This is a powerful result because it suggests that the update of magnet state information can be distributed amongst an array of dedicated “parallel” processors and accumulated by a second processor element to determine the collective magnet state and disturbance state. In neuromuscular control systems that are “communicating” with a large number of muscle-tendons—and so many implanted magnets—distributed processing and hencehardware and firmware scalability is easily achieved without sacrifice in sampling time or hardware complexity.

[0188] In magnetomicrometry applied to tracking of magnet state implanted in muscle- tendons which exhibit largely uncertain dynamics, magnet state information (precision) is not increasing without bound. Here, white-noise driven Markov dynamical models of the form in Equation 1.7 (see FIG.7) can be used for each component of magnet and disturbance state— these serving as an information-theoretic proxy for the stochastic muscle-tendon dynamics. Specifically, for the information-theoretic methods being deployed for use in the next generation of magnet tracker, a simple, two-state Markov model implemented as a series of low-pass filters each with a time-constant, ^ is employed. The time-constant is selected to be a “characteristic” time for the muscle-tendon movement. When the Markov time-constant is larger, the information filter can, in a sense, “average” more of these time-correlated measurement time samples, thereby increasing precision. The Markov model input noise variance is sized to create a probability distribution that matches the range of motion anticipated for the muscle-tendon dynamics, and hence implanted magnet, motion. By using this form of second-order model, the magnet velocity can easily be determined. Equation 1.7 presents the continuous-time and discrete-time second-order, Markov model formulation.

[0189] From an information-theoretic point of view the information deteriorates on every time step before a measurement update is performed as a result of the white-noise driven dynamics. So, Sk|k-1decreases before the measurement update is performed. The measurement update serves to increase the information driving the a posteriori information matrix to increase. The disturbance and measurement effects on the information ultimately balance and the information precision achieves steady-state—a fact that will be leveraged in the “rule of thumb” simulation results reported below and in the algorithms that are presented in Section 1.1 for both the Information Filter and Deep Learning algorithms. The Information Filter is a convenient way of implementing the Kalman Filter. See [1] for details.

[0190] 1.1. Model-based Magnet and Muscle-Dynamics State Estimation

[0191] The above learnings are applied to improve precision and scalability in both the application of distributed filtering and smoothing methods and Deep-Learning methods for Magnet and Muscle-Dynamics State estimation.

[0192] 1.1.1 Distributed Single Magnetometer Array Filter and Smoothing Methods

[0193] 1.1.1.1 Extended Kalman Filter State Estimation and Smoothing Methods

[0194] In the literature, the state estimator employing non-linear model linearization method proposed here is attributed to R. Kalman and is referenced as the “Extended Kalman Filter.”

[0195] 1.1.1.1.1. Magnet and Muscle-Dynamic State Estimation

[0196] FIG.9A shows a Magnetometer Array Information (MMAI) Processor Overview. FIG.9B shows Magnetic Field Information Vector Processor (MFIVP) Architecture Details. FIG.9C shows Information Aggregator Processor (IAP) Architecture Details.

[0197] FIGS.9A-9C illustrate the filter architecture that estimates (tracks) magnet state in real-time for a single MMA. When the neuromuscular controller initializes, a global optimization algorithm locates the magnets in quiescent conditions with some known degree of precision, S0, with an accuracy and precision equivalent to or better than that precomputed steady-state information matrix. Or, in other embodiments the nominal position of the magnets may be known and saved in memory for use at initialization by the tracking algorithm. From that point forward, the tracking algorithm computes the magnet and disturbance state using the linearized information filter model without the need for a global optimization algorithm step.

[0198] As shown, an array of magnetic field information vector processors (MFIVP) — which in some high-performance embedded processors are configured as concurrent or serialized threads and in others are instantiated as digital signal processors embedded in one or more Field- Programmable Gate Arrays (FPGA) or ASICs—to perform sensor data acquisition on a subset of the sensor population, where in low sensor density applications this subset may be the entire population. The MFIVP processors employ a dedicated shared memory that serves as the portal for bi-directional communication with an Information Aggregation Processor (IAP) that accumulates the information vectors from the distributed array of MFIVP. The MFIVP processors are digital signal processors (DSP) that are designed for efficient execution of the dot product computations—often of an integer type to reduce execution time—employed in executing measurement normalization and information vector computations. As shown in FIG. 9B, the first post-acquisition process step subtracts the sensor-specific bias as determined during a manufacturing calibration. The second process step multiplies the bias-free signal vector by a sensor-specific 3x3 matrix, βI, it also determined during a manufacturing calibration process stepthat corrects for scale, non-linearity (distortion), and coordinate frame irregularities. Here, for instance, the x-unit vector inside the between magnetometer IC may deviate from the x unit vector of the magnetometer array coordinate system, say, on a rigid printed-circuit board assembly (PCBA) substrate. Once the sensor signal has been normalized, the innovation vector, δzikis computed as difference between the measurement, bik, and the predicted measurement hi(E{xk}) as shown in FIG.9B.

[0199] Note: We will employ the term “innovation” and “information vector” synonymously for the remainder of this disclosure document.

[0200] The innovation is multiplied by a variance-scaled Jacobian matrix, Hi(E{xk} R,-1that we refer to as the information gain, Hik*, this evaluated at the expected value of xk|k-1, to yield the sensor-specific, information vector, iik|k. An Information Aggregator Processor (IAP) –this a processor or code thread within a processor— through information fusion, accumulates the information vectors from all of the sensors then multiplies the sum by the inverse of the steady- state Information Matrix, Sss-1–this equivalent to the steady-state state covariance—to yield the innovation vector update, E{xk|k}, that is added to the a priori state vector estimate, E{xk|k-1}. As shown in neuromuscular controllers where the muscle-dynamics is a linear function of the magnet target location and velocity (See 1.1.1.1.3), the IAP multiplies the magnet component of the state vector (the state vector with the disturbance states removed) by a pre-computed, application-specific, muscle-dynamics observation matrix itself the product of an output transformation matrix, C, that transforms the magnet state vector into a vector of magnet locations and velocity; and a matrix, ^, that transforms the resultant vector into a concatenated muscle dynamics vector, λ, comprising an application-specific concatenation of forces, torques, elongation and time-derivatives. The IAP then communicates the vector, ^, and the + / - n sigma confidence intervals based on the covariance of , to the robotic actuator, ormuscle tissue actuator. In the case of a wearable robotic controller, the robotic actuator uses this information to deliver a control response representative of the muscle-tendon dynamics of an intact limb or other tissue structure. In the case of a muscle controller, the muscle actuator uses this information to deliver a control response to modulate, in an updating manner, a muscle- tendon dynamic response. For instance, if the IAP communicates a muscle-tendon torque contribution, the actuator controller may deliver a motor torque in accordance. In orthotic orexoskeleton-type augmentation devices, a muscle-tendon velocity may be used to gate the delivery of torque to augment the torque delivered by the wearer’s tissue.

[0201] In summary, using the information-theoretic-inspired magnet state estimation method we create scalable, magnetic-field-to-magnet state processor architectures that can easily be distributed. The method builds upon and can coexist with the already-devised global optimization methods that may be used to determine the initial conditions—the state vector comprising the initial magnet state, including location and velocity—from which the tracker can begin execution. Using that initial state and application-specific knowledge relating to the Markov time-constant and driving noise variance the steady-state Information Matrix for the application can be computed along with the sensor-specific, information gain, Hi*, evaluated at the initial condition and loaded into the shared memory of the distributed MVIFP processors or processing threads. In applications that require adjustment arising from extreme geometry changes along the magnet-specific trajectory that invalidate the assumption of uniformity of Hi*at a nominal location, the elements of Hij*for a magnet, j, can be computed by the IAP on every time step to account for the expected value of changes in magnet position or can be approximated by the MFIVP as a tri-linear, bi-linear, linear or analytically derived Hessian-based non-linear function of the magnet location variation from nominal. In these Hij*approximations, the IAP can load the pre-computed coefficients of the approximation calculations into a shared memory along with the state estimate of the magnet location variation.

[0202] FIGS.10A-10C show numerical calculations associated with adjusting Hij* to account for large magnet excursions from the nominal operating point that invalidates use of a Jacobian excursion-invariant value. FIG.10A illustrates the changes in signal processing required for real-time adjustment of Hij*. As will be discussed in Section 3, the MFIVP can be implemented using special-purpose, integer multiply-add implemented in FPGAs, ASICs or in threaded software code—this to speed execution for applications involving large state vectors arising from a large number of implanted magnets in proximity to the MMAI, or a network of MMAI mounted on a limb.

[0203] 1.1.1.1.2. Smoothing of Magnet State and Muscle-Dynamics Estimates

[0204] Smoothing of the magnet state estimates and the consequent muscle-dynamics vectors can be accomplished by straightforward application of the commonly used Kalman filtersmoothing algorithms [1]. Here, we create a pipeline of smoothed estimates with a “fixed lag” of N samples. N is selected to minimize the phase delay in communicating magnet state and muscle-dynamics information. At the nominal 1 msec sampling time and a fixed lag of N=8, application of smoothing algorithm increases the information matrix precision by roughly sqrt(N)=2.8, and only introduces a phase delay of 8 msec. If smoothing is performed, the IAP performs these computations.

[0205] 1.1.1.1.3. Application-specific Muscle-Dynamics State Estimation Examples

[0206] 1.1.1.1.3.1. Muscle-Dynamics and Activation State Estimation on a Simple, Instrumented Muscle-Tendon Structure—An Example

[0207] The schematic in FIG.10A captures the magnet implantation and sensing paradigm suitable for determining muscle-tendon force and elongation. Here, one magnet (Magnet 1; M1) is implanted in a bone (Bone 1) that attaches to one end of the muscle-tendon complex. In one embodiment, the other two magnets (Magnet 2 and 3 respectively; M2, M3) are implanted into tendon tissue near the muscle-tendonous junction of each origin and insertion tendon associated with the muscle-tendon structure.

[0208] FIG.10B shows a Tendon Force, ftendon, and Muscle Length, ^lmuscle, derivation. Here the force of the muscle-tendon unit is computed from the measured elongation of a tendon about its nominal length using the positions of magnets 1 and 2, as well as information on tendon modulus. Muscle belly length is measured from the positions of magnets 2 and 3 implanted within or on the origin and insertion tendons, respectively, at the muscle belly output locations near the musculotendinous junctions.

[0209] Using the global optimization software, the homogeneous coordinate transformation,MMABoneT , between the bone-implanted magnet and the magnetometer array can be easily computed as an initial condition with an acceptable precision implied by its Information Matrix. Since there is no “driving noise” that introduces information matrix precision reduction, successive application of the magnet filter increases precision ofMMABoneT to any desired level. OnceMMABoneT is calculated, the Magnet 2 and 3 states are computed in the coordinate system of the bone (Magnet A). Equations 10 (see FIG.11) capture the Markov disturbance dynamics; the Information Matrix and Vector computation Sk|kand sk|k; the transformation from Information Matrix to the magnet state vector, yk|k; and the muscle-tendon dynamics to compute^ and the 1- σ confidence interval for λ. Here, the Information Matrix is assumed to have reached steady-state. Because the computations are referenced to a static magnet frame, the precision implied by Sk|k=Sssis equivalent to that of a two-magnet (2DOF) rather than a three- magnet (3DOF) application.

[0210] Extending the above to an agonist-antagonist application referenced to a first bone, only five magnets are needed in the implementation and the precision will match that of a 4-DOF rather than 5-DOF application.

[0211] 1.1.1.1.3.2. Integrated Muscle Activation Estimation via Information-theoretic Phase- Lock Loop—An Example

[0212] It is well known that muscle activation can excite lateral vibration (resonance) in the tendons attached to the muscle and that this vibration is related to muscle force generation and contraction. The vibration amplitude is on the order of 50-100 um and the frequencies range in the 30-50 Hz range. The vibration frequency drifts from low to high as the generated force increases.

[0213] (Note: The resonant frequency of a string under tension will vary as the square root of T / μ, where T is the tension and ^ is the mass density along the string.)

[0214] Equations 1.11 (see FIG.12) model a two-magnet Information Filter —essentially a synchronous demodulator--designed to estimate the vibratory state—comprising the frequency and amplitude of the muscle activation-induced lateral vibration. In one embodiment, one magnet (Magnet 1) is implanted in the bone, the other (Magnet 2) in the center of the tendon. Here, the quadrature amplitudes of the vibration on two axes, x and z, are each driven a stochastic models—in this case each by a single-state Markov process simulating Brownian motion—and the phase noise variation is driven by a two-state Markov process—simulating the integral of a Brownian motion noise process. The deflection in y is assumed to be a two-state Markov process with time-constant, τ. The magnetic field disturbances are ignored to simplify the equations. These of course could be incorporated without any substantive impact on the signal processing implementation.

[0215] Note: The lateral oscillations are assumed to lie only in the x-direction in the MMAI coordinate system—Information vector solution for the variable amplitude and frequency oscillator is derived here.

[0216] 1.1.1.1.3.3. Accommodating disturbance field variation while the Wearer is walking and changing directions

[0217] Here algorithmic improvements are considered that create a slowly-varying disturbance field as might manifest from geomagnetic fields, local fields generated by motors, generators and other equipment that might contain a residual magnetization. Even a static field will cause large disturbance transients on the order 0.5-2 Gauss when a Wearer is moving and changing direction. In this application, an inertial measurement unit is employed in the magnetometer PCBA hardware and used to update the orientation matrix embodied within the PCBA coordinate frame by Equations 1.13. By doing this, we reduce the amplitude and transient frequency thereby more nearly maintaining the precision implied by the information matrix during the state propagation step prior to the next measurement. Equations 1.13 (see FIG.13) capture an application with a single-magnet with disturbance field using two-state Markov noise models.

[0218] 1.1.1.1.4. Networked Magnetometer Array Filter and Smoothing Methods

[0219] The above magnet state estimation methods apply directly to networks of magnetometer arrays packaged in proximity to a biological limb. In such a network, the active, modulated-field, calibration infrastructure (See Section 2) creates a framework to determine a common coordinate system for the implanted magnets and for all nodes in the network to contribute, via their information vectors, to the information filter magnet state estimates. The Super-IAP node serves to integrate all the information vectors from all nodes to accomplish this. Figure 4 illustrates the network architecture employing this information fusion architecture.

[0220] FIG.14 shows an example embodiment of an algorithm deployment in a “Network of Magnetometer Arrays.” As shown, one of the magnetometer array assemblies serves as a container for the Super-IAP—the “aggregator of information aggregators.” The Ethernet network employs a publisher-subscriber paradigm where the nodes subscribe to the Super-IAP magnet / disturbance state estimates that are broadcast on the network. In this way, the state vector estimate is shared amongst all the nodes in the network. In turn, the Super-IAP subscribes to the partial information vector sums published by the nodes. This “time-sensitive” network (TSN) employs the Precision Time Protocol (PTP; IEEE 1588-2019) so that all published information is time-tagged with sub-microsecond resolution and synchronization. The equationsexecuted within this network architecture match exactly that which would be executed in the distributed architecture within a single magnetometer array. Smoothing works identically to that in a single node network. In this way the architecture scales to estimate the magnet and muscle- dynamics state for a large number of implanted magnets and takes advantage of the algorithm nuances afforded by the application-specific processing highlighted in Sections 1.1.1.1.3.

[0221] (Note: The network could also be implemented as an EtherCAT, Ethernet IP, or TSN network in which the Super IAP is the master with a network update rate of 4-8kHz to ensure that the information vector transfer occurs with sufficient time for all the vectors to be summed by the end of the 1 msec cycle. In other embodiments, the IAP could host a synchronous serial network to accomplish the above.

[0222] 1.1.1.1.5. Magnet State Tracking Using Kalman Filter on Magnetic Dipole Manifold

[0223] As a further improvement of magnet tracking method, we propose a novel theoretical framework that utilizes recent advancements in Riemannian geometry and manifold theory to enhance the magnet tracking technology stack for Magnetomicrometry and other relevant applications. By representing the true state of the magnetic dipole moment on a manifold and incorporating its dynamics, we aim to develop a more accurate and robust magnet tracking algorithm that can handle the complexities of real-world scenarios. [Yeon, Seong Ho. Enhancing Muscle Sensing Modalities for Advanced Bionics. Massachusetts Institute of Technology, 2023.]

[0224] By exploiting the geometric structure of the manifold and the physical properties of the magnet, we aim to develop a more accurate and robust tracking framework. The proposed algorithm consists of three key components: (1) representing the magnetic dipole state on the S2 manifold, (2) lifting the S2 manifold to the parallelizable SO3 manifold by incorporating the magnet's rotational dynamics, and (3) applying the Kalman Filter on the resulting manifold to estimate the magnet's pose and orientation.

[0225] A few methods to utilize manifolds structures on Kalman filters were recently developed and reported. Mathematical variants of Unscented Kalman Filter (UKF) and EKF methods with manifold structures in its state variables have recently developed. UKF is preferred to be used for magnet tracking with its ease of the implementation, while other variants including EKF with manifolds can be used with the proposed tracking framework[M. Brossard, S. Bonnabel, and J.-P. Condomines, “Unscented kalman filtering on lie groups,” in 2017 IEEE / RSJInternational Conference on Intelligent Robots and Systems (IROS), IEEE, 2017, pp.2485– 2491], [M. Brossard, A. Barrau, and S. Bonnabel, “A code for unscented kalman filtering on manifolds (ukf-m),” in 2020 IEEE International Conference on Robotics and Automation (ICRA), IEEE, 2020, pp.5701–5708.], [D. He, W. Xu, and F. Zhang, “Kalman filters on differentiable manifolds,” arXiv preprint arXiv:2102.03804, 2021]

[0226] Through this approach, we address the limitations of existing magnet tracking methods and provide a principled way to handle the complexities of real-world scenarios.

[0227] 1.1.1.1.5.A Magnetic Dipole Manifold Representation: Finding the magnetic dipole vector is analytically identical to finding a point on sphere with given radius; the solution space to describe the dipole moment vector thus lies on an analogous surface of a virtual sphere model. This solution space, which is a surface of a sphere, is called the S2 manifold [A. Alonso Izquierdo, G. Leon, J. Mateos Guilarte, and M. de la Torre Mayado, “On domain walls in a ginzburg-landau non-linear s2-sigma model,” Journal of High Energy Physics, vol.2010, no.8, pp.1–29, 2010.]. We propose to use this S2 manifold to represent the magnetic dipole state or magnet’s mechanical state.

[0228] 1.1.1.1.5.B Magnetic Dipole Manifold State Lifting and Parallelization: The S2 manifold is not parallelizable. Parallelizability is an essential property for the existence of a homogeneous transformation between the original manifold and the homogeneous tangent space [M. Brossard, A. Barrau, and S. Bonnabel, “A code for unscented kalman filtering on manifolds (ukf-m),” in 2020 IEEE International Conference on Robotics and Automation (ICRA), IEEE, 2020, pp.5701–5708.], [D. He, W. Xu, and F. Zhang, “Kalman filters on differentiable manifolds,” arXiv preprint arXiv:2102.03804, 2021]. Moreover, the homogenous transformation of manifold structures is the key function to augment the KF with the manifold. In order to utilize the manifold structure on tracking paradigm, we propose an approach to lift the S2 manifold representing the magnetic dipole moment to a parallelizable manifold structure.

[0229] By lifting the S2 manifold to the SO3 manifold, we can exploit the parallelizability of SO3 and apply Kalman filters effectively. This lifting is achieved by constructing the SO3 manifold using the direct product of the S2 and S1 manifolds [F. Ding and H. Geiges, “The diffeotopy group of s1× s2 via contact topology,” Composition Mathematica, vol.146, no.4, pp. 1096–1112, 2010] [M. Brossard, A. Barrau, and S. Bonnabel, “A code for unscented kalmanfiltering on manifolds (ukf-m),” in 2020 IEEE International Conference on Robotics and Automation (ICRA), IEEE, 2020, pp.5701–5708.] While the S2 x S1 manifolds is not strictly homogeneous to SO3 manifold, we can construct physically equivalent SO3 manifold using the right-hand rule. Here we propose to represent magnetic dipole using SO3 manifold by exploiting state lifting approach.

[0230] In physical world, we can define the S1 vector with velocity vector on the surface point of the magnetic dipole surface created by rotation of the magnet. With additional R3 manifold representing amplitude of rotation, this surface velocity vector from surface makes the system uniquely observable, and thus, we can uniquely define spatial pose of the magnet using SO3 manifold.

[0231] Implication of this approach suggests that we can define the unique magnetic dipole state by incorporating rotational movement of the magnetic dipole. Then, with rotational dynamics of the magnet, it is also feasible to define generic system dynamics for process modeling. Likewise, we also introduce a translational velocity state to model thorough dynamics of the magnet. The following table summarizes dynamic state lifting to enhance system state. For both dipole moment and position manifolds, I provided generic forward dynamics update method for proper update routine in KFs.

[0232] 1.1.1.1.5.C: Kalman Filter with Magnetic Dipole Manifold

[0233] To utilize the defined magnet manifold structure in the KF framework, we propose to establish a mapping between the manifold space and the tangent space. This mapping is achieved through the retraction and inverse retraction operations [M. Brossard, S. Bonnabel, and J.-P. Condomines, “Unscented kalman filtering on lie groups,” in 2017 IEEE / RSJ International Conference on Intelligent Robots and Systems (IROS), IEEE, 2017, pp.2485–2491], [M. Brossard, A. Barrau, and S. Bonnabel, “A code for unscented kalman filtering on manifolds (ukf- m),” in 2020 IEEE International Conference on Robotics and Automation (ICRA), IEEE, 2020, pp.5701–5708.].

[0234] Retraction operation maps a point from the homogeneous tangent space to the manifold, while inverse retraction maps a point from the manifold to the homogeneous tangent space. These operations allow us to perform the prediction and update steps of the KF in the homogeneous tangent space while ensuring that the estimated state remains on the manifold.

[0235] The following table presents the retraction and inverse retraction operations for the SO3 and the R3 manifolds used in our proposed algorithm. For the SO3 manifold, the retraction operation is performed using the matrix exponential (Exp), and the inverse retraction is performed using the matrix logarithm (Log) [D. He, W. Xu, and F. Zhang, “Kalman filters on differentiable manifolds,” arXiv preprint arXiv:2102.03804, 2021]. For the R3 manifold, the retraction and the inverse retraction operations are simply the addition and subtraction of vectors, respectively.

[0236] As one example of using the proposed algorithm, UKF on manifolds (UKF-M) algorithm can be used [M. Brossard, A. Barrau, and S. Bonnabel, “A code for unscented kalman filtering on manifolds (ukf-m),” in 2020 IEEE International Conference on Robotics and Automation (ICRA), IEEE, 2020, pp.5701–5708.]. The algorithm adapts the standard UKF to operate on the manifold by using the retraction and the inverse retraction transformations. The prediction and update steps are performed in the homogeneous tangent space, while the sigma points are generated and transformed using the manifold transformation operations.

[0237] By leveraging the manifold structure of the magnet state representation and the Kalman Filter variants algorithm incorporating the generic manifold, our proposed method effectively incorporates the dynamics of the magnet and provides a principled way to estimate its pose and orientation. The retraction and the inverse retraction transformations ensure that the estimated state remains consistent with the geometric constraints of the problem, leading to improved accuracy and robustness in magnet tracking.

[0238] 1.2. Model-based, Deep-Learning Magnet and Muscle-Dynamics State Estimation

[0239] As described earlier, the neuromuscular controller may be used as a key component in the efferent communications driving actuators realizing volitional intent. The magnetometer array serves as a magnetic field image acquisition and processing element that extracts real-time, high-fidelity muscle state information. This section describes the biophysical model-based design of the Deep Learning Array Information Processor (DLAIP)—a neural processing unit built upon deep learning design principles. The block diagram in FIG.15 illustrates the model- based architecture of the DLAIP.

[0240] As shown, a multi-level, convolutional neural network (CNN) extracts the muscle- dynamic state from normalized magnetic field images delivered by one or more MFIVP. The CNN is of a feedforward type and comprises multiple layers (typically four or more) employing leaky Rectified Linear Unit (ReLU) activation. The CNN is trained by Computer-Aided Engineering (CAE) process informed by a biophysical model that captures nominal implanted magnet array geometry, kinematics and kinetics as created via CAD models or other 3-D images of the Wearer muscle-tendon structures that include the implanted magnets. From this CAD data, a nominal magnetometer array placement is determined either from standard rigid-flex PCBA components or customized ones developed specifically for the Wearer. A Monte Carlo simulation model is built using the CAE-developed signature and randomized training and verification sequences are created—the randomization achieved via probabilistic-derived variation in: a) Nominal placement and fit of the MMAI in relation to the array of magnets; b) Nominal disturbance and magnet state; c) Disturbance variation driven by Markov or other stochastic muscle models that create randomized disturbance and disturbance rate-of-change representative of muscle dynamics; d) Magnet state variation driven by Markov or other stochastic models that create randomized magnet displacement and velocity vectors in accordance with the biophysical models of the muscle-tendon dynamics; and Note: In persons with amputation, the muscle agonist-antagonist structures that are created in the surgical process are modeled using biophysical principles and then transformed into Markov or other stochastic models.e) Biophysical input-output model-based variation—for instance tendon stiffness and elongation during muscle activation—that determine the dynamic vector outputs in response to the magnet state variation introduced in c).

[0241] The probability density functions (PDF) used to select a random sample from each variation distribution are typically uniform, Gaussian or mixed-Gaussian—these selected with the intent of training the CNN on “edge” cases that are possible to occur in multiple years of use across many use-case scenarios.

[0242] Millions of multi-second simulation sessions (typically 10 seconds each) are thereby created from selected operating points arising from probabilistic variations introduced by a)-e) above. Each session generates input / output sequences—typically in 1 msec sample intervals— the inputs comprising acquired noisy magnetic field images and the outputs comprising muscle- dynamic (truth model) vectors that that serve as a proxy for the signals that will be used to drive the Wearer-specific actuator used for neuromuscular augmentation.

[0243] In FIG.16A-16F simulation of a post-training DLAIP validation run is demonstrated for a “single” implanted magnet at a depth of 2.5 cm as might be used, for example, to estimate a tendon force achieved by tendon displacement relative to a magnet implanted in a bone. The x, y, and z tracking performance in relation to the magnet target are summarized in the graphs. The raw neural net output is plotted along with a 40 Hz IIR filter. The tracking error RMSE tr for the validation run illustrated is 75 ^m. As shown, the RMSE of the filtered output is only marginally better. In contrast, the corresponding tracking error RMSE for the information filter RMSE at this depth is 20 ^m.

[0244] FIGS.17A-17F simulate the x, y, and z tracking results for a post-training DLAIP validation run for a single implanted magnet trained at a nominal depth of 5 cm. The tracking error RMSE is 178 ^m. In comparison, the information filter prediction—an information floor— is 161 ^m.

[0245] 2. Modulated-Field Calibration Framework

[0246] The advances described in Section 1 are motivated by the need to create 1) an objective, information-theoretic measure of precision, and 2) a distributed and scalable Information filter framework that delivers optimal (least-squares) magnet and muscle-dynamicstate estimates in real-time and that can be adapted to use the DLAIP technology. Together these advancements close the gap in achieving precision of the state estimates.

[0247] This Section addresses how to close the gap in accuracy with an eye toward maximizing MMAI gauge repeatability and reproducibility (Gage R&R; Gage R&R is a term employed in Measurement Systems Analysis—a methodology used to define the amount of variation in measurement data due to the measurement system. www.qualityone.com) in manufacturing and in-use. The salient of this advancement lies in the use of programmable and addressable active coil circuitry employed in test fixtures and integrated within MMAI—the active coil circuitry projecting a magnetic field of known characteristics that can be used for manufacturing and in situ, real-time field calibration. This advancement leverages calibration methods detailed in Appendix A and captured in the two patents already filed—these methods employing the global optimization framework to process MMAI magnetic field imagery arising from manipulating one or more magnet dipoles in proximity of one or more MMA.

[0248] 2.1. Board-level Test and Calibration

[0249] Here we describe how the active calibration infrastructure is used for printed-circuit board assembly (PCBA) test in a manufacturing process. The block diagram in FIGS.18A-18C illustrates how a calibration test fixture employing the modulated field framework can be implemented.

[0250] The test fixture comprises: 1) MMAI Unit-under-Test (UUT) 2) NIST-traceable Current Monitoring instrument (Note: NIST stands for the National Institute of Standards and Technology; https: / / www.nist.gov). 3) Magnetic Field Generator PCBA A multi-layer PCBA comprising a rectilinear matrix of row / column addressable coils printed on multiple circuit layers that serves as the NIST-traceable reference field stimulus upon which the calibration is based. (Note: In a typical embodiment, the printed-circuit assembly might have as many at twenty or more layers to achieve the “amp-turn” field magnification.) 4) Programmable Current Source PCBA A programmable, calibrated current source that delivers the modulated coil current5) Test Controller PCBA A test controller assembly that: a. Measures the magnetic field bias on each of the three sensor outputs from each magnetometer IC (Sensor); b. Calibrates the coil current command using the NIST-traceable current monitoring instrument as the basis; c. Activates the row / column address switches (optically-coupled FET switches) in sequence to direct a specific coil magnetic field stimulus in a randomized, round- robin fashion until all coils have been stimulated and the calibration results confirmed (It is important to randomize the coil activation sequence to eliminate a possible signal correlation in time with a correlation in space); d. Modulates the coil current setpoint at a high-frequency (~125 Hz) that is not a multiple of AC line frequency (50 Hz or 60 Hz); e. Acquires the magnetic field image from the UUT; f. Demodulates, via information filter, the modulated field signals from each IC during each coil activation in the coil activation sequence; g. Processes the demodulated magnetic field image to determine for each magnetometer IC in the UUT: i. Centroid of the IC sensing circuitry (x0, y0 and z0) within the UUT coordinate-system; ii. 3x3 scale and magnetometer IC sensing axis normalization matrix that transforms the raw, digitized (quantized) magnetometer IC sensor outputs into a normalized, distortion-free (non-linearities removed) magnetic field vector—with bias and distortion removed—and aligned with the UUT coordinate-system; and iii. Measurement variance on each component of the magnetic field measurement Note: The bias cannot be measured by the demodulated image. It is the first calibration parameter to be measured prior to modulation.h. Loads the calibration parameter “signature” into a non-volatile memory on the UUT; i. Applies a spot verification of a randomized sequence of coil activations to confirm that the information vector is within tolerance; and j. In a separate step, checking that when a magnet is placed in a precise and accurate fixture proximal to the UUT that the magnet state (location, magnet strength) is within a specified tolerance.

[0251] In the preferred test fixture embodiment, the Magnetic Field Generator PCBA and the UUT are housed in a mu-metal box to substantially isolate it from stray magnetic fields in the manufacturing environment.

[0252] The active current modulation is the central idea that serves as the basis for some embodiments. The demodulated signal used in the calibration eliminates the impact of baseband signal disturbances on the calibration—these arising from unblocked geomagnetic fields, environmental variation in power, temperature fluctuations and EMI. By using information filters implemented in the Test Controller to demodulate the Sensor signals we achieve a high- level of precision in calculating the calibration parameters thereby ensuring repeatable and reproducible results.

[0253] Some details in the design of the Magnetic Field Generator PCBA are critical. It is important that the row and column PCBA traces are routed so as not to create near-field anomalies. To this end, the inbound and outbound traces to / from the coil are on different PCB planes—including in-bound and outbound traces to the Row / Column switches—lay exactly one atop the other to minimize loop area of the modulated field that would otherwise create innovation vector anomalies. Further, planarity and temperature stability in the PCB is paramount. A planar ceramic, aluminum-nitride (AlN) PCBA (AlN preferred because of its high thermal diffusivity and low coefficient of thermal expansion) ensure that geometric measurements are not confounded by a lack of planarity in the reference field generation.

[0254] 2.2. Self-calibrating Magnetometer Array

[0255] MMAIs are deployed in dynamic environments often wrapped around Wearer limbs, and often with rigid-flex or all-flex PCBAs to create “tiled structures.” (See Appendix A) While each MMAI is calibrated in a planar fashion per Section 2.1, in use these will be deformed inunpredictable ways. To this end, we embed one or more printed coils in the MMAI circuit board (either rigid, flex or rigid-flex) or in the MMAI assembly (comprising multiple rigid or rigid-flex PCBA subassembly), say, in PCBA mounting or alignment holes to save space. The IAP on the MMAI modulates calibrated current at high frequencies (125-200 Hz) to provide a magnetic field reference that can be demodulated by the MFMVP using an information filter. The mostly static fields generated by the magnets in the baseband are rejected by these filters. The demodulated fields are processed by information filters to compute sensor-to-coil the coordinate transformation,scT, that captures the sensor 9-DOF (DOF=Degree-of-Freedom referring to the number of matrix coefficients that need to be determined in a homogeneous transformation.) centroid and roll-pitch-and yaw in the coil coordinate system. The IAP computes an analytical Jacobian so that adjustments to the coordinate transformation can take advantage of the IAP- MFMVP architecture, including firmware and HDL code (for FPGA and ASIC) used in magnet tracking. Determining the transformation from a known field is just the “inverse” problem of determining magnet state (location and orientation) from observed magnetic field data. The only difference is that now the coil (dipole) position is known and that the precise position and orientation of the sensors is not. The integration times in the Markov models for the demodulated signals are set to large values (~ 1 sec) to ensure calibration precision. (Note: The deformations are assumed to not exhibit high frequency transients.) A quiescent state is enforced when donning the human augmentation device so that a high-precision, initial condition is achieved. Global optimization is used during initial fitting or even during routine donning to deliver high- precision magnet state estimates prior to information filter activation that signifies the start of tracking.

[0256] Inclusion of printed circuit coils with active, modulated field generation in the MMAI creates a valuable risk control measure for single MMAI or in heterogeneous networks comprising multiple MMAI to deliver a self-test function after donning or, via the modulation feature a real-time self-test to confirm integrity of the components and the network on a continuous basis—this to ensure compliance with international safety standards including IEC 60601 (Safety and Essential Performance of Medical Electrical Equipment), ISO 14971 (Reducing and Managing Risks in Medical Devices), and IEC 62304 (Software Lifecycle Processes for Medical Devices).

[0257] 2.3 Self-calibrating Networks of Magnetometer Arrays

[0258] The intra-MMA in situ calibration described in Section 2.1 is applied directly to MMAI networks using the state estimation architecture described in Section 1.1.1.1.4. Here, a plurality of unique, node-specific, modulated coil frequencies, including activation, in each MMAI, or in at least one, is coordinated by the Super-IAP. Configuration of the MFIVP demodulation filters is accomplished by the Super- IAP so that each node “knows” what frequencies must be employed in the demodulation and which node is the source. Now, the network of MFIVPs creates a separate set of information vectors aggregated by the Super-IAP to determine the coordinate frame transformations between the Super-IAP MMA frame and the other MMAI frames in the network using an analytic-Jacobian-based information filter. In this way magnet and disturbance state estimation are then referenced in the Super-IAP or any other biophysically-inspired reference frame that is in known relation to the Super-IAP frame. Here, global optimization can be used to configure the network prior to first use or at the time of donning at which time the initial condition—the a priori information--can then be passed to the information filter for continuous operation. In some applications the global optimization is computed using a steepest descent optimization in cases where the magnetic field structure is not amenable to use of the analytical Jacobian. FIGS.19A-19B capture the concept and the coordinate frame mathematics associated with the self-calibrating networks.

[0259] 3. Advanced Architecture and Packaging of Magnetometry-based Neuromuscular Control Applications

[0260] The advances described in Sections 2 and 3 address gaps in precision and accuracy in the prior art, but both build atop both the prior and current work captured in the filed patents and in Appendix A respectively. The third gap in performance relates to packaging characteristics of the neuromuscular solution—specifically in the areal and volumetric footprint and in the mass of the MMAI. Here, architecture and packaging are considered holistically. Architectural decisions align with how components are packaged and interconnected. Functional aggregation via architecture optimization have a profound impact on cost of packaging and weight. Indeed, the computational methods and framework outlined in the prior sections have been developed concurrently with the hardware architecture and packaging. This section describes how the methods and framework defined in Section 1 and 2 are packaged.

[0261] 3.1. Architecture for Integrating Rigid-Flex Packaging in Single and Multi- Magnetometer Networks

[0262] 3.1.1. High-density, Rigid-flex Packaging of Magnetometer Array and Signal Processing Electronics for High-precision Magnet State and Muscle-Dynamics State Estimation

[0263] The MMAI, including the IAP and MFIVP controller circuit functions, employs two rigid sections that are connected via embedded flexible circuitry that is sandwiched in the PCBA when it is fabricated. The flexible section serves as a hinge that enables the processor section to be positioned atop the magnetometer array to reduce the areal footprint. As shown, the magnetometer array employs sensors mounted on both sides of the rigid board enabling a 2x16x24 (768) sensors to be mounted in a smaller space than the 96 sensors in the prior MMAI packaging. The 8x increase in the number of sensors, per the rules of thumb established by information filter simulations and captured in Table 1, reduces the RMSE by a factor of sqrt(8)=2.8. Here, RMSE=Root-Mean-Squared Error, in this case represented as the square-root of the trace (sum) of the magnet location covariance matrix, Sk-1, divided by the number of magnets (to establish a mean). It represents what can be thought of as an effective, standard- deviation of the magnet location error. This increase in sensor density is accomplished through elimination of the high-frequency capacitor components—replacing these with a type of printed circuit fabrication called “buried capacitance” in which internal power and ground planes are separated by thin dielectrics. These layers create capacitance without taking up areal space on the PCBA. The footprint of the sensor array is reduced by folding-in flexibly attached wings that package the sensor multiplexor circuitry driven by the MFIVP. A thin metal shield is soldered into this rigid section to isolate it from EMI generated by the IAP-MFIVP section and to provide a means of attaching the multiplexor “wings” and second rigid section.

[0264] The second rigid section of circuitry packages a 9x15 mm FPGA that houses both the IAP and MFMVP functions, including the communication interface to the actuator—packaged as a personality module plug-in—and also custom circuitry to process internally generated EMG that is presented by the AMI—also packaged as a personality module plug-in. The FPGA employs a quad-core processor for IAP functions; a high-performance coprocessor for information vector accumulation and updating analytic Jacobian and innovation computation; and an array of application-specific MFIVP DSPs to perform the signal acquisition,normalization, and information vector computation. The FPGA employs power shedding capability in the form of shutting off internal power domains when the Wearer is not active.

[0265] 3.1.2. Integrated Magnet and Muscle-Dynamics State Estimation in a Heterogeneous Network of Rigid and Deformable Magnetometer Arrays

[0266] The integration of a network of MMAI rigid-flex assemblies and communications hub are placed onto or embedded in flexible and deformable substrates including “finger-braid” or soft-plastic sleaves. The interconnects between the various nodes and the hub are customized to achieve a “pigtail” that reduces the number of pluggable connectors to increase reliability. EMG interfaces are implemented as seamless rigid-flex assemblies that interface to AMI or internal or surface EMG. PoE is used to eliminate a separate power connection to each of the nodes. (PoE=Power-over-Ethernet which is a standard that embeds power connection in the Ethernet communications cable to reduce wiring complexity and mass.)

[0267] 3.2. Power Management

[0268] Power management is integral to closing the gap in the mass of the solution—where the mass of the solution is the sum of the mass of the: (1) MMAI, internal power supplies, and associated interconnects and (2) Battery energy storage (g / Joule * Joules needed to power the network) needed to power the MMA network.

[0269] To that end, the actuator processor monitors the PoE power consumption and Wearer activity and issues commands to the IAP to enter Standby and Sleep states as appropriate and to shut-off power to select magnetometer ICs and other power-hungry devices. Wearer activity, or lack thereof, is detected by multiple methods, including use of IMU, internal AMI EMG, or lack of significant change in magnet target or disturbance state over time.

[0270] Note: If a typical lower-extremity Wearer executes 1-2K gait cycles in a day, each lasting, say, two seconds, the “on-time” will be only 2+ hours a 4x reduction from that consumed in a 12-hour day without power saving measures applied. For upper extremity wearers different power saving measures can be applied, including direct-Wearer interventions.

[0271] 3.3. Advances in magnet implantation and Signal Integration

[0272] 3.3.1. Directed Magnet Implantation to Optimize Magnet and Muscle-Dynamic State Estimation

[0273] The Information Filter studies captured in Table 1 conclude that magnet orientation is critically important. Indeed, in a single-movable magnet application (for instance two magnets, one in bone and one in an attached tendon) the magnet location RMSE degrades by over 40% if the sub-optimal magnet orientations are applied. In multi-magnet deployments magnet orientation can have different and perhaps more signification impacts. For any specific wearer deployment in which the biophysically-inspired magnet implantation strategy is defined at the time of the implantation surgery, Information Filter simulations can be used prior to surgery to optimize magnet orientation for each in the population of deployed magnets. From an information filter standpoint, some variation is allowed. A + / -10 ^ variation is sufficient to gain at least 80% of the precision-optimization benefit. A method for ensuring magnet orientation optimization now follows.

[0274] Here, a magnetized, sterile pouch carries multiple pairs of magnets the pouch embedded with a planar magnet in its base that serves to align the sterilized magnet orientation when that magnet is placed by a human surgeon or robot. The package is then sealed for shipment and storage. At the time of magnet placement and deployment, the surgeon—guided by the optimized placement “floor plan”—uses a battery-operated instrument with a syringe tip that projects a programmable magnetic field orientation via a rigid-flex coil structure array activated by a microprocessor in the instrument. The orientation of the magnet is manipulated via magnetic field commutation using the coil structure array—in much the way a stepping motor is rotated by field commutation. Once acquired, the magnet can be implanted within a patient’s body. During extraction the syringe magnetic field is deactivated by the surgeon by a simple, poke-a-yoke HMI embedded in the instrument. The rigid-flex coil structure can be detached and disposed of at the end of the surgery. The floor planning described above can create magnet implant geometries that substantially eliminate forces between pairs of dipoles embedded in muscle tissue—specifically by orienting the dipoles to be orthogonal to each other and each orthogonal to the vector between the dipoles. This works to attenuate the dipole migration “forcing function” enabling higher implanted dipole density. By constraining the implanted dipole geometry in this way, the information content within the magnetic field created by the dipole pair is also enhanced thereby increasing tracking precision.

[0275] 3.3.2. EMG and Electrostimulation Packaging Integration

[0276] As discussed in Section 3.2 above, provisions are made in the customized sleeve to handle agonist-antagonist muscle-tendon extensions for connection to efferent internal or surface-mounted EMG and to afferent electrostimulation controllers.

[0277] 3.3.3. Post-Amputation Magnet Imaging for Computer-Aided-Engineering (CAE) of Deformable Magnetometer Array Packaging, Integration and Deployment

[0278] As described in Section 1 and 2, pre-Amputation procedures can include “floor planning” in magnet placement and orientation. That “floor planning” can be extended to placement of the MMAI via Information Filter simulations to design the MMAI attachment, via carbon fiber or other, to the bone structure exposed by the amputation surgery—this to reduce relative motion between the MMAI and the implanted array of magnets. This further defines the “initial conditions” for the magnet state vectors that in turn informs the global optimization if that is required to ensure convergence to the correct “minimum” or to inform the information filter with initial computation of analytic Jacobian to avoid using the global optimization. Also MFIVP and IAP firmware customization is determined—the customization comprising configuration of firmware or FPGA HDL for use with the information filter to incorporate customized Markov disturbance models that are application-specific. For instance, activation detection via lateral tendon vibration can be considered and configured in the firmware. As discussed in Section 1.2, supervised training and verification of the deep learning CNN will require the CAE component informed by post-AMI imaging.

[0279] 3.4 Preferred Magnet / Magnetometer Topology Embodiments Specific to Neuromuscular Prosthetic and Orthotic Applications

[0280] This section reviews three (3) specific magnet / magnetometer array topologies serving as preferred embodiments that promote dynamic state estimation each for specific prosthetic and orthotic neuromuscular control applications: (A) Magnets in tendon to measure muscle fascicle length, velocity and force; (B) Magnets implanted in bone to maximize precision of muscle- tendon dynamic estimates; and (C) Magnetometers structurally mounted onto osseointegrated implant abutments to achieve precision muscle dynamic state estimates without bone implanted magnets for prosthetic applications.

[0281] 3.4.1. Topology A: Magnets mounted in tendon and muscle belly to measure muscle fascicle length, velocity and force

[0282] 3.4.1.1. Upper Extremity Augmentation Application

[0283] FIG.21 shows muscle-tendon magnet deployment for estimation of muscle force, elongation and velocity of agonist muscle 204 and antagonist muscle 205 by magnetometer array 206 for thumb actuation of a limb 200. Here three magnets are implanted for each muscle-tendon unit 204 and 205. Two magnets are implanted in either an origin or an insertion tendon for the measurement of tendon elongation, which when combined with tendon stiffness, is used to compute muscle-tendon force. Still further, the relative separation distance between magnets positioned near the musculotendinous junction of each muscle-tendon unit 204 and 205 is used to estimate muscle belly length changes. More specifically, force in muscle 204 is estimated by extension between magnet 214 and magnet 202 multiplied by distal tendon 204 stiffness. Muscle 204 elongation is measured by relative displacement between magnets 202 and 210. Force in muscle 205 is estimated by extension between magnet 213 and 211 and distal tendon 205 stiffness. Muscle 205 elongation is measurement by the relative displacement of magnet 211 and 212.

[0284] 3.4.2. Topology B: Magnets implanted in bone to maximize precision of muscle- tendon dynamic estimates

[0285] 3.4.2.1. Upper Extremity Augmentation Applications

[0286] FIG.22 shows within arm 200, magnet 201 is implanted into or onto bone 208, and serves as a reference for measurement of muscle fascicle length changes using magnets 202 and 203 imbedded within muscles 204 and 205, respectively. Here an estimation of the relative distance between the bone-grounded magnet 201 and each muscle magnet 202 and 203 imbedded within muscles 204 and 205, respectively. Skin-mounted MMAI 206 and 207 of arm 200 are used to track the magnet muscle targets 202 and 203 each relative to bone-mounted reference magnet 201.

[0287] FIG.23 shows a cutaway of limb 100 shows reference bone magnets 123 and 121 implanted into bones 122 and 120, respectively. These bone-implanted magnets serve as a reference for muscle fascicle length change measurements. The relative distance between reference magnet 123 and muscle magnets 112, 111, and 110 are used to estimate muscle length changes for muscles 119, 118, and 117, respectively. Similarly, the relative distance between reference magnet 121 and muscle magnets 109, 107, 106, and 105 are used to estimate musclelength changes for muscles 116, 115, 114, and 113, respectively. These muscle magnets are imaged by magnetometer arrays 104, 103 and 102.

[0288] 3.4.2.2. Transradial Prosthetic Application

[0289] FIG.24 shows an amputated limb 220 that comprises reference bone magnet 201 implanted into bone 208. This bone attached magnet serves as a reference for muscle fascicle length change measurements. The relative distance between reference bone magnet 201 and muscle magnets 202 and 203 are used to estimate muscle length changes for muscles 204 and 205, respectively. These muscle magnets 202 and 203 are imaged by magnetometer arrays 222 and 221, respectively, that are attached to a prosthetic socket interface.

[0290] FIG.25 shows, for an amputated limb 140, a cutaway of the amputated residuum shows reference bone magnets 123 and 121 implanted into bones 122 and 120, respectively. These bone-attached magnets serve as a reference for muscle fascicle length change measurements. The relative distance between reference magnet 123 and muscle magnets 112, 111, and 110 are used to estimate muscle length changes for muscles 119, 118, and 117, respectively. Similarly, the relative distance between reference magnet 121 and muscle magnets 109, 115, 106, and 105 are used to estimate muscle length changes for muscles 116, 107, 114, and 113, respectively. These muscle magnets are imaged by magnetometer arrays 104, 103 and 102 attached to prosthetic socket interface 141.

[0291] 3.4.2.3. Transtibial Prosthetic Application

[0292] FIG.26 shows, for an amputated transtibial limb 320, reference bone magnet 306 is implanted into tibial bone 301. This bone-attached magnet 306 serves as a reference for muscle fascicle length change measurements. The relative distance between reference bone magnet 306 and muscle magnets 305 and 304 are used to estimate muscle length changes for muscles 303 and 302, respectively. These muscle magnets 305 and 304 are imaged by magnetometer arrays 322 and 321, respectively, attached to prosthetic socket interface 323.

[0293] FIG.27 shows, for transtibial amputated leg 350, a self-calibrating, tiled sensor array 352 integrated into socket 351 measures magnet displacements in muscle fascicles 302 and 303 relative to bone 301.

[0294] 3.4.2.4. Lower Extremity Augmentation Application

[0295] FIG.28 shows, for leg 300, reference bone magnet 306 is implanted into tibial bone 301. This bone-implanted magnet 306 serves as a reference for muscle fascicle length change measurements. The relative distance between reference bone magnet 306 and muscle magnets 305 and 304 are used to estimate muscle length changes for muscles 303 and 302, respectively. These muscle magnets 305 and 304 are imaged by magnetometer arrays 308 and 307, respectively.

[0296] 3.4.3. Topology C: Magnetometer arrays structurally mounted onto osseointegrated implant abutments to achieve precision muscle dynamic state estimates in prosthetic applications without bone-implanted magnets

[0297] 3.4.2.1. Transradial Application

[0298] FIG.29 shows magnetometer arrays 155, 156, 157, 158 and 159 are packaged onto rigid fiber-composite arms 160, 161, 162, 163 and 164 that are affixed to an osseointegrated implant 166, respectively. By attaching the magnetometer arrays 155, 156, 157, 158 and 159 to the osseointegrated implant 166, the arrays are effectively mechanically grounded to the residual bone 165. Given this bone reference, only one magnet is implanted into each residual muscle. Muscle magnets 151, 152, 153 and 154 are used to measure the muscle displacements in muscles 114, 116, 118, and 112, respectively.

[0299] FIG.30 shows, for transradial amputated arm 240, magnetometer arrays 241 and 242 are packaged onto rigid fiber-composite arms that are affixed to an osseointegrated implant 243. By attaching the magnetometer arrays 241 and 242 to the osseointegrated implant 243, the arrays are effectively mechanically grounded to the residual bone 208. Given this bone reference, only one magnet is implanted into each residual muscle. Muscle magnets 202 and 203 are used to measure the muscle displacements of muscles 204 and 205, respectively.

[0300] 3.4.3.2. Transtibial Application

[0301] FIG.31 shows, for transtibial amputated leg 340, magnetometer arrays 341 and 342 are attached onto rigid fiber-composite arms 343 and 344, respectively, that are affixed to an osseointegrated implant 345. By attaching the magnetometer arrays 341 and 342 to the osseointegrated implant 345, the arrays are effectively mechanically grounded to the residual tibial bone 301. Given this bone reference, only one magnet is implanted into each residualmuscle. Muscle magnets 304 and 100 are used to measure the muscle displacements of muscles 302 and 303, respectively.

[0302] References

[0303] 1. Information Filter and Kalman filter smoothing implementation references (See https: / / en.wikipedia.org / wiki / Kalman_filter)

[0304] 2. Applied Optimal Estimation, Technical Staff—The Analytic Sciences Corporation, Edited by Art Gelb, MIT Press, 1974.

[0305] 3. Trilinear Interpolation described( https: / / en.wikipedia.org / wiki / Trilinear_interpolation)

[0306] Appendix A. Further Magnetomicrometric Advances in Robotic Control

[0307] A.1 Background

[0308] We live in an era of wearable sensing, where our movement through the world is continuously monitored by devices. And yet, we lack a portable sensor that can continuously monitor muscle, tendon, and bone motion, allowing us to monitor performance, deliver targeted rehabilitation, and provide intuitive, reflexive control over prostheses and exoskeletons. Here, we introduce a sensing modality, magnetomicrometry, that uses the relative positions of implanted magnetic beads to enable wireless tracking of tissue length changes. We demonstrate real-time muscle length tracking in an in vivo turkey model via chronically implanted magnetic beads, while investigating accuracy, biocompatibility, and long-term implant stability. We anticipate that this tool will lay the groundwork for volitional control over wearable robots via real-time tracking of muscle lengths and speeds. Further, to inform future biomimetic control strategies, magnetomicrometry may also be used in the in vivo tracking of biological tissues to elucidate biomechanical principles of animal and human movement.

[0309] 1.1 Magnetomicrometry

[0310] Accurate, timely monitoring of user intent is necessary to provide volitional control over a prosthesis, exoskeleton, or other human-machine interface. As a result, substantial work has been undertaken towards developing approaches to measure intent by tracking the nervous, mechanical, and chemical signals generated by peripheral limbs (1–3). Amongst the mechanical parameters measured are muscle length and shortening speed, which must ideally be tracked on atimescale of tens of milliseconds with millimeter resolution to be useful for reflexive control of prostheses and exoskeletons (4, 5).

[0311] Non-invasive approaches to monitoring user intent, such as surface electromyography (EMG), ultrasound, and mechanomyography, reside outside the body but have poor, unstable signal quality (6, 7) or require substantial mass, power, and computation (5). For example, fluoromicrometry, which uses X-rays for high precision tissue position tracking, is wireless but is limited to short bursts due to ionizing radiation, requires an entire room, and involves substantial processing time (8). And whereas high-density surface EMG is portable and can be sufficiently accurate to decode spinal neural drives (9), signal drift and large artifacts due to skin- electrode impedance variations can be caused by changes in perspiration (10) or by dynamic pressure changes from, for instance, a prosthetic socket (11).

[0312] In contrast, highly-invasive approaches such as sonomicrometry, electrodes implanted in peripheral nerves, and EMG via implanted muscle electrodes provide improved signal quality but are expensive to implement, require delicate surgery, and are prone to damage or variable performance over time (6, 12). For instance, sonomicrometry uses implanted ultrasound crystals to yield high accuracy (13) but requires percutaneous wires and is difficult to miniaturize, precluding its use in humans. Additionally, EMG, whether invasive or not, only senses muscle activation, which without muscle length and velocity cannot be used to reliably observe, understand, or utilize muscle action (14). Despite the breadth of previous research, the field is missing a portable sensor that can perform accurate, minimally invasive, real-time measurement of muscle length to inform user peripheral intent.

[0313] This work introduces a low-footprint, minimally invasive device to measure the real- time length of tissues, including muscle tissues, that is accurate, easy to implement and provides high signal quality. It uses multiple implanted magnetic beads to wirelessly track tissue lengths via an array of magnetic field sensors, which senses the relative locations of the implanted magnetic beads. FIG.32 shows how this technique can be applied to tracking local muscle tissue lengths in the control of a prosthesis.

[0314] FIG.32 shows an example embodiment for Free-Space Control of a Robotic Prosthesis via Muscle Magnetomicrometry. Passive magnetic beads (indicated here by black circles) implanted in muscle can be used to wirelessly track muscle length via an array ofmagnetic field sensors (shown here as rectangular sensing arrays) mounted to the outside of the body. The pairs of magnetic beads shown here could each be placed in a single muscle in line with the muscle fiber orientation, or could be implanted in some other way. Muscle length data can be streamed to a control unit, which can in turn be used to stream commands to neuroprosthetic devices such as exoskeletons, muscle stimulators, or the robotic hand shown here. In a free-space control methodology, agonist and antagonist muscle states volitionally commanded by the user are mapped through a model of an intact biological limb to control joint angles by modulating motor torque. This control strategy can be extended beyond free-space control by incorporating muscle activation or direct force measurement.

[0315] Previously, magnets have been permanently implanted in humans alongside Hall sensors for joint tracking, successfully demonstrating the viability and safety of this approach (17). Because low-frequency magnetic fields are not affected by materials such as silicone, carbon fiber, or the human body, the magnetic field passes undisturbed from the muscles to the sensors as if these other materials are not present. This allows for accurate, transcutaneous, real- time tracking of the unpowered implants.

[0316] This real-time tracking of tissue length via magnetic beads is made possible by advances recently demonstrated in magnetic target tracking. Historically, magnet tracking methods have been slow, precluding real-time magnetic target tracking in high bandwidth applications. Further, traditional magnetic target tracking has suffered from inaccuracy due to ambient magnetic field disturbances, such as the geomagnetic field, restricting its use in a mobile context (15). Teachings herein demonstrate an improved method to track multiple magnets with high speed and accuracy while compensating for magnetic disturbances, enabling real-time, mobile use of magnetic target tracking in the control of human-machine interfaces (16).

[0317] A.1.2 Magnet Tracking

[0318] A commonly used method of tracking permanent magnets utilizes an optimization algorithm. At each step of the optimization, each of the magnet parameters (locations, orientations, and / or strengths) is estimated. Note that magnet strengths need not be estimated once they are known, but estimating the magnet strengths using the tracking algorithm can correct for errors in the factory specifications of magnets as well as their relationship to sensor gains. The estimate of magnet parameters is used to calculate a predicted magnetic field at eachknown sensor location in a sensor array. Comparing the magnetic field prediction with the actual magnetic field measurement at each sensor, a magnetic field prediction error is then computed corresponding to the current estimate of the magnet parameters. In the case of a gradient descent optimization, the derivative of the prediction error (i.e., the Jacobian matrix of the prediction error) is then determined with respect to each of the magnet parameter estimates, and these derivatives are used to update the magnet parameter estimates until the prediction error is minimized. The magnet parameters determined from the optimization solution are then used as the initial estimate to the subsequent tracking step.

[0319] The derivatives of the magnetic field prediction error are typically computed numerically. Computing these derivatives numerically is time intensive because it requires the prediction error to be computed at least once for every magnet parameter being tracked. The added computational time places limitations on real-time tracking bandwidth and, when the tracked magnets change position rapidly, can result in tracking instability.

[0320] A tracking algorithm is described below, implementing the use of analytic derivatives, to track spherical magnets via an optimization algorithm. The analytic derivatives in this tracking algorithm are implemented in a manner that has the benefits of numerical stability and allows for a significant decrease in latency compared with the latency inherent to other algorithms. Further, this tracking method is extended to the tracking of disturbance fields.

[0321] A.1.2.1 Cost Function

[0322] The magnetic field prediction error is used for an optimization cost function. At the ith sensor, the magnetic field prediction error, E_i, is the difference between the measured magnetic field B_tilde_i and the predicted magnetic field B_i,

[0323] To compute the predicted magnetic field B_i, we use our estimate of the magnet locations, orientations, and strengths.

[0324] If the estimated location of the jth magnet is (x_j, y_j, z_j) and the position of the ith sensor is (s_ix, s_iy, s_iz), it follows that a vector from the jth magnet to the ith sensor is given bywhere, for simplicity, we will define x_bar_ij, y_bar_ij, and z_bar_ij such that

[0325] Using the positive z-axis as an arbitrarily-chosen reference, the orientation of the jth magnet can be described bywhere theta_j and phi_j are the magnet's estimated orientation from vertical and around vertical, respectively, and m_j is the magnet's strength, or magnetic moment. By the definition of the magnetic moment, we havewhere B_rj is the jth magnet's residual flux density, R_j is its radius, and mu_0 is the permeability of free space. If we define the magnetic dipole weight of the jth magnet asthe strength of the magnet can be expressed more simply as

[0326] Using the equation for the magnetic field of a dipole, the magnetic field prediction B_i=(B_ix, B_iy, B_iz) at the ith sensor can then be expressed aswhere G=(G_x, G_y, G_z) is an estimate of the spatially uniform disturbance field and M is the number of magnets.

[0327] Substituting equations (6), (4), and (3) into (8), the three components of B_i are given by equations (8a), (8b), and (8c) below.

[0328] A.1.2.2 Analytic Jacobian

[0329] Having fully developed the cost function (equation 1) for an optimization, we now seek to form the Jacobian matrix, a matrix of the derivatives of the cost function elements with respect to each of the estimated magnet parameters (x_j, y_j, z_j, theta_j, phi_j and m_bar_j) for each of the magnets.

[0330] The measured magnetic field B_tilde_i does not vary with respect to our estimated magnet parameters (for example, the partial of E_i with respect to x_j is equivalent to the partial of B_i with respect to x_j), so the derivatives of the cost function can be written as a function of B_i alone. Thus, the Jacobian submatrix corresponding to the ith sensor and jth magnet can be calculated as

[0331] The first three columns of derivatives in equation (9) are written in terms of x_j, y_j, or z_j, but B_i in equations (8a), (8b), and (8c) is written in terms of x_bar_j, y_bar_j, andz_bar_j. From the definitions of x_bar_j, y_bar_j, and z_bar_j in equation (3) we see that x_bar_ij, y_bar_ij, and z_bar_ij are functions of x_j, y_j, and z_j, respectively. The chain rule can thus be used asfor x_j, and similarly for y_j and z_j, to get

[0332] These derivatives exist, and the analytic expressions for the elements of equation (11) are given by equations (11a)-(11r) (see FIG.33).

[0333] These derivatives were calculated by hand and verified using a symbolic equation solver. Because of the many repeated terms in equations (11a)-(11r), these elements are able to be efficiently calculated using common subexpression elimination.

[0334] The full 3N by 6M Jacobian matrix is composed of all of the Jacobian submatrices J_ij across the N sensors and M magnets, and is constructed as

[0335] A.1.2.3 Disturbance Field Compensation

[0336] Assuming that there is a uniform disturbance field seen across all magnetic field sensors (from, for example, the geomagnetic field), its Jacobian submatrix would be given bythe 3 by 3 identity matrix. The columns of the Jacobian submatrix corresponding to the disturbance field are thus given bya 3N by 3 matrix. Extending the Jacobian matrix equation (12) with the result in equation (14), we get the augmented Jacobian matrix

[0337] This augmented matrix assists in providing robust tracking information under sensor motion or far-field disturbance.

[0338] A.1.2.4 Calibration

[0339] The sensor array can be calibrated using rotation in a uniform ambient field to remove offsets and distortions, and then the sensors can be scaled relative to one another to achieve equivalent full-scale ranges.

[0340] A.1.2.5 Implementation

[0341] This magnet tracking algorithm can be written in a computer language such as C++, and can be implemented using an optimization algorithm such as an unconstrained Levenberg-Marquardt algorithm via C / C++ Minpack. The algorithm can be run in real-time on a computer, such as a single-board computer, a microcontroller, or an FPGA.

[0342] The three-axis magnetic field can be measured using magnetic field sensors (magnetometers), such as the LSM9DS1 iNEMO inertial modules (STMicroelectronics) or the LIS3MDL magnetometer (STMicroelectronics). The measurements can be communicated via SPI to a computer such as a microcontroller or FPGA and then relayed to a computer for the tracking, or the measurements can be directly relayed to the tracking computer. The tracking computer can be physically present or can be located remotely.

[0343] Magnets of various magnetizations, geometries, materials, and coatings could be used for tissue tracking, but as an example, uniformly-magnetized 3-mm-diameter spherical N48SH magnets coated in Nickel-Copper-Nickel, gold, and Parylene C could be used for tissue tracking.

[0344] A.1.4 Advantages of Magnetomicrometry

[0345] Single implanted magnets can be used to simultaneously monitor multiple muscles via external magnetic field sensors (18, 19). However, the single-magnet-per-muscle approach is limited in various ways. Muscle length can be passively cycled by the motion of a joint, such as when the elbow joint is engaged by a strong handshake from another person, or the muscle can be actively cycled when flexed, such as when holding a glass of water. In a controlled setting, a measurement of axial motion from a single point in the muscle could allow measurement of either the passive or active muscle length change (e.g., for free-space control or force control of a prosthesis), but these two sources of motion would confound one another when both are present. Further, single-magnet axial or radial displacement caused by muscle flexion (i.e., shortening and bulging of the muscle, which are roughly predictive of one another under the assumption of isovolume) would be challenging to measure due to movement of surrounding tissues or pressure from a prosthetic socket. These issues are solved by the use of multiple magnetic beads, such as a pair of magnetic beads in each muscle, allowing tissue length changes to be accurately measured regardless of tendon strain.

[0346] Using an approach we call magnetomicrometry, a pair of magnetic beads is implanted at a tissue, such as along the axis of each muscle or along the length of the muscle fascicle. Using externally mounted magnetic field sensor arrays, each magnetic bead pair is tracked wirelessly as outlined in previous work (16). The Euclidean distance between the three-dimensional positions of the beads is used to determine the length of the muscle, so the sensing of muscle length should remain unaffected by movement of the sensors or muscle relative to one another. The magnetic field sensors used for this tracking can be mounted to the skin, affixed to a prosthetic socket, or embedded in clothing, making this approach ideal for use in both stationary and mobile contexts.

[0347] As shown in FIG.32, one control strategy using magnetomicrometry maps muscle lengths to bionic joint angles through an intact biophysical limb model, providing the user intuitive volitional control over a robotic prosthesis or exoskeletal device. This strategy can be further extended beyond free-space control by incorporating muscle activation or direct musculotendon force measurement. For instance, muscle lengths and speeds from magnetomicrometry could be combined with EMG to calculate the force through a muscle model.

[0348] This system enables alternative device implementations for a variety of biomechanical applications.

[0349] The above background description is provided with more detail in the following references: h

[0350]

[0351] A.4 Tiled Sensor Arrays

[0352] Due to the nature of the field from magnetic dipoles, magnet tracking is inherently limited in how far the magnets can be from a magnetic field sensing away and still be tracked. This is because magnetic field sensing has limitations in both accuracy and precision due to sensor inaccuracy (such as inaccuracy due to nonlinearity) and sensor noise (such as Johnson noise).

[0353] In addition, when multiple magnet tracking systems (where a magnet tracking system is composed of at least one magnet and at least one subset of a magnetic field sensing array) are used concurrently in proximity to one another without information about one another, the presence of the neighboring magnet tracking systems can cause cross interference.

[0354] A.4.1 Sensor Array Tiles

[0355] These issues can be solved by using a smart array of sensor array tiles, where each sensor array tile is composed of one or more magnetic field sensors, and may additionally include one or more accelerometers, one or more gyros, one or more temperature sensors, as well as additional sensing components. Subsets of the multiple sensing components may share electronics components such as capacitors and may be temporally multiplexed for communications.

[0356] The use of additional magnetic field sensors in a magnet tracking system increases the accuracy of magnet tracking, in effect increasing the signal-to-noise ratio of tracking by providing additional information to the magnet tracking system.

[0357] The sensor array tiles may be one-, two-, or three-dimensional, may be multi-sided, and may be of any geometry. Embodiments including a two-dimensional hexagonal sensor array tile may maximize packing density.

[0358] Each tile can be powered via an integrated (on-board) power source, or can be powered via a power source shared between multiple boards, or can be powered via a power source shared with another sensing or robotic platform, or may be powered wirelessly or by harvested power, such as by motion of the user.

[0359] Each tile may be communicated with via an individual (parallel) wired connection, or wired for communication in series (daisy-chained), or communication may be performed wirelessly, such as via Bluetooth or WiFi communication.

[0360] The sensor array tiles may be affixed, either irreversibly (such as via epoxy) or reversibly (such as via Velcro, double-sided tape, 3M picture hanging strips, clips, or screws) to a rigid surface, affixed to one another via rigid connections or via flexible joints allowing flexibility while restricting any number of degrees of freedom, or may be stacked or interleaved. The sensor array tiles may have attachments on them allowing them to interlock with one another, and these attachments may also (in addition or alternatively) serve the purpose of creating electrical connections for power or communication. Some embodiments involve sensor array tiles rigidly affixed to a surface, such as a prosthetic socket as seen in the FIG.27.

[0361] A.4.2 Sensor Selection Strategy

[0362] In this context, “smart” refers to the use of the array of sensor array tiles to perform a trade-off attempting to optimally minimize the magnet tracking time delay of all magnet trackingsystems, the cross-interference between the magnet tracking systems, and the magnet tracking error of all magnet tracking systems, but choosing, previous or during real-time tracking, which subsets of the magnetic field sensors across all sensor array tiles should compose each magnet tracking system.

[0363] A.4.3 Calibration Strategy

[0364] To determine the orientations of the sensors relative to one another, the sensors may be first calibrated via the standard ellipsoid-to-sphere calibration followed by scaling, as described in briefly in 1.2.4 above (and in greater detail in section 6.1 of https: / / dspace.mit.edu / handle / 1721.1 / 130210), and then the relative orientations of the sensors may be either determined individually by rotation of the resultant spheres onto one another via optimization (as described in section 6.2 of https: / / dspace.mit.edu / handle / 1721.1 / 130210) or, as an improvement to this strategy, using singular value decomposition, as described by http: / / nghiaho.com / ?page_id=671 and http: / / nghiaho.com / ?page_id=846. The entire teachings at the links listed in this paragraph are incorporated by reference in their entirety.

[0365] Alternatively or in addition, the prior knowledge about the magnetic field sensor orientations on each tile (e.g., that they are all oriented in the same direction as one another within a sensor tile – as known from the manufacturing circuit board assembly process, or using information from a within-tile orientation calibration), may be used to orient the sensor arrays relative to one another. Some embodiments use singular value decomposition to rotate the resultant spherical data from all magnetic field sensors collectively from each given array onto the resultant spherical data from all magnetic field sensors collectively from a reference array.

[0366] Upon determining the relative orientations of the magnetic field sensors of the tiles relative to one another using one of the above strategies, a magnet (or multiple magnets) can be introduced in proximity to the magnetic field sensors, and the optimization of the pose (the multiple tracked positions and orientations over time) of the magnet over its path can be used to determine the positions of the tiled arrays relative to one another, as described in section 6.3 of https: / / dspace.mit.edu / handle / 1721.1 / 130210, but as an improvement using the known positions of the sensors, as known from the manufacturing circuit board assembly process (or from a within-tile position calibration). These one or more introduced magnets may be a magnet solely used for calibration (e.g., a large magnet, such as an 8-mm-diameter or 16-mm-diameterspherical magnet, or an electromagnet), or may be the magnets that are implanted for tissue tracking.

[0367] Alternatively, the relative positions and orientations of the tiles may be determined concurrently using optimization of the path of the introduced magnet seen by the sensor array tiles.

[0368] As a further alternative or for additional precision in this calibration, other sensors such as accelerometers and gyros may be used to perform this calibration using a similar calibration strategy.

[0369] As a further alternative, other strategies may be used to determine the positions of the sensor array tiles. For instance, an augmented reality headset equipped with position sensing, such as LIDAR, could be used to determine the relative positions and orientations of the arrays.

[0370] A.4.4 Other Uses

[0371] As is evident, these tiled sensing arrays could be used for other applications which involve the tracking of magnets or the sensing of magnetic fields for other purposes.

[0372] A.5 Error characterization

[0373] Error distributions from the magnetic field sensors propagate to resultant estimated poses of magnets. Due to the nature of modern magnet(s) pose tracking algorithms, quantitative evaluation of magnet pose error distribution is not trivial and has not been fully addressed. A novel algorithm is presented that enables estimation of this measurement error distribution. FIG. 34 shows an example generic pipeline of estimating error distribution of the estimated pose.

[0374] In the FIG.34, black rectangles represent typical pipeline for magnets poses tracking, and green rectangles represent the proposed algorithm. By re-evaluating the magnetic field errors from the magnetic field structure model of magnet, we can estimate variances of individual fields. Then, we can transform these to covariances of estimated magnets poses using the mathematical algorithm.

[0375] To illustrate, here is the generic magnets poses tracking algorithm using inverse optimization

[0376] The algorithm summarizes how an array of spatially distributed magnetic fields can be used to estimate magnets poses.

[0377] Based on the magnet structure models discussed above, the following figure is the rough summary of one example algorithm to transform magnetic field variances to magnets poses errors, while other algorithm can also be used under the same concept.

[0378] First of all, an effective whole sensor noise covariance matrix can be constructed by following equations.

[0379] Then, based on the ΣB , we can construct covariance matrix ΣP of the estimated state of magnet based on magnetomicrometry system using the magnetic field sensor-tracking result Jacobian (J) and pseudo inverse operation.

[0380] Applying singular vector decomposition or Eigen vector decomposition operations on ΣP, the error ellipsoid of the estimated magnets poses can be constructed.

[0381] Using this technique, we can evaluate how accurate the measurements of the magnets positions and orientations are in real-time.

[0382] A.5.1 Kalman Filtering

[0383] Based on the error characterization techniques, the measurement accuracy can be improved significantly using the technique called Extended Kalman filtering.

[0384] Therefore, as we now have a technique to estimate variance, we can apply the Extended Kalman Filter to increase accuracy. This technique enables iterative data updates based on the estimated model variances and past and current estimated data based on a mathematical tool to estimate magnet tracking accuracy.

[0385] Appendix B includes a document titled, “Enhancing Muscle Sensing Modalities for Advanced Bionics” by Seong Ho Yeon. The teachings in “Enhancing Muscle Sensing Modalities for Advanced Bionics” are incorporated by reference in their entirety. Any links and / or references in “Enhancing Muscle Sensing Modalities for Advanced Bionics” and / or Appendix B are incorporated by reference in their entirety.

[0386] The teachings of all patents, published applications, references, and web links cited herein are incorporated by reference in their entirety.

[0387] While example embodiments have been particularly shown and described, it will be understood by those skilled in the art that various changes in form and details may be madetherein without departing from the scope of the embodiments encompassed by the appended claims.APPENDIX B Enhancing Muscle Sensing Modalities for Advanced Bionics by Seong Ho Yeon B.S. Electrical and Engineering, Georgia Institute of Technolgy B.S. Electrical Engineering, Korea Advanced Institute of Science and Technology M.S. Program in Media Arts and Sciences, MIT Submitted to the Program in Media Arts and Sciences, School of Architecture and Planning in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY at the MASSACHUSETTS INSTITUTE OF TECHNOLOGY May 2024 © 2024 Seong Ho Yeon. All rights reserved. The author hereby grants to MIT a nonexclusive, worldwide, irrevocable, royalty-free license to exercise any and all rights under copyright, including to reproduce, preserve, distribute and publicly display copies of the thesis, or release the thesis under an open-access license. Authored by: Seong Ho Yeon Program in Media Arts and Sciences May 17, 2024 Certified by: Hugh M. Herr Professor, Program in Media Arts and Sciences, Thesis Supervisor Accepted by: Joseph Paradiso Academic Head Program in Media Arts and SciencesEnhancing Muscle Sensing Modalities for Advanced Bionics by Seong Ho Yeon Submitted to the Program in Media Arts and Sciences, School of Architecture and Planning on May 17, 2024 in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY ABSTRACT Muscle sensing technologies have significantly advanced our understanding of biomechan- ics and enhanced the efficacy of bionic devices. These technologies enable volitional control of prostheses and assistive devices by mapping the electrical and mechanical activities of mus- cles as control inputs. This dissertation presents novel paradigms and findings to improve the utility and efficacy of muscle sensing modalities for advanced bionic applications. In the first part, I introduce a comprehensive approach to improve the acquisition and processing of surface electromyography (sEMG) signals for bionic applications. This includes innovations in electrode materials and design to enhance user comfort and signal quality for long-term use within prosthetic sockets. Additionally, I propose a real-time impulse filtering algorithm to effectively suppress artifacts while preserving the underlying sEMG signal during dynamic movements. Furthermore, I demonstrate a synchronous sEMG and ultrasound acquisition method that enables simultaneous assessment of muscle electrical activity and mechanical deformation, providing valuable insights into muscle function and control. In the second part, I explore how Magnetomicrometry can serve as a new in-vivo and real-time mechanical muscle state tracking modality. Previous work has shown significant potential for Magnetomicrometry in muscle-state tracking via a tightly-controlled in situ setup. In this work, I demonstrate real-time tracking of muscle tissue length in freely-moving animals performing various motor activities, suggesting that Magnetomicrometry could be extended as a viable in-vivo and real-time muscle sensing modality. In the final part, I propose a novel theoretical framework leveraging Riemannian geom- etry and manifold theory to enhance the magnet tracking technology stack for Magnetomi- crometry. By representing the magnetic dipole state on a manifold and incorporating its dynamics, I develop a more accurate and robust magnet tracking algorithm that addresses the limitations of existing methods. Through simulations and real-world data evaluations, I demonstrate the superior performance of the proposed manifold-based tracking paradigm, showcasing its potential to improve the resolution and extend the observable depth of Mag- netomicrometry. The advancements presented in this dissertation have significant implications for the development of next-generation bionic devices, enabling more adaptive, versatile, and reliable myo-neural interfaces. Through this work, I hope to open up new possibilities for the designand control of advanced prostheses and assistive technologies with the advanced myo-neural control interface. Thesis supervisor: Hugh M. Herr Title: Professor, Program in Media Arts and SciencesEnhancing Muscle Sensing Modalities for Advanced Bionics by Seong Ho Yeon This dissertation has been reviewed and approved by the following committee members. Advisor Hugh M. Herr Professor of Media Arts and Sciences Co-Director of K. Lisa Yang Center for Bionics Associate Investigator, McGovern Institute MIT Media Lab Reader Joseph Paradiso Associate Academic Head Professor of Media Arts and Sciences MIT Media Lab Reader Thomas Roberts Vice Chair of the Department of Ecology and Evolutionary Biology Brown UniversityAcknowledgments First and foremost, I would like to express my deepest gratitude to my PhD advisor, Pro- fessor Hugh Herr, for his unwavering guidance, mentorship, trust, and support throughout my academic journey. Over the past seven years, during my Master’s and PhD studies, I have felt truly privileged to work alongside Hugh and learn from his invaluable insights and visionary perspective. His dedication to pushing the boundaries of bionics and his commit- ment to fostering my growth as a researcher have been instrumental in shaping my path and aspirations. I would also like to extend my sincere thanks to Professor Thomas Roberts for his men- torship and invaluable guidance during our collaboration on the Magnetomicrometry work. His expertise and dedications have been crucial in navigating the challenges and opportuni- ties of this innovative research, and I am grateful for the opportunity to have worked with him. Furthermore, I would like to express my appreciation to Professor Joseph Paradiso for serving on my PhD committee and providing valuable insights and feedback throughout my doctoral journey. His unique perspective and thought-provoking questions have consistently challenged me to think more critically about my research and its broader implications. I would like to thank to members of Biomechatronics group for being such a great team. This journey was only possible because we worked and collaborated as a team to tackle problems and challenges. A special mention goes to Dr. Cameron Taylor for his role in developing Magnetomicrometry technology; Rick Casler for his mentorship in information theory and control systems; Ellen Clarrissimeaux for her collaboration in Magnetomicrome- try validation efforts; Dr. Hyungeun Song, Tsung-Han Hsieh, Tony Shu, and Junqing Qiao for their collective work on our sEMG stack and mechatronics; and Guillermo Herrera-Arcos and Christian Landis for their engaging collaboration on animal experiments. Emily Rogers also deserves thanks for her contributions to our Mechatronics projects. I would also like to extend my thanks to Billy Clark and Mary-Kate O’Donnell from the Roberts Lab for their collaboration on the Magnetomicrometry validation work. I am thankful to all my wonderful friends in the US and Korea for their invaluable life and PhD advice, support, perspectives, and the joyful times we’ve shared over the past years and beyond. Last but not least, I want to express my sincere gratitude to my family. Thank you for believing in me unconditionally, supporting my decisions wholeheartedly, and providing boundless love. This journey was only possible because of your unwavering support.Contents Title page 1 Abstract 3 Acknowledgments 7 List of Figures 13 List of Tables 25 1 Introduction 27 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 1.2 Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.3 Specific Aims . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.4 Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2 Enhancing sEMG Utility for Bionics 31 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.2 Interface for Within Socket sEMG Acquisition: Flexible and Low-Profile Elec- trodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.2.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.2.2 Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.2.3 Quantitative Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.2.4 Experiment 1: sEMG Acquisition within Prosthetic Socket Systems . 44 2.2.5 Experiment 2: Comfort Evaluation . . . . . . . . . . . . . . . . . . . 52 2.2.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 2.3 sEMG Impulse Artifact Filter . . . . . . . . . . . . . . . . . . . . . . . . . . 57 2.3.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 2.3.2 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.3.3 Filtering Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 2.3.4 Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 2.3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 2.4 Synchronized sEMG and Ultrasonography Measurement . . . . . . . . . . . 75 2.4.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 2.4.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 2.4.3 Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 782.4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 2.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 3 In-vivo and Untethered Validation of Magnetomicrometry 85 3.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 3.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 3.2.1 Surgical Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 3.2.2 Magnetomicrometry . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 3.2.3 Accuracy Validation of Magnetomicrometry Against Fluoromicrometry 89 3.2.4 Untethered Muscle Tracking Across Various Activities . . . . . . . . . 90 3.2.5 Benchtop Magnetomicrometry Validation for the Variable Terrain Ac- tivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 3.2.6 Data Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 3.3 Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 3.3.1 Accuracy Validation of Magnetomicrometry Against Fluoromicrometry 93 3.3.2 Untethered Muscle Tracking Across Various Activities . . . . . . . . . 95 3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 3.4.1 Accuracy Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 3.4.2 Ambient Magnetic Field . . . . . . . . . . . . . . . . . . . . . . . . . 102 3.4.3 Range of Behaviors . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 3.4.4 Potential Applications . . . . . . . . . . . . . . . . . . . . . . . . . . 103 3.4.5 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 3.4.6 Sensing Improvements . . . . . . . . . . . . . . . . . . . . . . . . . . 105 4 Advancing Magnet Tracking Paradigm with Magnetic Dipole Manifold 107 4.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 4.2 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 4.2.1 Magnetic Dipole Model . . . . . . . . . . . . . . . . . . . . . . . . . . 109 4.2.2 Existing Magnet Localization and Tracking Methods . . . . . . . . . 112 4.2.3 Kalman Filter on Parallelizable Riemannian Manifolds . . . . . . . . 114 4.3 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 4.3.1 Magnetic Dipole True Manifold Representation . . . . . . . . . . . . 116 4.3.2 Dynamic State Lifting of Magnetic Dipole Manifold . . . . . . . . . . 117 4.3.3 Magnetic Dipole Manifold with Kalman Filter . . . . . . . . . . . . . 119 4.4 Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 4.4.1 Evaluation Metric . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 4.4.2 Simulation Setups . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 4.4.3 Static Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 4.4.4 Dynamic Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 4.4.5 Evaluation on Real-World Data . . . . . . . . . . . . . . . . . . . . . 148 4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1505 Conclusion 153 5.1 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 5.2 Future Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 5.3 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 References 157List of Figures 1.1 Schematics of myo-neural interface using sEMG and Magnetomicrometry for neural control of advanced Bionics devices. . . . . . . . . . . . . . . . . . . . 28 2.1 SLIP electrode design concept showing electrode worn against the skin of the residuum inside a standard silicone prosthetic liner. An integrated, flexible cable and connector allows the sensor to interface with an acquisition and processing system mounted externally on the socket. . . . . . . . . . . . . . . 36 2.2 Electrical model of the novel electrode. Each bipolar electrode has a diameter of 10 mm and inter-electrode distance of 16mm. Mechanical holes on the electrode were placed for potential anchoring and mounting inside a prosthetic liner or sock.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.3 Fabricated prototype of the novel electrode. (a) Size comparison of the elec- trode. (b,c) Demonstration of manufactured electrode flexibility. . . . . . . . 39 2.4 Example electrode configuration across the TA with Ag / AgCl electrodes in the medial position. An adhesive piece of tape is applied on top of the flexible electrode to ensure proper placement. Elastic exercise band applying uniform pressure to all electrodes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.5 Embedded 8-channel sEMG acquisition and processing platform used for eval- uation with example Ag / AgCl bipolar pair, SLIP electrode, and ground elec- trode attached. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.6 Bland-Altman plots of (a) relative muscle activity on-time and (b) relative muscle activity off-time. Positive differences indicate earlier detection by the SLIP electrode compared to commercial Ag / AgCl electrodes. Blue circles indicate differences from dorsiflexion trials (n = 22) while orange diamonds indicate differences from cocontraction trials (n = 21).. . . . . . . . . . . . . 45 2.7 (a-e) Representative TA muscle activations (α) processed from data collected during dorsiflexion trials for each of subjects 1-5. Vertical dashed red lines indicate on- and off-times calculated with data from the Ag / AgCl electrodes. Vertical solid green lines indicate on- and off-times calculated with data from the SLIP electrode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 452.8 Example use case of the SLIP electrodes within an LE prosthetic liner. The 16 SLIP electrodes are placed and attached on the residual limb. This figure does not represent the actual instrumentation setup of the Experiment B. (a,b) Preparation of the residual limb’s skin surface. Corresponding target musculature are labeled and bony anatomical landmarks are marked upon the skin’s surface. (c,d) Electrodes are placed on the surface of the residual limb using Kinesio tape. (e) Leads of the SLIP electrodes are routed through the prosthetic sock with modified holes. (f) Fully routed SLIP electrodes under the prosthetic liner. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 2.9 The developed 16-channel custom sEMG system. . . . . . . . . . . . . . . . 48 2.10 Example use case of the embedded 16-channel sEMG acquisition and pro- cessing platform. The sEMG system is connected to the 16 electrodes as configured in Fig. 2.8f. This setup does not represent the actual instrumen- tation setup of the Experiment B, which only instrumented the four muscles of the AMI muscle pairs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 2.11 Portable sEMG measurement setup used for Experiment B during both sta- tionary voluntary movement and ambulation. Four SLIP electrodes interface with the four muscles of the AMI muscle pairs (TA-LG and PL-TP). . . . . 50 2.12 Measured and normalized sEMG signals from the four SLIP electrodes over target AMI muscle pairs (TA-LG and PL-TP) during voluntary clockwise rotations of the phantom ankle and subtalar joints. . . . . . . . . . . . . . . 50 2.13 Signals measured using the 16-channel sEMG processing platform with four electrodes over target AMI muscle pairs during treadmill walking at a speed of 1.4 m / s. The measured raw sEMG signal, processed sEMG signal, and ground reaction force (GRF) data are shown. Cyclic and repetitive muscle activation patterns correlating with GRF are demonstrated. Stance phases are presented with a white background while swing phases are presented with a gray background. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.14 Skin of the residual limb after one hour of weight-bearing ambulation. Skin indentation and irritation are not visible where electrodes were placed. . . . 52 2.15 Electrode configurations for user comfort evaluation: (a) SLIP electrodes; (b) Commercial Ag / AgCl electrodes without connectors and wires; (c) Commer- cial Ag / AgCl electrodes with connectors and wires. . . . . . . . . . . . . . . 53 2.16 Normalized experimental probability density estimates of triceps muscle con- traction up to 50% maximum voluntary contraction, conducted as experiment II in

[0051] . Recreated from

[0051] with permission©2021 IEEE. . . . . . . . . . 60 2.17 Generic sEMG signal processing pipeline . . . . . . . . . . . . . . . . . . . . 60 2.18 Modeled impulse signal, its response after band-pass filtering, and correspond- ing features extracted. The impulse signal is modeled with unit amplitude and 10 ms duration. Peak of impulse onset occurs at 0 s. The response shown is derived in a non-causal manner with zero phase and group delay. . . . . . . . 61 2.19 Modified sEMG signal processing pipeline with the proposed CHF filtering method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 622.20 Graphic demonstration of Algorithm 1 and Algorithm 2 with filter coefficients hLB of 10% and hUB of 90%. (a) Example sEMG time series highlighting the short-time windowed signal for feature extraction (b) Histogram plot and cumulative histogram plot extracted from Fig. 2.20a. The filter coefficients, hLB and hUB, are highlighted with corresponding xLB of 4.4% and xUB of 56%. The distribution of relative amplitudes is observed to be skewed right.(c) Sorted and unsorted sEMG data based on Fig. 2.20a highlighting original data and filtered data with hLB, hUB, xLB, xUB indicated. . . . . . . . . . . . . 64 2.21 Reference sEMG data collected from TA and TP muscles during voluntary rotational movement of the residual phantom limb. The plot visualizes both RMS and MAV processing outputs. . . . . . . . . . . . . . . . . . . . . . . . 67 2.22 Partial results from all 289 filtered sEMG series using swept boundary coeffi- cients hUB and hLB, processed from the original reference sEMG of Fig. 2.21. Partial data from 9 to 16 seconds shown. . . . . . . . . . . . . . . . . . . . . 68 2.23 Effect of filtering on SNR of reference sEMG data from Fig. 2.21 with swept filter coefficients. (a) Relative SNR of the filtered TA sEMG data. Filtering preserves at least 95% relative SNR for both MAV and RMS features. (b) Relative SNR of the filtered TP sEMG data. As opposed to the TA muscle, filtering of the TP muscle data results in increased SNR. . . . . . . . . . . . 69 2.24 Synthesized artifact-affected sEMG data. (a) Artifact-affected and band- passed TA muscle sEMG. (b) Effect of impulse artifacts on output RMS and MAV features. Partial data from 9 to 16 seconds shown . . . . . . . . . . . . 70 2.25 Artifact suppression ratio analyses with swept hLB and hUB filter coeffi- cients. Artifact suppression ratio is largely correlated with hUB. (a) Artifact suppression ratio of the TA (b) Artifact suppression ratio of the TP . . . . . 71 2.26 Artifact suppression ratio of CHF. Synthesized sEMG data shown in Fig. 2.24 are processed with hLB of 10% and hUB of 60%, 70%, 80%, and 90%. RMS features show smaller impact artifacts with decreasing hUB. . . . . . . . . . 71 2.27 sEMG data collected during ambulation were processed using CHF with hUB at 70% and hLB at 10%. Impulse artifacts are synchronized with ground reaction force peaks, but suppressed in the resulting RMS feature. . . . . . . 72 2.28 Experimental setup for the validation of the proposed method through the phantom. (2.28a) The phantom used for the validation. (2.28b, 2.28c) Mea- surements of the phantom with and without the proposed electrode. (2.28d) The spatially synchronized electrode with the ultrasound probe. . . . . . . . 77 2.29 Experimental setup for validation of the proposed method through assess- ments of TA muscle dynamic movements. (2.29a) Sensor placements. The goniometer and ultrasound probe were placed on the medial aspect of the ankle joint and muscle belly of the TA, respectively. (2.29b) The flexible electrode was placed under the ultrasound probe, enabling spatiotemporally synchronized sEMG and ultrasonography recording. . . . . . . . . . . . . . . 78 2.30 Qualitative evaluation in ultrasonography of a phantom. (2.30a) Result with the flexible electrode. (2.30b) Result with the flexible electrode. . . . . . . . 792.31 Artifact in ultrasonography of the TA muscle due to the flexible electrode placed under the ultrasound probe. Only negligible artifact is introduced in the ultrasonography as a form of shade by the presence of the flexible electrode. (2.31a) Ultrasonogrpahy of the TA muscle without the flexible electrode. (2.31b) Ultrasonography of the TA muscle with the flexible electrode. 80 2.32 Experimental results of spatiotemporally synchronized electromyography and ultrasonography recording of the TA muscle. (2.32a) and (2.32b) are rep- resentative recordings of ankle joint angle and normalized fascicle length of the TA msucle during full dorsiflexion and plantarflexion with and without EMG recording from the muscle, respectively. (2.32c) Joint angle and normal- ized fascicle length relationships investigated without (w / o EMG) and with (w EMG) simultaneous EMG recordings during full dorsiflexion and plan- tarflexion (n=6 cycles). No significant differences were found between the two relationship curves. P values are reported. Wilcoxon signed-rank test at a significance level of α =0.05. . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 3.1 Magnetomicrometry Embedded System. We fabricated a custom sensor board (left) and a custom control board (right) for this study. The sensor board holds the magnetomicrometry sensing array, consisting of 96 magnetic field sensors arranged with a spacing of 5.08 mm. Digital multiplexers on the sensor board allow time-domain multiplexing, enabling a single microcontroller on the control board to communicate with and control all magnetic field sensors on the sensor board. The control board merges the data from the sensor board and streams the data wirelessly to the magnet tracking computer. The sensor board and control board weigh 24 g and 12 g, respectively. . . . . . . . . . . 89 3.2 Benchtop Magnetomicrometry Validation. (A) Two magnets were placed 40 mm apart in a 1x6 LEGO Technic block and centered under the sensing array. We used the tracked magnet z-position as a guide in setting up the minimum depth measurements, and we enforced the remaining depths by adding / re- moving 3.2-mm-thick 1x6 plates under the Technic block. (B) We manually swept the magnets out and back to center along the x and z axes. We set the sweep trajectory to sweep just beyond the volume within which the magnets were tracked during the variable terrain activities (see Supplementary Figure 6). For reference, the centroid of the two beads is labeled at the origin and at the extent of the sweep trajectories. (C) The vertical axis represents the magnetomicrometry error, and the horizontal axes represent the centroid x and y position for the two sweep trajectories at each depth. The submillime- ter error range is marked with a gray background. A maximum error of 1.463 mm was found at a test location just beyond the tracked bead position extent (bottom right plot). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 913.3 Validation of Untethered Muscle Tracking using Magnetomicrometry. (A) A magnetic field sensing array on the surface of the leg tracks the positions of two magnetic beads implanted into the muscle. A feather microcontroller (µC) in the turkey feathers wirelessly transmits the magnetic field data to a magnet tracking computer that calculates and displays the magnetomicrometry (MM) signal in real time. The turkeys walked and ran on a treadmill while x-ray video cameras recorded synchronized fluoromicrometry (FM) data for post- processing. (B) Comparison of MM (blue) with FM (red) to validate the MM accuracy. These representative results during running gait show the submillimeter accuracy of MM during untethered muscle length tracking. . . 93 3.4 Untethered Muscle Tracking During Treadmill Running: Magnetomicrometry Versus Fluoromicrometry. Changes in muscle tissue length measured by MM (blue) and FM (greenred) for three turkeys at five speeds (30 s shown for each speed). The column to the right of the plots gives the coefficients of determination (R2) between magnetomicrometry and fluoromicrometry cor- responding to each turkey and speed. Gaps in the fluoromicrometry data are due to researcher selection of full gait cycles during which both magnetic beads were visible in both x-ray images. Gaps in the magnetomicrometry data (gray) are due to packet drops during wireless transmission of the magnetic field signals to the tracking computer (gaps below 50 ms interpolated in gray, gaps above 50 ms highlighted in gray). The turkey gait diagram below the plots shows the corresponding gait phases over one gait cycle. . . . . . . . . 94 3.5 Coefficients of Determination (R2values) between Magnetomicrometry and Fluoromicrometry. We compared all magnetomicrometry measurements (hor- izontal axis) across 50 gait cycles of turkey running (10 gait cycles at each of 5 speeds, for each bird) against time-synchronized, interpolated fluoromi- crometry measurements (vertical axis). Data are plotted in blue, orange, and purple for Birds A, B, and C, respectively. Coefficients of determination (R2values, shown in corresponding colors) for each bird were 0.952, 0.860, and 0.967, respectively.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 3.6 Kernel Density Estimates of Differences between Magnetomicrometry and Flu- oromicrometry. We subtracted time-synchronized, interpolated fluoromicrom- etry measurements from magnetomicrometry measurements across all 50 gait cycles of turkey running (10 gait cycles at each of 5 speeds, for each bird). Kernel density estimates show the distribution (vertical axis) of these differ- ences (horizontal axis), with data plotted in blue, orange, and purple for Birds A, B, and C, respectively. Vertical lines indicate mean offsets (in correspond- ing colors), with the mean offsets and standard deviations labeled. Adjusted standard deviations that compensate for fluoromicrometry noise (0.098 mm, standard deviation) indicate an estimate of the magnetomicrometry measure- ment noise. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 963.7 Study-Specific Limitations to the Use of Fluoromicrometry. (A) Ten gait cy- cles of raw fluoromicrometry video data (Bird A, 3.5 m / s), independently manually re-labeled three times, are shown in light green, dark green, and brown. (B) The histogram shows the distribution (vertical axis) of the vari- ance (horizontal axis) between the three manual processed fluoromicrometry signals at each timepoint of the data, with a kernel density estimation curve overlay. A vertical line indicates the average variance, 0.010 mm2, and the corresponding square root of the variance is 0.098 mm. . . . . . . . . . . . . 97 3.8 Magnetomicrometry Tracking Time Delay. The magnet tracking computer recorded the times it received the magnetic field data and the times it com- pleted the magnet tracking algorithm. The distribution (vertical axis) of the difference between these times (horizontal axis) is the tracking time delay and indicates the bandwidth at which magnetomicrometry can track the muscle tissue length. The data are shown as a stacked histogram, with blue, orange, and purple data corresponding to Birds A, B, and C, respectively. Data are from all turkey gait cycles used to compare magnetomicrometry against fluo- romicrometry. The 99th percentile time delay (t99%) is labeled for each bird. The ninety-ninth percentile time delay for all birds was less than one millisecond. 98 3.9 Muscle Tissue Length During Non-Synchronous Ramp Ascent and Descent. We used magnetomicrometry to track muscle tissue length during ramp ascent and descent at two inclines for all three birds. Data for each bird and each slope are synchronized at right leg toe strike (indicated by the vertical gray line) and normalized from toe strike to toe strike. Variability between curves reflects gait cycle variability during untrained ramp navigation. Muscle tissue length is plotted in blue for right leg stance, in red purple for right leg swing, and in gray where video did not allow gait-phase labeling. We recorded at least three gait cycles of each activity for each bird. . . . . . . . . . . . . . 99 3.10 Muscle Tissue Length During Non-Synchronous Vertical Ascent and Descent. We used magnetomicrometry to track muscle tissue length during vertical as- cent and descent at three heights for all three birds. Data for each bird and each height are synchronized at right leg toe-off (start of the aerial phase, indicated by the vertical gray line). Variability between curves reflects move- ment variability during untrained vertical ascent and descent. Muscle tissue length during contact with the ground is plotted in blue, and muscle tissue length during the aerial phase is plotted in redpurple. All data are shown, including scenarios in which significant wing-flapping occurred during jump up or down. We captured at least three recordings of each activity for each bird. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 3.11 Extent of Bead Positions During Ramp Ascent and Descent and Vertical Ascent and Descent. Analysis of the magnetic bead position data (shown) recorded during all variable terrain activities revealed the extent of the mag- netic bead positions during the tracking. We used these bead position maxi- mums to design a benchtop test to verify the validity of our magnetomicrom- etry measurements during the variable terrain activities. . . . . . . . . . . . 1013.12 Muscle Tissue Length During Free Roaming Movement. Magnetomicrometry data was continuously collected for 150 seconds during free roaming activity. Muscle tissue length is plotted in blue during standing and walking and plotted in black purple during running. Blue highlighted regions indicate muscle tissue length during (a) feather ruffling, (b) jumping, and (c) balancing on one leg. Gray arrows indicate when the turkey was turning left (left arrows) or turning right (right arrows). Gaps due to wireless transmission packet drops are shown in gray, as described in Figure 2.. . . . . . . . . . . . . . . . . . . . . . . . . 101 4.1 Magnet State Representation. The magnetic dipole model describes the mag- net’s state using a position vector (p⃗) and a dipole moment vector (m⃗). The position vector represents the location of the magnet in 3-dimensional Eu- clidean space, while the dipole moment vector characterizes the strength and orientation of the magnet’s magnetic field. . . . . . . . . . . . . . . . . . . . 109 4.2 Magnetic field strength as a function of distance from the magnet (at a dipole moment vector aligned with the position vector). The magnetic field strength generated by a magnetic dipole decays exponentially with increasing distance between the sensor and the magnet. This inverse cubic relationship poses challenges for Magnetomicrometry in extending the observable sensing depth. 111 4.3 (a) Illustration of the Kalman filter on manifolds, adapted from He et al

[0119] . The Kalman filter is augmented with manifold operations, such as retraction and inverse retraction, to ensure that the state estimate remains consistent with the underlying system dynamics and preserves the inherent symmetries of the problem. (b) Schematic representation of the invariant Kalman filter on Lie groups, adapted from Barrau et al

[0120] . The Kalman filter pipeline is integrated with non-Euclidean manifolds and Lie groups using homogeneous transformations between the manifold space and the homogeneous space, with visualization of accurate innovation process on manifold space. . . . . . . . 115 4.4 Representation of the S2manifold. The dipole moment vector of a magnetic dipole model lies on the surface of a sphere, which is mathematically described as the S2manifold. Solutions for the dipole moment vector are represented as points on this spherical surface. Revised from Alonso et al

[0121] . . . . . . 116 4.5 Demonstration of vector discontinuity in the S2manifold, adapted from Brossard et al

[0117] . The S2manifold is not parallelizable, which means that it is im- possible to define a continuous and non-vanishing tangent vector field on its surface. This property poses challenges for the direct application of Kalman filters on the S2manifold. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 4.6 Illustration of the dipole moment vector and the rotational velocity vector. The S1vector, representing the rotational velocity on the surface of the mag- netic dipole, is introduced to lift the S2manifold to the parallelizable SO3 manifold. This lifting process enables the definition of the spatial pose of the magnet using the SO3 manifold. . . . . . . . . . . . . . . . . . . . . . . . . 118 4.7 Static grid with variations in dipole moment orientations and positions. The grid consists of 192 configurations with different combinations of (a) dipole moment orientations, (b) XY positions, and (c) z-axis positional heights. . . 1244.8 Illustration of variance caused by sensor noise for different sensing methods. The collective results of magnet tracking from both the fast-LM algorithm and the proposed tracking algorithm are shown. The volumes of the measured covariance of the proposed tracking algorithm are generally smaller than those of the fast-LM algorithm, indicating improved tracking performance. . . . . . 125 4.9 Single bead static measurement covariance. The covariance ellipsoids from both the fast-LM algorithm and the proposed tracking algorithm are shown for each static point in a subset of the grid with an upright dipole moment vector. The four subplots correspond to different z-axis heights ((a)2cm, (b)3cm, (c)4cm, (d)5cm). The covariance ellipsoids from the proposed track- ing method consistently demonstrate smaller volumes compared to those from the fast-LM algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 4.10 Position errors of the static simulation. The proposed method demonstrates superior performance in static single magnet localization, with lower mean errors and standard deviations compared to the fast-LM algorithm. . . . . . 127 4.11 Dipole moment errors of the static simulation. The proposed method exhibits better performance in estimating the dipole moment orientation, with lower average errors and standard deviations compared to the fast-LM algorithm. . 128 4.12 Real-time simulation results for the dynamic simulation case 1 (single magnet linear transverse). The proposed tracking algorithm demonstrates stable per- formance throughout the trajectory, while the error covariance of the fast-LM algorithm increases exponentially as the magnet moves further away from the center. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 4.13 Real-time position and dipole moment vector orientation errors for the dy- namic simulation case 1. The proposed tracking algorithm maintains lower errors compared to the fast-LM algorithm, particularly towards the end of the trajectory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 4.14 Real-time simulation results for the dynamic simulation case 2 (single mag- net linear transverse and rotation). The proposed tracking algorithm shows increased errors at the end of the trajectory compared to Case 1, but still outperforms the fast-LM algorithm, which exhibits significantly larger errors due to the added rotational movement. . . . . . . . . . . . . . . . . . . . . . 132 4.15 Real-time position and dipole moment vector orientation errors for the dy- namic simulation case 2. The performance of the fast-LM method significantly degrades during rotational movement, while the proposed method maintains its superiority, albeit with slightly increased errors compared to Case 1. . . . 133 4.16 Real-time simulation results for the dynamic simulation case 3 (two magnets linear transverse and rotation). The general tendencies of the tracked trajec- tories from the two methods are consistent with the previous observations, with the proposed tracking algorithm demonstrating better performance than the fast-LM algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 4.17 Real-time position and dipole moment vector orientation errors for the dy- namic simulation case 3. The proposed tracking algorithm maintains lower errors compared to the fast-LM algorithm in the presence of multiple magnets with linear transverse and rotational motion. . . . . . . . . . . . . . . . . . . 1354.18 Real-time simulation results for the dynamic simulation case 4 (two magnets rotational transverse and rotation). The proposed tracking algorithm demon- strates superior performance compared to the fast-LM method, with relatively stable tracking throughout the trajectory. The error covariance of the fast-LM method increases as the magnets move further away from the center. . . . . 136 4.19 Real-time position and dipole moment vector orientation errors for the dy- namic simulation case 4. The proposed tracking algorithm outperforms the fast-LM method, maintaining relatively stable tracking throughout the tra- jectory. Additionally, the error covariance of the fast-LM method increases as the magnets move further from the center. . . . . . . . . . . . . . . . . . . . 137 4.20 Real-time simulation results for the dynamic simulation case 5 (three mag- nets rotational and linear transverse and rotation). The proposed tracking algorithm maintains its superiority over the fast-LM method, even with the increased complexity of the scenario. The error covariance of the fast-LM method further increases compared to the previous cases. . . . . . . . . . . . 138 4.21 Real-time position and dipole moment vector orientation errors for the dy- namic simulation case 5. The proposed tracking algorithm demonstrates sig- nificantly lower errors compared to the fast-LM algorithm, confirming its abil- ity to handle challenging scenarios with multiple magnets and complex motion patterns. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 4.22 Real-time simulation results for the dynamic simulation case 6 (four magnets rotational transverse and rotation). The proposed tracking algorithm main- tains its stability and accuracy despite the increased number of magnets and the complexity of their motion. The fast-LM method suffers from even larger error covariances and higher position and dipole vector orientation errors. . . 140 4.23 Real-time position and dipole moment vector orientation errors for the dy- namic simulation case 6. The proposed tracking algorithm consistently out- performs the fast-LM algorithm, achieving significantly lower errors through- out the trajectory.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 4.24 Real-time simulation results for the dynamic simulation case 7 (two magnets random transverse and rotation). The proposed tracking algorithm demon- strates its robustness and superiority, maintaining its stability and accuracy despite the increased complexity and randomness of the motion. The fast-LM method exhibits significant error covariances and high position and dipole vector orientation errors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 4.25 Real-time position and dipole moment vector orientation errors for the dy- namic simulation case 7. The proposed tracking algorithm consistently pro- vides more accurate and reliable tracking results compared to the fast-LM algorithm, even in highly complex and unpredictable scenarios with random motion of the magnets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1434.26 Real-time simulation results for the dynamic simulation case 8 (three magnets random transverse and rotation). The proposed tracking algorithm demon- strates its robustness and superiority, maintaining its stability and accuracy despite the increased complexity and randomness of the motion. The fast-LM method exhibits significant error covariances and high position and dipole vector orientation errors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 4.27 Real-time position and dipole moment vector orientation errors for the dy- namic simulation case 8. The proposed tracking algorithm consistently pro- vides more accurate and reliable tracking results compared to the fast-LM algorithm, even in highly complex and unpredictable scenarios with random motion of the magnets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 4.28 Comparison of tracking algorithms on the case 5 (three magnets rotational and linear transverse and rotation). The traditional KF methods with 5-DoF state representation (EKF and UKF) demonstrate vulnerabilities in tracking the magnetic dipole moment vectors, resulting in larger biases in the estimated positions. The fast-LM algorithm and the proposed tracking algorithm show better performance, with the proposed method exhibiting the highest accuracy and robustness. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 4.29 Real-time position and dipole moment vector orientation errors for the com- parison of tracking algorithms on the case 5. The simulation results empha- size the importance of proper representation of the dipole moment vector and the appropriate process model for estimation and integration. The proposed tracking paradigm demonstrates its effectiveness and efficacy compared to the other tracking methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 4.30 Magnetomicrometry Sensing Error vs Implantation Depth of Two Implantable Magnetic Beads. (A) We affixed two magnetic beads (manufactured as de- scribed in the methods) into 1x1 LEGO round plates, and we attached these two round plates 40 mm apart to a 1x6 LEGO technic brick. We then cen- tered this LEGO brick over a 96-element (LIS3MDL) magnetic field sensing array parallel with the long axis of the array. We varied the depth in 3.2 mm increments from the initial depth of 11.5 mm by adding 1x6 LEGO plates. This centering presents a best-case scenario, suggesting maximum depth lim- its that could be achieved with proper placement of magnetic field sensors. (B) Seaborn strip plots show the error in mm (vertical axis) of the magne- tomicrometry signal for each depth. The horizontal axis shows each depth categorically. Note that the magnetomicrometry best-case accuracy with this sensing array and these magnets is sub-millimeter to a depth of about 33 mm. 148 4.31 Real-world benchmark data. (a) Tracking results for a magnet separation distance of 24 mm. (b) Tracking results for a magnet separation distance of 40 mm. The proposed tracking algorithm demonstrates reduced standard deviations of the measurements compared to the fast-LM algorithm, indicating improved tracking performance on real-world data. . . . . . . . . . . . . . . 1494.32 Practical implication demonstration. With the proposed tracking approach, it is hypothesized that the viable maximum sensing depth can be extended by approximately 1 cm compared to the current physical sensing frameworks (us- ing 3-mmmagnets and the sensor array from the previous Magnetomicrometry study). This extension in sensing depth would enable Magnetomicrometry to observe deeper tissues and be used with thicker prosthetic socket systems in bionic applications. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152List of Tables 2.1 Electrode Design Specifications . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.2 Electrode Material Specifications . . . . . . . . . . . . . . . . . . . . . . . . 38 2.3 Reported Comfort and Preference Scores . . . . . . . . . . . . . . . . . . . . 54 2.4 Generated Impulse Artifact Correlations . . . . . . . . . . . . . . . . . . . . 68 2.5 Generated Impulse Artifact Correlations . . . . . . . . . . . . . . . . . . . . 68 4.1 Static and Dynamic Magnet Manifold Representation . . . . . . . . . . . . . 119 4.2 Manifold Retraction and Inverse-Retraction . . . . . . . . . . . . . . . . . . 120Chapter 1 Introduction 1.1 Background Bionic devices and assistive technologies have significantly advanced in recent years, offering new possibilities for individuals with amputation or motor impairments[1], [2]. These devices aim to restore or enhance motor function by providing volitional control capabilities that mimic natural movements[1], [2]. However, the efficacy of these devices heavily relies on accurate muscle sensing modali- ties that can detect and interpret user intent in real-time. Muscle sensing technologies have not only deepened our understanding of biomechanics but have also played a crucial role in enhancing the performance of bionic devices. By mapping the electrical and mechanical activities of muscles, these technologies enable volitional control of prostheses and assistive devices, allowing users to interact with their environment more naturally[1]–[5]. As the de- mand for advanced bionic solutions continues to grow, there is a pressing need to explore and improve muscle sensing modalities that can provide reliable and intuitive control interfaces.Figure 1.1: Schematics of myo-neural interface using sEMG and Magnetomicrometry for neural control of advanced Bionics devices. 1.2 Objective The main objective of this dissertation is to present and evaluate myo-neural sensing paradigms that can improve volitional control capabilities of modern bionic devices. Specifically, this work focuses on two promising candidates for myo-neural interfaces: surface electromyogra- phy (sEMG) and Magnetomicrometry, as shown in Fig. 1.1[3], [6]. sEMG is a non-invasive technique that measures the electrical activity of muscles, while Magnetomicrometry is an emerging modality that tracks the mechanical states of muscles using implanted magnets. The dissertation aims to investigate methods to enhance the utility and efficacy of these sensing modalities in the context of bionic applications. 1.3 Specific Aims To achieve the objectives outlined above, this dissertation focuses on the following specific aims: • Presenting a new paradigm for sEMG acquisition that incorporates flexible and low- profile electrodes, advanced filtering algorithms, and synchronized ultrasound measure-ments to enhance its utility in dynamic bionic applications. • Demonstrating the in-vivo and real-time muscle tracking capabilities of Magnetomi- crometry in freely-moving animals performing various motor activities, validating its potential as a reliable myo-neural interface. • Introducing a novel theoretical framework for the magnet tracking problem, leveraging Riemannian geometry and manifold theory to improve the performance of magnet tracking system. 1.4 Thesis Organization The thesis is organized into five chapters. Chapter 2 focuses on enhancing utilities of sEMG sensing modality for bionic applications by presenting a novel electrode design, an impulse artifact filtering algorithm, and a method for synchronized sEMG and ultrasound acqui- sitions. Chapter 3 demonstrates the in-vivo validation of Magnetomicrometry, showcasing its real-time muscle tracking capabilities in untethered, freely-moving animals. Chapter 4 introduces a new analytic formulation method for the magnet tracking problem, aiming to improve the resolution and accuracy of magnet tracking framework and Magnetomicrom- etry. Finally, Chapter 5 summarizes the key contributions of this work and discusses the implications and future directions for myo-neural interfaces in bionic devices.Chapter 2 Enhancing sEMG Utility for Bionics This chapter discusses various methods to enhance and expand the utilities of sEMG sensing modality in Bionics applications. First, we introduce sEMG and its role in the control of pros- thetic devices. I first present a novel flexible and low-profile electrode design that improves user comfort and signal quality making it suitable for long-term use within prosthetic sock- ets. We then propose a real-time signal processing algorithm to filter out impulse artifacts that often contaminate sEMG signals during dynamic movements. Finally, we demonstrate a synchronous sEMG and ultrasound acquisition method that enables the simultaneous assess- ment of muscle electrical activity and mechanical deformation, providing valuable insights into muscle function and control. By addressing these challenges and presenting innovative solutions, this chapter aims to expand the capabilities of sEMG as a reliable, versatile, and effective myo-neural interface for bionic devices. 2.1 Introduction Electromyography (EMG) is a technique used to record and analyze the electrical activity of muscles, known as myoelectric signals[3]. Surface electromyography (sEMG), the least inva- sive EMG method, has been extensively researched for its potential in clinical applications, particularly in the control of prosthetic devices [1], [3], [7]Over the past few decades, advancements in sEMG have enabled the development of various applications for neural control in upper and lower extremity prosthetic devices[1], [7]. The introduction of novel surgical techniques, such as Targeted Muscle Reinnervation (TMR) and the Agonist-Antagonist Myoneural Interface (AMI), has further enhanced the efficacy and utility of sEMG-based control approaches in advanced neuroprosthetic systems[8], [9] As the utility of sEMG continues to grow, it is crucial to address the challenges associated with its implementation in real-world applications. These challenges include the need for reliable and comfortable electrode interfaces, effective signal processing techniques to miti- gate artifacts, and the integration of sEMG with other sensing modalities to provide a more comprehensive understanding of muscle function. In the following sections, we present novel methods and findings to address these challenges and enhance the utility of sEMG for bionic applications. 2.2 Interface for Within Socket sEMG Acquisition: Flex- ible and Low-Profile Electrodes This work was supported by the MIT Media Lab Consortia and Eunice Kennedy Shriver National Institute of Child Health & Human Development of the National Institutes of Health (NIH) under award number R01HD097135. This work was presented in Yeon et al

[0010] ,

[0011] . For persons with lower extremity (LE) amputation, acquisition of surface electromyog- raphy (sEMG) from within the prosthetic socket remains a significant challenge due to the dynamic loads experienced during the gait cycle. However, these signals are critical for both understanding the clinical effects of LE amputation and determining the desired control trajectories of active LE prostheses. Current solutions for collecting within-socket sEMG are generally (i) incompatible with a subject’s prescribed prosthetic socket and liners, (ii) uncomfortable, and (iii) expensive. This study presents an alternative within-socket sEMG acquisition paradigm using a novel flexible and low-profile electrode. First, the practical per-formance of this Sub-Liner Interface for Prosthetics (SLIP) electrode is compared to that of commercial Ag / AgCl electrodes within a cohort of subjects without amputation. Then, the corresponding SLIP electrode sEMG acquisition paradigm is implemented in a single subject with unilateral transtibial amputation performing unconstrained movements and walking on level ground. Finally, a quantitative questionnaire characterizes subjective comfort for SLIP electrode and commercial Ag / AgCl electrode instrumentation setups. Quantitative analyses suggest comparable signal qualities between SLIP and Ag / AgCl electrodes while qualita- tive analyses suggest the feasibility of using the SLIP electrode for real-time sEMG data collection from load-bearing, ambulatory subjects with LE amputation. 2.2.1 Motivation While sEMG has been widely applied to control upper extremity neuroprostheses, rela- tively fewer studies have explored sEMG control of lower extremity (LE) neuroprostheses, a disparity that Windrich et al. attribute to unique implementation challenges within the application domain

[0012] . This study specifically addresses sEMG acquisition within the LE prosthetic liner, a task whose difficulties are most pronounced during ambulatory activity when a residual limb experiences extreme dynamic pressure changes that cause structural deformation, electrode contact variation, pain, and additive sEMG signal noise. A solution for within-socket sEMG acquisition may inform subsequent solutions for downstream chal- lenges, including the need for a high-powered portable actuator design, a suitable neural control paradigm for ambulatory activity, and an optimal prosthetic socket design

[0013] . Several previous studies have successfully used sEMG signals from the residual limb to de- termine unconstrained neuroprosthetic movements in free space

[0014] –

[0017] . To further explore sEMG as a control signal during ambulation, others have attempted to place commercial pre-amplifiers inside the prosthetic liner within the socket

[0018] –

[0020] . While this approach has been shown to successfully acquire sEMG data, the bulkiness of the pre-amplifiers creates painful indentations on the skin which make it unsuitable for long-term use and translationbeyond clinical trials. Another alternative paradigm involves modifying existing or pro- totyping experimental prosthetic sockets and liners with embedded, minimally-protruding electrodes

[0021] –

[0026] . However, the necessity of fabricating a customized socket and liner for each individual subject is often both too time-consuming and prohibitively expensive for practical purposes. In consideration of these engineering trade-offs, we propose an efficient and cost-effective method for sEMG acquisition within prosthetic sockets using a novel, flexible, and low-profile electrode design. Importantly, the design is intended to be compatible with a user’s existing socket and liner without modification. In this study, we demonstrate the capability of these custom Sub-Liner Interface for Prosthetics (SLIP) electrodes to provide useful sEMG signals from the residual limb under dynamic load-bearing conditions. The flexible and thin nature of the custom electrodes minimizes skin indentation within the socket, especially compared to the methods with pre-amplifier dependencies from previous literature. Furthermore, SLIP electrodes can be readily fabricated with flexible printed-circuit-board (fPCB) technology at low cost. As a research tool, the SLIP electrode design may faciliate accelerated development of neuroprosthetic systems as a whole, and while not yet clinical grade, it introduces a general paradigm of sEMG collection which may have potential applications in translational neuroprosthetic systems due to the engineering advantages discussed. 2.2.2 Design Design Criteria The SLIP electrode is designed to overcome current challenges regarding utilization of sEMG- based control of Bionic prostheses in clinical trials. Specifically, the SLIP electrode design addresses three technical challenges: • Compatibility with standard sockets and liners : As discussed in Section 2.2.1, sEMG ac- quisition using custom subject-specific liners and sockets incurs high costs and requiressignificant lead times. Developing an electrode technology compatible with traditional sockets and liners would enable researchers to conduct experiments more efficiently from both of these perspectives. • Subject comfort : Prior research applications using commercial electrode technology to measure sEMG signals inside load-bearing sockets were primarily limited by subjects’ discomfort and pain. Subject comfort is critical in both research and translational applications

[0018] –

[0020] . • Manufacturability : Using standard fPCB electronics fabrication technologies to pro- duce the SLIP electrode itself minimizes cost and lead time when ordering in bulk from readily-available vendors. These design requirements for compatibility, comfort, and manufacturability are necessary for increasing research productivity and decreasing the cost of human subject experiments. Related Technologies Several research attempts have been made to develop flexible and dry surface electrodes for electroencephalogram (EEG), electrocardiogram (ECG) and EMG acquisition

[0027] –

[0035] . This push for flexible and dry electrode technology is driven by the desire for long-term biopoten- tial monitoring that extends beyond the lifespan of typical wet Ag / AgCl electrodes

[0027] –

[0032] . By using various substrate materials such as polydimethysiloaxane (PDMS) and parylene, these efforts successfully demonstrated that flexible and dry electrodes can obtain signal qualities comparable to Ag / AgCl electrodes. A couple of these attempts were also able to develop electrodes with minimal thickness (< 300um)

[0032] ,

[0033] . Nonetheless, most of these technologies remain too fragile or difficult to manufacture and are not directly translatable to dynamic Bionics prosthesis applications.Design Concept A conceptual diagram of the SLIP electrode design, including form factor and usage, is shown in Fig. 2.1 with additional design details being provided in the following subsections.Figure 2.1: SLIP electrode design concept showing electrode worn against the skin of the residuum inside a standard silicone prosthetic liner. An integrated, flexible cable and con- nector allows the sensor to interface with an acquisition and processing system mounted externally on the socket. The SLIP electrode was designed for insertion between a subject’s residuum and pros- thetic liner. As thinner profiles were hypothesized to provide better user comfort and expe- rience, electrode profile height was minimized. Given limitations in manufacturing capabili- ties, the profile height of the prototype electrode is slightly less than 100 µm. An insulated trace originating from the electrode contact heads leads to a terminal electrical connector that interfaces with a portable sEMG amplifier and processor mounted externally on the socket. This feature allows subjects to use their own prosthetic socket and liner without modifications.Electrical and Mechanical Design The shape of the electrode head was designed based on bipolar sEMG electrode guidelines from the European project Surface EMG for Non-Invasive Assessment of Muscles (SENIAM)

[0036] . The electrical CAD model of the electrode skin contact is shown in Fig. 2.2. Each bipolar electrode is circular with a diameter of 10 mm

[0036] ,

[0037] . An inter-electrode distance of 16 mm was chosen after considering the effects of inter-muscle crosstalk and differential sEMG signal power

[0038] . The electrical connections inside the electrode were routed according to differential signal routing standards for maximizing common mode rejection. Only a single layer of copper metal was utilized in order to minimize thickness and maximize flexibility of the electrode. Thread holes were placed on the electrode head and around the lead wire for potential mechanical integration with a prosthetic liner or sock by sewing. This approach would provide a semi-permanent installation for ease of donning and doffing the electrodes during repeated trials.Figure 2.2: Electrical model of the novel electrode. Each bipolar electrode has a diameter of 10 mm and inter-electrode distance of 16mm. Mechanical holes on the electrode were placed for potential anchoring and mounting inside a prosthetic liner or sock. Table 2.1 summarizes the design specifications of the SLIP electrode. Notably, the con- nector type and length of the electrode can be varied based on application. These design decisions are not necessarily optimal, but do represent efficient and convenient selections.Table 2.1: Electrode Design SpecificationsMaterials and Manufacturability Materials were chosen for compatibility with standard flexible PCB manufacturing tech- niques to increase accessibility for researchers. Details of the selected materials are shown in Table 2.2. Table 2.2: Electrode Material SpecificationsPolyimide (PI) film was used as a base substrate. PI film is one of the most commonly used insulating substrates in fPCB fabrication. Material safety data sheets for PI report that "no skin irritation is expected from handling film

[0039] ." The copper layer is gold-plated where exposed to improve sensor compatibility with the skin. While a thorough compatibility evaluation of the novel electrode with skin is beyond the scope of this paper, a growing body of literature suggests electrodes with PI substrates and gold contacts, including electrode implants, are compatible with long-term usage

[0040] ,

[0041] .Manufactured Prototype The prototype SLIP electrode was fabricated by an fPCB vendor (PCBWay, Xiacheng Dis- trict, Hangzhou, China) using generic fPCB fabrication techniques. The fabricated electrode is shown in Fig. 2.3. The manufactured flexible electrode costs less than $6.00 per sensor when purchased in quantities of 30. The per-electrode cost can be greatly reduced if pur- chased in larger quantities.Figure 2.3: Fabricated prototype of the novel electrode. (a) Size comparison of the electrode. (b,c) Demonstration of manufactured electrode flexibility.2.2.3 Quantitative Evaluation Subject Recruitment and Study Preparation To evaluate the measurement characteristics of the SLIP electrode, five subjects (18-22 years old; three male, two female) with no amputation and no self-reported neuromuscular patholo- gies were recruited. Subjects provided written informed consent through MIT Committee on the Use of Humans as Experimental Subjects (COUHES) protocol #1906898371. In a cross-validation study, a pair of commercially-available Ag / AgCl electrodes was placed alongside the SLIP electrode to measure bipolar sEMG from the left leg’s tibialis anterior (TA) during free space movements, as seen in Fig. 2.4. Ag / AgCl electrode chemistry was selected for comparison against the SLIP electrode due to its standardized use in sEMG acquisition

[0036] . The SLIP electrode was fixed to the skin with adhesive tape before applying uniform pressure to all electrodes using an elastic exercise band, as seen in Fig. 2.4b. The location of the muscle belly was estimated through palpation, and Ag / AgCl inter- electrode spacing was matched to the fixed spacing of the SLIP electrode. Because sEMG measurements are sensitive to the locations of electrodes relative to underlying muscle, the commercial wet electrodes’ and the SLIP electrode’s relative positioning on the limb was alternated between medial and lateral positions over the estimated location of the TA belly from subject to subject, as inspired by a study comparing two types of wet surface electrodes by Posada-Quintero et al.

[0042] . The TA was selected for analysis due to its proximity to the skin surface and corresponding ease of access. The low amount of subcutaneous fat superficial to the TA minimized the attenuation of sEMG signals for all subjects. sEMG Instrumentation Specification of the sEMG acquisition hardware is informed by qualitative assessment of generic electrical characteristics pertinent to the SLIP electrode. The SLIP electrode has a dry metal contact surface that results in higher skin-electrode impedance compared toFigure 2.4: Example electrode configuration across the TA with Ag / AgCl electrodes in the medial position. An adhesive piece of tape is applied on top of the flexible electrode to ensure proper placement. Elastic exercise band applying uniform pressure to all electrodes. that of Ag / AgCl electrodes

[0027] . Because of the dry interface between SLIP electrode and skin, the impedance of the system can vary with dynamic differential pressure changes that result in DC voltage swings an order-of-magnitude higher than the amplitude of the EMG signal. Moreover, due to the proposed acquisition architecture, raw sEMG signals need to be transferred over a relatively long-distance (> 40 cm) before interfacing with the analog-front- end (AFE) IC. Due to the single metal layer within the electrode, various methods to protect the sEMG signal through electrical shielding were not available for our applications

[0043] . Therefore, the AFE with a DC-coupled differential pre-amplifier (AD8422., Analog De- vices, Norwood, MA, USA) with high input impedance (Zin= 200GΩ ∥ 2pF ) and low gain (G = 10) was utilized in order to mitigate the high skin-electrode impedance and large dy- namic DC voltage swings anticipated. Along with the pre-amplifier, a 24-bit high resolution analog-to-digital-converter (ADC) IC (ADS1299., Texas Instrument, Dallas, TX, USA) was utilized for sEMG signal acquisition. Fig. 2.5 shows the custom embedded sEMG acquisitionplatform used in this experiment session

[0044] .Figure 2.5: Embedded 8-channel sEMG acquisition and processing platform used for evalu- ation with example Ag / AgCl bipolar pair, SLIP electrode, and ground electrode attached. Data Collection While seated and starting from a neutral ankle position, subjects performed five trials each of ankle dorsiflexion and ankle cocontraction in free space to generate raw bipolar sEMG data. Ankle dorsiflexion trials required subjects to fully dorsiflex their left ankle and hold the position. Ankle cocontraction trials required subjects to hold their left ankle in a neutral position while cocontracting associated muscles to increase the stiffness of the joint. In chronological order, each guided trial consisted of three seconds of rest, three seconds of activity, and three seconds of rest. Within these trials, raw sEMG signals were sampled and logged at a frequency of 1 kHZ. Data Analysis Data collected from commercial Ag / AgCl electrodes and the SLIP electrode were analyzed in the time domain to determine their interchangeability for the purposes of neuro-prosthesiscontrol. Namely, we measured muscle activation onset / offset times and Pearson’s correla- tion of discrete muscle activation envelopes between electrode types in a manner similar to Posada-Quintero et al.

[0042] . The practical control implication of identical onset / offset times and a unity Pearson’s correlation would be a neuroprosthesis that is identically controlled using either type of electrode. To obtain a linear envelope of sEMG roughly corresponding to muscle activation, all signals were first forward-backward filtered using a 4th-order Butterworth filter designed in MATLAB ver. R2019b (MathWorks, Inc., Natick, MA, USA) with a passband of 80 - 400 Hz. Passband thresholds were chosen to remove 60 Hz mains hum, motion artifacts, and high frequency noise

[0045] . Signals were then rectified and subjected to a moving average filter with a window length of 150 ms. Normalization was performed within each trial on each channel’s resulting waveform by dividing all voltages by the maximum voltage recorded on that channel. Applying a moving average filter to processed sEMG for an online estimate of muscle activation in this manner is widely used within the domain of neuromuscular modeling. • Muscle on- and off-times : The linear envelope of sEMG was used to estimate on- and off-times for muscle contraction during each trial. A single threshold was used to determine both time indices, specified by Threshold = µ+ 6σ (2.1) where µ and σ are the mean and standard deviation of the linear envelope during the first second of rest at the beginning of each trial. Due to the orders of magnitude difference between quiescent noise and true muscle activation, a 6σ threshold was chosen to conservatively define an absolute noise floor and avoid false detection of onset muscle activation. Onset and offset timings of the linear envelope measured by the SLIP electrode relative to the Ag / AgCl electrodes were compared in a Bland-Altman plot. Individual trials were systematically excluded from analysis if the linear sEMG envelope exceeded four seconds. These trials typically contained low amplitude sEMG peaks after the main sEMG envelope that made offset time indeterminate, indicative of the subject’s failure to fully relax according to instructions. • Pearson’s correlations : For each trial, Pearson’s correlations were calculated between the linear envelopes of sEMG measured by the SLIP electrode and Ag / AgCl electrodes between average on-time and average off-time. Quantitative Results Relative muscle on- and off-times suggest a high degree of agreement between the SLIP electrode and commercial Ag / AgCl electrodes, as seen visually in Fig. 2.6. Representative linear envelopes of muscle activity (α) measured by both commercially available Ag / AgCl electrodes and the SLIP electrode are shown in Fig. 2.7. Envelope correlations between Ag / AgCl and flexible electrodes also suggest high degree of agreement. The group average Pearson’s correlation coefficient for qualified dorsiflexion trials (n = 22) was 0.95 ± 0.92 (average ± std. dev.). The group average Pearson’s correlation coefficient for qualified cocontraction trials (n = 21) was 0.94 ± 0.04. 2.2.4 Experiment 1: sEMG Acquisition within Prosthetic Socket Systems Subject Recruitment To evaluate the efficacy of the novel electrodes within a load-bearing prosthetic socket, a single subject (43 years old, female) with unilateral left-side transtibial amputation was recruited. The subject’s amputation was performed under the novel Agonist-antagonist Myoneural Interface (AMI) surgical paradigm, forming two dynamically coupled agonist- antagonist muscle pairs—tibialis anterior (TA) paired with lateral gastrocnemius (LG), andtibialis posterior (TP) paired with peroneous longus (PL)—within the residuum

[0018] . The subject provided written informed consent through MIT COUHES protocol #1812634918. Physical Placement of the Electrode Proper physical placement of the electrode is required for user comfort as well as high signal quality. Specifically, because the SLIP electrode interfaces with the pre-amplifer outside of the socket after a relatively long transmission distance, it is important to isolate the signalfrom additional parasitic electrical elements. Because AMI muscles are readily palpable and precisely locatable, only four SLIP electrodes were required for this specific experimental section. However, for demonstrative purposes, Fig. 2.8 presents a generic skin preparation and fully-instrumented electrode placement process. Elastic Kinesio tape and hydrocolloidal matrix bandages were used to ensure good me- chanical anchoring of the electrode on the skin within the liner environment. The compliance of the bandages is hypothesized to support continuous contact between electrodes and skin surface during dynamic pressure changes so as to minimize skin-electrode impedance. Additionally, the wires of the electrodes were necessarily routed along the longitudinal plane of the residual limb toward the proximal end of the liner. In areas where the residual limb and liner experience bending, such as the knee joint, it is crucial to minimize friction between leads and skin. To mitigate this issue, a prosthetic sock was placed between the electrode wires and skin with electrode leads passing through small incisions in the sock. sEMG Instrumentation and Data Processing A custom embedded sEMG acquisition system was utilized to measure signals from within the LE prosthetic liner in this experiment session as shown in Fig. 2.9. The embedded sEMG system is an improved version of the sEMG system from Section 2.2.3, performing simultaneous measurement of 16 channels at an increased sampling rate of 2 kHz

[0044] . The sEMG system retains an identical AFE structure compared to the previous version. For demonstrative purposes, Fig. 2.10 presents a fully-instrumented example use case wherein electrode leads exiting the prosthetic liner are interfaced with the sEMG system through the FFC connector adapter. The sampled sEMG signals were band-pass filtered with a pass-band of 80 ∼ 340 Hz, rectified with a root-mean-square (RMS) operation, and summed within a 200 ms time window. When processing sEMG signals collected during ambulation, additional cumulative histogram filtering (CHF) process was performed to mitigate artifacts from dynamic groundFigure 2.8: Example use case of the SLIP electrodes within an LE prosthetic liner. The 16 SLIP electrodes are placed and attached on the residual limb. This figure does not represent the actual instrumentation setup of the Experiment B. (a,b) Preparation of the residual limb’s skin surface. Corresponding target musculature are labeled and bony anatomical landmarks are marked upon the skin’s surface. (c,d) Electrodes are placed on the surface of the residual limb using Kinesio tape. (e) Leads of the SLIP electrodes are routed through the prosthetic sock with modified holes. (f) Fully routed SLIP electrodes under the prosthetic liner. impacts

[0046] . A cumulative short-time histogram within a summation time window was first calculated, and a subsequent RMS operation was applied to the data points that fell withinExperiment Design and Data Collection For this experiment specifically, four of the flexible SLIP sEMG electrodes were applied to the residual limb to instrument the targeted AMI muscle pairs without channel redundancy. In the first session, the subject was asked to rotate her phantom limb in a clockwise direction while seated and actively contracting each AMI muscle pair. Then, the subject was asked to walk on an instrumented, split-belt treadmill (Bertec Corporation, Columbus OH) using a passive prosthetic ankle at a speed of 1.4 m / s. Two independent force plates under the split-belt treadmill recorded ground reaction force (GRF) data at a sampling rate of 1 kHz. The GRF data for the affected left side were subsequently processed with a 10 ms window median filter. The GRF for the intact right side are not reported here. Raw sEMG data were wirelessly logged during the experiment via a WIFI module (RN42., Microchip Technology, Chandler, AZ) at a frequency of 2kHz. Fig. 2.11 shows the sEMG acquisition setup for ambulation with the SLIP electrodes connected to the liner-mounted sEMG system. Qualitative Results Stationary Voluntary Movement The normalized sEMG signals measured from the residual limb during voluntary rotation of the phantom ankle and subtalar joints are shown in Fig. 2.12. A distinct sEMG activation pattern was clearly observed with high signal-to- noise ratio, as can be seen in the contrast between periods of rest and periods of activity. For each cycle of rotation, reciprocal activation patterns for each AMI muscle pair were observed which correlated with rotation of the phantom limb. The distinct patterns and high signal quality observed can be reasonably considered sufficient for various applications of sEMG, including closed-loop control of neuroprosthetic devices. Ambulation The normalized sEMG signals from the residual limb of the subject during ambulation are shown in Fig. 2.13. Measured ground reaction forces correlated with sEMG signals showed clear swing and stance phases within the gait cycle. During the gait cycle,repetitive muscle activation patterns for each muscle were clearly observed in the recorded data. Notably, sEMG activation patterns of an amputated residual limb are not identicalto those of an intact limb, and the subject was additionally amputated with the novel AMI amputation paradigm

[0018] ,

[0021] . In these results, the ground truth of the desired muscle activation patterns from AMI musculature is unknown, and direct comparison and evaluation of the sEMG signals based upon physiological musculature is challenging. However, cyclic and repetitive activation patterns synchronized with gait cycle ground reaction forces suggest the potential for their use in neuroprosthetic control.Figure 2.13: Signals measured using the 16-channel sEMG processing platform with four electrodes over target AMI muscle pairs during treadmill walking at a speed of 1.4 m / s. The measured raw sEMG signal, processed sEMG signal, and ground reaction force (GRF) data are shown. Cyclic and repetitive muscle activation patterns correlating with GRF are demonstrated. Stance phases are presented with a white background while swing phases are presented with a gray background.Post-Ambulation Skin Evaluation As stated in Section 2.2.2, subject comfort was one of the critical criteria when designing the SLIP electrode. Fig. 2.14 shows the appearance of the residuum after doffing the prosthetic socket and liner following over an hour of load-bearing activities. Purple marker outlines on the skin indicate the locations of the placed electrodes. Given the well-documented challenges to comfortably access muscles inside load-bearing sockets, this preliminary finding demonstrates minimal, if any, observable skin indentation due to the novel electrodes.Figure 2.14: Skin of the residual limb after one hour of weight-bearing ambulation. Skin indentation and irritation are not visible where electrodes were placed. 2.2.5 Experiment 2: Comfort Evaluation Subject Recruitment and Experiment Design To determine the practical feasibility of collecting sEMG data using the commercial Ag / AgCl electrodes within a load-bearing prosthetic socket, the same subject (43 years old, female) with unilateral transtibial AMI amputation was recruited for a third experimental sec- tion. This set of experiments was designed to provide a quantified comparison of relative comfort. The subject provided written informed consent through MIT COUHES protocol #1812634918.For each of the subsequent three experiments, four bipolar pairs of electrodes in various configurations were acutely placed superior the muscle bellies of the TA, LG, TP, and PL muscles. Three electrode configurations were investigated, involving: 1) the SLIP electrode, 2) the Ag / AgCl disposable electrode alone, and 3) the Ag / AgCl electrode with requisite connector and wires, respectively, as shown in Fig. 2.15. Configuration 2 represents an ideal use case involving commercial wet electrodes in which customized wire leads are fully embedded into the socket liner. Configuration 3 represents a typical use case involving commercial wet electrodes with commercial wire leads that are connected to an external sEMG recording system.Figure 2.15: Electrode configurations for user comfort evaluation: (a) SLIP electrodes; (b) Commercial Ag / AgCl electrodes without connectors and wires; (c) Commercial Ag / AgCl electrodes with connectors and wires. For each electrode configuration, a progressive series of tasks was presented to the user to test comfort in typical real-world scenarios. First, the patient was asked to don her prosthesis with the electrodes attached to the limb and remain seated for 15 mins. After 15 mins, the subject was asked to perform multiple motor tasks, including sitting, standing, walking (1.4 m / s), stair ascent, and stair descent. Finally, the subject was asked to report her subjective comfort and level of socket suspension, and preference for daily life usage on a 1-5 scale. If the subject indicated any significant discomfort at any point, the specific trial was stoppedimmediately and a ‘1 = not applicable’ score was given for that section. Reported Comfort and Preference Scores The comfort and preference scores reported by the subject are given in Table 2.3. The reported scores showed that the SLIP electrodes caused no noticeable differences in com- fort compared to wearing the prescribed prosthesis without electrodes over a wide range of movements. The subject reported that even detecting the presence of the electrodes was "extremely challenging", and she would be fully willing to use SLIP electrodes if they were able to provide volitional control of her prosthesis. Notably, Configuration 2, which involved only Ag / AgCl electrodes, elicited only a slight increase in discomfort. However, the subject reported the ability to feel each of the conventional Ag / AgCl electrodes within socket. Fur- thermore, the subject noted that these Ag / AgCl electrodes were applying pressure to the point of affecting circulation within the residual limb. During testing of Configuration 3, the subject reported severe discomfort during initial sock donning and was not able to proceed with any of given tasks. Table 2.3: Reported Comfort and Preference Scores2.2.6 Discussion This study investigated an effective and economical method of acquiring sEMG signals from the weight-bearing LE residual limb within the donned prosthetic liner and socket. The unique challenges presented by dynamic loading inside the prosthetic socket during gait have limited the applications of traditional sEMG acquisition techniques in the LE domain. Here, some of these challenges were addressed using flexible, low-profile SLIP electrodes. In a quantitative comparison of five subjects without amputation or LE pathology, SLIP electrodes produced muscle activation signals which were interchangeable with those recorded with commercial Ag / AgCl electrodes. When tested in a subject with an LE amputation during walking, the SLIP electrodes were found to be compatible with her prosthetic socket and liner, and did not cause any discomfort or skin indentation. With the caveat that the present walking study included only one single subject who had received an AMI transtibial amputation, it was remarkable to observe periodic activity of antagonistic residual limb muscle pairs (TA-LG and PL-TP) that appeared to be entrained to the gait cycle. In contrast, a previous study reported muscle activation patterns that were aperiodic despite subject entrainment into a gait cycle in subjects who had undergone standard transtibial amputations

[0019] . Additionally, due to their flexible and low-profile nature, the SLIP electrodes were shown to be compatible with the subject’s prescribed lower extremity prosthetic socket and liner without any discomfort and skin indentation. These characteristic of the SLIP electrodes allowed for confident measurement of signals which are potentially appropriate for neuroprosthetic control. However, additional extensive studies are still required to quantify their sensitivity to impact artifacts and understand the physiological basis for the observed sEMG patterns. While the proposed SLIP electrode sEMG acquisition paradigm tentatively demonstrates utility in LE prosthetic applications, future work should expand to consider the compatibil- ity of the paradigm in other,varied domains and application contexts. Due to the uniquerequirements of LE prosthetic systems, the SLIP electrode’s engineering concessions include its dry nature, potential susceptibility to mechanical perturbation, and lack of electrical shielding. Thus, as mentioned in Section 2.2.2, the SLIP electrodes require proper pairing with a high-end sEMG processor AFE. Although only four channels were presented for analysis during the gait portion of this study, the potential exists for extracting additional information from all sixteen channels. In particular, because measured signals are sensitive to environmental variables such as elec- trode placement and sweat, among other factors, machining learning techniques that do not require explicit labeling of muscles or other underlying assumptions may prove effective when applying the collected signals toward volitional control of a powered prosthesis during gait. The high level of channel redundancy provided by the flexible electrode system pro- posed in this manuscript, previously unavailable for practical reasons, presents the prosthesis researcher with new directions of investigation. For example, future investigations may ana- lyze the performance of combinations of multiple signal processing techniques and controller models in movement-based tasks.2.3 sEMG Impulse Artifact Filter This work was supported by the MIT Media Lab Consortia. This work was presented in S. H. Yeon, and H. M. Herr

[0046] . This work presents a cumulative histogram filtering (CHF) algorithm to filter impul- sive artifacts within surface electromyograhy (sEMG) signal for time-domain signal feature extraction. The proposed CHF algorithm filters sEMG signals by extracting a continuous subset of amplitude-sorted values within a real-time window of measured samples using infor- mation about the probabilistic distribution of sEMG amplitude. For real-time deployment of the proposed CHF algorithm on an embedded computing platform, I also present an efficient, iterative implementation of the proposed algorithm. The proposed CHF algorithm was eval- uated on synthetic impulse artifacts superimposed upon undisturbed sEMG recorded from a subject with transtibial amputation. Results suggest that the CHF algorithm effectively suppresses the simulated impulse artifacts while preserving a minimum signal-to-noise ra- tio of 95% and an average Pearson correlation of 0.99 compared to the undisturbed sEMG recordings. 2.3.1 Motivation Due to the non-stationary and chaotic nature of sEMG signals, numerous attempts have been made to understand and interpret sEMG signals in meaningful, physiological terms related to collective muscle unit action potentials (MUAPs), muscle state, muscle force, muscle synergies, and even the kinematics and kinetics of muscle[3], [4],

[0047] –

[0050] . Among these studies, early seminal results demonstrated the optimality of root-mean-square (RMS) and mean-absolute-value (MAV) processors for estimating muscle force output

[0051] . As a result, RMS and MAV signal processors have been extensively emphasized and utilized due to their practicality and effectiveness, especially in direct model-based control of prosthetic systems

[0015] ,

[0018] ,

[0026] ,

[0052] . These same sEMG signal features have also been quantitativelyevaluated to be some of the most meaningful features when estimating neural activities

[0053] ,

[0054] . However, Windrich et al. highlighted the relatively greater number of studies focused on neural control of upper limb prostheses using sEMG signals compared to studies involv- ing lower-extremity (LE) prostheses

[0012] . Indeed, there exist unique considerations for LE prostheses, including limited options for sEMG electrode interfacing and significant dynamic pressure changes on the residual limb that can cause structural deformation and electrode contact variation. In order to overcome these challenges, many studies have presented solu- tions within various electrical and mechanical engineering domains

[0010] ,

[0024] ,

[0025] ,

[0042] . In the signal processing domain, one of the dominant challenges in LE prosthesis appli- cations is mitigating sEMG impulse artifacts produced during ambulatory activities. Specif- ically, impulsive acceleration at heel-strike can generate drastic pressure changes at the skin- electrode interface as well as vibration in the mechanical connections of the sEMG platform. While De Luca et al. set an effective standard for filtering generic motion artifacts and baseline noise with band-pass filtering, this standard does not consider and handle impulse artifacts

[0055] ,

[0056] . Interestingly, sEMG denoising algorithms that suppress impulse artifacts have been developed within the domain of ECG artifact filtering. These algorithms, in- cluding Wavelet analysis, independent component analysis, empirical mode decomposition, and adaptive thresholding, are both insightful and easily applied. Nonetheless, due to these algorithms having been developed for filtering ECG signals, they are potentially unsuitable as real-time impulse artifact filters since they require multiple sEMG channels, demanding high computational costs, and require thorough a-priori knowledge of artifact type

[0057] ,

[0058] . Therefore, in this work, I present an effective sEMG signal processing method and algo- rithm for filtering impulse artifacts and extracting MAV and RMS sEMG features. Based upon a probabilistic understanding of simultaneously recorded sEMG signals and impulse artifact noise, the proposed algorithm leverages real-time cumulative histogram filtering with minimal computational demands. Individual performances of the proposed filtering methodsare evaluated on both synthetic sEMG data and real sEMG data collected during use of a transtibial (TT) prosthestic device. 2.3.2 Preliminaries Probabilistic Model of sEMG Signals Understanding the nature of sEMG signals is critical to utilizing them as control inputs. Though sEMG signals are naturally non-stationary and stochastic, there have still been several attempts made to understand and model their behavior

[0051] ,

[0054] ,

[0059] ,

[0060] . In the course of finding an optimal estimation of instantaneous sEMG amplitude, E. A. Clancy and N. Hogan determined that the amplitude distribution of a band-passed sEMG signal could be described as between those of Gaussian and Laplacian PDFs. Fig. 2.16, recreated from

[0051] , demonstrates this relationship. Correspondingly, because the RMS of Gaussian PDFs is equivalent to their maximum likelihood amplitude with optimal signal-to-noise ratio (SNR), and because the MAV of Laplacian probability density functions (PDFs) is equivalent to their maximum likelihood amplitude with optimal SNR, the researchers determined that both RMS and MAV could be effectively utilized to estimate the amplitude of band-passed sEMG. Real-time sEMG Feature Extraction Real-time feature extraction from sEMG follows a typical procedure shared among many related studies[4],

[0055] . Fig. 2.17 demonstrates a simplified generic sEMG signal processing pipeline. First, raw sEMG signals are required to be sampled with a frequency of at least 1 kHz to effectively capture the entire signal’s power spectrum. Then, sampled raw sEMG signals are processed with a band-pass filter to remove motion artifacts and power-line noise

[0055] ,

[0056] . A moving short-time window with a duration between 100 and 500 milliseconds can then be applied to extract a signal feature at each time step. In cases where normalizationFigure 2.16: Normalized experimental probability density estimates of triceps muscle con- traction up to 50% maximum voluntary contraction, conducted as experiment II in

[0051] . Recreated from

[0051] with permission©2021 IEEE. of the extracted signal features is required, signal features are divided by the maximum amplitude of those measured at maximum voluntary contraction (MVC). The normalized MAV or RMS features are often utilized as estimates of %MVC to provide a muscle activation ratio (α) for real-time control applications

[0015] ,

[0026] ,

[0052] ,

[0061] ,

[0062] .Modeling Impulse Artifacts Theoretically, the ideal delta function has unity frequency response along the entire frequency domain, a characteristic that renders band-pass filtering ineffective. When processed by the sEMG processing pipeline shown in Fig. 2.17, a distinct impulse artifact is produced. As ademonstration, we modeled an impulse as a unipolar sinusoidal signal with unit amplitude and 10 ms duration, as shown in Fig. 2.18. We then processed the sEMG using the pipeline shown in Fig. 2.17. The selected band-pass filter was designed as an FIR filter with a pass- band between 20 and 340 Hz and a stop-band attenuation of 80 dB. A short-time window with 200 ms duration was used to extract MAV and RMS features.Figure 2.18: Modeled impulse signal, its response after band-pass filtering, and corresponding features extracted. The impulse signal is modeled with unit amplitude and 10 ms duration. Peak of impulse onset occurs at 0 s. The response shown is derived in a non-causal manner with zero phase and group delay. After band-pass filtering, an impulse artifact presented itself in the form of a sinc signal with duration of around 60 milliseconds, approximately three times longer than that of the original peak (when counting up to the third dominant peak). MAV and RMS feature extraction processes on the band-passed artifact signal generated flat output features with amplitude around only 1% of the raw signal amplitude. However, in a situation where the impulse artifact has an order-of-magnitude larger amplitude than the original sEMG signal, it can clearly be seen how the artifact can significantly alter the results compared to processing the sEMG signal alone.2.3.3 Filtering Method Cumulative Histogram Filtering Fig. 2.19 introduces the proposed CHF to the generic sEMG processing pipeline from Fig. 2.17.Figure 2.19: Modified sEMG signal processing pipeline with the proposed CHF filtering method. Within the process demonstrated by Fig. 2.19, band-passed sEMG signals are sampled using a given time window for every computational time-step in the same way as the generic filtering scheme. However, the proposed CHF processes sEMG data in between the short- time signal windowing and feature extraction steps. Algorithm 1 describes the CHF process in detail. Fig. 2.20 visualizes Algorithm 1 in sequential order.The CHF is conducted over all sampled sEMG signals within the time window. First, samples are sorted by amplitude to generate a cumulative histogram. Then, data sampleswithin a specified continuous range of the histogram are utilized for subsequent feature extraction. Given the example sEMG dataset shown in Fig. 2.20a, Fig. 2.20b shows its histogram and cumulative histogram and Fig. 2.20c shows the filtered subset of the short- time windowed sEMG samples to be used in feature extraction. The proposed CHF method is only applicable for certain types of time-domain sEMG signal features. The types of sEMG features compatible with CHF include MAV, RMS, variance (VAR), log-detector, temporal moment (TM), and v-Order processors

[0063] . These signal features share a common characteristic, namely that the extracted signal features are correlated to overall signal power and energy of sampled sEMG data without considering sEMG signal dynamics. Because the CHF rearranges the signal with no regard to sampling order, time-domain sEMG feature extraction methods that depend on signal dynamics, such as Bayesian filtering, are not compatible with the proposed method. Implementation of the CHF Basic Implementation While Algorithm 1 demonstrates the procedure of CHF in se- quence, Algorithm 2 summarizes an intuitive implementation of the proposed CHF algo- rithm, as direct implementation of Algorithm 1 in real-time would inefficiently utilize mem- ory space and processing power. l i h i i l i fComputation routines in Algorithm 1 (computing cumulative histogram, finding lowerFigure 2.20: Graphic demonstration of Algorithm 1 and Algorithm 2 with filter coefficients hLB of 10% and hUB of 90%. (a) Example sEMG time series highlighting the short-time windowed signal for feature extraction (b) Histogram plot and cumulative histogram plot extracted from Fig. 2.20a. The filter coefficients, hLB and hUB, are highlighted with corre- sponding xLB of 4.4% and xUB of 56%. The distribution of relative amplitudes is observed to be skewed right.(c) Sorted and unsorted sEMG data based on Fig. 2.20a highlighting original data and filtered data with hLB, hUB, xLB, xUB indicated.and upper bound coefficients, and conditional signal extractions) are implemented in Al- gorithm 2 by sorting absolute sEMG data by amplitude and indexing a sub-array. The computational complexity of Algorithm 2 is bounded by the computational complexity of the sorting algorithm. Thus, the minimum computational complexity of Algorithm 2 is O(n log n) with fast sorting algorithms such as quick sort that make no assumptions on the underlying data. Real-time Iterative Implementation The computational complexity of Algorithm 2 can be significantly reduced and optimized by implementing CHF in an iterative scheme. Algorithm 3 provides an iterative CHF implementation.By inserting only newly sampled data into the running data buffers XABSand S, the sorting operation from Algorithm 2 is no longer required. Algorithm 3 can be considered to be a linked-list style data buffer update scheme. With Algorithm 3, computational complexity can be reduced to O(n) due to only requir- ing linear array traversals and O(1) operations. This scheme can easily be implemented in real-time on microcontrollers or microprocessors used with portable mechatronics systems.2.3.4 Evaluation Subject Recruitment and Data Collection A single subject (43 years old, female) with unilateral transtibial amputation participated in data collection for the evaluation of the proposed CHF method. The subject, a recipient of the novel agonist-antagonist myoneural interface(AMI) amputation surgery, possessed two dynamic AMI pairs (Tibialis Anterior (TA) coupled with Lateral Gastrocnemius (LG) and Tibialis Posterior (TP) coupled with Peroneous Longus (PL)) within her residuum

[0018] . The subject provided written informed consent through MIT COUHES protocol #1812634918. In the first experimental session, the subject was asked to voluntarily rotate her phantom limb while seated. Raw sEMG signals from TA and TP muscles were collected while the subject wore her prescribed prosthetic socket and a liner. A custom embedded sEMG acqui- sition system was utilized along with a flexible electrode designed for within-socket sEMG acquisition

[0010] ,

[0011] ,

[0044] . Raw sEMG signals were processed with a digital band-pass FIR filter with pass-band of 40 to 340 Hz and stop-band attenuation of 80 dB. A short-time window with duration of 200 ms was used for subsequent CHF processing. Fig. 2.21 shows the collected reference data processed using the generic technique illustrated in Fig. 2.17 to yield RMS and MAV features. In a second experimental session, sEMG data were acquired during ambulatory activity to evaluate the efficacy of the proposed CHF scheme illustrated in Fig. 2.19. Performance Analysis with Clean sEMG Signals This subsubsection utilizes the sEMG data from the first experimental session to quantita- tively compare the output of CHF to generic sEMG processing in the absence of impulse artifacts. In general, higher similarity between results is more desirable, implying that the proposed CHF does not suppress or distort original information from the sEMG signals. Using data shown in Fig. 2.21, we applied the proposed CHF while sweeping filterFigure 2.21: Reference sEMG data collected from TA and TP muscles during voluntary rotational movement of the residual phantom limb. The plot visualizes both RMS and MAV processing outputs. coefficient hLB from 0% to 40% and hUB from 60% to 100% with a step size of 2.5%, resulting in 289 sets of cutoff boundaries per channel and processor. We processed sEMG data from the subject’s TA and TP muscle channels with MAV and RMS processors, resulting in a total of 1156 processed sEMG time series. Pearson’s Correlation We analyzed Pearson’s correlation coefficients to quantify the correlation of the processed sEMG data. Table 2.4 shows the result of the correlation analysis. All of the correlation coefficients from both TA and TP muscles with RMS and MAV processors present mean values greater than 0.99, implying that CHF processed data are highly correlated with unfiltered sEMG data. Because the data from the TP muscle with RMS processor show the lowest minimum correlation values of 0.9837, individual correlation coefficients of these data are partially presented in Table 2.5. With a minimum observed Pearson’s correlation of 0.9837, all of the CHF processed dataTable 2.4: Generated Impulse Artifact Correlationsshow high correlation to their original signals. Fig. 2.22 presents a qualitative comparison of the CHF-processed RMS data from TA and TP muscles.Figure 2.22: Partial results from all 289 filtered sEMG series using swept boundary coeffi- cients hUB and hLB, processed from the original reference sEMG of Fig. 2.21. Partial data from 9 to 16 seconds shown. Signal-to-noise Ratio Signal-to-ambient-Noise-Ratio (SNR) of the CHF processed data are also compared with the SNR of the original data. Ambient noise is defined as the meanvalue of sEMG data collected and processed during the initial 5 seconds of the trial when the subject was instructed to rest, as shown in Fig. 2.21. SNR is subsequently defined as the ratio between the maximum measured value and the ambient noise. Fig. 2.23 shows the result of the SNR analysis.Figure 2.23: Effect of filtering on SNR of reference sEMG data from Fig. 2.21 with swept filter coefficients. (a) Relative SNR of the filtered TA sEMG data. Filtering preserves at least 95% relative SNR for both MAV and RMS features. (b) Relative SNR of the filtered TP sEMG data. As opposed to the TA muscle, filtering of the TP muscle data results in increased SNR. While the correlational relationship between CHF and SNR is shown to be non-linear, relative SNRs of the data remain higher than 95% compared to the SNRs of unfiltered data in all CHF processed data. Importantly, SNRs of the CHF-processed sEMG data from the TP muscle increase with decreasing hUB. Artifact Suppression Performance Analysis Synthetic Impulse Artifact Model A synthetic sEMG impulse was injected into the data collected from the first session to quantitatively evaluate the CHF’s ability to suppress artifacts. The impulse model from Section 2.3.2 was utilized to synthesize the artifact- affected sEMG dataset. A unipolar impulse with duration of 20 milliseconds and amplitude of 5 mV, approximately 25 times larger than the maximum of band-passed sEMG data was used. A periodic impulse train with impulses every 1 second was generated and added tothe raw sEMG data. Fig. 2.24 shows the synthesized reference data and its simulated signal output with RMS and MAV processors.Figure 2.24: Synthesized artifact-affected sEMG data. (a) Artifact-affected and band-passed TA muscle sEMG. (b) Effect of impulse artifacts on output RMS and MAV features. Partial data from 9 to 16 seconds shown Artifact Suppression Ratio Artifact suppression ratio is defined as the average ratio of outputs between filtered and non-filtered data sampled at the center of each artifact impulse. The synthesized sEMG data with impulse artifacts were processed using the CHF filter andswept boundary coefficients in a manner identical to the process of subsubsection IV-B. Fig. 2.25 shows the result of the artifact suppression ratio calculation.Figure 2.25: Artifact suppression ratio analyses with swept hLB and hUB filter coefficients. Artifact suppression ratio is largely correlated with hUB. (a) Artifact suppression ratio of the TA (b) Artifact suppression ratio of the TP From the data, artifact suppression ratio is inversely correlated with hUB. Fig. 2.26 visualizes the effect of varying hUB on noise suppression.Figure 2.26: Artifact suppression ratio of CHF. Synthesized sEMG data shown in Fig. 2.24 are processed with hLB of 10% and hUB of 60%, 70%, 80%, and 90%. RMS features show smaller impact artifacts with decreasing hUB.Case Study: Ambulation In this paragraph, I present the results from the second ex- perimental session involving ambulatory activity. Fig. 2.27 shows the sEMG data collected from the TA and TP muscles of the subject’s residuum within her socket during ambulation. The data are processed both with and without the proposed CHF method. The CHF was tuned with hUB at 70% and hLB at 10%. 70% was determined as the highest cutoff for hUB before relatively diminishing returns on artifact suppression percentage, as seen in Fig. 4.31a and Fig. 4.31b. Similarly, 10% was determined as the lowest cutoff for hLB, though the absolute differences produced by varying hLB were significantly lower.Figure 2.27: sEMG data collected during ambulation were processed using CHF with hUB at 70% and hLB at 10%. Impulse artifacts are synchronized with ground reaction force peaks, but suppressed in the resulting RMS feature. In Fig. 2.27, impulse artifacts synchronized with ground reaction force (GRF) peaks are clearly visible in the sEMG data collected from both TA and TP muscles. The results from non-CHF RMS and MAV processors demonstrate staircase-like plateaus in their output waveforms. With CHF, staircase-like plateaus are removed, yielding output waveforms whichmore qualitatively resemble natural activation patterns. Though there are no ground truth reference labels for sEMG data and impulse artifact sources within the given dataset, we feel the CHF’s results are compelling. 2.3.5 Discussion Through this work, I developed a filtering technique capable of robustly extracting time- domain features from sEMG signals in the presence of undesired impulse artifacts. The proposed CHF method was shown to effectively filter out synthetic impulse artifacts in a reference sEMG data set while preserving the signals’ underlying information. Intuitively, the CHF method exploits the different PDFs of the sEMG signal and noise models as described in Section 2.3.2. Because impulse artifacts under the tested condi- tions demonstrated a right-skewed PDF compared to the underlying sEMG’s PDF, the CHF method effectively filtered the artifacts out by discriminating via signal amplitude. The CHF method is able to leverage a limited range of information from the middle portion of a time window’s cumulative distribution, and as seen in Fig. 2.27, this omission of high-amplitude samples enabled extraction of sEMG features with comparable levels of SNR compared to an artifact-free reference data set. At the extreme limit, the CHF method becomes identical to the median feature extraction method with hLB and hUB at 50%. Admittedly, the optimality of the median processor as a muscle force estimator has not been compared to MAV or RMS processors. However, due to the stochastic nature of sEMG, features that utilize multiple samples when extracting information from an sEMG signal may vary less compared to the one produced by the median processor which depends on a single median sample. As CHF-filtered MAV and RMS features with reasonable hLB and hUB ranges demonstrated high correlations to non- filtered MAV and RMS features, it is suggested that the proposed CHF method does not disrupt the optimal nature of the MAV and RMS features discussed in Section 2.3.2. Though this preliminary study yielded promising results, the proposed filtering methodremains to be tested thoroughly with more data collected under additional testing conditions.2.4 Synchronized sEMG and Ultrasonography Measure- ment This work was supported by the MIT Media Lab Consortia and Eunice Kennedy Shriver National Institute of Child Health & Human Development of the National Institutes of Health (NIH) under award number R01HD097135. This work was presented in Yeon et al

[0064] . 2.4.1 Motivation Surface electromyography (sEMG) and ultrasonography provide easy-to-use, low-cost, and noninvasive modalities to assess muscle activity and fascicle dynamics, and have been widely used in both clinical and lab settings. However, due to size of these sensors and limited skin surface area, it is extremely challenging to collect data from a muscle of interest in a spatially as well as temporally synchronized manner. Here, we introduce a low-cost, non- invasive flexible electrode that provides high quality sEMG recording, while also enabling spatiotemporally synchronized ultrasonography recordings. The proposed method was veri- fied by comparing ultrasonography of a phantom and a tibialis anterior (TA) muscle during dorsiflexion and plantarflexion with and without the electrode acutely placed under an ul- trasound probe. The results show no significant artifact in ultrasonography from both the phantom and TA fascicle strains due to the presence of the electrode, demonstrating the capability of spatiotemporally synchronized sEMG and ultrasonography recording. 2.4.2 Method Flexible electrode for spatiotemporally synchronized sEMG and ultrasonography recording In this work, a custom-made flexible and low-profile sEMG electrode, the presented SLIP electrode in Section 2.2, was utilized to explore a novel sensor modality for spatiotemporallysynchonized sEMG and ultrasonogrpahy measurement. As stated, this recently developed electrode consists of a Polyimide film substrate with gold plated skin-contact surfaces with a thickness of 80 to 100 micron. The electrode was designed to ensure light-weight, easy-to- use, and comfort of subjects while enabling large scale manufacturing with generic flexible printed-circuit-board (fPCB) fabrication techniques at a low-cost

[0010] ,

[0011] . Evaluation of artifacts using ultrasonography of a phantom from a flexible elec- trode I first evaluated the quality of an ultrasound image of a phantom acquired from an ultrasonic probe spatially overlapped with the flexible sEMG electrode. A copper spring (length of 78.9 mm, outer diameter of 35.6 mm, inner diameter of 20.1 mm, and a pitch of 11.56 mm) was utilized as the phantom for the evaluation as shown in Fig. 2.28. The ultrasonograpy of the phantom was conducted by placing the phantom in the middle of a flat square plastic container filled with water while locating the ultrasound probe on the plastic container wall with ultrasound gel. Two-types of ultrasonography (B-mode images) of the phantom, with (w EMG) and without (w / o EMG) the overlapped sEMG electrode, were collected to evaluate artifact in ultrasonography due to the presence of the electrode. Quality of the collected images were then compared to analyze an ultrasonic artifact caused by the overlapped sEMG electrode with the probe. Evaluation on Spatiotemporally synchronized sEMG and ultrasonography record- ing of human muscle I further validated the proposed method through recording spatiotemporally synchronized sEMG and ultrasonography of the tibialis anterior (TA) muscle during full ankle dorsiflexion and plantarflexion. A subject with biologically intact limbs was asked to fully dorsiflex and plantarflex the ankle while joint kinematics were recorded through a goniometer attached to medial aspect of the ankle joint as shown in Fig. 2.29. The ultrasound probe was acutely(a) (b) (c) (d) Figure 2.28: Experimental setup for the validation of the proposed method through the phantom. (2.28a) The phantom used for the validation. (2.28b, 2.28c) Measurements of the phantom with and without the proposed electrode. (2.28d) The spatially synchronized electrode with the ultrasound probe. placed on the muscle belly the TA to assess muscle fascicle strains. Two types of trials were conducted; the ultrasonography of the TA was recorded with and without the flexible electrode placed under the probe (Fig. 2.29) and muscle fascicle strains were computed from the recorded ultrasound videos

[0065] . The subject was guided verbally to perform steady and consistent dorsiflexion and plantarflexion. Further, the placement of ultrasound probe was marked to ensure consistent ultrasound probe placements between trials. All experiments were carried out with informed consent at the Massachusetts Institute of Technology (MIT), under the approval of the Committee on the Use of Humans as Experimental Subjects (COUHES).Figure 2.29: Experimental setup for validation of the proposed method through assessments of TA muscle dynamic movements. (2.29a) Sensor placements. The goniometer and ultra- sound probe were placed on the medial aspect of the ankle joint and muscle belly of the TA, respectively. (2.29b) The flexible electrode was placed under the ultrasound probe, enabling spatiotemporally synchronized sEMG and ultrasonography recording. 2.4.3 Result Artifact in ultrasonography of the phantom from the flexible electrode The ultrasonography of the phantom with and without the flexible electrode placed on the ultrasound probe is shown in Fig. 2.30. To assess the quality of ultrasonography, the average ultrasonographic images with 600 frames of the phantom with and without the flexible electrode were computed. Based on the our qualitative assessments of the both images, we found no significant artifact in the image due to the presence of the electrode. Further, the ultrasonography was able to capture the features of the phantom, in which the outer diameter of spring and distance between rings were respectively 35.6 mm and 11.53 mm, while having the electrode under the ultrasoundFigure 2.30: Qualitative evaluation in ultrasonography of a phantom. (2.30a) Result with the flexible electrode. (2.30b) Result with the flexible electrode. probe. Spatiotemporally synchronized sEMG and ultrasonography recording of human muscle Experimental results of spatiotemporally synchronized sEMG and ultrasonography recording of the TA muscle are shown in Fig. 2.31. Based on the qualitative evaluation of ultrasonographic images with and without spatially synchronized electrodes, the results showed no significant artifact induced by the presence of the proposed electrode and only marginal ’shade-like’ artifact was introduced to the images. Estimated TA muscle fascicle strains from ultrasound videos were shown in Fig. 2.32 along with the temporally synchronized joint kinematics and sEMG. During steady full range of ankle dorsiflexion and plantarflexion, the recorded TA muscle fascicle strains with simultaneous sEMG recording showed comparable and consistent val- ues to the fascicle data recorded without simultaneous sEMG recording. To quantify theFigure 2.31: Artifact in ultrasonography of the TA muscle due to the flexible electrode placed under the ultrasound probe. Only negligible artifact is introduced in the ultrasonography as a form of shade by the presence of the flexible electrode. (2.31a) Ultrasonogrpahy of the TA muscle without the flexible electrode. (2.31b) Ultrasonography of the TA muscle with the flexible electrode. capability of simultaneous sEMG and ultrasonography, the relationships between ankle joint kinematics and TA muscle strains with (w EMG) and without (w / o EMG) simultaneous sEMG recording have been investigated (Fig. 2.32c). The similarity of the two identified relationships was evaluated by a Wilcoxon signed-rank test at a significance level of α = 0.05 across the range of motion (ROM) as shown in Fig. 2.32c. Our results indicated no significant differences between the ROM in the muscle fascicle relationship to joint kinematics identified from ultrasonography with stable simultaneous sEMG recording (Fig. 2.32b) and the one without sEMG recording. 2.4.4 Discussion Throughout this work, we have verified a novel sensor modality that enables spatiotempo- rally synchronized electromyography and ultrasonography recording. This is made possibleFigure 2.32: Experimental results of spatiotemporally synchronized electromyography and ultrasonography recording of the TA muscle. (2.32a) and (2.32b) are representative record- ings of ankle joint angle and normalized fascicle length of the TA msucle during full dorsi- flexion and plantarflexion with and without EMG recording from the muscle, respectively. (2.32c) Joint angle and normalized fascicle length relationships investigated without (w / o EMG) and with (w EMG) simultaneous EMG recordings during full dorsiflexion and plan- tarflexion (n=6 cycles). No significant differences were found between the two relationship curves. P values are reported. Wilcoxon signed-rank test at a significance level of α =0.05. through a flexible sEMG electrode that has only a 80-to-100 microns thickness which min- imizes the interference of mechanical impedance between an ultrasound probe and the skin surface, having the electrode invisible and only leaving marginal ’shade-like’ artifact in the images. For muscles with a slim-shape, small size, or that are squeezed between other muscles, it is extremely challenging to obtain temporally synchronized sEMG and ultrasonography recordings due to limited surface area allowed for the muscle of interest. Generally, plac- ing the sEMG electrode aside to ultrasound probe, or non-spatially synchronized electrode placement, is not preferred because of high cross-talk in sEMG recording from surrounding muscles, making it difficult to estimate the muscle activity level of interest. With con- ventional surface electrodes for sEMG recording such as a Ag / Cl electrode, it is arguably impossible to collect ultrasonography data while placing the electrode under the probe due to the volume of the electrode and its significant artifact in the images. Nevertheless, bysimply placing the flexible electrode under the ultrasound probe along with the conventional ultrasound gel, the proposed sensor modality was able to obtain temporally synchronized sEMG and ultrasonography with minimal surface area, offering the capability to study a wide rangse of musculature in a noninvasive manner. In developing a musculoskeletal model towards implementing biomimetic assistive tech- nologies and establishing effective rehabilitation, identifying muscle fascicle relationship against the joint kinematics is critical. The experiment results show that the proposed sensor modality provides non-inferior ultrasound imaging of muscle fascicle dynamics while also providing stable sEMG recording. This would allow the design of sophisticated musculoskeletal models with detailed descrip- tions of muscles based on the directly measured, temporally synchronized muscle activity and muscle fascicle data. Moreover, this would offer the opportunity to further verify and develop our understand- ing of afferent-efferent signaling and reflex, which is known to be mediated by both the fascicle dynamics and activity level of muscle, in an in vivo setting. Meanwhile, utilizing spatial sEMG data collected from high density grid (HD-sEMG) electrodes to estimate multi-modal physiological information has emerged as an active area of investigation

[0066] –

[0069] . By simultaneously considering the additional data from ultra- sound imaging along with HD-sEMG data collected from flexible HD-sEMG electrodes fab- ricated with our proposed approach, we envision the potential to capture high-resolution, synchronous mechanical and electrical muscle dynamics that provide new insights into muscle neurophysiology. The proposed sensor modality poses as a easy-to-use, low-cost, and noninvasive mea- surement of spatiotemporally synchronized sEMG and ultrasonography. Future works will include the development of a musculoskeletal model and investigation of neural circuitry behind complex motor control based on the spatiotemporally synchronized sEMG and ultra- sonography measurements. We hope that this work will enable more research on accuratemodeling of human biomechanics through the proposed multi-modal measurement method. 2.5 Summary In this chapter, I presented several novel methods and findings aimed at enhancing the utility and efficacy of surface electromyography (sEMG) for bionic applications. These ad- vancements address important challenges associated with the implementation of sEMG in real-world scenarios, including user comfort, signal quality, artifact mitigation, and the need for a more comprehensive understanding of muscle function. First, I introduced a flexible and low-profile SLIP electrode design that improves user comfort and signal quality, making it suitable for long-term use within prosthetic sockets. This innovative design overcomes the limitations of traditional sEMG electrodes, which can cause discomfort and skin irritation, leading to reduced user compliance and signal reliability. Next, I proposed a real-time signal processing algorithm to filter out impulse artifacts that often contaminate sEMG signals during dynamic movements. The developed cumulative histogram filtering (CHF) technique effectively suppresses impulse artifacts while preserving the underlying sEMG signal, enabling more accurate and reliable control of bionic devices during ambulatory activities. Finally, I demonstrated a synchronous sEMG and ultrasound acquisition method that enables the simultaneous assessment of muscle electrical activity and mechanical deforma- tion. This multi-modal approach provides valuable insights into muscle function and control, facilitating the development of more advanced and intuitive control strategies for bionic de- vices. The advancements presented in this chapter have significant implications for the devel- opment of next-generation bionic devices. By improving user comfort, signal quality, and artifact mitigation, the proposed methods can enhance the long-term usability and reliabil- ity of sEMG-based control systems. Moreover, the integration of sEMG with other sensingmodalities, such as ultrasound, opens up new possibilities for the design of more sophis- ticated control algorithms that can adapt to the user’s intent and dynamically changing environments. Indeed, these collective efforts in advancing sEMG modalities for bionics have played a critical role in enabling state-of-the-art research in bionic systems

[0023] ,

[0070] –

[0072] . Specifically, with the presented sEMG modality, we demonstrated the continuous neural control of a bionic ankle prosthesis without any intrinsic state machine, restoring natural biomimetic gait for the first time in the world

[0071] . I believe this is just the beginning and it will unlock new levels of versatility in modern bionics.Chapter 3 In-vivo and Untethered Validation of Magnetomicrometry This work was condcuted in collaboration with Professor Thomas Roberts at Brown Uni- versity. This work was supported by the Salah Foundation and Media Lab Consortia. This work was presented in Taylor., and Yeon., et al

[0073] . Magnetomicrometry is an emerging myo-neural sensing modality that has the potential to complement surface electromyography (sEMG) in providing a more comprehensive under- standing of muscle function and control for bionic applications. Previous work has demon- strated the potential of Magnetomicrometry in tracking muscle states using tightly-controlled in situ setups [6]. However, to establish Magnetomicrometry as a viable myo-neural inter- face for real-world applications, it is crucial to validate its performance in untethered, freely moving subjects. The main objective of this chapter is to demonstrate the capabilities of Magnetomicrometry in tracking muscle tissue length in real-time during various motor ac- tivities in freely moving animals. By validating its performance in untethered subjects, we aim to establish Magnetomicrometry as a reliable and effective myo-neural sensing modality for bionic applications.3.1 Motivation Muscle length measurements have driven important discoveries in movement biomechanics, informed models of motor control, and provided strategies for prosthetic and robotic de- sign.

[0074] –

[0077] . For decades, sonomicrometry (SM) has informed how muscles move, provid- ing high accuracy (70 m resolution) and high bandwidth (>250 Hz)

[0078] . Fluoromicrometry (FM) expanded the muscle tracking toolkit, enabling high accuracy (90 µm precision) and high bandwidth (>250 Hz) for high-marker-count tracking

[0079] ,

[0080] . Further, image-based ultrasound (U / S) added the capability to non-invasively track muscle geometries

[0081] –

[0083] . Yet, collecting direct muscle length measurements in natural environments remains in- feasible, and thus indirect muscle length estimation is still used for observing natural move- ments. For instance, muscle lengths are estimated using joint angles via biophysical mod- els

[0084] . These approximations are used due to the limitations of current muscle motion sensing techniques, all of which are tethered or bulky. SM and U / S both require tethered connections to bulky hardware for sensing, with SM requiring advanced surgery and per- cutaneous wires. And while FM does not require a tethered connection, it is limited to a volume approximately the size of a soccer ball, requires equipment the size of a small room, and is time-constrained due to thermal limitations and subject radiation exposure

[0079] ,

[0083] ,

[0085] . Present muscle length tracking technologies also require substantial post-processing time, hindering their use in longitudinal studies. SM requires accounting for and filtering out artifacts such as triggering errors, FM requires point labeling in stereo images, and U / S requires fascicle labeling, all of which require at least some manual processing

[0079] ,

[0086] ,

[0087] . While machine learning techniques have shown potential for automatic fascicle length tracking from ultrasound images, the current lack of reliability in tracking cross-activity measurements prevents such a strategy from being applicable toward sensing fascicle lengths during natural movement

[0088] .Researchers need a sensing platform that can operate untethered in natural environments, sensing the full dynamic range of muscle movement in context. To address this need, we developed magnetomicrometry (MM), a minimally-invasive strategy for portable, real-time muscle tracking. MM uses an array of magnetic field sensors to locate and calculate the dis- tance between two implanted magnetic beads with sub-millisecond time delay. This distance provides a measurement of the muscle tissue length between the implanted beads. MM al- lows continuous recording over an indefinite collection interval extending across hours, with the potential for continuous use across days, weeks, or years. Previously introduced magnetomicrometry, a method that uses magnetic beads to wire- lessly monitor muscle tissue length changes, demonstrated its capability and efficacy via tightly-controlled in situ testing[6]. In this work, we validate the accuracy of magnetomi- crometry against fluoromicrometry during untethered running in an in vivo turkey model. We demonstrate real-time muscle tissue length tracking of the freely-moving turkeys execut- ing various motor activities, including ramp ascent and descent, vertical ascent and descent, and free roaming movement. Given the demonstrated capacity of magnetomicrometry to track muscle movement in untethered animals, we feel that this technique will enable new scientific explorations and an improved understanding of muscle function. 3.2 Method All animal experiments were approved by the Institutional Animal Care and Use Committees at Brown University and the Massachusetts Institute of Technology. Wild turkeys (Meleagris gallopavo, adult female) were obtained from local breeders and maintained in the Animal Care Facility at Brown University on an ad libitum water and poultry feed diet. We used three animals in this study.3.2.1 Surgical Procedure One pair of 3-mm-diameter Parylene-coated magnetic beads (N48SH) were implanted into the right lateral gastrocnemius muscle of each turkey, with a target magnetic bead separation distance of 3.5 cm. For details on the surgical procedure and implants, see Taylor., et al. 2022

[0089] . A one-month recovery period was given before the start of the data collection. 3.2.2 Magnetomicrometry For this study, we designed a custom magnetic field sensing array as shown in Fig. 3.1. The sensing array was equipped with 96 of 3-axis magnetometers (LIS3MDL, STMicroelectronics) spaced 5.08 mm apart in an 8-by-12 grid. Each sensor was supplied with nonmagnetic capacitors (VJ1206Y105KCXAT and VJ0603Y104KCXAT, Vishay). Seven digital multiplexers on the sensing array allowed time-domain multiplexing (one 74HC138BQ,115 multiplexing into six 74HC154BQ,118, Nexperia) via a wired connection. The sensing array was connected through a custom adapter board to an off-the-shelf wireless microcontroller embedded system (Feather M0 WiFi microcontroller, Adafruit), which was powered by a lithium-ion polymer battery (3.7 V, 1800 mAh, 29 g). The microcontroller sampled the magnetic field signals at 155 Hz and wirelessly transmitted them to the magnet tracking computer via a WiFi router (Nighthawk R6900P, Netgear). The tracking algorithm ran in real-time on the magnet tracking computer, a Dell Precision 5550 laptop (Ubuntu 20.04 operating system) with 64 GB of random-access memory and an Intel i7 8-Core Processor, running at 2.30 GHz. The tracking algorithm used, including the strategy for disturbance compensation, is fully-detailed in the previous work

[0090] . We affixed the sensing array to the limb using Opsite Flexifix adhesive film (Smith Nephew). To secure the array, we first applied a base layer of the adhesive film to the skin. We then positioned the sensing array over the leg and wrapped the adhesive film around the sensing array and the leg. To maintain a sufficient minimum distance betweenFigure 3.1: Magnetomicrometry Embedded System. We fabricated a custom sensor board (left) and a custom control board (right) for this study. The sensor board holds the magne- tomicrometry sensing array, consisting of 96 magnetic field sensors arranged with a spacing of 5.08 mm. Digital multiplexers on the sensor board allow time-domain multiplexing, en- abling a single microcontroller on the control board to communicate with and control all magnetic field sensors on the sensor board. The control board merges the data from the sensor board and streams the data wirelessly to the magnet tracking computer. The sensor board and control board weigh 24 g and 12 g, respectively. the magnetic beads and the sensing array, we positioned layers of foam between the base adhesive and the array. We then secured the control board and battery within the back feathers of the turkey. 3.2.3 Accuracy Validation of Magnetomicrometry Against Fluo- romicrometry We used the W. M. Keck Foundation XROMM Facility at Brown University to perform FM

[0079] . We collected X-ray video from two intersecting X-ray beams oriented at 51 degrees relative to one another. We mounted a treadmill with a wooden base (TM145, Horizon Fitness) between the X-ray cameras and the X-ray sources and built a housing over thetreadmill with a movable wall to position the birds within the capture window. The turkeys walked and ran at five speeds (1.5, 2.0, 2.5, 3.0, and 3.5 m / s) in a randomized order until at least ten gait cycles were visible within the FM capture volume for each speed. We collected FM data at 155 Hz. Time syncing was performed via a coaxial cable connection from FM to an off-the-shelf microcontroller development board (Teensy 4.1, Adafruit). The time-syncing microcontroller relayed the time sync signal to the magnet tracking computer via a custom adapter board. 3.2.4 Untethered Muscle Tracking Across Various Activities We constructed a hallway to guide the turkeys through the variable terrain activities. We stacked plyometric boxes (Yes4All) in the hallway to heights of 20 cm, 41 cm, and 61 cm for vertical ascent and descent. Upon completion of the vertical ascent and descent tests, we placed ramps (Happy Ride Folding Dog Ramp, PetSafe) up to and down from the plyometric boxes at inclines of 10° and 18° for the turkeys to ascend and descend. Separate from the variable terrain activities, we then allowed one turkey (Bird A) to roam freely within its enclosure while we continued to record MM. 3.2.5 Benchtop Magnetomicrometry Validation for the Variable Ter- rain Activity For benchtop testing, we used super glue (Krazy Glue) to affix each of two N48SH magnetic beads into two 1x1 round LEGO plates. We attached these round LEGO plates to a 1x6 LEGO technic block, one at each end, to separate the pair of magnetic beads by a fixed distance of 40 mm as shown in Fig. 3.2 A. We imaged this pair of beads using FM to validate the 40 mm fixed distance between them. Specifically, after the MM accuracy validation of two of the birds, we collected FM data with the magnet pair statically in the volume and used the average of the last three seconds from each of these two FM collections to confirmthe distance between the two beads.Figure 3.2: Benchtop Magnetomicrometry Validation. (A) Two magnets were placed 40 mm apart in a 1x6 LEGO Technic block and centered under the sensing array. We used the tracked magnet z-position as a guide in setting up the minimum depth measurements, and we enforced the remaining depths by adding / removing 3.2-mm-thick 1x6 plates under the Technic block. (B) We manually swept the magnets out and back to center along the x and z axes. We set the sweep trajectory to sweep just beyond the volume within which the magnets were tracked during the variable terrain activities (see Supplementary Figure 6). For reference, the centroid of the two beads is labeled at the origin and at the extent of the sweep trajectories. (C) The vertical axis represents the magnetomicrometry error, and the horizontal axes represent the centroid x and y position for the two sweep trajectories at each depth. The submillimeter error range is marked with a gray background. A maximum error of 1.463 mm was found at a test location just beyond the tracked bead position extent (bottom right plot). To determine the volume to sweep this FM-validated 40-mm-distanced bead pair during benchtop testing, we analyzed the magnetic bead tracking data from all ramps and verti- cal elevation changes to determine the full range of the three-dimensional magnetic bead positions relative to the MM sensing array across all three turkeys. We first aligned and centered the magnet pair under the array, and we used the tracked magnet z-positions to place the magnets at the closest depth that was possible while stillwithin the full-scale sensing range of the sensors ( 1 cm). We then used 3.2-mm-thick 1x6 LEGO plates to enforce the remaining depths as shown in Fig. 3.2. At each depth, we manually swept from center out and back along the x and y axes to the point where the farthest magnet reached just beyond the test volume requirements derived from the variable terrain activity as shown in Fig. 3.2. 3.2.6 Data Analysis We post-processed the FM data using the XMA Lab software

[0091] . All FM and MM data were left unfiltered. We aligned the MM and FM data using the time sync signal and linearly interpolated the FM data at the MM measurement time points. Then, due to imprecision of the time sync signal from the X-ray system, we used local optimization to further align the MM and FM signals while iteratively interpolating FM. During one trial (one data collection for one turkey at one speed), where the tracking computer did not receive a time sync, we used global optimization to align the MM and FM signals. To validate the use of global optimization for synchronization, we tested this same global optimization on all other trials and found that the global optimization successfully located all time sync signals. We estimated the noise from manual FM processing by independently processing one set of ten gait cycles three times (manually re-processing the video data twice without reference to the previously processed data). We then calculated the variance at each time point and used the square root of the average variance as our estimate of the FM manual processing noise. We calculated the adjusted MM noise by subtracting the average variance of the FM manual processing noise from the variance of the difference between the MM and FM signals for each bird, then taking the square root. For gait cycles where only one toe strike was visible in the video, we normalized the gait cycle using the timing of the peak MM signals in the previous and current gait cycles.3.3 Result 3.3.1 Accuracy Validation of Magnetomicrometry Against Fluo- romicrometry To verify MM tracking accuracy during untethered activity, we tracked implanted magnetic bead pairs in turkey gastrocnemius muscles (right leg, three turkeys) using both MM and FM while the turkeys walked and ran at multiple speeds on a treadmill. The experiment setup overview and example data trajectory are demonstrated in Fig. 3.3.Figure 3.3: Validation of Untethered Muscle Tracking using Magnetomicrometry. (A) A magnetic field sensing array on the surface of the leg tracks the positions of two magnetic beads implanted into the muscle. A feather microcontroller (µC) in the turkey feathers wirelessly transmits the magnetic field data to a magnet tracking computer that calculates and displays the magnetomicrometry (MM) signal in real time. The turkeys walked and ran on a treadmill while x-ray video cameras recorded synchronized fluoromicrometry (FM) data for post-processing. (B) Comparison of MM (blue) with FM (red) to validate the MM accuracy. These representative results during running gait show the submillimeter accuracy of MM during untethered muscle length tracking.We compared the distances between the magnetic bead positions as measured by MM with their distances as measured by FM to evaluate accuracy during the treadmill activity as shown in Fig. 3.4. The coefficients of determination (R2values) between MM and FM were 0.952, 0.860, and 0.967 for Birds A, B, and C, respectively as shown in Fig. 3.5. The differences between MM and FM were -0.099 ± 0.186 mm, -0.526 ± 0.298 mm, and -0.546 ± 0.184 mm for Birds A, B, and C, respectively as shown in Fig. 3.6Figure 3.4: Untethered Muscle Tracking During Treadmill Running: Magnetomicrometry Versus Fluoromicrometry. Changes in muscle tissue length measured by MM (blue) and FM (greenred) for three turkeys at five speeds (30 s shown for each speed). The column to the right of the plots gives the coefficients of determination (R2) between magnetomicrometry and fluoromicrometry corresponding to each turkey and speed. Gaps in the fluoromicrometry data are due to researcher selection of full gait cycles during which both magnetic beads were visible in both x-ray images. Gaps in the magnetomicrometry data (gray) are due to packet drops during wireless transmission of the magnetic field signals to the tracking computer (gaps below 50 ms interpolated in gray, gaps above 50 ms highlighted in gray). The turkey gait diagram below the plots shows the corresponding gait phases over one gait cycle. To determine the study-specific reliability of the manual FM processing (marker position labeling in the X-ray video data), ten gait cycles of raw FM data were independently manually relabeled three times for one bird at one speed. Across these three labels for these ten gaitFigure 3.5: Coefficients of Determination (R2values) between Magnetomicrometry and Flu- oromicrometry. We compared all magnetomicrometry measurements (horizontal axis) across 50 gait cycles of turkey running (10 gait cycles at each of 5 speeds, for each bird) against time-synchronized, interpolated fluoromicrometry measurements (vertical axis). Data are plotted in blue, orange, and purple for Birds A, B, and C, respectively. Coefficients of de- termination (R2values, shown in corresponding colors) for each bird were 0.952, 0.860, and 0.967, respectively.. cycles, manual FM processing was consistent to a standard deviation of 0.098 mm as shown in Fig. 3.7. MM’s 99-th percentile tracking time delays were 0.698 ms, 0.690 ms, and 0.664 ms for Birds A, B, and C, respectively as shown in Fig. 3.8, and the MM data did not require any post-processing. In contrast, post-processing the FM data into marker-to-marker distances required approximately 84 manual processing hours spread across multiple months. 3.3.2 Untethered Muscle Tracking Across Various Activities To investigate the feasibility of using MM during dynamic, natural motion, we constructed a series of obstacles for the turkeys to navigate. Specifically, we provided the turkeys withFigure 3.6: Kernel Density Estimates of Differences between Magnetomicrometry and Flu- oromicrometry. We subtracted time-synchronized, interpolated fluoromicrometry measure- ments from magnetomicrometry measurements across all 50 gait cycles of turkey running (10 gait cycles at each of 5 speeds, for each bird). Kernel density estimates show the distribution (vertical axis) of these differences (horizontal axis), with data plotted in blue, orange, and purple for Birds A, B, and C, respectively. Vertical lines indicate mean offsets (in corre- sponding colors), with the mean offsets and standard deviations labeled. Adjusted standard deviations that compensate for fluoromicrometry noise (0.098 mm, standard deviation) in- dicate an estimate of the magnetomicrometry measurement noise. two ramp inclines (10◦and 18◦, as shown in Fig. 3.9) and three vertical elevation changes (20 cm, 41 cm, and 61 cm, as shown in Fig. 3.10). Because the purpose of these activities was to explore the range of dynamic motions that could be captured, we did not train the birds to navigate the ramps or vertical elevation changes repetitively, and thus variability is expected within the repeated tasks. To further validate the accuracy of MM used during navigation of ramps and vertical elevation changes, we analyzed the magnetic bead tracking data from these activities to find the range of the tracked three-dimensional magnetic bead positions as shown in Fig. 3.11. We then affixed two magnetic beads 40 mm apart, validated the distance between them using FM (40.000 ± 0.017 mm), and swept this FM-validated magnetic bead pair under the MM sensing array through a volume exceeding these ranges as shown in Fig. 3.2. We monitoredFigure 3.7: Study-Specific Limitations to the Use of Fluoromicrometry. (A) Ten gait cycles of raw fluoromicrometry video data (Bird A, 3.5 m / s), independently manually re-labeled three times, are shown in light green, dark green, and brown. (B) The histogram shows the distribution (vertical axis) of the variance (horizontal axis) between the three manual processed fluoromicrometry signals at each timepoint of the data, with a kernel density estimation curve overlay. A vertical line indicates the average variance, 0.010 mm2, and the corresponding square root of the variance is 0.098 mm. deviations from 40 mm in the MM signal during these benchtop tests and found a 99-th percentile error (e99%) of 1.000 mm (rounded up to the nearest micrometer). Finally, to explore whether untethered muscle tracking via MM is viable in a fully free roaming context, we tracked muscle tissue length while one turkey (Bird A) roamed freely about its enclosure. The results of this data collection are shown in Fig. 3.12. 3.4 Discussion We find that MM enables untethered muscle tissue length tracking with high correlation to FM (R2of 0.952, 0.860, and 0.967 for Birds A, B, and C, respectively) and submillimeter accuracy (AVG±SD of -0.099 ± 0.186 mm, -0.526 ± 0.298 mm, and -0.546 ± 0.184 mm forFigure 3.8: Magnetomicrometry Tracking Time Delay. The magnet tracking computer recorded the times it received the magnetic field data and the times it completed the magnet tracking algorithm. The distribution (vertical axis) of the difference between these times (horizontal axis) is the tracking time delay and indicates the bandwidth at which magne- tomicrometry can track the muscle tissue length. The data are shown as a stacked histogram, with blue, orange, and purple data corresponding to Birds A, B, and C, respectively. Data are from all turkey gait cycles used to compare magnetomicrometry against fluoromicrometry. The 99th percentile time delay (t99%) is labeled for each bird. The ninety-ninth percentile time delay for all birds was less than one millisecond. Birds A, B, and C, respectively). These findings enable tracking and investigation of muscle contractile behavior in settings previously inaccessible to biomechanics researchers. 3.4.1 Accuracy Validation The standard we used here to assess the accuracy of muscle length tracking using MM was FM. For magnets implanted superficially in muscles (at depths less than 2 cm), MM exhibits less noise than FM, but FM has the advantage of higher accuracy, especially at greater tissue depths (tracking depths in this study ranged from 11.2 mm to 26.6 mm). Indeed, our tests showed that for unobscured markers moving through the X-ray volume, FM was accurate to 0.030 mm. However, we note that marker tracking noise was a challenge for FMFigure 3.9: Muscle Tissue Length During Non-Synchronous Ramp Ascent and Descent. We used magnetomicrometry to track muscle tissue length during ramp ascent and descent at two inclines for all three birds. Data for each bird and each slope are synchronized at right leg toe strike (indicated by the vertical gray line) and normalized from toe strike to toe strike. Variability between curves reflects gait cycle variability during untrained ramp navigation. Muscle tissue length is plotted in blue for right leg stance, in red purple for right leg swing, and in gray where video did not allow gait-phase labeling. We recorded at least three gait cycles of each activity for each bird. in this particular study due to the use of a large animal and the presence of hardware (the MM sensing array) that regularly obscured the markers during the tracking. These factors resulted in substantial manual labeling noise in the FM signal of 0.098 mm, instead of the 0.030 mm noise we found in our FM accuracy test, affecting the accuracy standard deviations reported above. Accounting for this manual labeling noise gives adjusted accuracy standard deviations of 0.158, 0.281, and 0.156 mm, for Birds A, B, and C, respectively as shown in Fig. 3.6. Constraints to imaging volume make FM impractical during large-animal variable terrain activity, so we performed retrospective benchtop accuracy testing to further validate the MM data collected during navigation of the ramps and vertical elevation changes as shownFigure 3.10: Muscle Tissue Length During Non-Synchronous Vertical Ascent and Descent. We used magnetomicrometry to track muscle tissue length during vertical ascent and descent at three heights for all three birds. Data for each bird and each height are synchronized at right leg toe-off (start of the aerial phase, indicated by the vertical gray line). Variability between curves reflects movement variability during untrained vertical ascent and descent. Muscle tissue length during contact with the ground is plotted in blue, and muscle tissue length during the aerial phase is plotted in redpurple. All data are shown, including scenarios in which significant wing-flapping occurred during jump up or down. We captured at least three recordings of each activity for each bird. in Fig. 3.2. Soft tissue artifacts during dynamic movements, such as tissue deformation or movement, result in depth and position changes of the magnetic beads relative to the MM sensing array. The benchtop tests investigated the accuracy of the MM measurements across the range of depths and positions of the magnetic beads that we observed during those activities as shown in Fig. 3.6. The error we observed in the benchtop tests (e99% < 1 mm) was acceptable in comparison with the magnitude of the muscle contractions we observed during the variable terrain activity (average MM signal magnitude was 4.5 mm peak-to-peak). This suggests that MM robustly tracked the muscle tissue lengths during the variable terrain activities, despite any soft tissue artifacts that may have occurred duringFigure 3.11: Extent of Bead Positions During Ramp Ascent and Descent and Vertical Ascent and Descent. Analysis of the magnetic bead position data (shown) recorded during all variable terrain activities revealed the extent of the magnetic bead positions during the tracking. We used these bead position maximums to design a benchtop test to verify the validity of our magnetomicrometry measurements during the variable terrain activities.Figure 3.12: Muscle Tissue Length During Free Roaming Movement. Magnetomicrometry data was continuously collected for 150 seconds during free roaming activity. Muscle tissue length is plotted in blue during standing and walking and plotted in black purple during running. Blue highlighted regions indicate muscle tissue length during (a) feather ruffling, (b) jumping, and (c) balancing on one leg. Gray arrows indicate when the turkey was turning left (left arrows) or turning right (right arrows). Gaps due to wireless transmission packet drops are shown in gray, as described in Figure 2.the dynamic movements required by those activities. These tests, however, highlight the importance of sensor placement. Higher accuracy is achieved when the MM sensing array is properly placed – centered over the implanted beads. MM with perfect magnetic field sensing would, in theory, be unaffected by movement of the board relative to the implanted beads, but the errors we observed suggest that the sensors are nonlinear. Magnet tracking nonlinearity compensation (e.g., via sensor calibration or three-dimensional sensor geometries) is thus an important area for future research. Meanwhile, in future work, larger sensing arrays with broader coverage would be advantageous to mitigate the need for careful placement of the array. 3.4.2 Ambient Magnetic Field The software-based magnetic disturbance compensation we employed here was sufficient to compensate for ambient magnetic fields during untethered muscle tracking in the presence of large hydraulic ferromagnetic lift tables, a large, active treadmill motor, and a room full of active X-ray equipment

[0090] . However, our uniform disturbance compensation strategy may be insufficient for the exceptional situation where a large ferromagnetic object is immediately adjacent to (within a few centimeters of) the tracked muscle. Thus, software-based compen- sation for spatially-non-uniform ambient magnetic fields may still be a valuable direction for future work to extend the robustness of MM to that potential scenario. Alternatively, ferro- magnetic shielding could be used to physically perform disturbance compensation, but the shield would need to be sufficiently far away to prevent it from acting like a magnetic mirror, creating “image” magnets that would need to be tracked as well

[0092] ,

[0093] . Further, effective shielding would need to be thick enough to redirect most or all magnetic field disturbances, presenting a trade-off between the weight and the efficacy of the shielding.3.4.3 Range of Behaviors Fig. 3.9, 3.10, and 3.12 provide a sample of the range of behaviors that can be tracked using magnetomicrometry. Consistency in the curves was not strived for, expected, or desired. Rather, we intentionally preserved anomalous events in those data, such as single or multiple wing flaps during vertical ascent and descent and variable speed during ramp navigation, to explore the range of motor activities during which we could track the muscle activity. 3.4.4 Potential Applications MM has the potential to work across scales, from the ability to track both full-body and muscle movement of small organisms to the ability to track large magnetic beads implanted deep into large animal models. Mathematically, if the number of sensors is fixed and all system dimensions are scaled, the error as a percent of scaled magnetic bead excursion will remain unchanged

[0090] . However, larger sensing arrays can in principle be used when tracking very small or very large animals, resulting in an increase in tracking accuracy at those extremes. For instance, when tracking small animal muscle tissue, additional sensors can be embedded into the a...

Claims

CLAIMS What is claimed is:

1. A human augmentation device configured to estimate the state of at least one implanted magnetic target, the device comprising: an array of sensors, each being configured to measure one or more target field characteristics; a first processing element coupled to the sensors, the processing element being configured to: acquire measurements of the one or more target field characteristics; apply information-theoretic processing of the measurements to estimate in real-time at least one of the state and error covariance of at least one target state vector, the target state vector comprising at least one of: target position, velocity, orientation, angular rate, or target strength; and update the at least one target state employing a state update innovation method that incorporates at least one of: (1) a posteriori covariance; (2) an observation matrix relating the change in a sensor characteristic to the change in the state of the at least one target; (3) a measurement variance; or (4) the innovation representing a difference between the measurement and its expected value; and a second processing element configured to control an actuator of the human augmentation device.

2. The device of Claim 1, wherein the first processing element reduces the error covariance by employing at least one stochastic dynamic model of the at least one target state, the model taking the form of at least one of: a Markov model driven by white noise; a dynamic system described by differential or difference equations driven by noise;an oscillator comprising time varying quadrature amplitudes and frequency with the amplitude and frequency driven by a Markov process itself driven by white noise; or a disturbance field that is varying in an inertial frame that is driven by a Markov model, itself driven by white noise or by a dynamic system described by differential or difference equations driven by noise.

3. The device of Claim 1, wherein the information-theoretic processing of the measurements is accomplished by an Information Filter implemented in which the time-varying observation matrix used in the information gain calculation is derived through real-time computation of the Jacobian that relates the change of a signal characteristic to a change in target state.

4. The device of Claim 1, wherein the filter employs knowledge of the measurement noise characteristics, the characteristics comprising: measurement noise variance unique to at least one of the signal characteristics and to the at least one sensor component in the array; and measurement noise correlation in time.

5. The device of Claim 1, wherein the a posteriori state covariance is assumed to have substantially achieved steady-state.

6. The device of Claim 3, wherein the Jacobian is computed at a nominal location of the target.

7. The device of Claim 1, wherein the innovation contribution is a sum of Information Vectors computed by one of a plurality of multiple processors within the processing element where each processor is connected to a sub-population of sensors and an Information Vector unique to a sensor is computed.

8. The device of Claim 1, wherein the observation matrix is computed via interpolation in relation to its value at a nominal location of the target, the type of interpolation being selected from: a Linear interpolation, a Bilinear interpolation, a Trilinear interpolation, oran interpolation based on computation of the Hessian matrix evaluated at the nominal location.

9. The device of Claim 1, wherein the processing element normalizes and applies known corrections to each in the array of sensor outputs to correct for at least one of the following error sources: translation of the sensor with respect to its nominal position within the sensor coordinate frame comprising at least one of a component of a Cartesian translation; orientation of the sensor with respect to the nominal orientation within the sensor coordinate frame comprising at least one of a rotation about x, a rotation about y and a rotation about z; orientation offset of the unit vectors corresponding to the principal coordinate axes along which signal characteristics are acquired by the sensor, where the offsets are two-degree of freedom offsets represented by a rotation about two of the coordinate-frame axes; sensitivity and non-linearity of the measurement of the signal characteristics on at least one of principal coordinate axes; or bias offset in the measurement of at least one of the signal characteristics.

10. The device of Claim 1 where the sensor array and the processing element are packaged onto at least one of one side or both sides of at least one of: a rigid-flex PCBA; a rigid PCBA; or interconnected rigid or rigid-flex PCBAs.

11. The device of Claim 1 wherein the first and second processing elements are packaged in the same assembly.

12. The device of Claim 11, wherein the first and second processing elements are packaged into one or more FPGAs.

13. The device of Claim 1, wherein the communications between the second processing element and the human augmentation device is carried out on a network.

14. The device of Claim 13, wherein the network synchronizes time using at least one of a Network Time Protocol and Precision Time Protocol.

15. The device of Claim 13, wherein the network is an Ethernet, EtherCAT, or of a Time- Sensitive Network (TSN) type.

16. The device of Claim 1, wherein the communications is run on a multi-drop synchronous serial interface network.

17. The device of Claim 2, further comprising an inertial measurement unit (IMU) arranged to report the angular rate of the device; and wherein the first or second processing element is configured to: predict temporal changes in the magnetic field owing to the instantaneous angular rate; and cancels the predicted component of the magnetic field measurement changes so as to reduce the covariance of the target state estimates.

18. The device of Claim 1, wherein the initial vector of the at least one target state is informed by a global optimization algorithm that minimizes the prediction residuals.

19. The device of Claim 1, wherein an information-theoretic smoothing algorithm is applied to the state estimates retrospectively and those smoothed state estimates determine a target state with substantially reduced covariance.

20. The device of Claim 10, wherein the rigid-flex PCBA can be used to create a tile structure where each tile comprises at least one magnetometer.

21. The device of Claim 2, wherein the oscillatory response is detected through use of an adaptive information-theoretic tracker that applies synchronous demodulation and phase detection to track the evolution of amplitude and frequency of the oscillation.

22. The device in Claim 1, wherein the innovation update is computed as the product of an innovation gain matrix and the innovation.

23. The device of Claim 1 wherein the actuator comprises at least one of an electric motor, a hydraulic motor, or a pneumatic motor activated by the processing element.

24. The device of Claim 1, wherein the actuator is muscle tissue that is activated by artificial muscle stimulation delivered to the tissue by the processing element.

25. The device of Claim 1, wherein global optimization is used to create the a priori information that serves to initialize the information-theoretic tracking algorithm.

26. An apparatus employed as part of a neuromuscular controller of a device aimed at augmenting human biomechanical function, the apparatus comprising: an array of sensors each designed to measure one or more target field characteristics; a processing element coupled to the sensors, the processing element configured to: acquire measurements of the field characteristics; and employ a neural processing unit to determine at least one of a set of response parameters comprising a muscle-tendon force, a time derivative of a muscle-tendon force, an elongation or a time derivative of an elongation; the at least one of a set of response parameters being determined based upon a trained model of the relationship between the targets and the response parameters; and an actuator that communicates with the processing element,wherein the processing element commands the actuator to deliver the determined dynamic response.

27. The apparatus of Claim 26, wherein the neural processing unit is trained by a series of synthesized Monte Carlo runs representing a stochastic model of the muscle-tendon unit dynamics and the target dynamics therein, wherein the stochastic model captures the range of variation in at least one of signal characteristics, magnet strength, range of magnet motion and geometric relation between the magnets and the sensing array; or disturbance state in the target state that will present to the apparatus when in use.

28. The apparatus of Claim 26, wherein the neural processing unit employs a Convolutional Neural Network (CNN) to process at least one of sensed target field characteristic on at least one target.

29. The apparatus of Claim 26, wherein the actuator is muscle tissue that is activated by artificial muscle stimulation delivered to the tissue by the processing element.

30. The apparatus of Claim 26, wherein the actuator comprises at least one of an electric motor, a hydraulic motor, or a pneumatic motor activated by the processing element.

31. A self-calibrating apparatus for use in tracking the state of magnetic targets, the apparatus comprising: at least one addressable coil that can project a reference magnetic field with known characteristics; an array of sensors each designed to measure at least one component of the magnetic field characteristics arising from the presence of the magnetic targets and the coil, a processing element connected to the array and the coil, the processing element programmed to compensate for at least one of a sensor bias, sensitivity and non-linearity, or geometric alignment as a means of improving at leastone of precision and accuracy of the target state tracking, the target state tracking derived from the sensor measurements, where the target state comprises at least one of a target position, velocity, or orientation within an apparatus coordinate frame.

32. The apparatus of Claim 31, wherein the optimization is a global optimization.

33. The apparatus of Claim 31, wherein finding the coordinate transformation that minimizes the innovation employs an information-theoretic tracking algorithm.

34. The apparatus of Claim 31, wherein the global optimization is used to create the a priori information that serves to initialize the information-theoretic tracking algorithm.

35. The apparatus of Claim 31, wherein the tiled-sensor array and the modulated coil are integrated into the same rigid printed circuit assembly or into the same rigid printed circuit subassembly.

36. The apparatus of Claim 31, where the tiled-sensor array and the modulated coil are integrated into the same rigid-flex circuit assembly or into the same rigid-flex printed circuit subassembly.

37. The apparatus in Claim 31, where the global optimization is accomplished by gradient descent and not a Jacobian.

Images