Influence of averaging method on muscle deoxygenation interpretation during repeated‐sprint exercise

Near‐infrared spectroscopy (NIRS) is a common tool used to study oxygen availability and utilization during repeated‐sprint exercise. However, there are inconsistent methods of smoothing and determining peaks and nadirs from the NIRS signal, which make interpretation and comparisons between studies difficult. To examine the effects of averaging method on deoxyhaemoglobin concentration ([HHb]) trends, nine males performed ten 10‐s sprints, with 30 seconds of recovery, and six analysis methods were used for determining peaks and nadirs in the [HHb] signal. First, means were calculated over predetermined windows in the last 5 and 2 seconds of each sprint and recovery period. Second, moving 5‐seconds and 2‐seconds averages were also applied, and peaks/nadirs were determined for each 40‐seconds sprint/recovery cycle. Third, a Butterworth filter was used to smooth the signal, and the resulting signal output was used to determine peaks and nadirs from predetermined time points and a rolling approach. Correlation and residual analysis showed that the Butterworth filter attenuated the “noise” in the signal, while maintaining the integrity of the raw data (r = .9892; mean standardized residual −9.71 × 103 ± 3.80). Means derived from predetermined windows, irrespective of length and data smoothing, underestimated the magnitude of peak and nadir [HHb] compared to a rolling mean approach. Consequently, sprint‐induced metabolic changes (inferred from Δ[HHb]) were underestimated. Based on these results, we suggest using a digital filter to smooth NIRS data, rather than an arithmetic mean, and a rolling approach to determine peaks and nadirs for accurate interpretation of muscle oxygenation trends.


| INTRODUCTION
Near-infrared spectroscopy (NIRS) is a common tool to indirectly measure muscular oxygen availability and microvascular reactivity noninvasively. [1][2][3][4][5] Implementation of NIRS relies on the transparency of human tissue and the light absorbing characteristics of oxy-(O 2 Hb) and deoxyhaemoglobin (HHb) chromophores for the determination of their concentration ([O 2 Hb] and [HHb], respectively) in a localized tissue bed. 6 Changes in [O 2 Hb] and [HHb] reflect the dynamic balance between muscle oxygen (O 2 ) delivery and extraction in the underlying tissue. 7 In continuous exercise where NIRS responses are relatively stable, averages can be calculated over discrete and predetermined time points for identification of overall trends within the exercise bout. 1,2 When maximal sprint efforts are repeated, however, there is a rapid deoxygenation at exercise onset that slowly recovers at sprint cessation. [8][9][10] The evolution of peaks and nadirs across the NIRS signal is often used to describe the quality of metabolic recovery between sprint bouts. 9,11,12 Because of the rapid oxygenation adjustments and short duty cycle of repeated-sprint exercise, 10,12,13 accurate identification of peaks and nadirs in the NIRS signal is critical.
Analysis of NIRS data obtained during repeated-sprint exercise is often constrained to [HHb] 8,11,[13][14][15] due to Δ[O 2 Hb] being influenced by rapid blood volume and perfusion variations caused by forceful muscle contractions. 3,16 Additionally, the HHb signal is considered to be relatively independent of blood volume 2,3 and taken to reflect venous [HHb] which provides an estimate of muscular oxygen extraction. 1,2 However, across studies, there are differing methods used to smooth the NIRS signal and determine peak and nadir [HHb], which can potentially affect comparisons between studies and, therefore, interpretation.
To analyze a NIRS signal, single values for each sprint and recovery are typically determined for each peak and nadir. 8,10,13,17,18 A mean is calculated over a predetermined duration within the closing seconds of each sprint and recovery periods to smooth fluctuations in raw NIRS data during sprint exercise. [8][9][10][19][20][21][22] This method has been used on numerous occasions in acute settings, 9,13 varying inspired O 2 fraction, 11,[20][21][22] active vs passive rest, 8,23 after respiratory muscle warm-up, 24 and in response to training. 14,17,25 However, a possible drawback is that the true, physiological peak and/or nadir [HHb] may not fall within the predefined analysis window. It may be that [HHb] continues to rise if tissue O 2 consumption remains elevated post sprint and/or if O 2 delivery decreases. Additionally, the recovery nadir may be affected by limb activity when the athlete prepares for the next sprint (ie, leg movement to place the pedal in the right position and static contraction of the quadriceps). To overcome this, a rolling mean approach may be applied to smooth the data to determine the true peak and nadir of the NIRS signal. 10,12,15,18,19,26 But currently, there is no comparison of means calculated from predetermined time periods or a rolling mean approach. Additionally, there is no consistency of the moving average window duration, which may be constrained to sprint duration. 11,12 A digital filter is another typical technique used to attenuate noise and smooth raw data. 27 For example, when a low-pass filter is used, a cutoff frequency is chosen so that lower signal frequencies remain and higher frequencies are attenuated. 28 Such filters have been employed to smooth the NIRS signal during repeated-sprint exercise, 9,18,29 but again, the relevance of such technique has yet to be confirmed compared to more widely used averaging methods.
Therefore, the purpose of this study was to compare and evaluate the effect of different NIRS signal analysis methods (predetermined temporal window, rolling mean, and Butterworth filter) on muscle tissue oxygenation trends during a repeatedsprint protocol. We propose that the combination of a digital filter to smooth the NIRS signal and the identification of a local maximum and minimum for each sprint/recovery phase will improve our ability to detect changes in the signal.

| Participants
Nine males accustomed to high-intensity activity were recruited for this study (mean ± SD: age 25 ± 3 years; height, 183.2 ± 7.7 cm; body mass, 81.0 ± 8.7 kg; V ̇O 2peak , 54.6 ± 6.2 mL·min -1 ·Kg -1 ). All participants reported to be healthy and with no known neurological, cardiovascular, or respiratory diseases. After being fully informed of the requirements, benefits, and risks associated with participation, each participant gave written consent. Ethical approval for the study was obtained from the institutional Human Research Ethics Committee and conformed to the declaration of Helsinki.

| Experiment design
Participants were part of a larger project that required six separate laboratory visits. Data presented here were taken from the control trial that took place after familiarization. Testing was performed on an electronically braked cycle ergometer (Excalibur, Lode, Groningen, The Netherlands), set to "isokinetic" mode. In this mode, a variable resistance is applied to the flywheel proportional to the torque produced by the subjects to constrain their pedaling rate to 120 rpm. Below 120 rpm, no resistance is applied to the flywheel. This mode was chosen to avoid cadence-induced changes in mechanical power production, 30 and hemodynamics, 31 within and between sprints. After a 7-minute warm-up consisting of 5 min of unloaded cycling at 60-70 rpm and two 4-seconds sprints (separated by 1 minute each), participants rested for another 2.5 minute before the repeated-sprint protocol was initiated. The repeated-sprint protocol consisted of ten 10-seconds self-paced sprints separated by 30 seconds of passive rest. Participants were instructed to give an "all-out" effort for every sprint and verbally encouraged throughout to promote a maximal effort. Each sprint was performed in the seated position and initiated with the crank arm of the dominant leg at 45°. Before sprint one, subjects were instructed to accelerate the flywheel to 95 rpm over a 15-seconds period and assume the ready position 5 seconds before the commencement of the test. This ensured that each sprint was initiated with the flywheel rotating at ~90 rpm so that subjects could quickly reach 120 rpm. Five seconds prior to the initiation of each sprint, participants were asked to assume the ready position, followed by a verbal 3-2-1 countdown.

| NIRS
Participants were instrumented with NIRS probes (Oxymon MKIII, Artinis, the Netherland) fixed over the distal part of the vastus lateralis muscle belly of their dominant leg, approximately 10-15 cm above the proximal border of the patella and held in place with plastic spacers with an optode distance of 4 cm. Skinfold thickness was measured between the emitter and detector using a skinfold caliper (Harpenden Ltd.) to account for skin and adipose tissue thickness covering the muscle. The skinfold thickness (12.4 ± 6.9 mm) was less than half the distance between the emitter and the detector in each case. Probes were attached using double-sided stick disks, secured with tape, and shielded from light with black elastic bandages. Between sprints, participants were asked to minimize leg movement by remaining seated and relax their dominant leg in the extended position. A modified form of the Beer-Lambert law was used to calculate micromolar changes in tissue [HHb] across time using received optical density from one continuous wavelength of NIR light (763 nm). A differential pathlength factor of 4.95 was used. 20,21 Data were acquired at 10 Hz and exported to Excel for analysis. These data were expressed as a percentage so that resting baseline represented 0% and maximal [HHb] represented 100% (Δ%[HHb]). Maximal [HHb] was obtained with femoral arterial occlusion using a pneumatic tourniquet (inflated to 300-350 mm Hg) around the root of the thigh for 3-5 minutes until the [HHb] increase reached a plateau. Arterial occlusion was performed after the completion of the repeated-sprint protocol (within 10 minutes), while the subjects lay on an examination bed with the leg under examination at 90° knee flexion, and foot on the bed. During the trials, markers were placed in the NIRS software at sprint onset to demarcate the 40-s sprint/recovery windows for analysis.
The application of the 10th order zero-lag low-pass Butterworth filter was conducted in the R environment 32 using the signal package. 33 The filter order was determined based on previous research, 9,18 and the effects of filter order on the sharpness of filter response. The filters cutoff frequency (ƒ c ) was determined based on a combination of previous research, 18 residual analysis (of data from three subjects) of the effects of a range of different normalized ƒ c on HHb (Figure 1), and visual inspection with attention paid to local maxima and minima of filtered data compared to the raw signal. 34 Based on these, it was concluded that 0.1 was suitable ƒ c to be applied to the data for the remaining subjects. After the filter passed through the data, the resulting output was exported to Excel for standardization to occlusion values and determination of peaks and nadirs.

3.
Moving average with a window of 2 seconds applied to the data, followed by the identification peaks and nadirs within each 40-second exercise-recovery cycle: 2 MA . 4. Moving average with a window of 5 seconds applied to the data, followed by the identification peaks and nadirs within each 40-second exercise-recovery cycle: 5 MA 5. Application of a Butterworth filter to smooth the raw NIRS data, followed by the identification of peaks and nadirs from predetermined time points. A single value prior to each phase change (ie, end of exercise and end of recovery, 0.1 second): BWF PD. 6. The application of a Butterworth filter to smooth the data, followed by the identification of a peak and nadir using a rolling approach within each 40-seconds exercise-recovery cycle (0.1 second): BWF MA .
Tissue reoxygenation (ΔReoxy) was calculated as the difference between the peak and nadir for each analysis method.

| Statistical analyses
Data in text and figures are presented as mean ± SD. Relative changes (%) are expressed with 95% confidence limits (95% CL). Effects of the Butterworth filter on the NIRS signal was assessed by calculating Pearson's product-moment correlation (r), and standardized residuals of the raw vs filtered data in the R environment using the stats package. 32 The correlation between the raw and filtered NIRS signal was assessed by fitting a linear regression model to the pooled subject data. The following criteria were adopted to interpret the magnitude of the correlation between variables: ≥0.1, trivial; >0.1-0.3, small; >0.3-0.5, moderate; >0.5-0.7, large; >0.7-0.9, very large; and >0.9-1.0, almost perfect. 35 To determine the effects F I G U R E 1 A plot of the root-mean-square (RMS) residuals between filtered and unfiltered signals as a function of the filter cutoff frequency from the data of a representative subject. A line of best fit (ab) is projected to the Y-axis. At the intercept c, the horizontal line cd is drawn to intersect with the residuals. The chosen cutoff frequency ƒ c is at this point of intersection of analysis method, practical significance was also assessed by standardized effects and presented with 95% CL. 36 Effect sizes (ES) between0-< 0.2, >0.2-0.5, >0.5-0.8, and >0.8 were considered to as trivial, small, moderate, and large, respectively. Probabilities were also calculated to establish whether the chance the true (unknown) differences were lower, similar, or higher than the smallest worthwhile change (0.2 multiplied by the between-subject SD, based on Cohen's effect size principle). Quantitative probability of lower, similar, or higher differences was assessed qualitatively as follows: <1%, almost certainly not; >1%-5%, very unlikely; >5%-25%, unlikely; >25%-75%, possible; >75%-95%, likely; >95%-99%, very likely; and >99%, almost certainly. If the probability of having higher/lower values than the smallest worthwhile difference was both >5%, the true difference was assessed as unclear. 35,37 Data analysis was performed using a modified statistical Excel spreadsheet. 38 To examine the interaction effects between the method for identifying Δ%[HHb] peak, nadir and ΔReoxy (predetermined and moving), and the size of the analysis window (0.1, 2, and 5 seconds), two-way repeated measure ANOVAs were performed. Post hoc analysis was conducted using the Holm-Šídák method and adjusted for multiple comparisons. A threshold for significance was set at the P < .05 level. Analysis was performed GraphPad Prism 6.

| Application of the Butterworth filter
An example of the raw compared to the filtered data of a representative subject is presented in Figure 2. There was an almost perfect Pearson's correlation between the raw NIRS data and that after the Butterworth filter ( Figure 3A). The mean standardized residual of the raw data compared to the filtered was -9.71 × 10 3 ± 3.80 ( Figure 3B). When rectified, the mean residual was 2.51 ± 2.86 with a relative difference of 2.5% [CL: 1.7, 3.4].

| Peak muscle [HHb]
Mean results of the different analysis methods are presented in Figure 4A. Comparisons of analysis methods are shown in Table 1. There was a significant effect of the method for identify peaks on peak muscle Δ%[HHb] at the P < 0.  [11.7, 19.1]; P < .0001). There was also a likely small difference between 2 MA and 2 PD (8.2% [5.4, 11.0]; P < .0001). An almost certainly small effect was also observed when 2 PD was compared to 5 PD . Differences between 2 MA and 5 MA were almost certainly trivial. Means determined from BWF PD were almost certainly higher than 5 PD and almost certainly trivial compared to 2 PD . When the results from BWF MA were compared to other moving averages, BWF MA was almost certainly higher than 5 MA

| Nadir muscle [HHb]
Mean results of the different analysis methods are presented in Figure 4B. Comparisons of analysis methods are shown in Table 2. A significant effect of the method for identify nadirs on peak muscle Δ%

| DISCUSSION
Results indicate that using predefined averaging windows to analyze the NIRS signal within sprint and recovery periods, Δ%[HHb] peaks and nadirs are consistently underestimated compared to a moving average, regardless of the window length. Subsequently, muscle reoxygenation between efforts is underestimated using 5 PD and 2 PD compared to 5 MA and 2 MA . However, a drawback of using the arithmetic mean is that each data point contributes equally to the average, which allows outlying data points to bias the result. 39 To overcome this, we applied a 10th order zero-lag low-pass Butterworth filter to the NIRS signal, which incorporates a weighted mean from several data points across the signal. Correlation and residual analysis revealed the Butterworth filter attenuated the "noise" yet maintained the integrity of the raw data ( Figure 3). As NIRS responses are used as a surrogate for metabolic perturbations, detecting the magnitude of change is critical for assessing the influence of interventions and environmental factors. Thus, it appears that a digital filter combined with a rolling approach for determining peaks  and nadirs of the NIRS signal is the best method for accurate interpretation of oxygenation trends in repeated-sprint exercise. There were clear differences in peak Δ%[HHb] means between the predefined averaging methods. The difference between 2 PD and 5 PD can be attributed to the length of the averaging window. At sprint onset, there is a sharp rise in Δ%[HHb] from rest that peaks at sprint cessation or shortly thereafter (Figure 2 and Ref. [8,10,12,13,18]). In our representative data, Δ%[HHb] continued to increase ~40% in the final 5 s of the first sprint which led to a significant 10% reduction in peak Δ%[HHb] for 5 PD compared with 2 PD (Table 1). Consequently, an averaging duration that encapsulates a sharp rise within the window length will underestimate the amplitude of Δ%[HHb] change.
It is clear from Figure 2 that Δ%[HHb] continued to increase after the predefined averaging window. This would tend to underestimate the peak Δ%[HHb] response. To overcome this potential confounding factor, we applied both a 5-second and 2-second moving average to the NIRS signal. The 2 MA and 5 MA yielded 15% and 8% greater peaks compared to 2 PD and 5 PD , respectively. As maximal Δ%[HHb] may occur either at, or immediately after sprint cessation, identification of peaks using a moving average will always capture this maximal deoxygenation. Therefore, studies that only employed predetermined averaging windows 9,11,17,[20][21][22] may have underestimated the true magnitude of Δ%[HHb] change induced by repeated-sprint exercise. To more accurately represent muscle oxygenation changes in response to repeated-sprint exercise, values need to be determined from a moving average approach. 8,10,12,18,19,29 However, a moving average may not be appropriate when exercise protocols incorporate other muscular activity around repeatedsprint bouts. For example, some authors have used low-intensity exercise during recovery periods between sprints, 8,10,12 and others have used running-based protocols, which impose eccentric loading during the negative acceleration phase post-sprint. 8,14,17 Another limitation is that when using arithmetic means (ie, 5 PD and 2 PD ; 5 MA and 2 MA ), equal weight is given to all data points, which can lead to the result being distorted by outliers. 39 But in the field of repeated-sprint exercise, an arithmetic mean is commonly employed to smooth perturbations in NIRS data.
To address this methodological limitation, we used a Butterworth filter to smooth the data and obtained Δ%[HHb] peak from both predetermined time points within each sprint and from a rolling approach. The BWF PD approach yielded 13 and 3% greater peaks than both 5 PD and 2 PD methods, respectively. As BWF PD used the last value within each sprint (in the current instance, the 100th value of a 10-second sprint sampled at 10 Hz), this was the highest Δ%[HHb] value achieved within each 10-second sprint period. Similarly, peak Δ%[HHb] determined by BWF MA was 19% greater than 5 MA . However, only a trivial difference was found between BWF MA filtering and 2 MA . Various studies have applied digital filters to smooth biomechanical and biological data, 34,40,41 and yet few authors have used a filter to smooth a NIRS signal during repeated-sprint exercise. 9,18 When a low-pass filter is used, a ƒ c is chosen so that lower signal frequencies remain and higher frequencies (noise) are attenuated. 28 A low-pass Butterworth filter attenuates signal power above a specified ƒ c , but also included a weighted average across several data points, 27 which leads to lag in the signal output. This temporal shift can be removed by running the filter a second time in the reverse direction (zero-lag). Repeated-sprint exercise represents a particularly salient challenge to Δ%[HHb] signal due to the rapid changes in duty cycle. Our results suggest that an arithmetic mean under-represents [HHb] peak in most cases and that both BWF PD and BWF MA better reflect peak [HHb] after sprint exercise.
Differences in nadir Δ% [HHb] were less clear between the averaging methods. As means were calculated on the flatter portion of the NIRS signal during the late stage of recovery, differences between averaging methods were minimal. There was an 8% and 20% lower Δ%[HHb] nadir from both a 5 MA and 2 MA compared to 5 PD and 2 PD, respectively ( Table 2). While the difference between nadir 5 MA and 5 PD did not reach the typically adopted threshold for statistical significance of P < .05 ("trivial" standardized effect), we reported a small standardized effect (20% relative difference, P < .0001) between the 2 MA and 2 PD analysis methods. The 5-s averaging windows may not have the sensitivity to detect subtle differences between the nadir determined from predetermined and moving averages. It T A B L E 3 Comparison of smoothing method responses on ΔReoxy appears that a short moving average is better suited at detecting changes in nadir Δ% [HHb]. As the restoration of NIRS variable toward baseline has become a surrogate for metabolic recovery between sprints, 8,10,11,15,17 an accurate depiction of this variable is necessary to assess metabolic perturbations leading to greater peripheral fatigue. Additionally, the detection of the magnitude of change has important implications for assessing the potency of training programs and environmental factors. Unless the nadir of the signal is obtained from a rolling approach, the magnitude of [HHb] recovery will be underestimated. The ΔReoxy is determined from both peak and nadir Δ%[HHb] responses, and hence, the analysis method that yields the greatest peaks and nadirs will have the greatest ΔReoxy. Consequently, we observed clear and substantial differences between all ΔReoxy analysis methods apart from BWF MA vs 2 MA , with a range of relative differences from 5% to 31%. Studies reporting ΔReoxy from predetermined windows 11,17 have likely underestimated reoxygenation capacity. The most accurate representation of ΔReoxy would come from works that have chosen a moving average approach for the determination of peaks and nadirs. 12,18 Irrespective of the method, a narrow window for analysis (2 seconds) allows reporting greater magnitudes of response than a longer window (5 seconds). However, care should be taken when using a narrow analysis window. It contains less data from which the average is calculated, which increases the risk that outliers bias the calculated mean. 39 Such a methodological pitfall is displayed in the sprint phase of Figure 2A, where a small number of data points are far above the characteristics of the surrounding data. A narrow analysis window would place greater weight on these individual data points in the mean compared to a larger window. 39 Currently, there is no consistency on the length of the analysis window. In some cases, a window as large as 5 seconds 8,11,21,22 and as small as 1 seconds 10,12 has been used. When a Butterworth filter was employed, both peak and nadir values 18 and a 2-second average 9 have been used. However, using a Butterworth filter, a single data point from the resulting output can be used with assurance that it reflects the characteristics of the surrounding data. Although our choice to use a Butterworth filter was based on previous research, 9,18,29 other smoothing/filtering techniques which eliminate outliers may also yield similar results. Readers should also be aware that these data and analytical methods were collected during isokinetic sprints where cadence was constrained to 120 rpm, and, although muscle oxygenation patterns appear similar, one may exert caution when analyzing NIRS signal in nonisokinetic conditions where cadence is influenced by gear ratio and neuromuscular fatigue.

| PERSPECTIVE
NIRS-derived variables in sprint exercise are subject to rapid and large perturbations. Furthermore, during cyclic movements such as running or cycling, there are relatively large oscillations in the [HHb] signal response due to mechanical effects of muscle contraction on local blood flow. This requires appropriate smoothing of the signal to avoid either over or underestimation of peaks and nadirs. Sprint and recovery means calculated over a 1-5 seconds window in predetermined time frames are often reported, 8,9,11,13,14,17,[20][21][22][23][24][25] however, there remains little consistency between studies examining muscle oxygenation during repeated-sprint exercise. The current results reveal that moving averages derive greater changes in muscle oxygenation than means calculated from predetermined time points. Hence, this method is less prone to underestimation of the maximum rate of de-and reoxygenation. However, as these calculations are susceptible to signal bias from outliers, especially when a shorter averaging window is used, 39 we recommend using a digital filter 9,18,29 or other smoothing/filtering techniques prior to analysis. We also suggest that future studies should avoid predetermined analysis windows and focus on the determination of a single value for peaks and nadirs from a rolling approach.