An improved method for obtaining rotational accelerations from instrumented headforms

The following compares the effect of differentiation methods used to acquire angular acceleration from three types of un-helmeted headform impact tests. The differentiation methods considered were the commonly used 5-point stencil method and a total variation regularization method. Both methods were used to obtain angular acceleration by differentiating angular velocity measured by three angular rate sensors (gyroscopes), and a reference angular acceleration signal was obtained from an array of nine linear accelerometers (that do not require differentiation to obtain angular acceleration). For each impact, three injury criteria that use angular acceleration as an input were calculated from the three angular acceleration signals. The effect of the differentiation methods were considered by comparing the criteria values obtained from gyroscope data to those obtained from the reference signal. Agreement with reference values was observed to be greater for the TV method when a user-defined tuning parameter was optimized for the impact test and cutoff frequency of each condition, particularly at higher cutoff frequencies. In this case, mean absolute error of the five-point stencil ranged from 1.0 (the same) to 11.4 times larger than that associated with the TV method. When a constant tuning parameter value was used across all impacts and cutoff frequencies considered in this study, the TV method still provided a significant improvement over the 5-point stencil method, achieving mean absolute errors as low as one-tenth that observed for the five-point stencil method.


Introduction
Measures of linear and angular acceleration are central to head impact research. Anthropometric test device headforms are frequently used in a broad range of head impact research, ranging from the design and evaluation of helmets and headgear [1] to investigations of blast-related brain injury [2]. While several headforms have been used to investigate head impact in sport [3][4][5][6][7], the Hybrid III headform is one of the most frequently encountered in the literature [8,9]. The Hybrid III headform, like all others, features rigid elements to which sensors are mounted. Similar methods are also used in cadaveric studies of head impact biomechanics, where instrumented blocks are mounted to the skull [10]. Due to the frequency with which these methods feature in head impact research, accurate measures of rigid body acceleration are important.
Several instrumentation schemes have been developed to observe head or headform kinematics during impact, each with associated advantages and disadvantages. The two most common approaches involve instrumentation with either an array of linear accelerometers or a combination of linear accelerometers and gyroscopes (angular rate sensors). Linear accelerometer arrays frequently feature nine or more sensors in any one of several arrangements, with the most common being the 3-2-2-2 array also commonly referred to as the Nine Accelerometer Package (NAP) [11]. In this arrangement, three sensors are typically mounted at the center of gravity (CG) and the remaining sensors are mounted in biaxial pairs away from the CG (Fig. 1a). Linear accelerometer arrays like the NAP are widely used because of their robustness, reliability and computational stability, which has resulted in its broad application.
The most commonly used arrangement featuring angular rate sensors is the 3aω array, where all six sensing elements are typically mounted at the CG (Fig. 1b). The primary advantage of this instrumentation scheme is that it requires fewer channels, which decreases data acquisition costs and complexity. Further, the reduced complexity of the arrangement facilitates sensor packages with low mass and reduced volume form-factors which are desirable in wearable sensor applications. Additional benefits include better performance in capturing long duration impacts, insensitivity to mounting location and better accuracy in measuring angular velocity compared to linear accelerometer arrays [12,13]. Accuracy of angular velocity measurement is particularly important for wearable sensors where, through relationships induced by assuming rigid body dynamics, the quantity is needed to obtain accelerations at the head center of gravity [14]. Furthermore, increased accuracy in measuring angular velocity means gyroscopes have the benefit of being more accurate when considering criteria like the brain injury criterion (BrIC), which are computed directly from angular velocity [15].
While instrumentation with gyroscopes is simpler than with linear accelerometer arrays, the approach also has limitations. The principal limitation is the difficulty in obtaining accurate angular accelerations from observed angular velocities due to the noise amplification properties inherent in differentiation [12,16]. Several methods have been proposed to mitigate this shortcoming of the angular rate sensors. Filtering of angular velocity measures prior to differentiation improves accuracy generally [17,18], and optimization of channel-specific filtering cutoff frequencies has shown further benefit [13]. While these steps improve the accuracy of the 3aω system, they do not yield a noise-free angular acceleration signal.
Inaccuracies in angular accelerations obtained by differentiation are exacerbated at high sampling frequencies, which is of particular concern when computing head injury criteria sensitive to peak angular acceleration magnitudes or when considering head impacts of short duration. Analysis of cadaveric drop test data has indicated that the peak angular acceleration (PAA), the rotational injury criterion (RIC), the power rotational head injury criterion (PRHIC) and concussion correlate (CC) are the injury criteria most sensitive to the accuracy of angular acceleration measurements [16]. Of these, PAA, RIC and PRHIC were selected for consideration in the current study and computed for each impact because these measures depend exclusively on measures of angular acceleration [19]. As such, they are the most likely to benefit from improved differentiation methods. Likewise, helmets are designed to decrease head acceleration by extending contact durations. As a result, accuracy of kinematic measurements in un-helmeted impacts is likely more sensitive to the characteristics of differentiation methods.
The five-point stencil is a standard numerical technique for computation of derivatives [20] commonly used to compute angular acceleration from angular velocity signals. In the one-dimensional case, this technique may be thought of as an extension of the central difference method which largely preserves the computational simplicity and increases its ability to reduce noise amplification. More sophisticated numerical differentiation techniques have also been developed for analysis of biomechanical systems [21-30], but evaluation with regard to head impact studies is limited. One such method is the total variation (TV) regularization method of Chartrand [31] computes the derivative by solving an optimization problem with an objective function composed of a data fidelity and penalty term. The data fidelity term describes the difference between the original input signal to be differentiated and the signal obtained by integrating the computed derivative. The penalty term characterizes irregularity in the computed derivative. The relative influence of these two terms is determined by tuning a user-defined parameter α. This method has several potential advantages, including reduced noise amplification and the ability to capture discontinuities in the derivative, but requires the selection of a tuning parameter value.
The purpose of this study was to evaluate the performance of the TV differentiation method relative to the standard five-point stencil method in calculating angular acceleration from angular velocity. Un-helmeted impacts were examined in this study because of their shorter duration. Likewise, the PAA, RIC and PRHIC metrics described above were selected due to their sensitivity to angular acceleration signal accuracy.

Data collection
Impact tests were conducted using a 50th percentile male Hybrid III headform. The headform was instrumented with Headform kinematics were measured in two sports ball impact scenarios and one drop scenario (Fig. 2). A total of 36 impacts were conducted, 12 replicates for each scenario. In all impacts, the headform was bare (no helmet or headgear) and attached to a Hybrid III neck. In the sports ball impact scenarios, a solid ball (softball; Worth, Gold Dot 44/375) and inflated ball (soccer ball; Select, Club Size 5) were projected at the headform at speeds selected to be representative of play [32,33], though softball speeds were limited to 20 m/s to minimize risk of damage to the headform (Table 1). In the drop scenario, a NOCSAE compliant drop carriage was used to drop the head and neck assembly onto an 15.24 cm (6 in) thick MEP pad (Southern Impact Research Center, Rockford, TN, USA) from a height of 1.07 ± 0.003 m [34] (Table 1). A high speed video camera recording at 3000 fps was positioned to the side of the headform (normal to the sagittal plane) and video of each impact was used to qualitatively confirm consistent contact location for all impacts.

Data analysis
Accelerometer data were used to directly obtain linear acceleration of the headform at the center of gravity (CG). A "ground truth" angular acceleration signal was obtained from the NAP configuration (9 linear accelerometers) using a standard method [11,35]. This method for computing angular acceleration is not without error, but was selected for use as a reference because it has been extensively studied and is commonly used in investigations of sports-related head impact [14,[36][37][38][39].
Prior to analysis, all linear accelerometer data were filtered using a fourth-order Butterworth filter with a cutoff frequency f c = 1650 Hz in accordance with SAE J211 [40]. Angular velocity data are commonly filtered using lower cutoff frequencies when used to obtain angular acceleration to mitigate noise amplification during differentiation [12,14]. Higher cutoff frequencies are required, however, to obtain accurate head kinematic measures and capture some head injury criteria [16]. To investigate the noise rejection characteristics of the TV differentiation method, angular velocity data collected in this study were filtered using a fourth-order low-pass Butterworth filter and cutoff frequencies ranging from 1 to 2 kHz in increments of 25 Hz. The lower and upper bounds of this range were selected based on the duration of the shortest contacts in the study and the manufacturer-specified frequency response limit of the angular rate sensors, respectively. Kinematic data were processed using MATLAB (Math-Works, Inc., Natick, Massachusetts).
To investigate the performance of the differentiation methods, PAA, RIC and PRHIC were computed for each impact replicate from all three angular acceleration measures (One from the NAP and two from differentiation of angular velocity data). Angular acceleration was obtained from the gyroscopes using two differentiation methods, the five-point stencil and TV method described above, and performance of each differentiation technique was assessed using the mean absolute percentage error (MAPE): where n is the number of impact replicates in each impact scenario, and ŷ i and y i are the injury criterion value for the ith replicate computed from differentiated angular velocity and NAP data, respectively. For both differentiation methods, MAPE was calculated from all replicates in each impact scenario. These calculations were repeated at every cutoff frequency.
Output of the TV differentiation method also depends on α, so for this method MAPE was calculated for 100 logarithmically spaced values of α ranging from 5 × 10 −8 to 500. Thus, it was possible to find α values which yield minimum MAPE for each injury criterion, impact scenario and cutoff frequency. To characterize the maximum benefit that might be achieved through the TV method these locally optimal α values were identified and the corresponding MAPE values reported. To characterize the performance of the TV method when impact types are not known a priori, a globally optimal α value was also identified by computing the mean MAPE across criteria, impact scenarios and cutoff frequencies for each α in the range investigated and identifying the value which achieved minimum overall error.

Results
Softball impacts produced the shortest duration contacts (μ = 4.4 ms), followed by drop impacts (μ = 7.7 ms) and soccer ball impacts (μ = 17.8 ms; Fig. 2). It is not surprising then that the greatest disagreement between angular accelerations obtained by differentiation of angular velocity data were observed at low cutoff frequencies in the softball impact scenario. High-frequency noise was observed, however, in five-point stencil measures of angular acceleration for soccer ball and drop impacts at higher cutoff frequencies.
With the use of a locally optimal α, the TV method output contained substantively less noise at high cutoff frequencies (Fig. 3).
With respect to PAA, the TV method and five-point stencil produced similar results at low cutoff frequencies. The five-point stencil error (compared to ground-truth signal) was sensitive to cutoff frequency, with MAPE increasing at high cutoff frequencies, while the TV method error remained nearly constant across cutoff frequencies (Fig. 4). This was observed when both locally and globally optimal α were used. Utilization of the globally optimal α substantively increased error of the TV method for softball impacts, but had a smaller effect for soccer ball and drop impacts (Fig. 4).
For both RIC and PRHIC, the same trends were observed across all impact scenarios: five-point stencil error increased with cutoff frequency, and TV method yielded lower error and remained constant across cutoff frequencies. In these cases, only marginal increases in error were observed for the TV method when the globally optimal α was used, even in the case of softball impacts (Figs. 5, 6).
Across cutoff frequencies, for the locally optimal case, 5 × 10 −8 < α < 2.997, 1.182 < α < 314.0 and 1.883 < α < 396.2 were found for PAA, RIC and PRHIC, respectively. Contour plots of MAPE for combinations of α and cutoff frequency indicated a "corridor" corresponding to near-optimal performance which extends across cutoff frequencies (Figs. 4, 5,  6). This feature is evident in plots for all criteria and impact scenarios with the exception of PAA for softball impacts. In that case, the locally optimal α appears to lie at the end of the range investigated and MAPE appeared more sensitive to α than cutoff frequency (Figs. 4, 5, 6). The globally optimal α identified from these results, including PAA for softball impacts, was 2.997. Contours of MAPE as a function of angular velocity filter cutoff frequency and α for each impact scenario PRHIC computed using angular velocity data differentiated using the TV method was less sensitive to filter cutoff frequency than the five-point stencil.
Error of the TV method, compared to the ground truth reference signal, was the same or lower than the five-point stencil method for all injury criteria examined in this study (Right) Contours of MAPE as a function of angular velocity filter cutoff frequency and α for each impact scenario when a locally optimal α was used. In particular, the TV method error was substantively lower for the RIC and PRHIC criteria, while gains in accuracy for PAA were more dependent on impact scenario and cutoff frequency. Substantive decreases in error obtained by the TV method for RIC and PRHIC is likely due to the exponential weighting of their respective integrals. In the process of computing these criteria, noise in the gyroscope angular velocity signal is first amplified by differentiation. Then, because these criteria consider angular acceleration magnitudes, the noise is subsequently summed through integration and raised to a power of 2.5. Therefore, these criteria are particularly sensitive to the processing of angular velocity data, as observed here (Figs. 5, 6).
The exception to the trends observed in these results was PAA for softball impacts, where the five-point stencil and TV methods yielded identical results. This is not surprising, given that the signal of interest for the short duration softball impacts are near or exceed the sensor bandwidth. While these results indicate that the TV method has the potential to out-perform the five-point stencil method, they were obtained through the selection of an optimal α for each injury criterion, impact scenario and cutoff frequency. Identification of an α which yields good performance across impact scenarios and injury criteria is particularly important for application in wearable sensors where, compared to laboratory tests, a wider range of impact types may occur and less information is available to inform analysis.
The sensitivity analyses conducted as part of this investigation suggest it may be possible to identify an alpha value which yields good performance. Contour plots of the TV method error against α and cutoff frequency indicate a "corridor" of good performance for many impact scenarios and injury criteria. Moreover, these corridors of near-optimal α extend across the range of cutoff frequencies examined and there appears to be substantive overlap across injury metrics and impact scenarios. Notably, the TV method error remained relatively insensitive to cutoff frequency when the globally optimal α was used. This suggests the method may have application in cases where the impact modality (e.g., ball-to-head, head-to-head, etc.) is unknown a priori, or in cases where multiple impact modalities are anticipated.

Limitations
While this study lays the ground work for selecting an appropriate α, only a limited number of impact types were investigated. Future applications may require similar effort to identify appropriate α for different impact scenarios. Such studies will need to consider both the sensor bandwidth and the metrics of interest. It should be noted that the same signal could be processed multiple times with differing α, one for each metric or impact scenario, to increase accuracy when there are multiple metrics of interest or the impact scenario is known.
The method investigated here possesses several positive qualities, but it does require substantively more computational effort than the five-point stencil method to which it was compared. This limitation is particularly pertinent if multiple α are necessary to increase accuracy for a given application. This may prohibit on-board computation in lowpower wearable devices, but could potentially be performed offline on more robust hardware. It should also be noted that the present study considered cutoff frequencies from 1 to 2 kHz and examined data from sensors rigidly mounted to the headform. Wearable devices would require gyroscope sampling rates in this range to benefit from the findings and the (non-rigid) wearable device attachment interface may affect performance gains.

Conclusion
The five-point stencil and TV methods produced similar results for low cutoff frequencies, but diverged as cutoff frequency increased. Performance of the five-point stencil decreased as the cutoff frequency increased. This result is not surprising, as differentiation amplifies noise and this is the reason that angular rate sensors are not commonly used to obtain angular acceleration in these environments. In contrast, the TV method error was less sensitive to filter cutoff frequency. This suggests the method may be better suited to applications where the full bandwidth of angular rate sensors is required to characterize angular accelerations. Improved accuracy of gyroscope-based measures of angular acceleration has several potential benefits, including the practical benefit of reducing the cost of instrumentation for scientific investigations of head impact. Angular rate sensors also have the benefit of being more accurate when considering criteria, like BrIC, which are computed directly from angular velocity. Furthermore, accurate estimation of angular accelerations from angular velocity data may have application in cases, as with wearable sensors, where instrumentation with linear accelerometer arrays is not practicable.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.