Wavelet Canonical Coherence for Nonstationary Signals

Assessing the dependence structure between node clusters in a network is one of the most critical aspects of network time series analysis. Many models and frameworks have been developed to capture between-clusters association (e.g., correlation, coherence, and causality). Most existing methods characterize the dependence between two clusters through the dependence between (many) node pairs. However, in many scenarios, the primary interest lies in understanding the dependence structure between two groups of multivariate time series rather than individual processes. Figure 1 illustrates this perspective using brain activity signals. In this experiment, local field potential (LFP) activity was recorded from multiple electrodes implanted in different subregions of the hippocampus of rodents (rats) as they performed a complex sequence memory task. Instead of focusing on coherence between individual channels (electrodes), the main goal is to quantify time-varying functional interactions between groups of electrodes in order to understand how information processing differs in these two subregions. Similar challenges arise in other domains. For instance, in finance, understanding the dependence between entire market sectors (e.g., technology and energy) can be more informative than analyzing associations between individual stocks. These scenarios require a framework capable of capturing dynamic coherence between sets of nonstationary multivariate signals. In this paper, we propose a novel framework called scale-specific wavelet canonical coherence (WaveCanCoh) to characterize time-localized and scale-specific coherence between two clusters of multivariate time series. By leveraging the time-frequency localization properties of wavelets, WaveCanCoh is well-suited for analyzing nonstationary multivariate signals in neuroscience, finance, and other fields where transient, cross-group interactions are of scientific interest.

Canonical variate analysis (CVA) (Hotelling (1936)) provides a method for measuring the correlation between two vector variables, its application to time series data having started in the 1950s and 1960s, primarily in the econometrics and signal processing fields. Box and Tiao (1977) and Geweke (1982) extend CVA to time series for forecasting and causality detection. Brillinger (2001) provides a spectral domain formulation of canonical correlation, useful for frequency-domain time series analysis. This approach enables the analysis of canonical coherence between two sets of time series across different frequency bands. Many modern studies have been developed within this framework (e.g., Talento et al. (2024) and Viduarre et al. (2019)), and related methods have been widely applied across various fields, including neuroscience, finance, speech processing, and machine learning.

However, previous methods rely on the assumption that time series are (weakly) stationary, meaning their statistical properties (e.g., expectation and covariance, spectral signature) remain constant over time. In practice, this assumption often does not hold for time series that arise in practice. Moreover, two sets of time series commonly exhibit time-varying global coherence, which can sometimes be crucial for analysis. Thus, a method capable of handling nonstationary time series is necessary. Wavelet analysis is a widely used tool for studying nonstationary time series, as its localization property allows for the examination of localized correlations between two time series across both time and frequency domains. Wavelets are particularly effective for capturing transient properties of nonstationary signals due to their compact support, which can be compressed or stretched to adapt to the dynamic characteristics of the signal. Wavelet coherence has been well-defined and extensively studied in previous research, with applications spanning various fields (Embleton et al., 2022; Grinsted et al., 2004). To the best of our knowledge, no existing work has extended classical canonical coherence to the wavelet domain to measure the time-evolving canonical coherence between two groups of multivariate time series.

The key novelty of this paper lies in the development of a comprehensive and rigorous framework based on wavelets for measuring canonical coherence between two sets of nonstationary multivariate time series. Specifically, our main contributions include: (1) we define scale-specific wavelet canonical coherence (WaveCanCoh) and introduce its use as a tool to quantify the coherence between two sets of multivariate time series, (2) we provide a complete and theoretically justified algorithm for its estimation, and (3) we apply our method to LFP activity data from multiple electrodes to quantify dynamic interaction patterns among different subregions in the hippocampus. Multivariate locally stationary wavelet processes (MvLSW) underpin our construction, and the reader is directed to Nason et al. (2000) and Park et al. (2014) for details on their construction. Our framework not only captures the time-varying coherence between two sets of signals but also determines the contribution of each individual channel to the global coherence. Compared to previous models, the proposed approach provides a detailed, localized characterization of interactions within the multivariate time series. Our findings on the LFP activity data offer new insights into the functional relationship between hippocampal subregions, demonstrating the potential of our method to advance the study of functional brain connectivity.

The format of the paper is as follows. Section 2 overviews the current methodology for assessing time series canonical coherence. Section 3 provides a brief overview of MvLSW processes, supporting a detailed introduction to our proposed WaveCanCoh framework. Its estimation, and that of related parameters, is tackled in Section 4. Section 5 validates the proposed framework and demonstrates its performance through simulation. In Section 6 we apply the WaveCanCoh method on LFP activity data collected from rats to investigate the dynamic interactions of different subregions in the hippocampus during memory tasks. Section 7 concludes the paper.

First, we provide a brief introduction to classical canonical correlation analysis for time series, with the primary goal to characterize dependence between two clusters of time series, where each cluster features several nodes. In particular, we consider two multivariate time series with dimensions $p$ and $q$ respectively, denoted by $\mathbf{X}_{t}=\left(X_{t}^{(1)},\ldots,X_{t}^{(p)}\right)^{\top}$ and $\mathbf{Y}_{t}=\left(Y_{t}^{(1)},\ldots,Y_{t}^{(q)}\right)^{\top}$, for $t=\{1,\ldots,T\}$. Typically, $\{\mathbf{X}_{t}\}$ and $\{\mathbf{Y}_{t}\}$ are assumed to be zero-mean weakly stationary time series. Letting $\mathbf{Z}_{t}$ denote the concatenated ($p+q$) dimension time series, $\mathbf{Z}_{t}=\left(X_{t}^{(1)},\ldots,X_{t}^{(p)},Y_{t}^{(1)},\ldots,Y_{t}^{(q)}\right)^{\top}$, its covariance matrix at lag $\tau$ is

where $\boldsymbol{\Sigma}_{\mathbf{XX}}(\cdot)$, $\boldsymbol{\Sigma}_{\mathbf{YY}}(\cdot)$ are the autocovariance matrices of $\{\mathbf{X}_{t}\}$ and $\{\mathbf{Y}_{t}\}$ respectively, and $\boldsymbol{\Sigma}_{\mathbf{XY}}(\cdot)$, $\boldsymbol{\Sigma}_{\mathbf{YX}}(\cdot)$ are their cross-covariances. The canonical correlation between $\{\mathbf{X}_{t}\}$ and $\{\mathbf{Y}_{t}\}$ at lag $\tau$, $\boldsymbol{\rho}(\tau)$, is defined as

where $\mathbf{a}\in\mathbb{R}^{p}$ and $\mathbf{b}\in\mathbb{R}^{q}$ are called canonical correlation vectors, subject to standardized constraints $\mathbf{a}^{\top}\boldsymbol{\Sigma}_{\mathbf{XX}}\mathbf{a}=1$, $\mathbf{b}^{\top}\boldsymbol{\Sigma}_{\mathbf{YY}}\mathbf{b}=1$ (Brillinger, 2001). Thus, the canonical correlation can be rewritten as

The solution for $\mathbf{a}$ and $\mathbf{b}$ in equation (1) can be obtained from the eigenvectors of the following matrices, respectively (in what follows, $\tau$ is dropped for brevity),

Here, $\mathbf{a}$ and $\mathbf{b}$ are the eigenvectors corresponding to the largest eigenvalue, $\lambda$, of matrices (3) and (4), respectively, and for largest canonical correlation coefficient we have $\boldsymbol{\rho}=\sqrt{\lambda}$ (Mardia et al., 1979). The framework above yields the canonical correlation between $\{\mathbf{X}_{t}\}$ and $\{\mathbf{Y}_{t}\}$. However, in many practical cases, such as EEG analysis, canonical correlation in the spectral domain is more meaningful than in the time domain, as components at different frequencies (or scales) reveal crucial information about neural dynamics and functional connectivity (Ombao and Pinto, 2024) , Brillinger (2001) extends the time-domain canonical correlation into the spectral domain). Namely, suppose the spectral matrix of $\mathbf{Z}_{t}=\left(\mathbf{X}_{t}^{\top},\mathbf{Y}_{t}^{\top}\right)^{\top}$ is

where $\mathbf{f_{XX}}(\boldsymbol{\omega})$ is the $p\times p$ autospectral matrix of $\{\mathbf{X}_{t}\}$, $\mathbf{f_{YY}}(\boldsymbol{\omega})$ is the $q\times q$ autospectral matrix of $\{\mathbf{Y}_{t}\}$, and $\mathbf{f_{XY}}(\boldsymbol{\omega})$ is the $p\times q$ cross-spectral matrix between $\{\mathbf{X}_{t}\}$ and $\{\mathbf{Y}_{t}\}$. Given vectors $\mathbf{a}\in\mathbb{C}^{p}$ and $\mathbf{b}\in\mathbb{C}^{q}$, such that $\mathbf{a}^{\top}\mathbf{f_{XX}}(\boldsymbol{\omega})\mathbf{a}=\mathbf{b}^{\top}\mathbf{f_{YY}}(\boldsymbol{\omega})\mathbf{b}=1$, the canonical coherence at frequency $\boldsymbol{\omega}$ is

By solving the maximization problem in equation (5), the canonical coherence vectors $\mathbf{a}$ and $\mathbf{b}$ are determined, leading to the quantification of the canonical coherence at frequency $\boldsymbol{\omega}$.

Note the classic canonical coherence in (5) completely ignores temporal dynamics, a consequence of the stationarity assumption where dependence between clusters is imposed to remain constant over time, whilst most real-world data, such as EEG, exhibit nonstationarity (Huang et al., 2004; Knight et al., 2024; West et al., 1999). Hence the lack of time-localization information in the above approach may result in misleading results and a novel method capable of capturing time-varying canonical coherence is needed.

Our WaveCanCoh framework is built upon the multivariate locally stationary wavelet (MvLSW) process (Nason et al. (2000), Park et al. (2014)), which is a model based on wavelet analysis for time series. A brief overview of wavelets and highlight of their differences from Fourier-based methods are provided in Appendix A.

A new representation for discretely sampled nonstationary time series based on discrete non-decimated wavelets is the locally stationary wavelet process introduced by Nason et al. (2000), later extended to a multivariate framework in Park et al. (2014). A $P$-variate stochastic process with time evolving second-order structure, $\mathbf{X}_{t}=\left(X_{t}^{(1)},X_{t}^{(2)},\ldots,X_{t}^{(P)}\right)^{\top}$, where $t=1,\ldots,T$, can be represented with the MvLSW formulation

where $\mathbf{V}_{j}(k/T)$ is a $P\times P$ transfer function matrix assumed to have a lower-triangular form; $\left\{\psi_{j,k}\right\}_{j,k}$ is a set of discrete non-decimated wavelets; $\left\{\mathbf{z}_{j,k}\right\}_{j,k}$ is a set of $P\times 1$ uncorrelated random vectors with (column) mean vector $\mathbf{0}$ and $P\times P$ identity covariance matrix. Since the wavelet basis $\psi_{j,k}(t)$ is localized in both time and frequency, the transfer matrix $\mathbf{V}_{j}(k/T)$ provides a measure of the time-varying contribution to the variance among channels at a specific scale $j$ and rescaled time $u=k/T$, thus enabling the statistical properties of the process $\{\mathbf{X}_{t}\}$ to change smoothly over time.

The time-varying statistical properties of $\{\mathbf{X}_{t}\}$ can be captured through the localized, scale-specific local wavelet spectral matrix (LWS, Park et al. (2014)), $\mathbf{S}_{j}(u)$, defined at scale $j$ and rescaled time $u\in(0,1)$, as

Note $\mathbf{S}_{j}(u)$ is a $P\times P$ symmetric, positive semi-definite matrix and its $(p,q)$ entry, $S_{j}^{(p,q)}(u)$, denotes the cross-spectrum between channels $p$ and $q$. We now extend the LWS matrix construction from a single set of multivariate time series to a cross-group LSW matrix, between $\mathbf{X}_{t}=\left(X_{t}^{(1)},\ldots,X_{t}^{(P)}\right)^{\top}$ and $\mathbf{Y}_{t}=\left(Y_{t}^{(1)},\ldots,Y_{t}^{(Q)}\right)^{\top}$. Denoting $\mathbf{Z}_{t}=\left(\mathbf{X}_{t}^{\top},\mathbf{Y}_{t}^{\top}\right)^{\top}$, the LWS matrix of $\{\mathbf{Z}_{t}\}$ at scale $j$ and rescaled time $u$, $\mathbf{S}_{j;\mathbf{ZZ}}(u)$, is

In equation (9), $\mathbf{V}_{j;\mathbf{Z}}(u)$ denotes the $(P+Q)\times(P+Q)$ transfer function matrix of the MvLSW process $\{\mathbf{Z}_{t}\}$, and $\mathbf{S}_{j;\mathbf{ZZ}}(u)$ is its corresponding LWS matrix. The main diagonal blocks $\mathbf{S}_{j;\mathbf{XX}}(u)$ $(P\times P)$ and $\mathbf{S}_{j;\mathbf{YY}}(u)$ $(Q\times Q)$ denote the auto-LWS matrices of the $\{\mathbf{X}_{t}\}$ and $\{\mathbf{Y}_{t}\}$ processes, respectively, while $\mathbf{S}_{j;\mathbf{XY}}(u)$ and $\mathbf{S}_{j;\mathbf{YX}}(u)$ denote their cross-LWS matrices.

The LSW matrix quantifies the localized contributions to the process variance for individual and cross-channels, which motivates us to next define the localized canonical coherence between two sets of locally stationary time series at a specific scale (corresponding to a determined frequency band).

The canonical coherence vectors $\mathbf{a}_{j}(\cdotp)$, $\mathbf{b}_{j}(\cdotp)$ can be obtained by maximizing (10) and the solution can be found by solving the eigenvalue and eigenvector problem associated with the following matrices

Denote by $\Lambda_{j;\mathbf{a}}^{(k)}(u)$ the $k$-th largest eigenvalue of matrix $\mathbf{M}_{j;\mathbf{a}}(u)$ in equation (11), and by $\Lambda_{j;\mathbf{b}}^{(l)}(u)$ the $l$-th largest eigenvalue of matrix $\mathbf{M}_{j;\mathbf{b}}(u)$ in equation (12), for $k,l=1,\ldots,\min(P,Q)$. An important observation is that $\mathbf{M}_{j;\mathbf{a}}(u)$ and $\mathbf{M}_{j;\mathbf{b}}(u)$ share the same eigenvalues (see Appendix B for details), hence denoting by $\Lambda_{j}^{(1)}(u)=\Lambda_{j;\mathbf{a}}^{(1)}(u)=\Lambda_{j;\mathbf{b}}^{(1)}(u)$ their largest eigenvalue, the canonical coherence between $\{\mathbf{X}_{t}\}$ and $\{\mathbf{Y}_{t}\}$ at rescaled time $u$, as defined in equation (10), becomes

The eigenvectors of $\mathbf{M}_{j;\mathbf{a}}(u)$ and $\mathbf{M}_{j;\mathbf{b}}(u)$ corresponding to $\Lambda_{j}^{(1)}(u)$ provide the solutions to the canonical directions of $\{\mathbf{X}_{t}\}$ and $\{\mathbf{Y}_{t}\}$, respectively. Proof: see Appendix B.

An important extension of our framework is the incorporation of leading-lag relationships into the scale-specific wavelet canonical coherence (WaveCanCoh). To account for potential causal effects, we define a lagged joint process, $\mathbf{Z}_{t}(h)=\left(\mathbf{X}^{\top}_{t},\mathbf{Y}^{\top}_{t+h}\right)^{\top}$, where $h$ is the value of lag, and $t=1,\ldots,T-h$ for $h>0$. We define the LWS matrix of $\{\mathbf{Z}_{t}(h)\}$ at scale $j$, as

where $\mathbf{S}_{j;\mathbf{XY}}(u,h)$ denotes the cross-LWS matrix between $\mathbf{X}_{[uT]}$ and $\mathbf{Y}_{[uT]+h}$, capturing the interaction between current values of $\mathbf{X}$ and future values of $\mathbf{Y}$. Based on this construction, we can define and estimate the lagged version of WaveCanCoh, enabling us to infer potential causal relationships between two groups of time series.

The framework above allows us to capture the time-varying overall association between two sets of multivariate time series, as well as the time-varying contributions from each individual channel within these sets. However, a natural consideration is how to project this time-varying coherence at each scale $j$ into the frequency domain in a manner consistent with the Fourier-based method described in Section 2, as a key concern in many analyses is to interpret the results in the frequency domain. As mentioned earlier, each scale in the wavelet analysis corresponds approximately, but not exactly, to a specific frequency band. This correspondence is governed by the unique filtering mechanism of wavelets, an explanation for this corresponding relationship is provided in Appendix A.

In Section 3, we developed a rigorous framework that allowed us to introduce the localized, scale-specific wavelet canonical coherence. In this section, we propose a well-behaved estimation procedure for quantifying the canonical coherence and corresponding canonical vectors. We start by estimating the local wavelet spectrum (LWS) matrices in equations (11) and (12) in the spirit of Park et al. (2014), given by

$k$ represents the shift of the wavelet function and is equivalent to time $k=[uT]$ in our context, and $M$ is the half-width of the rectangular smoothing kernel, controlling the amount of temporal smoothing. The matrix $\mathbf{I}_{l,k}$ is the raw periodogram at scale $l$ and time $k$, obtained as

This multistep procedure yields consistent estimators of the LWS matrices under the asymptotic conditions $T,M\rightarrow\infty$ and $M/T\rightarrow 0$ (Park et al., 2014). We propose the following estimator for the scale-specific wavelet canonical coherence (WaveCanCoh)

where $\widehat{\Lambda}_{j}^{(1)}(u)$ is the largest eigenvalue of $\widehat{\mathbf{M}}_{j;\mathbf{a}}(u)$ and $\widehat{\mathbf{M}}_{j;\mathbf{b}}(u)$, defined as

The estimated localized, scale-specific canonical direction vectors $\widehat{\mathbf{a}}_{j}(u)$ and $\widehat{\mathbf{b}}_{j}(u)$ are the eigenvectors of $\widehat{\mathbf{M}}_{j;\mathbf{a}}(u)$ and $\widehat{\mathbf{M}}_{j;\mathbf{b}}(u)$ respectively, associated with $\widehat{\Lambda}_{j}^{(1)}(u)$. These quantities provide estimates of the time-varying global coherence and the channel-specific contributions at scale $j$. The proposed estimators are consistent with the true quantities they aim to approximate, provided certain asymptotic conditions are met. These include increasing sample size and appropriate smoothing bandwidth, ensuring reliable estimation in the limit. Proof: See Appendix B.

Algorithm 1 summarizes the estimation procedure for WaveCanCoh and its results can be further used to investigate the temporal channel contributions to the global association, at a particular scale.

In this section, we implement the proposed framework using simulated data under two distinct scenarios, one adhering to the MvLSW assumptions underpinning our method, while the other introduces nonstationarity without strictly satisfying the MvLSW assumptions. These setups allow us to validate both the theoretical soundness and empirical performance of the proposed approach, as well as to assess its robustness and practical applicability in real-world scenarios where model assumptions may be violated. To further evaluate performance, we also compare the results with those obtained from the classical Fourier-based canonical coherence approach.

To demonstrate the practical utility of our proposed WaveCanCoh framework, we analyze LFP activity recorded from the hippocampus of rats engaged in a sequence memory task (Allen et al., 2016; Shahbaba et al., 2022). The data were recorded using a 22-electrode microdrive implanted in the CA1 subregion to capture high-resolution LFP signals across all channels at a sampling rate of $1000Hz$. In this task, rats were tested on their memory of a sequence of five odors (odors ABCDE). Each odor was presented for $\sim$1.2 second ($s$) and a variable delay of $\sim$5$s$ separated each odor (see Figure 1). For each trial (i.e., each odor presentation), the rat had to judge whether the odor was presented "in sequence" (e.g., ABC…) or "out of sequence" (e.g., ABD…) and indicate their decision by holding their nosepoke response until a tone signal (at 1.2$s$) or withdrawing before the signal, respectively. Correct-response trials (i.e., correct "in sequence" or "out of sequence" decisions) were rewarded. LFP activity data were recorded over a 4$s$ period ($T=4000$ time points) per trial, with $t=0$ marking the moment the rat initiated a nosepoke to receive the odor stimulus. This paradigm provides a well-controlled setting to investigate dynamic, time-varying functional interactions in the hippocampus during memory-guided decisions. We employ WaveCanCoh framework to quantify frequency-specific functional coherence between two groups of hippocampal electrodes (T1, T2, T4, T5 and T13–T17), and to examine how coherence patterns differ between correct- and incorrect-response trials ("in sequence" trials only). Specifically, we analyze LFP data from the rat Mitt, which included 40 correct-response trials and 32 incorrect-response trials. Figure 4 presents the estimated wavelet canonical coherence at scale $j=5$, corresponding to the $15.625-31.25Hz$ frequency band. The results, averaged across trials for each condition, reveal dynamic changes in inter-regional coherence, with a pronounced peak around the time of odor stimulus delivery ($t=0$). Notably, distinct patterns emerge between correct and incorrect trials, suggesting that coherent activity in this frequency band may play a role in supporting successful memory retrieval and decision making. More results for several other scales can be found in Figure 10 in Appendix D.1.

To further interpret the coherence patterns, Figure 5 provides a spatial summary of the canonical coherence between the two electrode groups at several selected time points. The double-headed arrows represent the magnitude of estimated coherence between the two regions, while the numbers in the circles reflect the individual channel contributions to the global coherence, derived from the elements of the canonical vectors $\mathbf{a}_{5}(\cdotp)$ and $\mathbf{b}_{5}(\cdotp)$.

The results highlight that both the strength and structure of inter-regional coherence vary dynamically over time and differ across trial outcomes, reflecting the nonstationary nature of neural interactions during task performance. Specifically, incorrect-response trials exhibit lower coherence at time points immediately following the odor stimulus delivery and are often driven by a few dominant channels, while correct trials show higher coherence values and relatively balanced contributions across channels. These findings illustrate the necessity of a framework like WaveCanCoh, which simultaneously identifies time-varying and scale-specific dependence structures between multivariate time series. Traditional stationary or pairwise approaches would fail to capture such nuanced dynamics, highlighting the need for a method like WaveCanCoh to extract such complex relationships in brain activity. The directed interactions between brain regions are also explored in Appendix E with proposed Causal-WaveCanCoh framework.

To further validate the existence of significant differences in the activity between correct- and incorrect-response trials, we propose a time-specific detection procedure based on the permutation test to determine temporally localized differences in the wavelet canonical coherence at a given scale between conditions, while maintaining the nonparametric nature of the statistical inference. The detailed inference procedure steps are shown in Algorithm 2 (Appendix D.2) and the results of the time-specific permutation test for WaveCanCoh across five wavelet scales and four selected time points are summarized in Table 1 (Appendix D.3). Statistically significant differences in canonical coherence between correct and incorrect trials emerge at time points following odor sampling ($t=0$), predominantly at intermediate wavelet scales corresponding approximately to the $8-62Hz$ frequency range. On the other hand, no significant differences between the two groups are observed prior to stimulus onset, indicating that the coherence patterns distinguishing inter-regional communication among trial types are tightly linked to task engagement. These findings demonstrate the effectiveness of the proposed WaveCanCoh framework and associated permutation test in capturing localized, frequency-specific differences in neural coordination between behavioral conditions, thus offering a powerful tool for analyzing complex brain interactions.

In this paper we introduced a novel methodological framework, scale-specific wavelet canonical coherence (WaveCanCoh), designed to quantify the dynamic multiscale coherence between two sets of nonstationary multivariate time series. Our primary contributions include the rigorous definition of WaveCanCoh within the multivariate locally stationary wavelet (MvLSW) framework and the development of a comprehensive estimation and theoretically-backed inference procedure based on wavelet analysis. We validated our proposed methodology through simulation studies, demonstrating its accuracy in tracking true coherence structures. The application to local field potential (LFP) activity recorded from subregions of the hippocampus effectively showcased the capability of WaveCanCoh to identify nuanced spatio-temporal coherence patterns associated with cognitive performance, while our permutation-based inference procedure provided a robust, nonparametric approach for detecting significant coherence differences between conditions. Compared to existing stationary and Fourier-based canonical coherence methods, WaveCanCoh is shown to offer significant advantages, particularly through its ability to adaptively capture transient and time-localized interactions, hence it is highly suitable for analyzing signals in neuroscience and other fields with data exhibiting dynamic cross-group interactions.