## 2 Materials and Methods
### fNIRS Data Analysis  
Signal Processing was performed using Homer3 software [62], and it consisted of four main steps: (1) filtering out channels with excessive noise according to their scale-coupling index (SCI); (2) removal of motion artifacts and unwanted physiological signals, (3) conversion of data into oxygenated hemoglobin (HbO) and deoxygenated hemoglobin (HbR) concentration changes, and (4) making inferences on brain activity from the hemodynamic response.  
Step 1--The quality of the data was determined based on the SCI, which is a measure of how well the source and detector optodes for a given channel are coupled to the skin [30; 63]. To determine the SCI, a bandpass filter (0.7-1.5 Hz) was used to detect the frequency range corresponding to blood volume pulsations related to the heartbeat for each of the two NIR wavelengths. Then, the correlation between the two filtered NIR wavelengths in each channel was calculated based on the assumption that good optode-skin contact will lead to good correlation of the heartbeat signals that would be present in both NIR wavelengths [64]. Channels with absolute correlation values less than 0.6 [32] were excluded from further analysis. Table 3 details the number of participants excluded for each channel, divided by group. Note that more channels were excluded from the CI group than from the two NH groups because some optodes were above the processor coil (see, for example, the location of the striped circle in Figure 3).  
Step 2--After channel exclusion, preprocessing of the raw data included: (a) conversion of the data from each channel (light intensity) to optical density [62], (b) removal of high-frequency physiological noises, including the cardiac response and respiration in the hemodynamic response using a low-pass filter of 0.5 Hz, (c) exclusion of motion artifacts in each channel by performing a cubic spline correction of the artifacts using the algorithm described by [PERSON] et al. [65], followed by wavelet analysis [47]. The latter analysis performs wavelet transformation of the optical density data and computes the distribution of wavelet coefficients using an algorithm that follows the procedure described by [PERSON] and [PERSON] [66]. This method has been shown to effectively diminish motion artifacts during experiments with speech tasks [67].  
Step 3--The processed optical density data were converted to the concentration changes of HbO and HbR, i.e., the hemodynamic response, using the Modified Beer-Lambert Law [68].  
Step 4--To compare the quality of fit between the predicted and actual hemodynamic response curves, an ordinary least squares GLM was applied to the processed fNIRS signal [69]. The GLM matrix included a modified gamma function applied to both HbO and HbR: \([\gamma(t,\tau,\sigma)=(t-\gamma^{2}/\sigma^{2})^{*}\exp{[1-(t-\tau)^{2}/\sigma^ {2}]}]\), convolved with a square-wave of  
\begin{table}
\begin{tabular}{c c c c} \hline \hline
**Group/Channel** & **NH** & **NHV** & **CI** \\ \hline  
\hline \end{tabular}
\end{table}
Table 3: The number of participants who had SCI \(<\) 0.6 [32], and thus were excluded for each channel, divided by group. SCI = scalp coupling index. CI = cochlear implant (n = 14), NH = normal hearing (n = 15), NHV = Normal hearing participants who listened via a vocoder (n = 15).  
duration (T). Base values for \(\tau\) and \(\sigma\) were obtained from SPM (e.g., [70]). Simultaneous drift regression using 3 rd order polynomial drift correction and regression of: (a) the most correlated short-distance separation channels (SSC) and (b) auxiliary measurements from the accelerometer was conducted to further remove systemic physiological noise and motion artifacts. The beta weights of the canonical HRF term were extracted for each stimulation condition, measurement channel, and participant.