Phase coherence—A time-localized approach to studying interactions

Oscillations and waves are ubiquitous in nature. They occur in mechanical and dynamical systems in virtually all areas of science: many physiological processes are oscillatory, such as the beating of the heart, breathing, or neuronal oscillations in the brain; the ecology abounds with seasonal cycles; most dynamical phenomena in astrophysics and space science are oscillatory, as are geological and hydrodynamics phenomena, such as ocean waves or earthquakes; there are business cycles in economy; strings in musical instruments produce vibrations, as do many man-made devices. Most electronic devices, the Internet, TV signals, communication systems, and medical imaging, use electromagnetic waves. The study of oscillations and waves is, therefore, essential for understanding the universe, as stated by Tesla in the quote: “If you want to find the hidden secrets of the universe, you must think in terms of energy, frequency, and vibration.”

While the underlying dynamical system may be very different in distinct cases, oscillatory processes share two key time-dependent features: amplitude (associated with the energy of the oscillation) and phase (associated with the time evolution of the oscillation). To identify interactions between different parts of a system, we can calculate the similarity of these features using the physical property known as coherence.

In this paper, we provide a review of coherence, beginning in its conceptualization in physics and subsequently evaluating relevant numerical methods used to measure coherence. In particular, we improve current understanding by both establishing the fundamental importance of taking a time-localized approach to coherence and comparing a method based on amplitude and phase to one only using phase information.

In Sec. II, we provide an overview of the development of coherence in physics and its adoption in time-series analysis. We also provide a definition of coherence based on the Fourier transform and explain the differences between coherence and the related concept of synchronization.

In Sec. III, we provide a model for a dynamical system, which is used to numerically illustrate the differences between the phase-only and amplitude-weighted methods of measuring coherence when the system is perturbed by different forms of noise.

In Sec. IV, we introduce wavelet-based coherence and explain the consequences of moving to the time–frequency domain that arise from the uncertainty principle. In this section, we also specify the alternate definitions of coherence in amplitude and phase and, based on results found using the illustrative model, argue that phase coherence is more resistant to the effects of noise and particularly movement artifacts.

In Sec. V, we provide an in-depth guide to the application of coherence in time-series analysis, including how to identify significant coherence.

In Sec. VI, we apply this knowledge and evaluate the two methods considered by considering four real-world problems, including the cardio-respiratory interaction, noisy electroencephalography (EEG) and functional near infrared spectroscopy (fNIRS) data, and electron dynamics on the surface of liquid helium.

We conclude in Sec. VII with a discussion of the time-localized approach to coherence and the impact of using methods based on only phase to those that rely also on amplitude information.

The theory of waves was initially developed by Young, Huygens, and Fresnel.¹ Along with providing explanations for phenomena, such as diffraction and refraction, they also studied wave interference. In this latter case, multiple waves combine to produce a characteristic pattern of spatially and time-localized maxima and minima. However, this effect is only seen clearly when the change in the phase of the waves is the same. It is this property of the waves that we term coherence.

The study of interference and wave coherence has already led to many well-known discoveries. These include the Michelson–Morley experiment, which disproved the existence of the luminiferous ether.² Variations of Young’s double-slit experiment have also played an important role in the understanding of wave–particle duality.^3–5 In addition, the drive to develop a coherent source of light led to the invention of the laser.⁶ Subsequent to the development of lasers, larger-scale interference experiments have been possible, which resulted in the discovery of gravitational waves.⁷ Coherence is now studied across a broad spectrum of domains. This includes solid state and quantum physics,^8–11 remote sensing,¹² electrophysiology,^13–15 communications,¹⁶ and space science.¹⁷ 

With the advent of computers, the study of coherence is no longer restricted to physical experiments. Numerical methods allow for the analysis of oscillations in recorded data. Using this recorded data, coherence can be investigated.¹⁸ Coherence between different parts of a dynamical system can result from either synchronization or from modulation by a common process. While one can separately analyze two variables and qualitatively assess the common features present in each, interactions are often nonlinear in nature and, hence, difficult to discern. Coherence, therefore, provides a useful quantitative measure to identify these interactions.

An important aspect of coherence is that it is a time-localized phenomenon. This makes it particularly useful for analyzing dynamics comprised of oscillations with time-dependent quantitative characteristics. Such dynamics has been modeled using chaotic, stochastic, and non-autonomous systems.^19,20 Time-series analysis methods that give a non-time-dependent representation of a time-series, such as its histogram or Fourier transform, may yield some insight into the amplitudes of oscillations present. However, these methods will generally provide little understanding of phase dynamics if the quantitative characteristics of the oscillations, or of their interactions with each other, are being modulated over time. In contrast to this, the time evolution of phases carries a great wealth of information about the underlying system when such time modulation exists.²¹ 

Time-evolving time-localized analysis is typically performed in the time–frequency domain. This type of analysis was originally developed in quantum mechanics, with the distribution proposed by Wigner providing the highest possible frequency resolution that is mathematically possible within the limitations of the uncertainty principle.²² Ville later applied this function in the context of time–frequency analysis more generally.²³ At the same time, the windowed Fourier transform was also developed,²⁴ and the field has since been advanced with the introduction of the continuous wavelet transform.^25,26 Time–frequency analysis has been applied most commonly to deal with simple forms of nonstationary data, with applications in communications, radar, sonar, and acoustics.²⁷ Recently, it has also been invaluable in the analysis problems, such as turbulence,²⁸ brain signals,²⁹ blood flow,³⁰ and excited electron oscillations on liquid helium.³¹ These systems involve multiple potentially mutually interacting oscillatory processes that take place simultaneously across a range of timescales; we refer to such systems as multi-scale systems.

One specific advantage of the time–frequency methods is that they, to various degrees, allow for the time-localized extraction of instantaneous phases over time (see, e.g., Ref. 32). These phases can be studied further to give insight into the system. This can be seen in phase synchronization methods, which have been applied to the cardiorespiratory system.³³ Phase differences can also be observed and point to delays in coupled networks of oscillators, such as those seen in biology.³⁴ Beyond this, we can estimate coupling functions and infer the directionality of coupling (see Ref. 35 and the references therein). In the case of weakly coupled oscillator networks, connectivity can be inferred directly from the phases.³⁶ There are also phase stability methods, which have been used to find stable oscillations in the heart rate variability.³⁷ 

In the case of coherence represented in the time–frequency domain, the initial development of the methods was motivated by applications to biomedical data. Specifically, it has been of great importance to the mapping of functional connectivity and study of synchronization in the brain.^38–45 At the same time, the development of time–frequency coherence has spearheaded investigations into microvasculature dynamics.^46–51 It has since been used in other biomedical studies and found use as a marker for ageing of the cardiorespiratory system,⁵² as well as revealing the relation between the width of the subarachnoid space and blood pressure.⁵³ Moreover, the generality of time–frequency coherence means that it has found applicability elsewhere. In particular, these methods have also been used extensively in the analysis of solar, geophysical, and meteorological time-series to determine the Earth–Sun dynamical relationship.^54–57 Coherence has also found use in the analysis of economic time-series, where it has been used to identify instability and risk in specific markets as well as the relation between the monetary policy and the macroeconomic activity.^58–62 It has also been applied in the case of cyclo-nonstationarity, where it has been used to analyze mechanical systems, such as engines and wind turbines.⁶³ Further examples include the evaluation of electron dynamics,³¹ behavioral rhythms in mice,⁶⁴ and social networks.⁶⁵ 

It is worth noting that coherence should not be confused with synchronization. In terms of dimensionality, synchronization is defined specifically in the time dimension and, therefore, applies to the dynamics of oscillations in time. In contrast, coherence refers to a more general phenomenon, which extends to waves that are defined across space as well as time.

There are also important differences in the context of time-series generated by dynamical systems. While many types of synchronization exist, they all result from an interaction between two or more oscillations.^33,67 As such, synchronization refers to a process of adjustment of rhythms caused by interactions. In contrast, coherence implies that two oscillations are observed to have the same frequency and frequency modulation, but this does not necessarily imply that they are coupled.

As examples, consider two linear oscillators with the same frequencies or two autonomous nonlinear oscillators with the same parameters and initial conditions. In both of these cases, the oscillations produced by the two systems will be coherent. However, since the state of one oscillator does not depend on the state of the other, they are not coupled.

Despite this difference, there is still a strong connection between coherence and specific types of synchronization. The states of complete 1:1 synchronization or 1:1 phase synchronization are more or less the same as coherence as the strength of the interaction reduces to a small value when two oscillators are completely synchronized. One can also consider indirect synchronization, such as two non-autonomous oscillators becoming synchronized via the same time-dependent modulation. In each of these cases, the effect can be measured directly using coherence.³⁹ 

In order to illustrate the factors affecting the measurement of coherence, we consider a pair of time-series, which contain common oscillations generated by non-autonomous systems with independent perturbations. To ensure that we are not biased toward perturbations in amplitude or phase, we consider a system with a separable amplitude and phase dynamics.

The important feature of this system is that $r$ and $θ$ vary independently. This means that the amplitude of the oscillator can be perturbed without affecting the phase and vice versa. However, comparing the effect of amplitude and phase perturbations this way using the current system would not be a fair comparison since $r$ has a stable point attractor while $θ$ does not. This leads to the perturbations to $r$ being suppressed over time, while perturbations to $θ$ are integrated over time.³⁷ 

To generate each time-series, the amplitudes ${ r 1, r 2}$ and phases ${ θ 1, θ 2}$ of two modified Poincaré oscillators were numerically modeled and summed together in a time-series $X(t)= r 1(t)cos⁡( θ 1(t))+ r 2(t)cos⁡( θ 2(t))$⁠. However, even with perturbations, this time-series would appear as two noisy sinusoids with approximately stationary dynamics. To simulate more realistic nonstationary time-series, the system was made non-autonomous by modulating the oscillator frequencies with $ω(t)=2π ω 0+Asin⁡(2π ω mt)$⁠. To investigate the effect of phase differences, the phase offset of the oscillations was also adjusted by changing the initial value of $ϕ$⁠.

In the numerically modeled examples used in Sec. IV, we considered a high-frequency mode with $ω 0=1$⁠, $ω m=0.008$⁠, $A=0.8$ and a low-frequency mode with $ω 0=0.5$⁠, $ω m=0.0055$⁠, $A=0.2$⁠. For the other oscillator parameters, we used $a=1$ and $α=5$ in each case.

Noise plays a significant role in the evaluation of coherence. Consider two time-series with a single, identical sinusoidal oscillation with frequency $f sin$⁠. By analyzing Eq. (3), we can see that $S a b( f sin)≡ S a a( f sin)$ and $S a a( f sin)≡ S b b( f sin)$⁠, which results in the expected value $C( f sin)=1$⁠. However, since the time-series contain no other oscillations, this relation holds true not just for $f sin$ but for all values of $f$⁠. This means that we might mistakenly believe that coherent oscillations exist at all frequencies.

Similar behavior is apparent whenever dominant oscillations are present in both time-series. Without independent fluctuations at adjacent frequencies, significant coherence will be observed at values far from the frequencies of the corresponding oscillations.

In most real data, this is not an issue as they are usually influenced by both system noise and measurement noise. We must, therefore, take care to approximate real-world examples in our analysis by including noise in the numerical model.

For the analysis of coherence of phases of oscillations in time-series, the Fourier-based definition of coherence is perfectly valid when the time-series are stationary. However, for multi-scale, nonstationary time-series, the dynamics cannot be approximated by assuming a constant time-averaged phase and amplitude, as is assumed in the Fourier transform. As discussed in Rowland Adams et al.,²¹ such time-series must not be analyzed from the infinite-time, non-time-evolving framework of analysis that is designed for stationary time-series—which is precisely the framework within which Fourier coherence exists—but rather, such time-series need to be analyzed from within the framework of time-evolving time-localized analysis of oscillatory characteristics.

Accordingly, it is natural to seek a way to compute coherence from time–frequency representations of the data. As already mentioned, we can compute a time–frequency representation using an ordinary Fourier transform with a moving window, which is also known as a short-time Fourier transform. However, as soon as we do this, we must ask what size of window? A large window gives us excellent frequency resolution, but then it is more difficult to determine the time at which oscillation frequencies change. Similarly, while a small window enables us to track the change in frequency more precisely, the frequency resolution is lower and makes it difficult to determine the exact frequencies of oscillations. These characteristics of the measurement of waves are well known in quantum mechanics and famously summarized in the Heisenberg uncertainty principle.

The main limiting factor in the choice of window size is the lowest-frequency oscillation that we wish to observe. It is necessary to choose a window that contains enough cycles of this oscillation to determine its frequency to reasonable precision. However, this window size is larger than the window needed to have the same frequency resolution for higher-frequency oscillations. For higher-frequency oscillations, this window size will represent a slower timescale than the timescale of these oscillations, making the analysis effectively equivalent to the kind of long-time-averaging associated with the classical non-time-evolving, long-time-asymptotic-statistics framework designed for stationary time-series described above.

Therefore, to achieve a time-localized analysis of multi-scale time-series, we would need to use an adaptive window size to increase the time resolution at high frequencies while maintaining an optimal frequency resolution overall.

The difference between this time-localized approach and the slow-timescale averaging that takes place in the fixed-window-size approach is illustrated in Fig. 1. Here, time–frequency analysis is performed on a time-series from the illustrative Poincaré oscillator model. In this case, the oscillators were not perturbed with phase noise, and only minimal amplitude noise, $ξ r=0.005$⁠, was introduced. In addition, background fluctuations were numerically modeled by adding independent $1/f$ noise to each time-series. In Fig. 1(a), depicting the fixed-window approach, the idea is to characterize all aspects of the dynamics at a given time using the data in a given window. This means that all of the analysis for every frequency is performed within the same window (note that this window is shown as rectangular for illustrative purposes only—a Gaussian window was used in the short-time Fourier transform to enable a fairer comparison of the two approaches).

By contrast, as depicted in Fig. 1(b), the time-localized approach uses a variable-sized window depending on which frequency is being analyzed. For the former approach, where at each time a full-frequency-spectrum Fourier transform is performed inside a pre-specified window, the result is that the time–frequency analysis can be optimized around one frequency only. However, in the time-localized approach, the analysis is centered around each frequency under analysis, much like adjusting an optical focus. This means that the time–frequency plots for this latter approach provide much greater detail across time at high frequencies, as well as much greater detail across frequency at low frequencies. An alternative version, with a Fourier transform presented with a logarithmic scale, is provided in Fig. 1 of the supplementary material. When comparing the two figures, it is obvious that a logarithmic scale is disadvantageous for the Fourier transform, which is calculated with linear frequency resolution.

It is worth noting an issue that arises from the fact that the integral in Eq. (7) is bounded between 0 and $T$⁠. This means that when $t$ is close to one end of $x(t)$⁠, a significant part of the amplitude of the wavelet function extends beyond the bounds of the integral. This bounded integral is also equivalent to an unbounded integral where the ends of $x(t)$ are padded with infinite zeros. This problem is common among methods using a moving window and other strategies include using reflected data or predicted data equal to half the length of the window. However, each of these methods causes boundary effects that result in errors in the time–frequency representation.⁶⁸ The other alternative is to not include these regions in the plot. This results in a cone of influence, which is larger in size at lower frequencies due to the larger-sized wavelets reaching the ends sooner than smaller wavelets.

From the wavelet transform, one can extract an instantaneous amplitude and phase associated to each frequency $f$ at each time $t$ by expressing $W(s,t)= |W(s,t) | e i θ ( s , t )$ and taking $|W(s,t) |$ as the amplitude and $θ(s,t)$ as the phase.

The application of a wavelet-based approach vs a Fourier-based approach has a significant effect on the information provided by coherence analysis. This can be seen by comparing the studies of Karavaev et al.⁷¹ and Mizeva et al.,⁷² both of which consider cardiovascular time-series recorded over similar timescales (15 and 20 min, respectively). In the former study, the macroscopic autonomic control is characterized by dividing the Fourier coherence into a “high-frequency” (0.15–0.4 Hz) and “low-frequency” (0.05–0.15 Hz) band. In the latter study, the wavelet coherence is divided into five separate frequency bands with ranges 0.6–2, 0.145–0.6, 0.052–0.145, 0.021–0.052, and 0.0095–0.021 Hz, which allows for the characterization of both the macroscopic and microscopic dynamics. The logarithmic scale provided by the wavelet coherence, therefore, acts much like a telescope or microscope, allowing us to zoom in and out of all frequencies of interest at every moment in time.

This definition of coherence was developed independently by Lachaux et al.³⁸ (where it was termed single-trial phase coherence) and Bandrivskyy et al.⁴⁶ While Eq. (13) defines PC for a pair of time-series, it has since been extended to groups of three or more time-series.^73,74

Like Fourier coherence, both PC and AWPC take values between $0$ and $1$⁠. Note, however, that for oscillations with time-dependent characteristics, strong coherence will not typically manifest as a coherence value of $1$⁠, but often as distinctly less than $1$⁠.

In the examples shown in this paper, PC was calculated using MODA—an interactive MATLAB toolbox.⁷⁵ We also encourage readers to consult the MODA user guide, which contains practical information for performing PC and other time–frequency analyses.⁷⁶ 

The differences between PC and AWPC are shown in Figs. 2–4 using the previously defined illustrative Poincaré model. In each case, the two time-series, their corresponding WT, and the PC and AWPC plots are shown. The methods were applied using three different time–frequency resolutions by changing the central frequency $f 0$ of the Morlet wavelet. The effect of adjusting $f 0$ can be seen in the WT, where the frequency width of the bands corresponding to the oscillatory modes is decreased with increasing $f 0$⁠. This effect is also seen for the coherence plots. Here, the darker bands of coherence reveal the common frequency modulation of the two modes, which becomes more localized in frequency as $f 0$ is increased.

An additional effect seen when increasing the frequency resolution is that the background coherence between the modes also increases. The reason for this effect is due to the fact that larger wavelets average over more cycles, leading to extracted wavelet components that are more stationary in frequency. These components, therefore, appear coherent, but only because the rate of change in frequency converges to the same value (i.e., 0) for all oscillations as $f 0$ is increased.

Figure 2 shows the effect of amplitude perturbations on the modes following the two coherence measures. The coherence bands associated with the modes are lighter and less well-defined in the case of AWPC, with the effect being greatest for the lowest frequency resolution. The explanation for this can be found in the independent fluctuations seen in the amplitude of the WT. As highlighted in Eq. (11), AWPC is dependent on the wavelet amplitude, which means that the amplitude perturbations result in lower coherence. In contrast, PC is not dependent on the wavelet amplitude and is, therefore, resistant to such perturbations.

As one might expect, the effect is similar for both approaches when the perturbations are instead applied to the phase of oscillations. Figure 3 illustrates the effect of phase perturbations, where PC and AWPC are affected similarly by the noise due to both methods being dependent on the phase of the wavelet components.

A significant difference between PC and AWPC can be seen in the additive noise case shown in Fig. 4. Here, the common dichotomous noise results in time-localized spikes in the time domain. These can be seen as large cones of amplitude permeating into the lower frequencies in the WT. In the coherence plots, this effect has the most significant impact on the low frequencies, as larger wavelets have a lower time resolution and span across a greater period. Furthermore, it can be seen that the case for $f 0=5$ is most affected by the amplitude perturbations due to the increased temporal width of the wavelets. Generally speaking, therefore, smaller values of $f 0$ should be used in cases where extremely time-localized noise features are present, such as movement artifacts in biomedical measurements.

Also worth noting in the additive noise example is that even though the added dichotomous and $1/f$ noise affect both the phase and amplitude of the wavelet components, the coherence bands of the modes are more strongly defined in the PC plots and the low-frequency coherence is reduced. This is caused by the time-localized properties of the dichotomous noise, which only affect a relatively small number of cycles at each wavelet scale. Since the window used to calculate the coherence averages the phase difference over a relatively large number of cycles, the effect on PC is reduced. In contrast, as shown in Eq. (11), the phase difference in AWPC is weighted by the amplitude. This means that even though the noise spikes last only a small number of cycles, the relative weight to the calculation of the coherence is increased due to the large associated amplitude.

Beyond coherence, it is often useful to extract the instantaneous wavelet phase difference (⁠ $θ a(f,t)− θ b(f,t)$⁠) and analyze this directly. This has been done in many studies to investigate deterministic phase differences in oscillations from two time-series.^77–80 While phase is technically a time-independent measure, the direction and magnitude of the phase difference are still a valuable measure that can be used to determine time lags, which provide weight to statements of causality.

In the studies cited above, analysis of the phase difference involves extracting individual pairs of instantaneous phases and examining the change in the phase difference over time. However, in time-series containing many modes, it is often useful to analyze the phase differences in the frequency domain. Doing this reveals the phase relationships present across different timescales of the dynamics.

Note that while this definition of the time-averaged phase difference correctly identifies the phase differences of the coherent modes, it does not necessarily provide a meaningful value for areas of zero coherence. This is because the result will be the argument of the sum of random phasors. While the amplitude of this sum correctly gives a value of the time-averaged PC at the background level, the argument will be a random angle between 0 and 2 $π$⁠. It is, therefore, important to assess such a measure of the phase difference in conjunction with the actual coherence and only to evaluate its values where the coherence is significant.

In all four cases, we have given formulas according to the assumption that the WT is defined over the whole of $[0,T]$⁠. Once again, this requires that padding has been added to the signal before time $0$ and after time $T$⁠; if, instead, the WT has been computed only over the cone of influence, then the integrals $∫ 0 T$ or averages $1 T ∫ 0 T$ taken over the time-interval $[0,T]$ in the above formulas need to be taken instead over the $f$-section of the cone of influence.

In this paper, we work with the former of the two approaches, namely, Eqs. (15) and (16).

Coherence analysis is restricted by the properties of the measured data. Each dataset is likely to contain idiosyncrasies that require specific attention to avoid false representation of the results. By unlocking the temporal dimension with time-resolved analysis methods, one may properly view and assess the type of data under investigation, and once this step is completed, choose and perform the appropriate analysis. The multi-scale nature of the present analysis also enables simultaneous observation of the behavior across a number of frequencies, which in many cases are representative of various independent behaviors in the system. A review of the statistical properties of wavelet coherence is provided in Cohen and Walden.⁸¹ However, here, we focus on the practical implementation and application of these methods.

To demonstrate the nuance required when selecting parameters for analysis, we consider two sets of time-series containing two common modes. As before, the modes are generated using the modified Poincaré system and have independent perturbations. The key difference is that the first set of time-series has modes with frequencies $ω 0=1$ and $ω 0=0.2$ that are stationary in time, with $ω m=0$ (the leftmost set of Fig. 5). In contrast, in the second set of time-series, the frequency of the modes varies with $ω m=0.016π$ for the high-frequency mode and $ω m=0.010π$ for the low-frequency mode (the rightmost set of Fig. 5).

Importantly, when considering coherence between simultaneously measured time-series, one may use two sets of apparatus with varying sampling frequencies, $f s$⁠. For the calculation of coherence, a common $f s$ must be established. While it is theoretically possible to up-sample the data series with the smaller sampling frequency, this is not recommended as it will not recover information regarding higher-frequency oscillations. Instead, the solution is to downsample the larger time-series so that a common $f s$ is established.

The value of $f s$ determines the maximum observable frequency, $f max$⁠, because we need at least two points in each cycle to capture an oscillation. Consequently, the upper-frequency limit, or the Nyquist frequency, is defined as $f N= f s/2$⁠. A low value of $f max$ can introduce problems when assessing data, as seen in Fig. 5(j). In this case, the system was simulated with $f s$ = 4 Hz, which means that $f max$ = 2 Hz is selected. The coherent mode seemingly passes above $f max$⁠, illustrating the need for a higher $f s$⁠.

The lowest attainable frequency, $f min$⁠, is determined by the length of the time-series. In the examples demonstrated in this work, AWPC and PC are evaluated across ten cycles of oscillation at a given frequency. It follows that the length of the time-series restricts $f min$ and that the length must be at least ten times the length of the minimum frequency of interest. If the interaction is time-varying, then more cycles are needed to account for the modulation present, dependent upon the frequency of the modulation. The time-varying example shown in Figs. 5(c), 5(d), 5(g), 5(h), 5(j), and 5(l) demonstrates a situation where the simulated mode may be interpreted as being centered upon a greater frequency (0.25 Hz) than it really is. Specifically, the mode should be centered upon 0.2 Hz. Due to the shortness of the recording, the cone of influence contains only the upper half of the modulation cycle, resulting in an apparently higher value. In the non-time-varying frequency case, there is no issue, and the peak coherence is centered around 0.2 Hz.

The presence of oscillatory dynamics can be confirmed by first considering the time–frequency representation of the data. In addition, this step will provide information on the frequency range of interest if this is not known beforehand. Limiting the coherence analysis to this range will reduce the burden on computational capacity and save time. The WT will guide the choice of the resolution parameter. However, one must consider that this is always a trade-off, as discussed in Sec. IV C and seen in Figs. 2–4.

The considerations outlined above will help to reduce false conclusions regarding the data. However, to further reduce the chance of falsely representing spurious coherence as significant, a further step must be performed.

Even with the existence of independent fluctuations in both time-series, the interpretation of coherence is not straightforward, as illustrated by Holm.⁸² This is because even two completely independent noise time-series will contain fluctuations that appear at the same time and frequency, resulting in a non-zero value of coherence.

We must, therefore, determine whether observed coherence is significant. This is necessary both for being able to make physical inferences from the observation of coherence values and for being able to make physical inferences from phase-shift values $ψ(f)$ associated with high coherence. Consideration of significance of coherence values can be divided into two aspects: First, the coherence values themselves need to be statistically significant in terms of exceeding some critical threshold, i.e., some baseline coherence value. Second, when one computes the time-averaged phase difference $ψ(f)$ as a function of $f$⁠, where there is significant coherence, one should observe a plateau—i.e., an approximately constant phase difference—over the frequency range in which the phase-coherent oscillations manifest in the time–frequency representation. One should only regard coherence as significant if it is found to satisfy both of these aspects of testing for significance.

In regard to the first aspect, defining the baseline coherence value for significance is not trivial, as it is dependent on the nature of the background dynamics generated by the system under investigation. For example, in the system described above, the independent fluctuations generated from perturbations to the phase and amplitude will result in a different level of background coherence to the case of independent $1/f$ additive noise. Furthermore, in real systems, the deterministic dynamics cannot be separated from the noise perturbations, which increases the difficulty of defining a coherence baseline.

A more formulaic approach is to use a hypothesis test. Specifically, we would like to test a null hypothesis that two time-series are not coherent at a specific frequency. Such a hypothesis can be tested through the use of surrogate data.⁸³ Surrogate data are numerically modeled time-series that are designed to preserve all features of the measured time-series apart from the feature under investigation. In this method, a set of surrogate time-series is first randomly generated. The same analysis that is performed on the real time-series is then performed on the surrogates, with the end result being the discriminating statistic corresponding to the factor of interest. This results in a distribution of values for these statistics, which can then be used to define a specific confidence interval (i.e., the value of a percentile) for discerning significance and rejection of the null hypothesis.

The optimal percentile to use in the test varies from case to case. This can be due to a number of factors. For example, a high intensity of the difference between the noises affecting the two time-series will decrease the coherence between the two time-series to a greater extent than it would decrease the coherence between surrogates, making a lower percentile for the surrogate threshold more appropriate. In this paper, we will adopt a 95th percentile threshold for most cases. However, in some cases, due to factors like the one we have just mentioned, we will use a lower threshold.

One of the most common uses of surrogate data is to test for nonlinearity, where it is possible to apply methods, such as amplitude-adjusted Fourier transform surrogates, that preserve only the linear statistical properties of the time-series (see Ref. 84 for a review of surrogate data methods). However, in testing for significant coherence, we must also preserve the effects of nonlinearity in the surrogate data. Otherwise, even if the surrogates preserve the linear statistical properties, such as the amplitude probability distribution and the frequency spectrum, the null hypothesis may still be spuriously rejected due to increased coherence resulting from nonlinearity.

Mismatched surrogates, also known as intersubject surrogates in the context of biomedical data, are one of the simplest ways to preserve potential nonlinearity in the surrogate data. With this method, pairs of real measurements of the same system (such as the human body, measured across different subjects) are separated and then re-paired with the corresponding time-series from an independent measurement (i.e., another subject). This has the advantage of preserving all properties of the time-series apart from the time-specific information. However, coherence is not preserved as the oscillations are no longer ordered in time.

While mismatched surrogates usually apply only to measured data from real systems, it is still possible to generate time-series approximating mismatched surrogates with the illustrative model defined in Sec. III. In this case, we can simply modify the frequency modulation of the two modes, $ω(t)=2π ω 0+Asin⁡(2π ω mt+ψ)$⁠, where $ψ$ is a phase offset of the modulation. Each pair of surrogate time-series is then generated using different values of $ψ$ for each mode, which are uniformly sampled on the interval $[0,2π]$⁠.

It is also worth noting that surrogate testing is not the only method for determining significance thresholds for coherence values. The method proposed by Sheppard et al.⁸⁵ provide analytically derived significance thresholds based on higher-order statistics, which was shown to give better performance than amplitude-adjusted Fourier transform surrogates.

The effect of time-averaged surrogates is illustrated in Figs. 5(k) and 5(l), which show the 95th percentile of 99 mismatch surrogates. These surrogate thresholds give a much clearer indication of the coherence values that are present in the system vs the spurious coherence. One may also choose to illustrate the time-localized effective coherence. This is demonstrated in Fig. 6, with parameters identical to those in Figs. 2(e), 2(f), 2(k), and 2(l). The threshold here was chosen as the 75th percentile of 99 mismatch surrogates. One can now discriminate the coherence due to the modes vs the background fluctuations in the time-averaged coherence. However, many areas of significant coherence still remain in the time-localized plot distributed away from the modes. This illustrates the fact that it is easier for spurious significant coherence to occur in the time–frequency domain, where the testing area is essentially squared.

The other effect of surrogates can be seen on the effective coherence of the low-frequency mode, which is much reduced compared to the high-frequency mode. This is due to the fact that spurious coherence between random fluctuations is more likely to be found since the average coherence is calculated over fewer cycles. This essentially reduces the observable frequency range, adding to the effects already caused by the size of the wavelets (parameterized by $f 0$⁠) and the window size used for the coherence calculation (parameterized by $n c y$⁠). Taking into account these cumulative effects, we generally recommend that effective coherence can only be assessed if a minimum of 30 cycles can be observed, giving the lowest observable frequency of $30/T$⁠.

Now, to illustrate the second aspect of considering significance of coherence: In the two frequency bands where Fig. 6(b) shows coherence values exceeding the surrogate threshold, Fig. 6(c) shows the phase difference plateauing at about 0.75 $π$⁠. These plateaus in conjunction with the statistical significance of the coherence values suggest that the coherence in these two frequency bands is significant. Moreover, as a consequence, we can conclude that the value 0.75 $π$ around which the phase difference plateaus is the amount by which the first time-series leads the second, consistent with the numerically modeled input values.

Therefore, we have seen that the surrogate threshold and the phase difference are invaluable tools when interpreting coherence; this will be demonstrated in Sec. VI via a series of examples.

The heart rate is modulated through several processes, with respiration being an important factor. During inhalation, the heart tends to beat quicker, and during exhalation, it tends to slow down. This interaction is known as respiratory sinus arrhythmia.⁸⁶ Cardio-respiratory interactions are perhaps one of the most widely studied interactions. Several methods have been employed,^87,88 including coherence analysis based on the Fourier and wavelet transforms.^89,90 Utilizing PC to study cardio-respiratory interactions has also proven valuable,⁵⁰ for example, in the context of ageing,⁵² malaria,⁹¹ and hypoxia.⁹² 

In this example, we evaluate cardio-respiratory interactions based on the simultaneously recorded respiratory effort and the electric activity of the heart. The 1400 s recordings are taken from a 28-year healthy male participating in the study of ageing,⁹³ where the sensor/electrode placements are described. A time-insert of respiration is shown in Fig. 7(a) and the ECG in Fig. 7(b). The instantaneous frequencies of respiration [IRR, Fig. 7(c)] and beating of the heart [IHR, Fig. 7(d)] are extracted by ridge extraction³² after the WT was obtained. Two types of interactions are investigated: (a) between the original respiratory time-series and the IHR and (b) between both instantaneous rates, IRR and IHR. The PC and AWPC for both cases are shown in Figs. 7(e), 7(f), 7(i), and 7(j). The surrogate threshold was set to the 95th percentage of 140 intersubject surrogates, as used in the original study.⁹³ The time-averaged values of PC and AWPC from the entire 1400 s recordings are shown in Figs. 7(g) and 7(k) for the cases (a) and (b), respectively. The phase differences, as a function of frequency, obtained for case (a) and (b), are shown in Figs. 7(h) and 7(l).

It is clear that both PC and AWPC are much higher for the respiration-IHR case, compared to IRR-IHR case, and that the highest values of coherence are at the frequency of respiration (around 0.2–0.3 Hz), consistent with earlier studies. This indicates that, in the resting state, the heart rate is strongly modulated by the amplitude of respiration and to a much lesser extent by the frequency of respiration. In Fig. 7(g), one can see that the PC and AWPC are similar. The phase difference at the respiration frequency is around 0 rad.

Coherence analysis is often applied to find common oscillatory behavior between brain signals from different locations. This can elucidate the functional connectivity of the brain, which is known to change in various conditions.^94,95 Spontaneous activity in the brain can be measured noninvasively at a relatively low cost using EEG or fNIRS with minimal discomfort to the subjects. However, both methods are susceptible to movement artifacts.⁹⁶ Several approaches exist to remove these artifacts from the data, although they often compromise the quality of the data and may additionally remove information of interest.^97,98 As seen in Secs. IV, phase-based approaches may be more resilient against movement artifacts and noise and, as such, can circumvent some of the more draconian preprocessing requirements. In this section, we investigate two examples of movement artifacts, one using EEG and the second using fNIRS.

Non-invasive brain activity measurements in children are fraught with artifacts due to difficulties in keeping younger subjects still for extended periods. Analysis of signals derived from younger cohorts, therefore, necessitates methods that are robust to movement artifacts. In addition, when considering the presence of interactions between time-series, it can be important to assess how the nature of these interactions changes over time. Time-localized methods can reveal temporal dependencies in this mutual behavior. In a wide array of neurological conditions, it is not only the intensity of interaction between brain regions but the duration of interaction that is altered.^99,100 By observing the time-localized coherence, one may deduce the regularity and strength of time-varying interactions.

We consider a resting-state measurement with eyes open of two simultaneously recorded EEG time-series. These data were measured in a cohort of male children aged 3–5 years old with a diagnosis of autism spectrum disorder (ASD). The time-series were captured using a Nicolet cEEG instrument (Viasys Healthcare, USA) at a sampling rate of 256 Hz. A 20-min recording period was used to collect the data, and a 180-second interval was analyzed, with the central 60 s illustrated in Fig. 8 as it contained a clear artifact. Measurement sites corresponding to F3 and F4 in the international 10–20 system were chosen, as the initial objective of the investigation was to assess reports of reduced frontal connectivity in children with ASD.^101–104 

The effect of the movement artifact is clearly seen in both the time domain, Figs. 8(a) and 8(b), and the WT, Figs. 8(c) and 8(d), where at the instance of the movement, all frequencies are present (around 93 s) in the spectrum. The effect on the coherence is much stronger and can be seen in Figs. 8(e)–8(g) for the AWPC compared to the PC. A threshold of the 75th percentile of 156 intersubject surrogates was used, leaving only the significant coherence. The time-localized coherence [Fig. 8(e)] shows that the magnitude and presence of the interactions vary over time. Both the time-localized, Figs. 8(e) and 8(f), and the time-average, Fig. 8(g), coherence are elevated for the AWPC compared to the PC.

Now, we consider two time-series recorded from the temporal brain areas, in a study that investigated coherence between neuronal and vascular function.¹⁰⁵ These locations often have artifacts due to movement of the jaw. The data are from a participant with a positive genetic test for Huntington’s disease (HD), who has not yet developed the movement disorder known as chorea. Still, as chorea is a hallmark of the disease, HD research would benefit from methods that are resistant to movement artifacts.

We compared PC and AWPC of two resting-state oxygenated hemoglobin (oxyHb) time-series measured using a fNIRS device (NIRScout, NIRx, Germany) with a sampling frequency of 31.25 Hz over 20 min (for further details on measurements, see Ref. 93). The measurement sites correspond to T7 and T8 in the international 10–20 system (left and right temporal locations). The resolution parameter $f 0$⁠, Eq. (8), was set to 1, as to minimize the spread of an artifact.

The results are shown in Fig. 9. The time-series contain two movement artifacts, which appear as high-amplitude cones in the WT and have the greatest impact at low frequencies. The artifacts have a very significant impact on the AWPC plot and affect an even wider area of time and frequency than is visible in the WT plots. This is a consequence of the moving window used to calculate wavelet coherence. In the plot of PC. the effect of the artifacts is not obvious. This illustrates how any simultaneous increase in amplitude, even if not phase coherent, results in AWPC appearing significant over large areas of the time–frequency domain. This can also be seen in the time-averaged coherence plot, where the AWPC (orange line) is much higher than the PC (black line). The two dashed lines show the 95th percentile of the 136 intersubject surrogates.

This example illustrates that PC is relatively resistant to artifacts, which is beneficial when analyzing time-series from various non-invasive measurement techniques.

Time–frequency and coherence analysis can provide valuable information about the dynamics of a system. In addition, the phase difference between oscillations can give information about the direction of influence. We consider the movement of electrons on the surface of liquid helium at very low temperatures, as discussed in Siddiq et al.³¹ At very low temperatures, the helium will be a superfluid. Since such a system can be used for constructing the qubits that are needed for quantum computers, increasing the understanding of its dynamics is important.

In the experiments, the electrons were just above the liquid helium, trapped between the helium and a vacuum. They were in a perpendicular magnetic field and subjected to microwave radiation and varying pressing voltage. Current oscillations were recorded from five electrodes for 60 s at 100 kHz. The full experimental setup is described in Ref. 31. We chose an example with low electron density and 4.18 V pressing voltage. Currents measured from electrodes E4 and C in the time-interval 30–31.4 s were selected for analysis, and high coherence was obtained as in the original paper.³¹ 

Figure 10 shows the PC and AWPC between current oscillations at the two electrodes. The current signals were first downsampled to 20 kHz, as in this example, we will focus on oscillations around 0.5 kHz. The resolution parameter was set to 3, in line with the original paper.³¹ 100 iterated amplitude-adjusted Fourier transform (IAAFT) surrogates were used to calculate the surrogate thresholds.⁸⁴ 

Both PC and AWPC methods pick up a time-varying coherence following the dominant mode in the WT plots, which resembles a non-autonomous influence on the system. The time-averaged PC (black line) and AWPC (orange lines) are similar, with the AWPC having a slightly higher value at the higher frequencies. This could indicate that there is some amplitude covariance. The surrogate thresholds are very similar for both methods. The time-averaged phase difference is positive, meaning that the oscillation at E4 is preceding that at C.

The existence of coherence indicates that the electrons are moving, and the phase difference suggests that they are moving toward the C electrode from the E4 electrode. This is consistent with the microwave radiation being applied closer to E4. Furthermore, by studying the time–frequency representations, we see a clear mode with a time-varying frequency. Using ridge extraction,³² which essentially tracks the maximum amplitude within a frequency range, we can extract a time-series of the instantaneous frequency. The WT of this time-series shows a clear amplitude peak at around 5.2 Hz, indicating modulation of the electron movement at this frequency. This was shown to be caused by slow gravity waves on the liquid helium.³¹ It is important to note that in the WT of the original current data, there is also a peak at around 5.2 Hz. However, this peak is relatively weak compared to the rest of the spectrum, and, in particular, compared with the dominant oscillatory component, making it challenging to observe and identify directly from the frequency spectrum. This illustrates that time-localized, time–frequency methods can uncover a great deal of physically meaningful information.

The study of coherence has its foundations in physics, where methods were first developed to measure the coherence between the phases of waves. It has then been extended to considering coherence between the phases of more general oscillatory processes occurring in a wide variety of scientific disciplines; for this, one of the most fundamental issues is the quantification of such coherence from measured data. Accordingly, it is a subject particularly treated by harmonic analysis in mathematics and by signal-processing theory. We have approached this same question again from a physics perspective, but still with this greater generality than the kinds of setup that initiated the study of coherence—namely, from the perspective of multi-scale time-dependent oscillatory dynamics.

We have seen that for time-series data recorded from systems involving interacting oscillations, key information about the interactions is contained in the time evolution of the phases of the oscillations. Moreover, we have seen that for the analysis of systems involving oscillations with temporally modulated quantitative characteristics, such as frequency and amplitude, time-series analysis methods that are fundamentally designed for time-series with stationary statistics are inappropriate. For example, the measure of coherence of phases intended to be revealed by Fourier coherence will have little meaning for systems with frequency-modulated oscillations. Rather, tools designed to extract time-evolving, time-localized information about systems exhibiting time-dependent far-from-equilibrium dynamics are needed. In particular, phase information needs to be extracted in a suitably time-localized manner.

Such time-localization inherently needs to be understood relative to the timescale of the dynamical process under investigation, and therefore, for multi-scale time-series involving oscillations of a range of frequencies, this time-localization needs to be adaptive to the range of timescales involved. This has been illustrated in Fig. 1.

In the setting of time–frequency analysis, where the Heisenberg uncertainty principle requires a trade-off between precision in measurements of frequency and precision in location in time, this multi-scale adaptivity corresponds to a logarithmic frequency resolution. This is achieved by the continuous wavelet transform, where the scale variable is taken as inversely proportional to the frequency under investigation. Accordingly, we have seen that the wavelet transform is the appropriate tool for extracting phase information from multi-scale nonstationary time-series, and in particular, WT-based coherence analysis is the appropriate approach to investigating coherence of phases manifesting in such time-series.

In this paper, we have provided an introduction to wavelet-based coherence analysis and evaluated several related issues—some already established and others that had not previously been addressed.

Existing discussions of different approaches to quantifying coherence, and more generally of different approaches to time–frequency analysis, have mainly treated the different approaches as if on essentially equal footing, and practical choices, such as the use of WT over STFT, or of PC over AWPC, are often treated merely as a matter of quantitative optimization, without theoretically reasoned or experimentally explored consideration of the qualitative impact of such choices on the resulting analysis.

To address this issue, we have provided a systematic explanation of the practicalities and pitfalls of how to carry out wavelet coherence analysis in practice. In particular, we have provided a detailed review of the methodology for reliably testing for and detecting significant phase coherence from measured data.

Since the wavelet transform provides extractable phases and amplitudes, one can calculate^38,46 a measure of coherence independent of changes in the amplitude, namely, PC, as well as a measure of coherence that is weighted in time by amplitude, namely, AWPC. Prior to this work, an in-depth comparison between AWPC and PC had not been performed. Perhaps counterintuitively, we found a consistent difference in the performance of the two definitions of coherence when applied to noisy time-series. PC is, in general, more robust to noise and particularly to time-localized perturbations, meaning that it is affected to a much lesser degree by phenomena, such as movement artifacts.

Along with the definition of PC, one can also analyze phase differences in the oscillations present in the pair of time-series under investigation. This is first needed as one of the aspects of determining significance of coherence, along with statistical significance of the coherence values themselves: the time-averaged phase difference as a function of frequency needs to have a plateau in the frequency band where coherent oscillations manifest in time–frequency representation. Second, where there is significant coherence, the phase difference can provide indications of which oscillation is leading.¹⁰⁶ However, it is important to note that causality (i.e., which process is the origin of the common oscillations) is not always possible to infer from the phase shift. This can be because the phase shift is wrapped on the interval $[−π,π]$ or because of the existence of higher degree interactions, such as node triples.¹⁰⁷ 

This investigation of coherence has revealed the wealth of information provided by the phase. Part of the utility of phase over amplitude comes from the fact that phase dynamics is constrained by the frequency interval within which an oscillation lies. For example, each wavelet has a defined frequency response, which limits the rate at which the phase can change. However, in contrast to the phase, the amplitude is not bounded to frequency in such a manner, making the separation of amplitude dynamics from noise harder to satisfactorily achieve. This is analogous to the advantages of frequency modulation over amplitude modulation in radio communications.¹⁰⁸ 

The practical implications of the main points highlighted by our work are effectively illustrated in our analyses of real data in Sec. VI. In the examples shown in Sec. VI B, the presence of movement artifacts has a significant effect on the results of any analysis. Such artifacts usually need to be removed prior to analysis, which often requires subjective identification. The techniques used to remove identified artifacts may also introduce artificial manipulations in the data.¹⁰⁹ The resistance of amplitude-independent phase-based methods to these sorts of artifacts allows for the analysis of noisy data without the need for preprocessing or constraints on the measurement setup, leading to better research into conditions, such as HD and ASD.

In the example of electron dynamics on the surface of liquid helium, using phase coherence analysis, we identified the existence of gravity waves. Without the time-localized approach, these waves might not have been detected. This illustrates how the application of coherence methods and time–frequency analysis can be used to identify specific properties of a physical system. Thus, we see the importance of using a time-localized approach instead of an asymptotic approach, i.e., infinite-time, non-time-evolving approach. The wider adoption of explicitly finite-time and time-localized methods should, therefore, lead to similar discoveries in systems characterized by non-autonomous dynamics involving nonstationary amplitudes and frequencies.

We review the current understanding of coherence, a universal phenomenon that can appear between oscillations or waves, irrespective of their origin. We start from its definition in physics and review numerical methods for analyzing coherence from modeled or real data. We focus particularly on coherence between non-autonomous oscillatory processes whose oscillations have deterministically time-varying frequencies. For this reason, we apply wavelet-based, time-resolved coherence analysis. We discuss differences between amplitude-weighted phase coherence and phase coherence. In the former case, time-resolved information includes both the amplitude and the phase; in the latter case, only the information about phase is considered. We illustrate that the amplitudes are more readily perturbed than phases by noise or movement artifacts, and consequently, that phase coherence provides more robust information about interacting oscillatory systems. We illustrate this in relation to several real-world examples.

An alternative to Fig. 1 in the main paper is presented in the supplementary material. It shows the Fourier transform in a logarithmic scale. By comparing Figs. 1(c), 1(d), 1(g), and 1(h) in the main text and supplementary Figs. 1(c), 1(d), 1(g), and 1(h), it is clear that the logarithmic scale is disadvantageous to the short-time Fourier transform, which is obtained with linear frequency resolution.

We are grateful to Peter McClintock and Kostya Nasyedkin for useful comments on the paper. The work of P.T.C. and A.S. is funded through the Sony Research Award Program. A.S. is also funded by the Engineering and Physical Sciences Research Council, UK (Grant No. EP/X004597/1). J.B. is supported by the Sir John Fisher Foundation. S.J.K.B. is jointly supported by the EPSRC, UK and the MyMind GmbH—Brain Hero, Vienna, Austria (Grant No. EP/T518037/1). The High End Computing facility at Lancaster University was used for some of the computations. The experimental part of the ASD study was supported by the Action Medical Research (UK) MASDA Project (GN1963) and partly by the Slovenian Research Agency (Program No. P20232). The HD and ageing studies were funded by the Engineering and Physical Sciences Research Council, UK (Grant No. EP/M006298/1) and the Slovenian Research Agency (ARRS) (Program No. P20232). The experimental data used for the example on electron dynamics on the surface of liquid helium were obtained by Kostyantyn Nasyedkin in the Quantum Condensed Phases Research Team, RIKEN CEMS, Japan, headed by Kimitoshi Kono. The development of the MODA toolbox used for analyses has been supported by the Engineering and Physical Sciences Research Council, UK (Grant Nos. EP/100999X1 and EP/M006298/1), the EU projects BRACCIA (517133) and COSMOS (642563), the Action Medical Research (UK) MASDA Project (GN1963), and the Slovene Research Agency (Program No. P20232).

The authors have no conflicts to disclose.

The authors have been listed alphabetically due to their approximately equal contributions.

S. J. K. Barnes: Conceptualization (supporting); Data curation (equal); Formal analysis (equal); Investigation (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). J. Bjerkan: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). P. T. Clemson: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). J. Newman: Conceptualization (supporting); Data curation (supporting); Formal analysis (supporting); Funding acquisition (supporting); Investigation (equal); Methodology (supporting); Validation (supporting); Visualization (supporting); Writing – review & editing (equal). A. Stefanovska: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal).

The data used in this work are available in the Publications and Research (Pure) portal on Lancaster University’s research information management system. The data used in the cardio-respiratory example can be found at doi.org/10.17635/lancaster/researchdata/630. The data used in the ASD example can be found at doi.org/10.17635/lancaster/researchdata/604. The data used in the HD example can be found at doi.org/10.17635/lancaster/researchdata/631. The data used for the electron dynamics on the surface of a liquid helium example can be found at doi.org/10.17635/lancaster/researchdata/655.