A multi-imu and plantar pressure synchronous acquisition system and method based on master-slave architecture

Abstract

This invention discloses a lower limb motion and plantar pressure synchronous acquisition system and joint angle measurement method based on a master-slave architecture. System: The main control board acquires quaternions from five IMUs (Imported Mutual Arrays) for the waist and legs via I2C; the left and right foot slave boards scan a 45-point plantar pressure array and acquire ankle joint IMU data via SPI interface and DMA state machine, packaging the data into frames with CRC16 checksums and uploading them via serial port. The main control board aggregates the data into 118-dimensional data packets and uploads them at 100Hz. Measurement method: The quaternions of the seven IMUs are normalized and converted into Euler angles; a dynamic time warping algorithm is used to align the reference signal; noise is removed using a fourth-order bandpass filter; mixed time-domain and frequency-domain features are extracted; and joint angles are predicted using a Stacking ensemble learning model. This invention achieves high real-time and high-precision synchronous acquisition and analysis of lower limb multimodal motion information, applicable to gait analysis and rehabilitation training.

Description

Technical Field

[0001] This invention relates to wearable sensor systems, and more particularly to a multi-IMU and plantar pressure synchronous acquisition system based on a master-slave architecture for lower limb joint angle measurement and gait analysis. Background Technology

[0002] Lower limb joint angle measurement is a fundamental technique for gait analysis, rehabilitation assessment, and sports injury prevention. Clinically, doctors need to use the flexion and extension angles of the hip and knee joints to determine whether a patient's gait is abnormal and to develop targeted rehabilitation training programs. In recent years, with the development of exoskeleton robots, real-time and accurate joint angle data has also become a key input for human-machine collaborative control.

[0003] The most mature measurement solution currently is the optical motion capture system. This type of system uses multiple high-speed cameras to track reflective markers on the human body surface, reconstructing a three-dimensional motion trajectory, with angular measurement errors controlled within 1°. However, optical systems are expensive, must be used in specialized laboratories, and the markers are easily obstructed, making them unsuitable for outdoor measurements or home rehabilitation needs.

[0004] Existing wearable solutions based on inertial measurement units (IMUs) have the following significant shortcomings. While nine-axis sensors worn on multiple body parts can calculate the relative posture of each joint, these products are highly portable but still have significant drawbacks: magnetometers are susceptible to interference from indoor metal structures, and gyroscope integration drift can lead to a decrease in measurement accuracy over long periods; wireless transmission latency is generally 60-120 milliseconds, making it difficult to use in scenarios requiring real-time feedback, such as exoskeleton control which typically requires latency of less than 20 milliseconds.

[0005] Existing solutions also have problems with multi-sensor fusion. When simultaneously acquiring IMU and plantar pressure data, most systems use independent acquisition and post-processing timestamp alignment, which easily leads to data asynchrony. Some high-end systems use hardware synchronization lines for triggering, but this increases wiring complexity and is difficult to implement in wireless scenarios. In addition, traditional algorithms that directly use Euler angles are at risk of gimbal lock, while machine learning models are mostly single-network structures with limited generalization ability across different individuals and movement patterns. Plantar pressure data is closely related to gait phases, but existing IMU systems generally do not utilize this complementary information.

[0006] Therefore, developing a low-cost, high-precision, real-time response lower limb measurement system is of practical significance. An ideal solution should integrate multiple IMUs and pressure sensors, simplify the communication topology through a master-slave architecture, employ DMA and a circular buffer to ensure data synchronization, utilize dynamic time warping algorithms to handle time offsets, and combine this with an ensemble learning model to improve prediction robustness. Ultimately, this would achieve a 100Hz sampling frequency, end-to-end latency of less than 20 milliseconds, and system costs controlled within several thousand yuan, thus promoting the application of this technology in community healthcare and home rehabilitation scenarios. Summary of the Invention

[0007] Purpose of the invention: To overcome the shortcomings of existing technologies, this invention provides a multi-IMU and plantar pressure synchronous acquisition system based on a master-slave architecture. The technical problems to be solved are: how to achieve quasi-synchronous acquisition of multiple IMUs and plantar pressure sensors under low cost, how to ensure the reliability and real-time performance of high-speed data transmission, how to eliminate the impact of sampling time offset on data analysis, and how to improve the accuracy and robustness of joint angle prediction through algorithm optimization, ultimately realizing real-time measurement of lower limb joint angles applicable to gait analysis, rehabilitation assessment, exoskeleton control and other application scenarios.

[0008] Technical solution: The multi-IMU and plantar pressure synchronous acquisition system based on a master-slave architecture provided by the present invention includes: The main control board is connected to five first inertial measurement units via a first bus. The five first inertial measurement units are respectively located at the waist, left thigh, right thigh, left calf, and right calf. The main control board also receives data from the left foot slave plate and the right foot slave plate via a second bus. The left foot follower plate and the right foot follower plate are used to collect quaternion data from the second inertial measurement unit of the ankle joint and pressure data from the plantar pressure sensor array, respectively, and package the collected data into data frames containing frame header, frame sequence number, data content and cyclic redundancy check code and upload them to the main control board. The main control board and the slave board operate independently at a preset sampling frequency. The main control board summarizes the quaternion data and pressure data received from all inertial measurement units into a data packet of a predetermined dimension and uploads it.

[0009] Furthermore, the slave board acquires pressure data from 45 measuring points in the plantar pressure sensor array via multiplexing scanning; the main control board reads data from the first inertial measurement unit via direct memory access and receives and parses data frames from the slave board via a circular buffer; the preset sampling frequency is 100Hz, and the predetermined dimension data packet is 118-dimensional.

[0010] The method for measuring lower limb joint angles based on the aforementioned system includes the following steps: S1: Collect quaternion data and plantar pressure data from all inertial measurement units; S2: Normalize each quaternion q = (q0, q1, q2, q3): Calculate the modulus. ,like Replace with a unit quaternion; otherwise, perform normalization. ; S3: Calculate the relative rotation quaternions of each joint in the lower limb based on the normalized quaternions, where the conjugate of the quaternion is... And convert the relative rotation quaternion into Euler angles; S4: The Euler angle sequence is aligned with the reference signal using a dynamic time warping algorithm, and the aligned sequence is then bandpass filtered. S5: Extract time-domain features and frequency-domain features based on Welch power spectrum estimation from the filtered signal; S6: Input the extracted features into the Stacking ensemble learning model to predict the joint angles. .

[0011] Further, in step S3, the relative rotation quaternions (q0, q1, q2, q3) are converted into Euler angles in ZYX order. The formulas for rotation α, abduction β, and flexion / extension γ are as follows:

[0012] Furthermore, in step S4, the dynamic time warping algorithm constructs a cost matrix. And backtrack to achieve the optimal path.

[0013] Furthermore, in step S4, the bandpass filter adopts a 4th-order Butterworth bandpass filter with a passband frequency range of 0.5Hz to 5Hz, and a zero-phase filtering method is used to eliminate phase delay.

[0014] Further, in step S5, the frequency domain features are obtained based on Welch power spectrum estimation. The Welch power spectrum estimation uses a Hanning window with a window length of 256 points, an overlap rate of 50%, and 512 FFT points. The power spectral density is calculated using the following formula: .

[0015] Further, in step S6, the first base learner of the Stacking ensemble learning model includes a random forest, a gradient boosting tree, and a ridge regression model, and uses 5-fold cross-validation to generate meta-features; the parameters of the random forest are n_estimators=400, max_depth=9; the parameters of the gradient boosting tree are n_estimators=400, max_depth=7.

[0016] Furthermore, the second-layer meta-learner of the Stacking ensemble learning model is a Bayesian ridge regression model, which updates hyperparameters iteratively. Training was conducted, among which .

[0017] Furthermore, in step S6, the predicted joint angles Calculated by the following formula: in, , , These are the predicted outputs of the first-layer random forest, gradient boosting tree, and ridge regression model for the input feature x, respectively, where w0, w1, w2, and w3 are the weight coefficients learned by the Bayesian ridge regression model.

[0018] Furthermore, the main control board is equipped with an STM32H7R7 microcontroller and connects to five HWT606 inertial measurement units via an I2C bus. These units are fixed at the waist, both thighs, and both calves, with I2C addresses configured from 0x50 to 0x54. The main control board receives data from the left and right foot slave boards via two UARTs and uploads the aggregated 118-dimensional data to the host computer in real time via a third UART. The slave boards are equipped with STM32G431 microcontrollers, connected to the ankle joint IMU via a serial port and to a 45-point plantar pressure sensor array via an SPI interface.

[0019] Both the master and slave boards operate independently with a 100Hz timer, requiring no hardware synchronization line. The slave board packages the 4D quaternion and 45D pressure data into a 206-byte frame, including a frame header of 0xAA55, a frame sequence number, and a CRC16 checksum, ensuring reliable transmission. The master board uses DMA cyclic reception and a 512-byte circular buffer to extract complete data frames, summarizes the 118-dimensional data, and uploads it at a 100Hz frequency, with an end-to-end latency of less than 15 milliseconds.

[0020] Both the left and right foot slave boards are STM32G431, including: a second processor; a second IMU, fixed to the ankle joint and connected via a serial port; a plantar pressure sensor array containing 45 measuring points, connected via SPI+ADC; and a slave communication interface. Synchronization mechanism: The main control board and the slave board both operate independently at 100Hz. The main control board reads five first IMU quaternions sequentially via DMA every 10ms. The slave board actively uploads a 206-byte data frame (including frame header 0xAA55, frame number, 4-dimensional quaternion, 45-dimensional pressure, and CRC16 checksum) after acquiring the second IMU and pressure data every 10ms. Data output: The main control board aggregates the 20-dimensional first IMU quaternion, the 8-dimensional second IMU quaternion, and the 90-dimensional stress data into a 118-dimensional data packet, which is then uploaded to the host computer at a frequency of 100Hz.

[0021] The five primary IMUs are configured with I2C addresses 0x50 (waist), 0x51 (right thigh), 0x52 (right calf), 0x53 (left thigh), and 0x54 (left calf) via hardware address pins. The primary processor reads the data sequentially using DMA and a state machine, with each round of data acquisition taking 3-5ms.

[0022] The plantar pressure sensor array contains 6 X channels and 10 Y channels, and scans 45 measurement points through multiplexing.

[0023] The second processor uses a 9-state finite state machine to control SPI+DMA to complete the data acquisition. The states include: initialization, two warm-ups, two pre-settings, cyclic data acquisition, two terminations, and switching to the next set. After completing 6 sets, it outputs a 45-dimensional pressure vector.

[0024] The main control board receives data from the slave board through two independent serial port DMA loops, using a 512-byte circular buffer. The receiving process includes: calculating the write pointer position using the DMA counter; scanning the 0xAA55 frame header from the read pointer to the write pointer; copying the complete 206-byte frame across boundaries; calculating the CRC16 checksum to verify data integrity; and detecting packet loss by using the frame sequence number difference.

[0025] The third communication interface uses UART+DMA double buffering for transmission, which includes two 1024-byte buffers. During the DMA transmission of the current buffer, the next buffer is simultaneously formatted. Data formatting uses a custom fast floating-point to string function that directly manipulates pointers, with a fixed 3 decimal places. A single floating-point conversion takes about 10μs, and a complete 118-dimensional data format takes about 1.2ms, meeting the requirements for 100Hz output.

[0026] The CRC16 check uses the CRC-16 / MODBUS standard, with an initial value of 0xFFFF and a polynomial of 0xA001. It calculates the first 204 bytes of the data frame and compares them with the last 2 bytes of the frame. If they do not match, the frame is discarded.

[0027] This invention also provides a joint angle measurement method, which includes four processing steps: data acquisition and transmission, quaternion solution and time alignment, feature extraction and filtering, and machine learning regression.

[0028] In the data acquisition stage, the main control board sequentially reads 5 IMU quaternions through I2C DMA, and the slave board scans 45 pressure measurement points through the SPI state machine and packs and uploads them. In the quaternion calculation stage, the quaternions are normalized and checked, the relative rotation quaternion of the joint is calculated, and it is converted into Euler angles in the ZYX order, and the sampling time offset between this system and the reference system is eliminated through the dynamic time warping algorithm. In the feature extraction stage, noise is removed through a 4th-order Butterworth band-pass filter (0.5 - 5 Hz), time-domain features such as angles, angular velocities, and angular accelerations are extracted, and frequency-domain features such as average power and main frequency are extracted through Welch power spectrum estimation to form a 7-dimensional feature vector. In the machine learning stage, a Stacking ensemble learning architecture is adopted. The first layer contains three base learners, RandomForest, GradientBoosting, and Ridge, and meta-features are generated through 5-fold cross-validation. The second-layer BayesianRidge meta-learner adaptively fuses and outputs the final joint angle.

[0029] The measurement method is preferably: S1: Data acquisition The main control board reads 5 first IMU quaternions through I2C + DMA; the slave board collects the second IMU and pressure data, packs them into a 206-byte frame and uploads them through the serial port; the main control board analyzes the data frame, aggregates them into a 118-dimensional data packet and uploads it at 100 Hz; S2: Data preprocessing Perform normalization check on 7 quaternions and calculate the norm , if then replace it with the unit quaternion; perform normalization ; S3: Kinematic calculation Calculate the relative rotation quaternion of the joint, for example, the hip joint: ; Convert to ZYX Euler angles (rotation - abduction - flexion and extension):

[0030] S4: Dynamic time warping (DTW) alignment Construct the cost matrix:

[0031] Trace back from the lower right corner to the upper left corner, and resample and query the signal according to the optimal path; S5: Signal filtering Design a 4th-order Butterworth band-pass filter (0.5 - 5 Hz), and use zero-phase filtering to eliminate phase delay; S6: Feature engineering Calculate time-domain features (angle, angular velocity, angular acceleration) and frequency-domain features (average power, standard deviation, dominant frequency, and number of spectral lines of the Welch power spectrum) to construct a 7-dimensional feature vector; S7: Ensemble Learning Regression Training Level-0 base learners: RandomForest (n_estimators=400, max_depth=9), GradientBoosting (n_estimators=400, max_depth=7), Ridge (α-grid search); generating meta-features using 5-fold cross-validation; training Level-1 meta-learners: BayesianRidge, final prediction:

[0032] S8: Performance Evaluation Calculate the Pearson correlation coefficient and standard deviation, requiring flexion and extension... Outreach Rotation , radian.

[0033] The fast floating-point to string conversion algorithm includes: handling negative numbers, rounding to 3 decimal places (value += 0.0005f), separating the integer and decimal parts, reversing the integer to a character, writing the decimal point, and outputting the decimal with a fixed 3 decimal places. It is about 20 times faster than sprintf.

[0034] The quaternion conjugate is Quaternion multiplication is:

[0035] The DTW algorithm employs dynamic programming, with a time complexity of [missing information]. Optimize space using a rolling array When backtracking, prioritize the diagonal (synchronization point), and then choose vertical or horizontal (compression / stretch).

[0036] The Butterworth filter has an order of n=4 and a cutoff frequency of... rad / s, rad / s; Zero-phase filtering is achieved through forward filtering, reverse filtering, and then reverse filtering again, ultimately resulting in a phase delay of 0.

[0037] The Welch power spectrum estimation uses a Hanning window with a window length of 256 points, an overlap rate of 50%, 512 FFT points, and a frequency resolution of 0.195 Hz; power spectral density: .

[0038] The Stacking ensemble uses 5-fold cross-validation to generate meta-features, and BayesianRidge automatically adjusts hyperparameters.

[0039] in .

[0040] Compared with the prior art, the present invention has the following significant advantages: The cost is significantly reduced. The material cost of the complete system is approximately 2000 RMB, only 4% of that of commercial systems, making it possible for applications in primary healthcare institutions and home rehabilitation. Real-time performance is excellent, with end-to-end latency of less than 15 milliseconds and stable 100Hz output, meeting the real-time application requirements of exoskeleton control and other applications. Measurement accuracy is high, with a Pearson correlation coefficient of 0.96 for hip flexion-extension and 0.94 for knee flexion-extension, approaching the level of optical systems. Data synchronization is reliable; the quasi-synchronous design eliminates the need for hardware synchronization lines, and CRC16 checksum and frame sequence number mechanism ensure a packet loss rate of less than 0.1%. The algorithm is robust; the dynamic time warping algorithm adapts to different paces and individual differences, and Stacking ensemble learning performs stably under different subjects and movement modes. The system has good scalability; the master-slave architecture easily expands to upper limb monitoring, and the 118-dimensional raw data supports various post-processing analyses such as gait phase recognition and balance assessment. In summary, this invention achieves low-cost, high-precision, real-time, and reliable lower limb joint angle measurement, expanding the application scope of wearable sensor technology in rehabilitation medicine, sports science, exoskeleton robotics, and other fields, and has good application prospects and social value. Attached Figure Description

[0041] Figure 1. Schematic diagram of the overall system architecture; Figure 2. System control flowchart; Figure 3. Block diagram of hierarchical control system; Figure 4. Schematic diagram of the control module; Figure 5. Communication and data flow architecture diagram of master and slave boards. Detailed Implementation

[0042] The following description, with reference to the accompanying drawings, illustrates several preferred embodiments of the present invention to make its technical content clearer and easier to understand. The present invention can be embodied in many different forms, and the scope of protection of the present invention is not limited to the embodiments mentioned herein.

[0043] like Figure 1As shown in this embodiment, the lower limb motion and plantar pressure synchronous acquisition and joint angle measurement system based on a master-slave architecture adopts a hierarchical processing architecture, forming a complete data processing chain from the bottom-level hardware acquisition to the top-level algorithm prediction. The overall system architecture includes four core modules: the master-slave board synchronous data acquisition module 1 is responsible for the real-time acquisition and aggregation of raw sensor data, outputting a 118-dimensional raw data vector; the quaternion solving module 2 converts the quaternion data into relative joint angles, outputting an 18-dimensional angle sequence; the feature extraction and filtering module 3 performs signal processing and feature engineering on the angle sequence, outputting a 7-dimensional hybrid feature; and the machine learning regression module 4 predicts the final joint angle through ensemble learning.

[0044] like Figure 2 As shown, the system data flow starts from the physical layer sensors, undergoes analog-to-digital conversion at the hardware acquisition layer, and is aggregated via master-slave board communication at the data transmission layer. Quaternion solving, time alignment, and signal filtering are performed at the data processing layer, and finally, the predicted joint angle values ​​are output at the algorithm layer. The entire process runs in real-time at a frequency of 100 Hz, with an end-to-end latency of less than 15 ms, meeting the real-time requirements of rehabilitation assessment and gait analysis.

[0045] like Figure 5 As shown, this invention employs a master-slave architecture to achieve synchronous data acquisition from multiple sensors. The main control board is equipped with an STM32H7R7 microcontroller (ARM Cortex-M7 core, 600 MHz clock speed, built-in single-precision floating-point arithmetic unit), which connects to five HWT606 inertial measurement units in DMA mode via an I2C bus. These units are fixed at the waist, right thigh, right calf, left thigh, and left calf, respectively. The I2C bus is configured in fast mode (400 kHz clock frequency), and the device address is configured from 0x50 to 0x54 via the hardware pin AD0.

[0046] The main control board triggers the data acquisition task with a 00 Hz timer interrupt. The timer auto-reload value is precisely configured using the following formula:

[0047] in, MHz is the system clock frequency. Hz is the sampling frequency. This is the prescaler coefficient. This configuration ensures that the timer triggers an interrupt precisely every 10ms, with a measured period standard deviation of less than 50μs, meeting the accuracy requirements for synchronous acquisition.

[0048] Each foot slave board uses an STM32G431 microcontroller (ARM Cortex-M4 core, 170 MHz). Each slave board connects to an HWT606 ankle joint IMU via a UART serial port and to the plantar pressure sensor acquisition circuit via an SPI interface. The plantar pressure sensor array uses 45-point (9 rows × 5 columns) resistive thin-film sensors, with a single point size of 10 mm × 10 mm and a resistance range of 1-100 kΩ corresponding to a pressure of 0-500 kPa. The slave board sequentially selects the 45 sensors using three 8-channel multiplexers 74HC4051, and uses a 12-bit SAR-type ADC chip ADS7883 to complete the analog-to-digital conversion.

[0049] The ADC sampled values ​​are converted into physical pressure values ​​using calibration coefficients:

[0050] in, V is the reference voltage. For the first Zero-point offset voltage of each sensor, These are the calibration coefficients. The calibration coefficients are obtained by applying known weights (100 g, 200 g, 500 g) and fitting the data using the least squares method. They are stored in the flash memory of the slave board for automatic loading.

[0051] Seven IMU sensors were secured to specific locations on the body using elastic straps, ensuring a rigid connection between the sensors and the limb segments. The lumbar IMU was fixed at the level of the L4-L5 vertebral bodies, with the X-axis pointing forward (sagittal plane), the Y-axis pointing to the left (coronal plane), and the Z-axis vertically upward (horizontal plane). The thigh IMU was fixed on the lateral aspect of the mid-shaft of the femur, approximately one-third of the way from the hip joint. The lower leg IMU was fixed on the lateral aspect of the mid-shaft of the tibia, approximately one-third of the way from the knee joint. The ankle IMU was fixed 5 cm above the lateral malleolus. All IMUs used a right-handed coordinate system, with the positive Z-axis pointing opposite to the direction of gravity. During initial calibration, the subject maintained an upright standing posture to determine the zero-point reference.

[0052] like Figure 5 As shown, the master and slave boards use UART serial communication for data transmission at a baud rate of 921600 bps and a data format of 8N1. The slave board packages the 4D ankle joint quaternion and 45D pressure data into a 206-byte data frame. The frame format includes: frame header (2 bytes, fixed value 0xAA55 for frame synchronization), data length (2 bytes, fixed value 200), frame sequence number (4 bytes, for packet loss detection), quaternion (16 bytes), pressure data (180 bytes), and CRC16 checksum (2 bytes).

[0053] The CRC16 checksum uses the CRC-16 / MODBUS standard, with an initial value of 0xFFFF, a polynomial of 0xA001, and an iteration formula as follows:

[0054] This verification algorithm can detect single-bit errors, double-bit errors, and burst errors of less than 16 bits in length during data transmission, ensuring data integrity.

[0055] The main control board connects to the left and right slave boards via UART1 and UART2 respectively, and uses DMA circular receive mode to listen to serial port data. The receive buffer size is 512 bytes, using a circular buffer management mechanism. The readable data length is calculated using the following formula:

[0056] in, Write pointer to DMA. The software reads the pointer. The main control board searches for the frame header 0xAA55 using a sliding window algorithm, extracts the complete data frame, verifies the CRC check, and parses the frame sequence number to detect packet loss. The number of lost packets is calculated using the difference in frame sequence numbers.

[0057] like Accumulate packet loss counters and record anomaly flags to ensure data quality is traceable.

[0058] After the main control board completes data aggregation, it packages a total of 118 dimensions of data, including quaternions from 5 waist and leg IMUs (20-dimensional), 2 ankle joint IMUs (8-dimensional), and 2 pressure array data points (90-dimensional). This data is then uploaded to the host computer via UART8 at 921600 bps. The total data frame length is 482 bytes (472 bytes of data + 8 bytes of timestamp + 2 bytes of CRC16), with a transmission time of approximately 4.2 ms. The system end-to-end latency can be broken down as follows: 6.0 ms for sensor acquisition from the board + 2.2 ms for UART transmission + 1.5 ms for main control board data parsing + 0.8 ms for data aggregation + 4.2 ms for UART transmission = 14.7 ms, meeting the 100 Hz real-time output requirement.

[0059] like Figure 1 , Figure 2 , Figure 4As shown, the quaternion calculation module receives 118-dimensional data uploaded from the master and slave boards. First, it extracts 28-dimensional quaternion data from seven IMUs (each IMU outputs a 4-dimensional quaternion). Data validity is ensured through normalization checks. Then, it calculates the relative rotation quaternions for six pairs of joints, finally converting them into ZYX Euler angles to output an 18-dimensional joint angle sequence (three degrees of freedom for each joint pair: flexion / extension, abduction, and rotation). The output of this module provides the foundational data for subsequent filtering and feature extraction.

[0060] Quaternion Representing rotation in three-dimensional space, the modulus of the unit quaternion is 1:

[0061] The received quaternion data is first subjected to a normalization check; if... (threshold) If any of the above conditions are found, it is determined to be a sensor malfunction or data transmission error. For abnormal data, if there are 3 consecutive abnormal frames, a sensor fault flag is recorded; if the abnormality is intermittent (<3 frames), linear interpolation is used to fill in the gaps.

[0062] This interpolation method maintains the continuity of the angle sequence and avoids drastic jumps that could affect subsequent filtering results.

[0063] like Figure 4 As shown, the relative angles of the joints are calculated using the proximal coordinate system. and remote coordinate system The relative rotation quaternion is obtained. The conjugate (i.e., inverse operation) of a quaternion is defined as:

[0064] Relative rotation quaternions are calculated using quaternion multiplication:

[0065] Quaternion multiplication The expanded form is:

[0066] For the six pairs of joints, the relative rotation calculation follows the anatomical kinetic chain relationship: right hip joint (Pelvic rotation relative to the thigh); Right knee joint (Rotation of the thigh relative to the lower leg); right ankle joint (Rotation of the lower leg relative to the foot); the calculation method for the three pairs of joints on the left side is symmetrical. This calculation process ensures that the joint angles reflect the relative motion between adjacent limb segments, rather than the absolute posture relative to the global coordinate system.

[0067] This invention uses the ZYX Euler angle sequence, and the conversion sequence is as follows: first rotate the yaw angle around the Z-axis. (Joint rotation), then rotate the pitch angle around the new Y-axis. (Joint abduction), and finally, rotation angle around the new X-axis. (Joint flexion and extension). From relative rotation quaternions The conversion formula to ZYX Euler angles is:

[0068] To eliminate sensor installation errors, the subject is required to maintain an upright standing posture when the system is started, and the initial quaternions of all IMUs are recorded. As the zero-point reference, the formula for calculating the relative quaternion is modified during subsequent motion as follows:

[0069] This formula automatically cancels out the initial posture offset, bringing all joint angles in the upright standing position close to 0°. At this point, the quaternion solving module converts the 28-dimensional quaternion data into an 18-dimensional joint angle sequence, providing input data for subsequent dynamic time warping.

[0070] like Figure 1 and Figure 2 As shown, this invention solves the sampling time offset problem between the current system and the reference system using the Dynamic Time Warping (DTW) algorithm. In practical applications, due to the independent clock sources, different sampling triggering mechanisms, and inconsistent data transmission delays of the two systems, even with synchronization via hardware trigger signals, there is still an uncertain time offset (typically 0-100 ms). Directly aligning and comparing based on timestamps will lead to significant errors due to sampling time mismatch, especially during rapid joint angle changes (such as the landing moment and acceleration phase of the gait cycle). The DTW algorithm constructs a cost matrix to find the optimal time alignment path, achieving the best morphological match between the two sequences, thereby eliminating the influence of time offset.

[0071] Given the joint angle sequence measured by this system and reference system sequence The DTW algorithm first constructs a local cost matrix and defines the local distance between two time points as the Euclidean distance:

[0072] Then, the cumulative cost matrix is ​​calculated recursively using dynamic programming. :

[0073] Boundary conditions are , ( Vertical movement indicates Skipping a frame in the sequence (the sampling rate of this system is temporarily higher than that of the reference system), horizontal movement indicates The sequence skips one frame (the reference system's sampling rate is temporarily higher than this system's), and diagonal movement indicates that the two sequences are moving synchronously. To avoid the path deviating excessively from the diagonal, a Sakoe-Chiba band constraint is introduced:

[0074] in, This is the maximum allowed time offset. This constraint reduces the time complexity from... Reduce to This significantly improves computational efficiency.

[0075] from Backtracking To obtain the optimal alignment path The backtracking rules are as follows:

[0076] in, This is the set of predecessor nodes at the current position.

[0077] The system sequence is obtained through linear interpolation. Resampling to the reference system timeline Aligned sequence Calculated using the following formula:

[0078] in, and This aligns adjacent points in the alignment path. Alignment quality is assessed using the Pearson correlation coefficient, which is typically [value missing] before alignment. After alignment, it is promoted to The root mean square error was reduced from 8.3° to 3.2°. At this point, the time alignment module eliminated the influence of sampling time offset, outputting an angle sequence synchronized with the reference system's time axis, providing high-quality input data for subsequent filtering and feature extraction.

[0079] like Figure 1 , Figure 3 , Figure 4As shown, the feature extraction and filtering module receives the 18-dimensional joint angle sequence aligned by DTW. First, it removes low-frequency drift and high-frequency noise using a Butterworth bandpass filter. Then, it extracts time-domain and frequency-domain features, and finally, it uses principal component analysis to reduce the dimensionality and obtain a 7-dimensional hybrid feature vector. This module is designed based on the spectral characteristics of human gait: the dominant gait frequency is concentrated in the 0.5-5 Hz range (step frequency approximately 1-2 Hz, second harmonic approximately 2-4 Hz, third harmonic approximately 3-5 Hz). Components below 0.5 Hz mainly consist of low-frequency sensor drift and baseline offset, while components above 5 Hz mainly consist of high-frequency noise and muscle tremors. Therefore, bandpass filtering can retain effective motion signals while suppressing noise interference.

[0080] This system uses a 4th-order Butterworth bandpass filter, with the following design parameters: lower passband frequency. Hz, upper limit of passband frequency Hz, sampling frequency Hz. The formula for calculating the normalized cutoff frequency is:

[0081] A fourth-order Butterworth bandpass filter consists of a cascaded second-order high-pass filter and a second-order low-pass filter. The transfer function of the high-pass filter is:

[0082] The transfer function of the low-pass filter is:

[0083] in, The quality factor for the Butterworth filter ensures the flattest amplitude-frequency response within the passband. The analog filter is converted to a digital filter using a bilinear transform:

[0084] in, s is the sampling period. The transformed result is the difference equation of the second-order digital filter:

[0085] To avoid phase delay introduced by the filter, a zero-phase filtering method is used: the signal is processed twice, forward and reverse filtered, to automatically cancel out phase distortion. The filtered 18-dimensional angle sequence then enters the feature extraction stage.

[0086] For each joint angle sequence Extract the following features. Temporal features include: the mean angle reflects the average pose of the joint.

[0087] The standard deviation of the angle reflects the range of motion of the joint.

[0088] The angular velocity is calculated using the first-order difference method, and its mean is extracted. and standard deviation :

[0089] in, ms is the sampling interval.

[0090] Frequency domain features were obtained through Welch power spectrum estimation. The sequence was divided into 256-point windows with 50% overlap. A Hanning window was applied to each window, followed by a 512-point FFT. The power spectrum of all windows was then averaged. The power spectral density was calculated using the following formula:

[0091] Extract the average power (reflecting the total signal energy), power standard deviation (reflecting the power spectrum flatness), and dominant frequency (corresponding to step frequency or harmonic frequency):

[0092] Seven-dimensional features (mean angle, standard deviation angle, mean angular velocity, standard deviation angular velocity, average power, standard deviation power, and dominant frequency) were extracted from each of the 18 joint angles, resulting in a total of 126-dimensional feature vectors. To reduce computational complexity and avoid overfitting, principal component analysis (PCA) was used to retain the top 20 principal components with a cumulative contribution rate >95%, which were then used as input features for the machine learning module.

[0093] like Figure 1 , Figure 4 As shown, the machine learning regression module employs a two-layer Stacking ensemble learning architecture. The input is 20-dimensional PCA features, and the output is 18-dimensional joint angle predictions. Stacking overcomes the limitations of a single model by integrating the advantages of multiple heterogeneous base learners. The three base learners in the first layer capture different mapping patterns between features and angles: RandomForest handles non-linear relationships through decision tree ensemble, GradientBoosting optimizes detail fitting through residual iteration, and Ridge Regression provides a global linear benchmark. The second-layer BayesianRidge meta-learner learns how to optimally combine the predictions of the three base learners, adaptively adjusting the weights of each model to achieve higher prediction accuracy and robustness than a single model.

[0094] The first layer contains three heterogeneous base learners. RandomForest ensembles 400 decision trees with a maximum depth of 9 and a minimum number of split samples of 5. A single tree minimizes Gini impurity through recursive splitting:

[0095] in, For the first Class sample proportion. GradientBoosting integrates 200 decision trees with a maximum depth of 7 and a learning rate of 0.1. Gradient boosting progressively fits the residuals through an additive model:

[0096] in, For the first Wheel model pseudo residuals fitted to the current decision tree , The learning rate is used, and the loss function is the mean squared error. Ridge regression suppresses overfitting by minimizing the following objective function:

[0097] Wherein, regularization coefficient The method is determined through grid search. The three base learners are trained in parallel on the same training set, each outputting a predicted value, forming a three-dimensional prediction vector that serves as the input to the meta-learner.

[0098] The first layer of base learners is generated using 5-fold cross-validation. The training set is randomly divided into 5 mutually exclusive subsets, and for each fold... Use the remaining 4 folds to train three base learners, and use the trained model to predict the first... The matrix is ​​folded five times and then concatenated to obtain the complete meta-feature matrix. The meta-learner uses Bayesian Ridge regression, and the regularization coefficients are adaptively determined through maximum a posteriori estimation. The iterative formula is as follows:

[0099] in, , These are the singular values ​​of the characteristic matrix.

[0100] During the testing phase, new samples Predictions are obtained by sequentially passing the data through three base learners. Then, input the result into the meta-learner to obtain the final joint angle:

[0101] This two-layer architecture enables end-to-end prediction from 118-dimensional raw sensor data to 18-dimensional joint angles, with Pearson correlation coefficients. The root mean square error (RMSE) is < 5°, which meets the accuracy requirements for rehabilitation assessment and gait analysis.

[0102] like Figure 2 and Figure 3 As shown, the system implements a complete data flow from the physical layer sensors to the algorithm layer prediction through a hierarchical control architecture. The physical layer uses seven IMUs and two pressure sensors (90 pressure change points) to acquire raw motion data in real time. The master-slave board architecture of the hardware acquisition layer completes analog-to-digital conversion and data packaging. The data transmission layer ensures data integrity through UART and CRC checksums. The data processing layer sequentially performs quaternion solving, DTW time alignment, Butterworth filtering, and feature extraction. The algorithm layer outputs the final joint angle prediction through stacking ensemble learning. The modules are loosely coupled through standardized data interfaces, ensuring both real-time performance and good scalability and maintainability.

[0103] The above are merely preferred embodiments of the present invention and do not constitute any limitation on the present invention. Any equivalent substitutions or modifications made by those skilled in the art to the technical solutions and content disclosed in the present invention without departing from the scope of the present invention shall be deemed to have remained within the protection scope of the present invention.

Claims

1. A multi-IMU and plantar pressure synchronous acquisition system based on a master-slave architecture, characterized in that, include: The main control board is connected to five first inertial measurement units via a first bus. The five first inertial measurement units are respectively located at the waist, left thigh, right thigh, left calf, and right calf. The main control board also receives data from the left foot slave plate and the right foot slave plate via a second bus. The left foot follower plate and the right foot follower plate are used to collect quaternion data from the second inertial measurement unit of the ankle joint and pressure data from the plantar pressure sensor array, respectively, and package the collected data into data frames containing frame header, frame sequence number, data content and cyclic redundancy check code and upload them to the main control board. The main control board and the slave board operate independently at a preset sampling frequency. The main control board summarizes the quaternion data and pressure data received from all inertial measurement units into a data packet of a predetermined dimension and uploads it.

2. The system according to claim 1, characterized in that: The slave board acquires pressure data from 45 measuring points in the plantar pressure sensor array via multiplexing scanning; the main control board reads data from the first inertial measurement unit via direct memory access and receives and parses data frames from the slave board via a circular buffer; the preset sampling frequency is 100Hz, and the predetermined data packet has 118 dimensions.

3. A method for measuring lower limb joint angles based on the system described in claim 1 or 2, characterized in that, Includes the following steps: S1: Collect quaternion data and plantar pressure data from all inertial measurement units; S2: Normalize each quaternion q = (q0, q1, q2, q3): Calculate the modulus. ,like Replace with a unit quaternion; otherwise, perform normalization. ; S3: Calculate the relative rotation quaternions of each joint in the lower limb based on the normalized quaternions, where the conjugate of the quaternion is... And convert the relative rotation quaternion into Euler angles; S4: The Euler angle sequence is aligned with the reference signal using a dynamic time warping algorithm, and the aligned sequence is then bandpass filtered. S5: Extract time-domain features and frequency-domain features based on Welch power spectrum estimation from the filtered signal; S6: Input the extracted features into the Stacking ensemble learning model to predict the joint angles. .

4. The method according to claim 3, characterized in that: In step S3, the relative rotation quaternions (q0, q1, q2, q3) are converted into Euler angles in ZYX order. The formulas for rotation α, abduction β, and flexion / extension γ are as follows: 。 5. The method according to claim 3, characterized in that: In step S4, the dynamic time warping algorithm constructs a cost matrix. And backtrack to achieve the optimal path.

6. The method according to claim 3, characterized in that: In step S4, the bandpass filter uses a 4th-order Butterworth bandpass filter with a passband frequency range of 0.5Hz to 5Hz, and a zero-phase filtering method is used to eliminate phase delay.

7. The method according to claim 3, characterized in that: In step S5, the frequency domain features are obtained based on Welch power spectrum estimation. The Welch power spectrum estimation uses a Hanning window with a window length of 256 points, an overlap rate of 50%, and 512 FFT points. The power spectral density is calculated using the following formula: .

8. The method according to claim 3, characterized in that: In step S6, the first base learner of the Stacking ensemble learning model includes a random forest, a gradient boosting tree, and a ridge regression model, and uses 5-fold cross-validation to generate meta-features; the parameters of the random forest are n_estimators=400 and max_depth=9; the parameters of the gradient boosting tree are n_estimators=400 and max_depth=7.

9. The method according to claim 8, characterized in that: The second-layer meta-learner of the Stacking ensemble learning model is a Bayesian ridge regression model, which updates hyperparameters iteratively. Training was conducted, among which .

10. The method according to claim 3, 8 or 9, characterized in that: In step S6, the predicted joint angles Calculated by the following formula: in, , , These are the predicted outputs of the first-layer random forest, gradient boosting tree, and ridge regression model for the input feature x, respectively, where w0, w1, w2, and w3 are the weight coefficients learned by the Bayesian ridge regression model.

Images