A Novel Computationally Efficient Approach for Exploring Neural Entrainment to Continuous Speech Stimuli Incorporating Cross-Correlation

Abstract

A linear system identification technique has been widely used to track neural entrainment in response to continuous speech stimuli. Although the approach of the standard regularization method using ridge regression provides a straightforward solution to estimate and interpret neural responses to continuous speech stimuli, inconsistent results and costly computational processes can arise due to the need for parameter tuning. We developed a novel approach to the system identification method called the detrended cross-correlation function, which aims to map stimulus features to neural responses using the reverse correlation and derivative of convolution. This non-parametric (i.e., no need for parametric tuning) approach can maintain consistent results. Moreover, it provides a computationally efficient training process compared to the conventional method of ridge regression. The detrended cross-correlation function correctly captures the temporal response function to speech envelope and the spectral–temporal receptive field to speech spectrogram in univariate and multivariate forward models, respectively. The suggested model also provides more efficient computation compared to the ridge regression to process electroencephalography (EEG) signals. In conclusion, we suggest that the detrended cross-correlation function can be comparably used to investigate continuous speech- (or sound-) evoked EEG signals.

1. Introduction

Brain processes in response to continuous speech have been investigated through the acquisition of brain electrical activity. Electroencephalography (EEG) has facilitated the understanding of underlying brain activity with high temporal resolution [1,2,3]. Event-related potential (ERP), which is generally computed by averaging EEG signals over trials, has been widely used to extract time-locked neural activities evoked by short-duration speech tokens (e.g., words, syllables, and phonemes) [4,5]. However, using short tokens as stimuli limits our understanding of how the brain processes continuous speech [6].

System identification (SI) is a more direct method to investigate neural entrainment to continuous stimuli by modeling how stimulus features are mapped to neural responses [7]. The SI model, which considers the human brain a linear time-invariant (LTI) system, has been used to understand brain functions in response to specific visual and auditory stimuli [8,9,10,11]. In these studies, reverse correlation was commonly used to characterize the temporal response function (TRF) in an LTI system [12]. This approach, which implements reverse correlation by computing the cross-correlation between the input and output of the system, has been employed to investigate neural entrainment to continuous speech stimuli [13,14,15] as well as to continuous visual stimuli [16,17]. Although this approach has proven able to investigate neural encoding in response to continuous stimuli, it is more appropriate for stochastic process stimuli (e.g., Gaussian white noise) than for slow-modulating and non-stationary stimuli, such as continuous speech [18].

However, natural stimuli for humans rarely exhibit the characteristics of stochastic processes. To address cases involving modeling of the TRF to slow-modulating and non-stationary stimuli, normalized reverse correlation has been considered, which effectively estimates the spectral–temporal receptive field (STRF) of brain neurons for more complex stimuli [19]. Although this normalized reverse correlation can map slow-modulating and non-stationary stimuli to neural responses, problems of matrix inversion and overfitting due to noise along the low-variance dimensions of datasets remain [20]. Regularized linear regression (i.e., ridge regression), developed initially to investigate neural responses to pure-tone stimuli in the primary auditory cortex, has been employed to address these problems [21]. However, incorrect selection of the regularization parameter of ridge regression can result in an inaccurate model [22]. Thus, the intensive optimization of the regularization parameter should be performed with ridge regression. To overcome this inefficiency, in this study, we proposed a non-parametric and computationally efficient approach to model linearization between stimuli and neural responses. The estimated model was developed using a detrended cross-correlation function. The computational efficiency while analyzing EEG data in response to continuous speech stimuli was compared between our approach and the ridge regression method.

2. Materials and Methods

2.1. Experiments

Twenty Korean native speakers (mean age 23.4 ± 3.0 years; 11 males and 9 females) with normal hearing participated in this study. Ten continuous Korean sentences from the Korean Standard Sentence List [23] were used as speech stimuli; the mean sentence duration was 1.9 s, with a standard deviation of 0.2. Each sentence was randomly presented 100 times via a loudspeaker placed 1 m away from the participant at a level of 65 dB SPL. Participants gave written consent before the experiment. All the experimental procedures were approved by the Institutional Review Board of the University of Ulsan.

2.2. EEG Acquisition and Preprocessing

Brain activity in response to continuous speech sentences was recorded using a 64-channel EEG system at a sampling rate of 2048 Hz (Biosemi Co., Amsterdam, The Netherlands). The EEG data were preprocessed using average-referencing and 1–57 Hz band-pass filtering and then resampled at a rate of 256 Hz. Each epoch had a −0.5 to 2.5 s interval relative to the onset of the spoken sentence. Epochs and channels with large-scale artifacts were rejected, and specific channel data were spherically interpolated using nearby channels. Independent component analysis was performed to remove eye blink-related artifacts using EEGLAB [24]. Finally, the data were filtered using a zero-phase anti-aliasing 10-Hz low-pass filter after baseline correction.

2.3. Speech Feature Extraction

The broadband speech envelope and speech spectrogram were used as stimulus features to compute their correlation with the EEG data. The speech envelope extracted using the Hilbert transform was downsampled at a sampling rate of EEG data, at 256 Hz. The spectrogram was computed using a filter bank with a series of 64 band-pass filters with logarithmically spaced bandwidths between 250 and 8000 Hz [25].

2.4. Detrended Cross-Correlation and EEG Analysis

2.4.1. Detrended Cross-Correlation Computation

The detrended cross-correlation approach was originally developed based on the reverse correlation technique and the convolution of the derivatives [12]. The derivative is most commonly used as a pre-whitening method to reduce the effect of spurious correlation caused by trends and internal oscillations in time-series inputs [26,27,28]. Detrended cross-correlation analysis was used to estimate the impulse response function of the LTI system. Here, it assumes that a brain can be modeled as an LTI system, which is characterized by the impulse response of $h\left(t\right)$. Then, the output $y\left(t\right)$ of EEG signals in response to speech envelope $x\left(t\right)$ can be described as follows:

$$x\left(t\right) \ast h\left(t\right)=\sum _{\tau }x\left(t-\tau \right)h\left(\tau \right)=y\left(t\right),$$

(1)

where $x\left(t\right)$ is the system input, and $h\left(t\right)$ is the impulse response function. Equation (1) yields the derivative of the convolution as

$$\nabla x\left(t\right) \ast h\left(t\right)=\sum _{\tau }\nabla x\left(t-\tau \right)h\left(\tau \right)=\nabla y\left(t\right),$$

(2)

where $\nabla x\left(t\right)=x\left(t\right)-x\left(t-1\right)$,$\nabla y\left(t\right)=\mathrm{y}\left(t\right)-y\left(t-1\right)$ (Appendix A). From Equation (1), if the input $x\left(t\right)$ is the white noise, the impulse response function $h\left(t\right)$ can be measured by cross-correlating the output $y\left(t\right)$ with the input $x\left(t\right)$ [12]. The impulse response function $h\left(t\right)$ can be estimated using the detrended signals of $x\left(t\right)$ and $y\left(t\right)$, i.e.,$\nabla x\left(t\right)$ and $\nabla \mathrm{y}\left(t\right)$, respectively, according to Equation (2). The detrended cross-correlation function, $C\left(\tau \right)$, between $x\left(t\right)$ and $\mathrm{y}\left(t\right)$ is obtained as

$$C\left(\tau \right)=<\nabla x\left(t\right)\nabla y(t-\tau )>.$$

(3)

We examined the autocorrelation of the original and differentiated speech envelopes to assess the similarity between the autocorrelation and impulse functions. When the autocorrelation of the differentiated input signal converges on the impulse function, the detrended cross-correlation function can be computed as

$$C\left(\tau \right)=<\nabla x\left(t\right)\nabla y(t-\tau )>=Kh\left(t\right),$$

(4)

where $K$ represents the power spectral density, which is computed by dividing the root mean square of the input $\nabla x\left(t\right)$ by the spectral bandwidth [29]. Although the derivation action can generally cause noise in the signals, the cross-correlation between the input and output derivatives effectively mitigate noise in the signals. Figure 1 shows examples of speech envelopes of sentence stimuli, the autocorrelation of speech envelopes, and the autocorrelation of the corresponding differentiated speech envelopes. The shape of the autocorrelation of the differentiated speech envelope more closely resembles an impulse function than that of the original speech envelope, indicating that the differentiated speech envelope is more suitable for identifying the brain system.

2.4.2. Forward Model Computation

The detrended cross-correlation function was computed between the brain response and the differentiated speech envelope to predict the brain response and investigate brain mechanisms under speech stimuli processing. This study used a detrended cross-correlation function to predict EEG signals in response to sentence speech stimuli. The measured EEG data, denoted as $r(t,n)$, were acquired at a specific channel $n$ in response to continuous speech represented by $s\left(t\right)$. Similar to Equation (1), the measured EEG data can be modeled as follows:

$$r(t,n) \propto h(t,n) \ast s(t)+ϵ(n)={\sum }_{\tau }h(\tau ,n)s(t-\tau ) +ϵ(n).$$

(5)

where $h(t,n)$ represents the TRF of the brain, $ϵ\left(n\right)$ is the noise process, and $\tau$ corresponds to the lag between the stimulation and neural response. The TRF of the brain was estimated using the following formula:

$$h\left(\tau ,n\right)\propto <\nabla r\left(\tau ,n\right)\nabla s\left(t-\tau \right).$$

(6)

In this study, the trained TRF was computed over a lag range of −100 ms to 500 ms, which is relevant to the stimulus onset. The forward model (i.e., the encoding model) was obtained through ridge regression using the mTRF toolbox implemented in MATLAB as the reference model [18]. The optimal alpha value for ridge regression was selected within the range of 10⁻⁶ to 10⁶. The baseline-normalized TRF was computed as follows:

$$h\left(\tau ,n\right)=\frac{h\left(\tau ,n\right)-\overline{h_b\left(n\right)}}{\sigma \left(h_b\left(n\right)\right)},$$

(7)

where $\overline{h_b\left(n\right)}$ and $\sigma (h_b\left(n\right))$ represent the mean and standard deviation of TRF during the baseline period of −250 to 0 ms before stimulus onset.

The univariate model uses a speech envelope as the system input to estimate the TRF of the brain in response to speech stimuli. The multivariate model was computed using the TRFs of 64 frequency bands to approximate the STRF of the brain-to-speech stimuli. The multivariate forward model aims to determine the STRFs for a speech spectrogram ranging from 250 to 6000 Hz.

2.4.3. Model Training and Performance Evaluation

An individual model was constructed by averaging across $n$ sweeps for each subject, in which $n-1$ subsets were used for training. The remaining subset, independent of the training set, was used to test the trained model. The model was validated using a leave-one-out cross-validation approach [30]. The Pearson correlation between the measured EEG data and the predicted EEG data over 64 channels was used to evaluate the model performance. Subsequently, the grand-averaged model was computed by averaging the models’ weights for all 20 subjects.

2.4.4. Model Significance Evaluation

The models were computed using lag intervals of −100 to 500 ms and −250 to 0 ms, relevant to speech onset, as the trained and baseline models, respectively. A one-tailed paired $t$-test was used to evaluate the statistical significance of the trained models of the detrended cross-correlation function and ridge regression against the baseline models.

2.4.5. Comparisons of Computational Efficiency and Cross-Validation Performance

In addition to comparing the model performance of the two methods, we also compared the two methods for computational efficiency. The model performance was evaluated using the Pearson correlation coefficient between the predicted EEG and the original EEG signals. The computational efficiency was measured as the elapsed time required to train each model. A two-tailed paired $t$-test was used to evaluate the differences between the detrended cross-correlation function and ridge regression in terms of computational efficiency and model performance.

3. Results

3.1. Performance of the Univariate Forward Model

Figure 2A illustrates the TRFs estimated using the detrended cross-correlation function (left panel) and ridge regression (right panel) for a typical Cz channel. The P1-N1-P2 peaks in the TRFs were observed at 75, 120, and 210 ms, respectively. The Pearson correlation coefficient between the TRFs of the two methods is 0.941, indicating a high level of similarity between the results obtained from the detrended cross-correlation function and ridge regression. Figure 2B shows an example of EEG prediction at Cz and the prediction accuracy of the 64 channels in the topology format for both approaches. Both the detrended cross-correlation model and the ridge regression model performed well in the central and temporal areas, with Pearson’s r values of 0.308 and 0.300, respectively, for the Cz channel (Figure 2C).

3.2. Performance of the Multivariate Forward Model

The top portion of Figure 3A depicts the speech STRF for channel Cz, calculated using detrended cross-correlation and ridge regression. These two methods produced highly correlated STRFs in the lower frequency range (up to approximately 3300 Hz) with a mean ± standard deviation Pearson correlation coefficient of 0.723 ± 0.109. The correlation decreased in the higher frequency range, with a Pearson correlation coefficient of 0.564 ± 0.116 (Figure 3A, far-right panel).

Figure 3B shows the plots of predicted EEG and original EEG with their similarity measurement in terms of the Pearson correlation coefficient at Cz. The EEG prediction accuracies at Cz were comparable between the detrended cross-correlation function (0.382) and ridge regression (0.389). The topographical distribution of the Pearson correlation coefficients indicated high performance in the central, right frontal, and left temporal lobes for both the detrended cross-correlation function and ridge regression (Figure 3C). However, the topography of the ridge regression revealed more detailed information regarding the performance of individual channels than the detrended cross-correlation function.

3.3. Significance of the Detrended Cross-Correlation

Figure 4 depicts the computational time and model performance regarding the cross-validation Pearson correlation of the 20 individual models, computed over the lag range of −250 ms to 0 ms (baseline model) and −100 ms to 500 ms (trained model), for the detrended cross-correlation and ridge regression.

In comparison with ridge regression, the detrended cross-correlation function required significantly less computational time to train both the univariate (p < 0.05) and multivariate (p < 0.01) forward models (Figure 4A). In terms of cross-validation performance, the trained models of both the detrended cross-correlation and ridge regression significantly outperformed the corresponding baseline models (p < 0.01; Figure 4B). No significant difference was observed in the cross-validation Pearson correlation between the two models.

Figure 4A shows that the detrended cross-correlation outperformed ridge regression as the complexity increased.

Table 1. Summary of computational time for training each model and Pearson correlation coefficient for validating EEG prediction using the detrended cross-correlation function and ridge regression. The notations of **, *, and n.s (not significant) denote p < 0.01, p < 0.05, and p > 0.05, respectively.

	Computational Time (s)			Validation Pearson Correlation
	Detrended Cross- Correlation	Ridge Regression	p-Value	Detrended Cross- Correlation	Ridge Regression	p-Value
Univariate forward model	1.09 ± 0.27	1.29 ± 0.33	*	0.024 ± 0.010	0.027 ± 0.009	n.s
Multivariate forward model	3.73 ± 1.56	21.69 ± 2.12	**	0.038 ± 0.012	0.038 ± 0.011	n.s

summarizes the mean and standard deviation of the computational time (s) and validation of the Pearson’s correlation coefficients for the trained models of the detrended cross-correlation function and ridge regression, as seen in Figure 4B. In the univariate forward model, the computational time for ridge regression was 1.18 times that required for detrended cross-correlation, while that number increased to 5.82 times in the case of the multivariate forward model. Regarding prediction accuracy, Figure 4 and Table 1 show that both approaches converged to a similar performance in EEG prediction for both the univariate and multivariate models (p > 0.05). Figure 4C illustrates the correlation (using the Pearson correlation coefficient) between performance of individual models (n = 20) computed by the detrended cross-correlation (ordinate) and ridge regression (abscissa) functions. The Pearson correlation coefficients denoted as 0.98 and 0.95 are correlations between model performance computed by the detrended cross-correlation and ridge regression for the univariate and multivariate model, respectively.

4. Discussion

In this study, we aimed to demonstrate the potential of the detrended cross-correlation function for capturing neural entrainment in continuous speech. The EEG signals were predicted from the speech envelope (univariate forward model) and speech spectrogram (multivariate forward model). The detrended cross-correlation function performed significantly better within a meaningful lag range than the baseline model computed from the baseline period. Furthermore, the proposed model has computational advantages while achieving the same level of prediction performance as the widely used ridge regression model.

Ridge regression implemented in the mTRF toolbox [18] served as a reference for our study. It is widely applied in sensory neuroscience to model brain responses to continuous stimuli [31,32]. Our results showed that the detrended cross-correlation function resulted in TRFs or speech STRFs comparable to those of ridge regression. It also exhibited a similar performance to individual models in terms of the Pearson correlation coefficient between the predicted EEG and original EEG. Moreover, the detrended cross-correlation function can more efficiently interpret neurophysiological data. The detrended cross-correlation function has a significant advantage in terms of computational time for both the univariate and multivariate forward models (Figure 4A). Moreover, ridge regression requires cross-validation to determine the optimal value of the regularization parameter, which can be computationally expensive if a wide range of parameter values are considered. Specifically, the TRFs estimated by ridge regression are sensitive to the regularization parameter, which can lead to unstable results, depending on the range of values chosen for regularization. In contrast, the detrended cross-correlation function is a parameter-free model that does not require a parameter optimization procedure.

Previous studies applied cross-correlation analysis to investigate how the brain perceives attended and unattended speech [33,34]. Although cross-correlation analysis is more straightforward than the ridge regression method, it is better suited for inputs modulated by a stochastic process. Because of this limitation, ridge regression has been widely applied in sensory neuroscience for continuous stimuli [31,33,35,36]. In this study, we proposed a detrended cross-correlation function to investigate neural entrainment to speech envelopes. This method overcomes the limitations of conventional cross-correlation analysis in analyzing nonstationary or slow-modulating signals while maintaining its straightforwardness. Additionally, the estimations obtained from the detrended cross-correlation function were comparable to those of ridge regression, demonstrating the utility of the detrended cross-correlation function in auditory scene analysis research.

The short trial duration may lead to overfitting, particularly in the large response to speech onset. Thus, even though the results of this study revealed a significant degree of similarity in the TRFs generated by the detrended cross-correlation function and ridge regression, the influence of trial duration should also be considered. Therefore, this study was extended to analyze a publicly available dataset featuring longer EEG trials in response to natural narrative speech in this section [37]. The open dataset is briefly summarized as follows. The 128-channel EEG data in response to narrative speech were recorded from 17 subjects, and the EEG data were segmented to have a duration of 10 and 60 s. The performance of the detrended cross-correlation function to generate long-speech TRFs was evaluated. The detrended cross-correlation function generated TRFs comparable to those generated by ridge regression, with an averaged correlation coefficient across 128 channels of 0.91 as seen in Figure 5. Moreover, although there was no difference in reconstructing EEG data, the computation of the detrended cross-correlation function was more efficient. These results indicated that the detrended cross-correlation can also effectively and comparably capture neural responses to long continuous speech.

Recently, SI methods have been applied to investigate the neural source dynamics in response to continuous stimuli [38,39] and the neural response to the temporal fine structure of continuous stimuli [40]. However, the computational demands of these methods increase due to the significant large quantity of cortices at the source level or the necessity of a high sampling rate for temporal fine structure analysis. Thus, the detrended cross-correlation function, which inherently provides better computational efficiency, may be properly employed. The validation of source level analysis and temporal fine structure will be performed in the near future.

Although the detrended cross-correlation function exhibited performance comparable to ridge regression in exploring neural responses to continuous features (such as speech envelope and speech spectrogram), further investigation to examine neural responses to discrete features such as phonemic onset and semantic dissimilarity should be considered [31,37]. As the detrended cross-correlation function may generate a more zero-centered TRF (i.e., dampening trends) compared to those generated by the ridge regression function, careful consideration is needed to employ the discrete features. This will be systematically investigated in future studies.

5. Conclusions

The detrended cross-correlation function offers advantages in terms of straightforwardness and computational simplicity. By leveraging these advantages, applying the detrended cross-correlation function to SI methods could potentially provide a more efficient and computationally feasible approach for investigating neural responses to continuous stimuli at the source level or analyzing the temporal fine structure of stimuli. This research direction could offer insights into the neural mechanisms underlying the processing of continuous stimuli while addressing the computational challenges associated with existing SI methods.