I Introduction
The need to identify the existence (or lack) of second order nonlinear interactions in dynamical systems is a widespread problem in numerous disciplines. Examples can be found in population dynamics [1], geophysics [2], oceanology [3] or in plasma physics [4]. The basic idea is to investigate the threewave coupling in a signal, or more accurately, to characterize the fraction of the signalenergy of two waves that is quadratically phase coupled to a third wave at the sum frequency. The data analysis methods for stationary (on the scale of a single investigated time slice [5]) processes are well developed [6, 7]. However, more care is necessary when working with nonstationary processes and transient signals. There have been a number of previous approaches to extend the validity of the bicoherence method to rapidly changing systems. However, all of these required quasistationary behaviour on some short timescale. Good examples are the waveletbicoherence [8]
, the shorttime Fourier transform based bicoherence
[9, 10], or the method of appying adaptive windowing to select an optimal window length with maximized local coherence [11]. All of these, however, require quasistationarity on a short time scale, which is not always the case in various reallife systems.This paper introduces an approach to characterize nonlinearity of nonstationary systems based on the “classical” bicoherence technique [6, 7]
, complemented by statistical analysis. The goal is to estimate the significance level of the measured bicoherence for each frequencyfrequency point, in order to filter likely false positives (type I error), i.e. values of high bicoherence despite the lack of actual phase coupling. The method is based on the estimation of the
random bicoherence probability density function
. A similar Monte Carlo based method has been used earlier to estimate the significance level of wavelet spectra by Torrence et.al. [12].Section II describes bicoherence calculation for stationary processes based on Kim et.al. [6], and exposes the problem with nonstationary systems. In section III a model for understanding the statistical properties of bicoherence calculations is introduced, and we present a method to detect false positives by estimating the significance level of the measured bicoherence. A validation using simple model systems follows in Section IV. Finally, implications for lower and higher order coherence calculations are discussed in Section V.
Ii Bicoherence for stationary processes
Bicoherence calculation is a standard method to investigate second order (i.e.: quadratic) nonlinearities, which appear as phase coupling between different Fourier components of a signal. We define the coupling condition with the following (1)(2) equations [6]:
(1)  
(2) 
where () is the frequency and is the phase of the corresponding Fourier components of the signal. Quadratic nonlinear coupling occurs between and , generating a third component at the sum frequency . Throughout the paper we are going to work with Fourier components of real (i.e. not complex) signals, and note these Fourier components with capital letters. We approximate the analytic Fourier transformation with the following expression:
(3) 
where
is the imaginary unit vector, in the second line
defines the locations of independent time windows. This definition assumes an infinite long signal from which infinite number of independent realizations of the spectrum can be calculated. For infinite long stationary signals (originating from ergodic processes) the time average is identical to the average of the independent realizations. In reality, infinite long signals do not exist, therefore we approximate the transform with a ”sufficient” number of averages:(4) 
This means that we divide the measured signal into pieces of long time slices, and carry out Fourier transformation separately, before averaging the spectra. We will use ensemble average throughout this work. To emphasise this, we introduce the averaging operator:
(5) 
where marks a derived quantity (e.g.: Fouriertransform) from the th time slice.
In the limit of ergodic stationary processes will lead to the usual definition of the expected value, which we used in equation (II). For the analysis of nonlinear quadratic coupling, first we define the bispectrum [6]:
(6) 
where marks complex conjugation. Bispectrum is defined on a frequencyfrequency plane. The procedure by which it is estimated can be visualized on the complex plane, illustrated in figure 1.
In the case of phase coupling between components , the absolute value of the bispectrum will be large; otherwise its value will tend to zero with increasing number of averages. Bicoherence is then defined as the bispectrum normalized in the following way:
(7) 
With the definition in eq. (7), , and can be interpreted similarly to classic coherence [6]. If components are coupled (the conditions (1)(2) are met), , else it will tend to zero as the number of averaging tends to infinity.
In the case of signals with finite length the frequency domain of the bicoherence is bounded with the Nyquistfrequency [13], which is shown as the red hexagonal frame in figure 2. Due to the symmetries of the bispectrum (detailed in appendix A) in the case of real signals it is enough to plot an even smaller region, which contains all the independent information (marked with “P” in figure 2). For further analysis we will only present results in this area.
Iia Numerical implementation
Reallife digital signals sampled at a given sampling frequency have finite length (with discrete elements). The calculation of the averages in (6) & (7), based on the definition in eq. (5) requires the signal to be divided it into parts with equal length. The length of a single part and the sampling frequency determine the Nyquist frequency and frequency resolution of the calculations:
(8)  
(9) 
It means that increasing the number of averages can only be achieved at the cost of degrading the frequency resolution, making a tradeoff necessary for any given phenomenon under investigation. From a practical point of view, first we should consider the minimal frequency resolution required for the investigated phenomena, which will then set the maximum number of averages. Splitting the signal into parts (practically windowing with a boxcar window function) leads to sidelobes and other undesired features in the Discrete Fourier Transform. For the suppression of these effects it is beneficial to apply window functions to the smaller parts. Throughout the paper we utilize the Hann window with overlap [14], where the window function is defined as
(10) 
extending for a single period between two zero function values.
Let us now introduce the Discreet Fourier Transformed (DFT) complex vector of the signal’s th part, which has elements, and its frequency resolution is defined by (8). A single element of this vector will be marked as , where and the th element corresponds to the frequency component. (Note that, as far as we are analysing real signals, we consider only the positive frequencies of the transform). Most programming languages support optimized matrix operations, therefore it is beneficial to rearrange our calculations in a matrix form to make calculations faster. First we define the crossfrequency matrix with the following equation:
(11) 
which is the dyadic product of with itself. Then we define the cyclically shifted frequency matrix of the complex conjugate vector elements :
(12)  
where the th row of the matrix is generated by shifting the vector cyclically by , and finally zeroing out all elements outside of the upper quadrant defined by and . Taking the elementwise product of the above defined matrices will produce the bispectrum matrix, and these matrices can be used for the normalization to produce the bicoherence. These statements are summarized in the following equations:
(13)  
(14) 
where denotes componentwise multiplication, and the absolute value, division and square of and are taken componentwise as well at the nonzero elements. Utilizing the matrix implementation on a modern computer can significantly speed up the calculation (the exact speedup depends on problem size, numerical library, and architecture).
Iii Bicoherence for nonstationary processes
The calculations introduced so far were developed to investigate the nonlinear properties of stationary processes [6, 7]. Problems in interpretation may arise when basic assumptions – such as the requirement of stationary input signal – are not fulfilled. First, we demonstrate the characteristic problem of calculating bicoherence by definition (7) for nonstationary processes. In figure 3 we return to the visualization introduced in figure 1, which illustrates the bispectrum calculation as averaging vectors on the complex plane. In the case of nonstationary processes, the length of the steps taken –which depend on the Fourier amplitudes – may vary. In the case of phase coupling, the vectors align in the same direction. Therefore the bicoherence value is high, similarly as for stationary signals, as illustrated in figure 3(a). In the case of a nonstationary signal, the amplitudes in the blocks to be averaged vary significantly. This introduces a bias to the estimate, which may result in high bicoherence values even without phase coupling, as illustrated in 3(b). Therefore it is important to understand the effect of the distribution of Fourier amplitudes on the classical bicoherence calculation.
Iiia Statistical properties of bicoherence
Let us now analyse the statistical properties of bicoherence calculated for a real life (measurable) signal with arbitrary Fourier amplitude distribution. As a simple demonstration we model two processes – one stationary with constant Fourier amplitude in time, the other nonstationary with timevarying amplitudes –, both of these without phase coupling. We then generate the windowed Fourier components (at a selected frequency) of the signals in the following way. For the stationary process the amplitude is constant in all time windows, while the nonstationary process has a uniform random distributed amplitude in the
range. In each time window the phase of both processes is a uniform random variable
. (This realization is analogous to the random walk picture illustrated in figures 1 & 3).In this particular example (shown in figure 4), for both processes blocks are generated for a single realization (to simulate short reallife signals) and then bicoherence is calculated using (7), that corresponds to the specific realization of these random processes. Repeating this process a large number of times (in this example: ) with newly generated phases represents further different realizations of the random processes. With this Monte Carlo method it is now possible to numerically estimate (using a histogram) the random bicoherence probability density function, hereinafter abbreviated as . A different distribution corresponds to each frequency pair on the frequencyfrequency plane (see e.g. figure 2). The describes the bicoherence probability density of a random process with given Fourier amplitudes in the absence of phase coupling. For these two specific test signals described above, the estimated is shown in figure 4.
As we show in figure 4 – for a fixed number of averages – there is significant difference between the characteristics of the two different distributions. For example, in the stationary case the expectation value is much lower than in the case of the nonstationary process. The generated creates the basis for comparing the measured value of bicoherence with the bicoherence distribution of a process with similar amplitude, but completely random phase. In the following we will discuss, how the probability density function can be used to identify possible false positives.
IiiB Application of confidence filtering
Numerically estimating for a given amplitude distrubution gives us the possibility to analyse the bicoherence distribution in the absence of phase coupling. In practice, will act as a reference to which the measured value can be compared to. With this in mind, the bicoherence analysis of the components of an arbitrary (even nonstationary) signal goes as the following:

We calculate the Fourier amplitudes, for the given , components.

Calculate the corresponding measured bicoherence value, using definition (7).

Generate the phaserandomized bicoherence density function using the real life Fourier amplitudes of the the given frequency components with a random phase, using a sufficiently large (e.g.: ) number of realizations.

Define an confidence level (e.g.: ).

Calculate critical bicoherence value:
(15) 
If , then it can be stated that at confidence level does NOT originate from a random process.

If , we assume that at confidence level the bicoherence in that point is not significant.
Naturally, we can carry out the above process for the entire frequency range of interest. The procedure to estimate can be trivially paralellized for all couples.
Steps 3)–7) of the bicoherence analysis procedure are visualized in figure 5. We chose a global confidence level, which defines a critical bicoherence value through the probability density function. In the case shown in figure 5 the measured bicoherence value is higher than the critical one, therefore we accept it at the chosen confidence level as a sign (and measure) of phase coupling.
In section IV we will use different model systems to test the method introduced above, and to demonstrate the effect of different transients on the bicoherence calculation.
Iv Analysis of simple model systems
We chose to simulate a simple physics model system (illustrated in figure 6) in order to validate the suggested method. The model consists of two masses and attached to each other and to the walls with springs, where two springs are ideal and one of the springs was chosen to have nonlinear characteristics described by the parameter .
The equations of motion for the model system are:
(16) 
where . Setting and we can calculate the eigenfrequencies of the investigated system:
(17)  
(18) 
The time signal used in the analysis is the displacement of the first body. We can choose and as input parameters (), and we can calculate the corresponding and constants, using e.g. kg. For convenience we have parametrized the system such that the eigenfrequencies are far enough from eachother at Hz and Hz. Nonlinearities are controlled by varying the parameter. Initial values for the calculations were , , , . The differential equation system was solved numerically with 4th order RungeKutta method, implemented in IDL language. s with a sampling frequency of
kHz was simulated. Additive white noise
, with signal to noise ratio , was mixed to the simulated signal to eliminate possible type divisions when evaluating the (7) bicoherence. Finally, broadband perturbations (bursts) were mixed to the signal. The bursts were created as a sum of Gaussian envelope functions multiplying independent white noise signals, which can be written in the following form:(19) 
The (signal to noise ratio of the white noise multiplier) was chosen to be , and the () characteristic width of a single perturbation was chosen as ms. The signal for the bicoherence analysis is the sum of these three components:
(20) 
We applied Short Time Fourier Transformation [15] to visualize the timefrequency evolution of the generated signal and its components. In the following subsectionsections we will demonstrate the effects of transients on the calculated bicoherence, and discuss the results of confidence filtering.
Iva Linear case with bursts
The first natural test case is one without nonlinearities (), and broadband perturbations. We will use this case to demonstrate the effects of transients on bicoherence calculation. On the spectrogram in figure 7 we can observe the two eigenfrequencies at Hz and Hz, and the broadband perturbations that appear around sec.
The results of bicoherence calculation, using blocks is shown in figure 8 (Only the relevant part of the matrix is plotted, indicated with the yellow dashed triangle in figure 2). We can observe, that even in the lack of phasecoupling there are several frequency points with high bicoherence values. Since we know that our physical system is linear, the main source of these high values are the amplitude perturbations. We now apply the procedure outlined in section III in order to attempt the identification of false positives.
We have calculated realizations with the real signal amplitudes to generate the distribution functions for the entire frequency plane.
Using a confidence level filtering with , we can eliminate any points as false positives, and as shown in figure 9, most of the points disappear, as expected in the lack of phase coupling. (The selected level of is considered high, but as we will show later in this section, “real” bicoherence remains unfiltered even at a high value.)
Some unfiltered points are still left after filtering at any confidence level. The origin of these is the statistical nature of the procedure: at a filtering level of we will make mistakes with roughly probability. A bicoherence matrix for samples will contain elements, therefore we expect false positives at to be unaffected by the filtering procedure in the entire bicoherence matrix. Due to symmetries and plotting limits, in figure 9, we expect false positives, which is in good agreement with the number of points observed in the figure. Furthermore, we have signals of finite length, which will lead to the broadening of the peaks of the Fourier spectrum. Therefore, in the case of phase coupling, we expect to have a broader, contiguous area with high bicoherence, instead of randomly spread points.
IvB Nonlinear case with bursts
The second test case used the same eigenfrequencies ( Hz, Hz) and the same 4 broadband perturbations, but with the added nonlinearity. The parameter was chosen to cause deviation in the force at maximum displacement compared to the linear case. The corresponding spectrogram is shown in figure 10, where again we can observe the broadband perturbations, the basic harmonics and at the sum frequency of these ( Hz) a weak frequency component appears due to the introduced nonlinear coupling.
The result of bicoherence calculation is presented in figure 11. Although in the model system we only define nonlinearity between the basic frequency components, high bicoherence appears at a large number of frequencies. Based on this picture alone, we cannot distinguish false high values from the ones caused by the actual nonlinear coupling.
The filtering process was carried out, as described above, with the result shown in figure 12. Apparently, based filtering eliminated almost all bicoherence values, except a small number of residual points (as explained earlier), and the “true” bicoherence corresponding to the nonlinearity of the system. The Hz frequency component is strongly nonlinear, as evidenced by high bicoherence appearing around the Hz region. We can also see the indication of significant interaction between the the Hz componentand the Hz component, as high bicoherence appears around Hz, at a confidence level of . This example illustrates that the filtering method based on the estimated distribution can distinguish between bicoherence caused by actual phase coupling, and false positives introduced by the nonstationarity of the signal, thereby greatly aiding the interpretation of the results.
V Discussion
The evaluation method of bicoherence for nonstationary signals presented in this paper can be generalized to arbitrary nth order coherence, and to any signal processing technique where timeaveraging of complex spectra plays a role. Particularly, the coherence between two signals characterizing the linear coupling can be calculated this way by substituting the bispectrum to the crossspectrum, and calculating the random coherence probability density function. The main idea is, that if we can formulate ’randomness’ in a system in the quantity in which we are looking for systematic behaviour, then this random probability density can be compared to the actually measured quantity, and draw conclusions on its significance.
It is important to note, that the bicoherence calculation presented in the paper is ,,blind” to higher than second order nonlinearities. In case the second order nonlinearity is forbidden due to symmetry reasons, it might be necessary to study the third order nonlinearity by tricoherence [16]. Also one could aim for a complete description of the system by studying many levels on nonlinear terms at the same time. All of these techniques are based on detecting different higher order phase couplings, and so they are suitable for the type of significance filtering presented in the paper.
Finally, a word on the confidence level. It is straightforward to select high confidence levels in our analysis. However, one should be aware that while this helps to filter out more and more false positives, – which corresponds to eliminating errors of first kind while testing the null hypothesis of the calculated bicoherence being the result of a random process, – this also can cause the elimination of points exhibiting real phase coupling, – which are errors of the second kind. Selecting the right level of confidence level requires the balancing of this tradeoff, and may necessitate detailed MonteCarlo simulations of the partially coupled systems. Choosing a moderate confidence level will cause some false positives to be sustained, which can be discriminated based on their scattered spatial distribution.
Vi Conclusion
The lowest order nonlinearity in physical systems is frequently the secondorder nonlinearity, which leads to phase coupling of signal components of different frequency. The classical bicoherence calculation was originally developed to investigate quadratic, secondorder coupling in nonlinear systems, which are stationary from the time scale of the selected time windows upto the time scale of a large enough number of time windows to provide a convergent unbiased estimate. In real world applications however, the requirement for long stationary signals for the data processing cannot always be fulfilled, and applying bicoherence calculation in such cases can lead to the appearance of false positives: high bicoherence values even in the lack of nonlinear coupling.
This paper introduces a possible way to identify false positives in the estimated bicoherence, caused by nonstationary signals. The approach is based on a Monte Carlo method, where test signals are constructed with the measured Fourier amplitudes but with random phases, and a sufficiently large ensemble of these signals is used to generate the random bicoherence probability density function () for each point of the entire frequencyfrequency plane. Comparing the value of bicoherence estimated from the real signal to the critical bicoherence level calculated from at a given confidence level for each frequencyfrequency point can help decide if the estimated bicoherence is really due to second order phase coupling, or is it just an artefact caused by the changes of Fourier amplitudes in the investigated time window.
The method was tested using numerical simulations of physical test systems. We demonstrated that when the requirements for stationary signals are not met false positives emerge throughout the calculation domain. The method presented in the paper helped identify these false positives at a high confidence level (in the presented examples ) while retaining actual physical nonlinearities which were deliberately introduced in the model system for the purpose of the tests. The filtering process therefore provides an opportunity to make the difference between actual phase coupling and false positives.
Appendix A Symmetries of the bispectrum
From the (6) definition it follows that the bispectrum has several symmetries when the signal is real.
(21)  
(22)  
(23)  
(24) 
These symmetries reduce the range in which the bispectrum (or bicoherence) has to be plotted in. Let us now refer to figure 2. (21) expresses the symmetry vs the mirroring across the line. (22) expresses the symmetry vs the mirroring across the line. (22) & (23) leads to the equivalence of regions P and Q in figure 2, while the equivalence of regions P and R can be shown by applying (21), (22) & (24).
Acknowledgment
This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 20142018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.
References
 [1] B. Cazelles and L. Stone. Detection of imperfect population synchrony in an uncertain world. Journal of Animal Ecology, 72(6):953–968, 2003.
 [2] L. Dyrud, B. Krane, M. Oppenheim, H. L. Pécseli, K. Schlegel, J. Trulsen, and A. W. Wernik. Lowfrequency electrostatic waves in the ionospheric eregion: a comparison of rocket observations and numerical simulations. Annales Geophysicae, 24(11):2959–2979, 2006.
 [3] X. Xie, X. Shang, and G. Chen. Nonlinear interactions among internal tidal waves in the northeastern South China Sea. Chinese Journal of Oceanology and Limnology, 28:996–1001, September 2010.
 [4] P. Manz, P. Lauber, V. E. Nikolaeva, T. Happel, F. Ryter, G. Birkenmeier, A. Bogomolov, G. D. Conway, M. E. Manso, M. Maraschek, D. Prisiazhniuk, and E. Viezzer. Geodesic oscillations and the weakly coherent mode in the imode of asdex upgrade. Nuclear Fusion, 55(8):083004, 2015.
 [5] Julius S. Bendat and Allan G. Piersol. Random data: analysis and measurement procedures, volume 729. John Wiley & Sons, 2011.
 [6] Y. C. Kim and E. J. Powers. Digital bispectral analysis and its applications to nonlinear wave interactions. IEEE Transactions on Plasma Science, 7(2):120–131, June 1979.
 [7] Y. C. Kim, J. M. Beall, E. J. Powers, and R. W. Miksad. Bispectrum and nonlinear wave coupling. The Physics of Fluids, 23(2):258–263, 1980.
 [8] B. Ph. van Milligen, C. Hidalgo, and E. Sánchez. Nonlinear phenomena and intermittency in plasma turbulence. Phys. Rev. Lett., 74:395–398, Jan 1995.
 [9] Vinod Chandran. Timevarying bispectral analysis of visually evoked multichannel eeg. EURASIP Journal on Advances in Signal Processing, 2012(1):140, 2012.
 [10] G. I. Pokol, L. Horvath, N. Lazanyi, G. Papp, G. Por, V. Igochine, et al. Continuous linear timefrequency transforms in the analysis of fusion plasma transients. In Proc. of the 40th EPS Plasma Physics Conf., 2013.
 [11] H. Ombao and S. Van Bellegem. Evolutionary coherence of nonstationary signals. IEEE Transactions on Signal Processing, 56(6):2259–2266, June 2008.
 [12] Christopher Torrence and Gilbert P Compo. A practical guide to wavelet analysis. Bulletin of the American Meteorological society, 79(1):61–78, 1998.
 [13] Harry Nyquist. Certain topics in telegraph transmission theory. Transactions of the American Institute of Electrical Engineers, 47(2):617–644, 1928.
 [14] Peter Welch. The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms. IEEE Transactions on audio and electroacoustics, 15(2):70–73, 1967.
 [15] Stéphane Mallat. A wavelet tour of signal processing. Academic press, 1999.
 [16] Vinod Chandran, Steve Elgar, and Barry Vanhoff. Statistics of tricoherence. IEEE transactions on signal processing, 42(12):3430–3440, 1994.
Comments
There are no comments yet.