Скачать 30.97 Kb.
|Instantaneous-phase wavelet bispectral analysis|
Janez Jamšek1 and Aneta Stefanovska2
1 Department of Physics and Technology,Uni. of Ljubljana, Faculty of Education, Kardeljeva pl. 16, 1000 Ljubljana
2 Department of Physics, University of Lancaster, Lancaster, LA1 4YB, UK
Abstract - Bispectral analysis, recently introduced as a technique for revealing time-phase relationships, is extended to make use of wavelets rather than Fourier analysis. It is thus able to encompass instantaneous phase-time dependence for the case of two or more coupled nonlinear oscillators. The method is demonstrated and evaluated by use of test signals from a pair of coupled Poincaré oscillators. It promises to be useful in a wide range of scientific contexts for studies of interacting oscillators whose basic frequencies are significantly time-variable.
The phase relationships between a pair of interacting oscillators can be obtained from bivariate data (i.e.\ where the coordinate of each oscillator can be measured separately) by use of the methods recently developed for analysis of synchronization, or generalized synchronization, between chaotic and/or noisy systems. Not only can the interactions be detected , but their strength and direction can also be determined [2-5]. The next logical step in studying the interactions among coupled oscillators was to determine the nature of the couplings  as the methods developed for synchronization analysis are not capable of answering this question.
Systems are usually taken to be stationary. For real systems, however, the mutual interaction among subsystems often results in time-variability of their characteristic frequencies. Frequency and phase couplings can occur temporally, and the strength of coupling between pairs of individual oscillators may change with time. In studying such systems, bispectral analysis for stationary signals, based on time averages, is no longer sufficient. Rather, the time evolution of the bispectral estimates is needed. Schack et al.  have recently introduced a time-varying spectral method. They assume, however, that the interacting oscillators are harmonic. Millingen et al.  introduced the wavelet bicoherence and were the first to demonstrate the use of bispectra for studying interactions among nonlinear oscillators. They used the method to detect periodic and chaotic interactions between two coupled van der Pol oscillators, but without concentrating on time-phase relationships in particular.  succeeded in extending bispectral analysis to encompass time dependence, and demonstrated the potential of the extended technique to determine the type of couplings among interacting nonlinear oscillators. The method has the advantage that it allows an arbitrary number of interacting oscillatory processes to be studied. It is applicable both to univariate data, and to multivariate data. It yields results that are applicable quite generally to any system of coupled nonlinear oscillators.
Our principal motivation has been to develop a technique for studying the human cardiovascular system , including the interactions among its subsystems, and the nature of these interactions. The technique was applied to univariate cardiovascular blood flow signals . In the present paper, we extend the technique to wavelets. Here, however, we are concerned with basic principles, and in demonstrating/testing the technique on a well-characterized simple model
Bispectral analysis belongs to a group of techniques based on high order statistics that may be used to analyze non-Gaussian signals, to obtain phase information, to suppress Gaussian noise of unknown spectral form, and to detect and characterize signal nonlinearities .
The bispectrum involves third-order statistics. Spectral estimation is based on the conventional Fourier type direct approach through computation of the 3rd order moments, which in the case of 3rd order statistics are equivalent to 3rd order cumulants .
The classical bispectrum estimate is obtained as an average of estimated 3rd order moments (cumulants)
where the 3rd order moment estimate is performed by a triple product of discrete Fourier transforms at discrete frequencies k, l and k+l
with i = 1,…K segments into which the signal is divided. The bispectrum B(k, l) is a complex quantity, defined by magnitude A = |B(k, l)| and phase =B(k, l). Consequently, for each (k, l), its value can be represented as a point in a complex space, Re[B(k ,l)] versus Im[B(k ,l)], thus defining a vector. Its magnitude (length) is known as the biamplitude. The phase, which for the bispectrum is called the biphase, is determined by the angle between the vector and the positive real axis.
Generalization of bispectrum based on Fourier transform to wavelets can be seen as a generalization of the Fourier analysis  by adding time resolution - in a more fundamental way than is permitted by the Short-Time Fourier Transform (STFT) .
Morlet first introduced wavelet analysis . Within this transform, the window length is adjusted to the frequency currently being analysed. It is a scale independent method. Window function is called a mother wavelet or basic wavelet. It can be any function that satisfies the wavelet admissibility condition . This function introduces a scale s (its width) into the analyses. The mother wavelet is also translated along the signal to achieve time localization. Morlet mother wavelet was chosen to be the most suitable one . The frequency resolution changes with frequency; at low frequencies (large scales), the resolution is better than at the high frequencies (small scales).
The definitions are completely analogous to the definitions used in bispectral analysis based on Fourier transform [5, 8]. The Wavelet Bispectrum (WB) is given by
where 1/s=1/s1+1/s2. The WB measures the amount of phase coupling in the interval T that occurs between wavelet components of scale lengths s1 and s2 and s of signal g(t), in such a way that the frequency sum-rule is satisfied. It is a complex quantity, defined by magnitude A and phase .
The instantaneous biphase and biamplitude are calculated: from (3) and (4)
If two scale components s1 and s2 are scale and phase coupled, s = s1 + s2, it holds that the biphase is 0 2 radians. For our purposes, the phase coupling is less strict because dependent scale components can be phase-delayed. We consider phase coupling to exist if the biphase is constant (but not necessarily = 0 radians) for at least several periods of the highest scale component. Simultaneously, we observe the instantaneous biamplitude from which it is possible to infer the relative strength of the interaction.
For ease of interpretation, the WB is plotted in the (f1, f2)-plane, rather than in the (s1, s2)-plane. It has the same symmetries in frequency domain as in the case of Fourier based Bispectrum (FB) . The non-redundant region is the principal domain of the wavelet bispectrum. Similarly, the principal domain can be divided into two triangular regions in which the wavelet bispectrum has different properties: the inner triangle (IT), and the outer one. The IT is of our interest defined in .
To illustrate the essence of the method, and to test it, we use a generic model of two interacting systems whose basic unit is the Poincaré oscillator
The activity of each subsystem is described by the two state variables, xi and yi, where i = 1, 2 denotes the subsystem and i, ai and i are constants and 2 is the coupling amplitude. The parameters of the model are set to 1 = 1, a1 = 0,5 and 2, a 2= 1. Here (t) is zero-mean white Gaussian noise (t) = 0, (t), (0) =D(t) and D=0,08 is the noise intensity. In this way we obtain a test signal x1A(t), Fig. 1 (a).
The test signal x1A(t) is the variable x1 of the first oscillator, recorded as a continuous time series. Prior analyzing the signal was first normalized between 0 and 1 and than removed its mean value. For the first 400 s, there was quadratic coupling with amplitude 2 = 0,2. After a further 400 s, there was no forcing, i.e., the amplitude was set to zero. In the last 400 s the coupling
amplitude was increased again to 0,2. In this case the characteristic frequency of the second oscillator is linearly increasing from f2 = 0,24 to f2 = 0,34. The first 15 s for each coupling strength are shown in Fig. 1 (a).
A quadratic nonlinear interaction generates higher harmonic components in addition to the characteristic frequencies . In order to demonstrate the changes in spectral content and behaviour caused by the coupling power spectrum is showed Fig. 1 (b). The peak labelled as f1 = 1,1 Hz represents the first oscillator and the f2 = 0,24 Hz represents the second oscillator. These frequencies are deliberately chosen to have an irrational ratio. Clearly, the test signal x1A has a richer harmonic structure when nonlinear coupling is present. In addition to the characteristic frequency of the first oscillator it contains components with frequencies 2f1, 2f2, f1+f2 and f1-f2. As well as having a particular harmonic structure, the components of the signal x1A also have related phases, 21, 22, 1+2 and 1-2. For bispectral analysis the whole signal is analyzed as a single entity, but the transients caused by the changes in coupling strength are removed prior to processing. First the wavelet bispectrum is estimated Fig. 2 (a) and (b). Its close inspection resolves that all the necessary peaks arise as a result of a nonlinear interaction between the two oscillators f1 and f2 are present. The quadratical coupling was already discussed in detail in  and is not a subject of this paper.
Bifrequencies where peaks provide evidence of possible frequency interactions are then further studied by calculation of the biphase/biamplitude as functions of time. Bifrequency (f1, f2) is of our primary interest. Time evolution for its biphase and biamplitude are shown in Fig. 2 (c) and (d). The results for non-zero coupling are quite different from those where coupling is absent (second 400 s). The biphase is constant in the presence of quadratic coupling (first 400 s) and the biamplitude is above zero. In the third 400 s one of the bifrequency is varying, f2. In the bispectrum we concentrate only on peaks, where possible phase and/or frequency couplings occur, by estimating time-biphase evolution for the peak bifrequency, i.e., where the average biamplitude has the highest value. This is why while the characteristic frequency of the second oscillator f2 is close approximately to 0,24 Hz we obtain constant biphase and high biamplitude at approximately 800 s. Further when the f2 continues to grow linearly to 0,34 Hz the biamplitude drops to 0 and the biphase has phase slips false indicating that there is no coupling Fig. 2 (c) and (d).
Instantaneous phase of the first f1(t) and the second f2(t) oscillator activity was obtained by characteristic or marker events (zero-crossing as the marker event, as we proposed for the cardiovascular ECG and respiration signals [10, 17]). In synchrogram (not shown), for m = 1 a short-lasting horizontal structure that resolves 1:4 phase locking during approximately 940 s can be seen, otherwise no synchronization is evident in the synchrogram. Biphase and biamplitude for instantaneous bifrequency are shown in Fig. 3 (a) and (b). During the last 800 s the biamplitude is above zero and the biphase is constant during the whole time of observation. In such a way we obtain correct results.
Although there is no coupling, synchronization can onset in synchrogram, due to constant frequency ratio. In case of weak synchronization. bispectrum has advantage over synchrogram. It still clearly detects the coupling at bifrequency of our primary interest (f1, f2), there can be more phase slips.
Using bispectral analysis to obtain the time-dependent biphase/biamplitude estimate based on STFT means using a window of constant length. The optimal window length depends, however, on the frequency being studied. The effective length of the window used for each frequency can be varied by applying the wavelet transform. With broader frequency content where natural frequencies lie, the STFT in not sufficient for good time and phase/frequency localization and the wavelet transform needs to be applied.
Using wavelet bispectrum intermittent phase couplings can be detected, whereas Fourier bispectrum (FB) averages out most of the time relevant information. There are many other advantages of wavelet bispectrum over the FB. It allows an arbitrary frequency step to be chosen, achieving the optimum time and frequency resolution, there is a simple relation between scale and frequency; it has smaller statistical error; and is computationally less demanding. The only drawback is that it has to be normalized to obtain signal energy, and it is not orthogonal. However normalization can be preformed, whereas we are not concerned with the inverse wavelet transform.
Mutual interaction among cardiovascular subsystems often results in time-variability of their characteristic frequencies. In studying such systems, we wish to be able to trace the time-variations of characteristics frequencies. This can be achieved by implementing instantaneous bifrequency in the WB and thus obtaining instantaneous biphase and bifrequency. In such a way we obtain better information about the cardiovascular subsystems coupling. In general the proposed wavelet bispectral analysis provides a promising tool for studying the nature of coupling between two or more nonlinear oscillators whose basic frequencies considerately change in time.
 A.S. Pikovsky, M.G. Rosenblum, and J. Kurths, Synchronization; A Universal Concept in Nonlinear Sciences, Cambridge, Cambridge University Press, 2001.
 T. Schreiber, Phys. Rev. Lett. vol. 85, p. 461, 2000; M.G.Rosenblum, and A.S.Pikovsky, "Detecting direction of coupling in interacting oscillators", Phys. Rev. E, vol. 64 , p. 045202, 2001.
 M.G. Rosenblum, L. Cimponeriu, A. Bezerianos, A. Patzak, and A. Mrowka, "Identification of coupling direction: Application to cardiovascular interaction", Phys. Rev. E, vol. 65, p. 041909, 2002.
 M. Paluš, V. Komárek, Z. Hrnčiř and K. Svebrová, "Synchronization as adjustment of information rates: Detection from bivariate time series", Phys. Rev. E, vol. 63, p. 046211, 2001.
 N.B. Janson, A.G. Balanov, V.S. Anishchenko and P.V.E. McClintock, " Detecting synchronization of self-sustained oscillators using wavelet analysis of univariate data for variable external drive frequency ", Phys. Rev. E, vol. 65, p. 036211, 2002.
 J. Jamšek, A. Stefanovska, P.V.E. McClintock and I.A. Khovanov, " Time-phase bispectral analysis ", Phys. Rev. E vol. 68, p. 01620, 2003.
 B. Schack et al., " Time-variant non linear phase coupling analysis of EEG burst patterns in sedated patiens during electroencephalic burst suppression period", Clinical neurophysiology, vol. 112, p. 1388, 2001.
 B.Ph. van Milligen, C. Hidalgo, and E. Sánchez, " Nonlinear phenomena and intermittency in plasma turbulence", Phys. Rev. Lett., vol. 74, p. 395, 1995.
 A. Stefanovska and M. Bračič, "Physics of the human cardiovascular systems", Contemporary Phys., vol. 40, p. 31, 1999.
 J. Jamšek, A. Stefanovska and P.V.E. McClintock, "Nonlinear cardio-respiratory interactions revealed by time-phase bispectral analysis", Physics in Medicine and Biology, vol. 49, p. 4407, 2004.
 C.L. Nikias and A.P. Petropulu, Higher-Order Spectra Anlysis: A Nonlinear Signal Proc. Framework, Englewood Cliffs, Pren.-Hall, 1993.
 J.G. Proakis and D.G. Manolakis, Digital Signal Processing, New Jersey, Prentice-Hall, 1996.
 G. Kaiser, A Friendly Guide to Wavelets, Boston, Birkhäuser, 1994.
 A. Grossmann and J. Morlet, "Decomposition of Hardy functions into square integrable wavelets of constant shape", SIAM J. Math. Anal., vol. 15, p. 723, 1984.
 M. Bračič and A. Stefanovska, "Wavelet-based analysis of human blood-flow dynamics", Bull. Math. Biol., vol. 60, p. 919, 1998.
 L.A. Pflug, G.E. Ioup, and J.W. Ioup, "Sampling requirements for nth-order correlations", J. Acoust. Soc. Am., vol. 95, p. 2762, 1994.
 A. Stefanovska and M. Bračič, "Reconstructing cardiovascular dynamics", Control Eng. Pract., vol. 7, p. 161, 1999.