Identifying finite dimensional behaviour from broadband spectra

4 April 1994 PHYSICS LETTERS A

ELSEVIER

Physics Letters A 187 (1994) 59-66

Identifying finite dimensional behaviour from broadband spectra J . J . Healey Department of Engineering, University of Cambridge, Trumpington Street, Cambridge CB2 1PZ, UK Received 27 November 1993; revised manuscript received 11 February 1994; accepted for publication 15 February 1994 Communicated by A.P. Fordy

Abstract

A broadband Fourier spectrum can arise either from a large number of linear modes, or from complex nonlinear interactions among a few. A method for differentiating between these two cases is presented. Nonlinear dependences amongst the bases of the reconstructed phase space are tested for using radial basis functions.

1. Introduction

The decomposition of a signal into a set of Fourier modes remains the pre-eminent data analysis technique in many scientific and engineering applications. This approach is the natural one when there are physical grounds for believing that the signal consists of a number of harmonic components. However, the great majority of physical problems are nonlinear to some extent, and so it is useful to know how this nonlinearity manifests itself in the power spectrum. In general, weak nonlinearity results in the coupling of linear modes. This generates additional spikes in the power spectrum at frequencies which are integer combinations of the fundamental components. For example, two linearly independent modes with frequencies jq and 3~, under the influence of nonlinear coupling, will have spectral components at mJq + nJ~ for integer pairs (m, n ) (the values of (m, n) may be constrained by symmetry). Bispectral techniques aim to identify these nonlinear interactions, see for example Ref. [ 1 ]. However, the rapidly growing field of dynamical systems theory has produced many examples of low dimensional nonlinear equations which have non-

periodic, nonquasiperiodic "chaotic" states (see e.g. Refs. [2,3] ). Such states have sensitive dependence on initial conditions i.e. infinitesimal perturbations grow exponentially, and their power spectra contain broadband components. Low dimensional chaotic behaviour has also been observed in physical systems (see e.g. Refs. [4,5]) and sometimes a nonlinear mechanism can be found that explains the observed dynamics (see e.g. Ref. [6] ). Before properties of nonlinear dynamical systems became widely known, there was the danger that an experimentalist, upon unexpectedly observing a broadband spectrum, might conclude that high levels of noise were present. Thus the possibility of relatively simple nonlinear behaviour could have been overlooked. On the other hand, it is now tempting to see finite dimensional chaos when, in fact, only stochastic noise is present. There are other techniques for interpreting broadband spectra in terms of low-dimensional systems (see e.g. Ref. [7] ). The approach there is to model the evolution on the attractor and then to classify the behaviour in terms of invariants of the fitted models like Lyapunov exponents and fractal dimensions. In order to estimate these invariants it is first necessary

0375-9601/94/$07.00 (~ 1994 Elsevier Science B.V. All rights reserved SSDI0375-9601 (94)00135-C

60

J.J. Healey / Physics Letters A 187 (1994) 59-66

to find an embedding dimension for the phase portrait reconstruction. This is done by calculating the correlation dimension, D (see Ref. [8] ) for a number of embedding dimensions (this involves a s s u m i n g finite-dimensional dynamics to start with). For a sufficiently large embedding dimension, D becomes independent of the embedding. Firstly, these techniques tend to be adversely affected by even low noise levels and secondly, Ruelle [9] has shown that for a finite time series the results always give a finite dimension irrespective of the true dimension. Thus there exists the danger that if the (rather strict) data requirements are not met then the conclusion of finite dimensional dynamics comes to rest on a circular argument. Further, a n y noisy signal will give nonzero Lyapunov exponents and noninteger fractal dimension. It would seem somewhat imprudent to assert that a small positive Lyapunov exponent can be used to prove that a physical system is chaotic, that it is low dimensional and that the broadband components in the power spectrum derive from nonlinear interactions among these few modes. In this paper a technique is presented for estimating the number of independent modes present in the signal (the harmonics discussed above, although linearly independent, are nonlinearly dependent on the fundamental frequencies). The method does not attempt to identify chaos. It provides a means for determining which components of a (possibly broadband) spectrum can be attributed to nonlinearity. It also defines a "model dimension" which is the smallest number of independent variables necessary to construct a model which re-creates the dynamics. This dimension applies to the whole attractor and may be thought of as a global generalisation of the local model dimension calculated using the method of Broomhead et al. [10]. The idea behind the technique is discussed in Section 2 and the results of a number of test cases are given in Section 3. An application to "real" data from a wind tunnel experiment is given in Section 4 and the conclusions are discussed in Section 5.

ear dependences among the bases of the reconstructed space. This is easiest to understand when the phase space is constructed using successive derivatives of the time series (x, x, 3¢. . . . ). Consider the self-excited van der Pol oscillator: - e ( 1 - - X 2 ) X + X = 0.

Arbitrarily many derivatives will all be linearly independent of one another (if e # 0), but the second derivative can be expressed as a nonlinear graph of x and ) . The trivial solution is stable for e < 0, it undergoes a Hopf bifurcation at e = 0 and a stable limit cycle exists for e > 0. The behaviour becomes increasingly nonlinear as e is increased but no further bifurcations occur. In this example it can be seen that a cubic curve (surface) of the form x = f ( x , Jc) could be fitted to the data points. Such a curve would pass exactly through the points indicating that the second derivative is dependent on the first two (even though it is linearly independent). Given a single time series, the function f could be estimated by taking successive derivatives. This corresponds to reconstructing a phase portrait and testing for nonlinear dependences among the bases of the embedding space. The basic idea of the approach is to determine whether such a surface can be found for a given time series. In the presence of experimental noise, a more robust phase portrait reconstruction can be obtained based on the "method of delays" (see Ref. [ 11 ] ). Delay vectors are constructed using (Xn, Xn+z, Xn+2z, •. •, X n + ( m _ l ) z ) T where xn is the nth element of the time series, m is the embedding dimension and z is the embedding delay. The delay coordinates can be thought of as a discrete analogue to the derivatives ( x , ~ , 3~. . . . ). To illustrate the approach when using delay vectors, consider a simple Euler integration scheme for ( 1): xn -

-

The approach is based on reconstructing the phase portrait from a time series and then looking for nonlin-

x~ - xn-l h

(2)

and 3~

2. M e t h o d

(1)

:~ - k~-i h

-

x~ - 2x~-t + x~-2 h2 ,

(3)

where h is the time step. Substitution into ( 1 ) gives xn - 2 x n _ l + x n _ 2 - h e ( 1 - x 2) (Xn - x n - 1 ) + h2xn

= 0

(4)

J.J. Healey / Physics Letters A 187 (1994) 59-66

and hence Xn-2 m f (Xn,Xn-l),

(5)

where f is the cubic function defined by (4). Therefore, in the method of delays, the data projected onto one basis would be a nonlinear function of the data projected onto two other bases• Broomhead and King [ 12 ] have shown that better signal to noise ratios are obtained if the delay vectors are projected onto the optimal basis calculated using singular value decomposition (SVD). The set of N delay vectors (with ~ = 1 ) are combined to form the trajectory matrix:

X =

X1

X2

X3

X2

X3

X4

Xn Xn+l

"•"

Xm Xm + I

Xn+2

Xn+m-I

XN XN+I XN+2

XN+m-1

J .

f (Ynl,

Yn2 . . . . .

Ynk).

we had attempted to fit a model of the form xn = f ( X n - l , xn-2) then a more complicated functional form would be required. A class of functions known to span a wide range in function space are the "radial basis functions". They can be thought of as splines generalised to arbitrary dimension (see Refs. [ 14,15 ] ) and have been successfully applied to reconstructed phase portrait data (see e.g. Refs. [ 16,17 ] ). An important difference is that here the model does not describe the evolution on the attractor, but instead, it describes the graph of the attractor in the embedding space• This point is made clearer in the next section. In order to test for the nonlinear dependence of the (k + 1 )th coordinate, a sum of functions of the form

P

(6)

' ~ Unk

(8)

Y ~ tot~(llUnk -- Utkll) 1=1

!

The SVD of X is performed (see Ref. [13] ) giving a set of singular vectors that span X and the singular spectrum that describes the relative energies of the singular vectors. However, the singular spectrum only identifies the number of significant linearly independent modes i.e. the rank of X, and so, like the Fourier spectrum, may contain many modes which are nonlinearly dependent on one another. A similar procedure can be used here as with the straightforward method of delays. Let Ynk be the projection of the nth delay vector onto the kth singular vector, then, analagous to (5), Yn(k+l) =

61

(7)

In addition to improved signal to noise ratios, an advantage of using (7) rather than (4) is that it provides a systematic method for deciding which coordinate should be treated as a function of the other coordinates. It is assumed that singular values dominated by nonlinear interactions have lower energy than the independent modes that give rise to them. For a partitular value of k, (7) gives a relationship for the nonlinearity. The problem now is to choose a function basis for f . If a quadratic model had been chosen for (5) then the cubic interaction would have been missed, or if

is chosen where Unk .m Y n ( k + l ) , Unk = ( Y n l , Yn2, . . . , Ynk) T and I[ II denotes the Euclidean norm. Once the eentres {Utk} have been chosen, the weights tot are varied to minimise

o~2+, = ~N n=l

Unk

p - )-~ to,~(llunk - U,kll)

2.

(9)

1=1

The c e n t r e s {Ulk} were chosen at random from the data set itself and ~(r) = r 3. The least squares minimisation problem posed by (9) is solved using a method based on SVD (see e.g. Ref. [18]). It is not necessary to calculate tot explicitly, only a~+~ is required. Let ak+x be the (k + 1)th singular value of X. If a~:+l ,~ ak+~ then the fitted model is uncorrelated with the data and no nonlinear dependence has been uncovered (although there could be some more complicated nonlinear structure that would require more, or different, basis functions in order to be detected). However, if a~+~ << Ok+~ then most of the energy in the (k + 1 )th singular vector of X can be attributed to nonlinear interactions between the first k singular vectors. In the following examples, the nonlinear dependences for k = 1,2 . . . . . 8 are tested for with p = 50 and N = 1000.

62

J.J. Healey / Physics Letters A 187 (1994,) 59-66

loglo (p.s.d.)

o:] .,o

i

-~

-;

i

-,

;

*

o

i

;

,

;

(b) 1

-0.5 -6

2

f i

-4

l

-2

|

0

l

2

l

4

Fig. 1. Projections of the phase portrait reconstructed from (1) with e = 0.5. (a) First and second singular vectors, (b) second and third singular vectors.

(~) lo' 0"~0"! 1o0

3. Validation lo' i 3.1. Self-excited van der Pol oscillator

Eq. (1) has been integrated for E = 0.5 to give a time series (the x coordinate). The phase portrait has been calculated using the optimal basis method [9] and two projections are shown in Fig. 1. Fig. la shows the limit cycle projected onto the two singular vectors corresponding to the two largest singular values and Fig. lb shows the limit cycle projected onto the second and third singular vectors. If the behaviour were linear then the limit cycle would lie on a plane in the embedding space: Fig. 1a would show an ellipse and Fig. lb would give the side view i.e. a straight line. The singular spectrum would contain only two significant singular values. In the presence of nonlinearity, the ellipse is distotted and there is now a thickness in higher dimensions i.e. many significant singular values. In this case the thickness of the attractor in higher dimensions can be interpreted as nonlinear graphs of the attractor. The particular shape of Fig. lb is clearly a consequence of the cubic nonlinearity in ( 1 ). Fig. 2a shows the Fourier power spectral density. The fundamental frequency has the greatest energy and a number o f harmonics are visible at odd multiples of the fundamental (only odd multiples due to the inversion symmetry of ( 1 ) ). The crosses in Fig. 2b are the singular values (or) and the ten most energetic sin-

,o'

;

;o

(i,)

(~)

Fig. 2. Same data as Fig. 1. (a) Fourier spectrum. (b) Singular spectrum (crosses) and energy remaining after nonlinearity has been accounted for (circles). (c) The ten strongest singular vectors. gular vectors are shown in Fig. 2c. The singular spectrum conceals the fact that the model equation contains only two independent variables. The additional singular vectors are analogous to the harmonics in the Fourier spectrum as can be seen from the increasing frequencies of the corresponding singular vectors. The circles in Fig. 2b show the behaviour of the standard deviations, a', i.e. the energy remaining after nonlinearity has been accounted for. The sharp drop in ~' associated with the third vector (over three orders of magnitude smaller than the corresponding a ) indicates that the radial basis functions have modelled the dynamics in the third dimension as a nonlinear function of the dynamics on the first two dimensions. It is this sharp cut-offin the a ' spectrum that enables the model dimension of the system to be estimated (two in this case).

J.J. Healey ~Physics Letters A 187 (1994) 59--66 3.3. Gaussian noise

10°

tO

63

.5

10

'15

(b)

I0 °

Here the method is tested on data which has structure but still no nonlinear interactions. The time series is generated by performing the inverse Fourier transform of a Gaussian power spectrum with randomised phases. The singular spectrum and singular vectors are shown in Figs. 3c and 3d. There are now six significant singular values representing a finite band-width signal. The singular vectors resemble polynomials of increasing order (a warning against a too literal interpretation of the singular vectors as "shapes" or "coherent structures" appearing in the time series). Again a ~ a ' showing that no nonlinear structure has been detected (even though chance multiples of any given frequency will exist).

0"~ 0 "I 10"

3.4. Chaotic data

10"~

The Rossler model [ 19] is a simple example of a nonlinear chaotic system with few independent variables:

tO"s

tO"

Jc = - y 104

tO

~ - : - - : : ~ r:~ ~

5,

tO

z,

= b + xz15

(d) Fig. 3. (a) and (b) are for white noise, (c) and (d) are for Gaussian noise. In each ease the meaning is the same as for Figs. 2b and 2e. 3.2. White noise

Here the method is tested on white noise which is modelled by a superposition of many independent linear modes. The time series is constructed by performing the inverse Fourier transform on a fiat power spectrum with randomised phases. The singular values (crosses) and standard deviations (circles) are shown in Fig. 3a and the first ten singular vectors in Fig. 3b. Both spectra are essentially fiat and the singlar vectors show no coherent patterns (the slight slope in the spectrum is a consequence of using a finite length data set). For all the dimensions tested a ~ a ' showing that no nonlinear interactions have been detected.

~ = x + ay, cz,

(10)

where a = b = 0.2 and c = 5.0. The method will be tested on the x time series. The reconstructed phase portrait is shown in Fig. 4a and the power spectrum is shown in Fig. 4b. The power spectrum contains broadband components and gives little indication of the simplicity of the model. The singular spectrum, see Fig. 4c, also gives many significant singular values. However, when nonlinearity is accounted for there is almost an order of magnitude drop in the fourth dimension, i.e. nearly 90% of the energy in the fourth singular vector can be attributed to nonlinearity. Hence the three dimensionality of the model has been recovered. 3.5. Chaotic data plus noise

The x time series from (10) is taken and random numbers with a uniform distribution added. The random numbers have a peak-to-peak amplitude that is 10% as large as the peak-to-peak amplitude of the original x data set. They are generated using an algorithm from Ref. [13] and model the effects of measure-

64

J.J. Healey / Physics Letters A 187 (1994) 59-66

(~)

1

loglo (p.s.d.)

loglo (p.s.d.) 2

/

(~) to ~ 0"~ 0";

i

i

1

2

f io o

0"~ 0 "1

Io

~,

1'0

(b)

~0 a

(~)

Fig. 5. Data from (10) with 10% noise. (a) Fourier spectrum. (b) and (c) have the same meaning as Figs. 2b and 2c.

Io

m

Io

5

10

(,]) Fig. 4. Data from (10). (a) Reconstructed phase portrait. (b) Fourier spectrum. (c) and (d) have the same meaning as Figs. 2b and 2c. merit noise in an experiment. The Fourier spectrum is shown in Fig. 5a. Nonlinearity is evident from the harmonics, but it is not clear whether the broadband component can be associated with noise or chaos. The singular spectrum and singular vectors are shown in Figs. 5b and 5c. The first two singular vectors, with almost equal singular values, form a sine and cosine pair and suggest a fundamental periodic component analogous to the strongest peak in the Fourier spectrum. However, the rest of the singular spectrum and singu-

lar vectors are indistinguishable from that of Figs. 3c and 3d which shows only Gaussian noise. However, the graph of the standard deviations, tr', shows a sharp cut-off with only three significant independent directions. Many physcial experiments can achieve cleaner signals, but this test demonstrates the robustness of the method.

4. Wind tunnel data A series of wind tunnel experiments on the evolution of disturbances in a growing flat plate boundary layer have been performed in the low-turbulence research facility in the Cambridge University Engineering Department. The experiments were motivated by the work of Gaster [20]. A disturbance is introduced into the boundary layer via a small loudspeaker em-

J.J. Healey /Physics Letters A 187 (1994) 59-66

(~)

O ioglo (p.s.d.)

65

spectrum is shown in Fig. 6c and contains five significant singular values. However, when the nonlinear component is subtracted away there are only three remaining independent directions i.e. the fourth and fifth singular values correspond to the harmonic in the Fourier spectrum. Therefore, in order to construct a nonlinear model for the dynamics it is only necessary to use three independent variables. However, this is n o t to say that the motion is chaotic, only that it is three dimensional. In fact to understand the dynamics it is necessary to study the full spatio-temporal behaviour and not a single time series [ 18 ].

5. Conclusions

f(Hz) O'~ iT t IO'

io ~

io ~

IO"

IO"

5

to

(~)

15

(d)

Fig. 6. Wind tunnel data. (a) Reconstructed phase portrait. (b) Fourier spectrum. (c) and (d) have the same meaning as Figs. 2b and 2c. bedded near the leading edge of the plate. A hot-wire probe measures the velocity fluctuations of the disturbance at a point downstream. The disturbance studied was a 200 Hz sine wave. The phase portrait for one of the hot-wire signals is shown in Fig. 6a. The basic circling motion corresponds to the driving frequency and the thickness of the band is associated with irregular modulation of the fundamental. The power spectrum is shown in Fig. 6b. In addition to the fundamental frequency, there is a harmonic and a broadband component. The singular

In this paper a new method for estimating the number of independent modes in a time series has been presented. The method is based on the principle that in a reconstructed phase portrait, the signal projected onto the (k + 1 )th coordinate, although linearly independent from the other k coordinates, may be nonlinearly dependent on them. The natural way of constructing a set of linearly independent coordinates is via singular value decomposition, and in this basis nonlinear dependences are tested for explicitly using radial basis functions. The approach has been demonstrated on a number of test data sets and has also been applied to a data set taken from a wind-tunnel experiment. In all cases it was possible to estimate the number of independent variables needed to model the dynamics. The method appears to be robust in the presence of noise. It should be emphasised that the identification of finite-dimensional dynamics does not depend on trying to prove that the behaviour is chaotic. In spatially developing flows the notion of temporal chaos may not even be appropriate.

Acknowledgement The author was supported by the SERC (UK).

66

ZJ. Healey / Physics Letters A 187 (I 994) 59-66

References [ 1] Y.C. Kim and E.J. Powers, IEEE Trans. Plasma Sci. 7 (1979) 120. [2 ] P. Berg6, Y. Pomeau and C. Vidal, Order within chaos (Hermann/Wiley, Paris/New York, 1984). [3] J. Guckenheimer and P. Holmes, Nonlinear oscillators, dynamical systems and bifurcations of vector fields, 2nd Ed. (Springer, Berlin, 1986). [4] A. Libchaber and J. Maurer, in: Nonlinear phenomena at phase transitions and instabilities, ed. T. Triste (Plenum, New York, 1982) p. 259. [5] T. Mullin and T.J. Price, Nature 340 (1989) 294. [6] J.J. Healey, D.S. Broomhead, K.A. Cliffe, R. Jones and T. Mullin, Physica D 48 (1991) 322. [7] H.D.I. Abarbanel, R. Brown and J.B. Kadtke, Phys. Lett. A 138 (1989) 401. [8 ] P. Grassberger and I. Procaccia, Physica D 9 (1983) 189. [9] D. Ruelle, Proc. R. Soc. A 427 (1989) 241. [10] D.S. Broomhead, R. Jones and G.P. King, J. Phys. A 20 (1987) L563.

[ 11 ] F. Takens, in: Lecture notes in mathematics, Vol. 898, eds. D.A. Rand and L.-S. Young (Springer, Berlin, 1981) p. 366. [ 12 ] D.S. Broomhead and G.P. King, Physica D 20 (1986) 217. [13] W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling, Numerical recipes (Cambridge Univ. Press, Cambridge, 1986). [14] M.J.D. Powell, Proc. IMA Conf. on Algorithms for the approximation of functions and data, RMCS Shrivenham (1985). [ 15 ] M.J.D. Powell, The theory of radial basis function approximation, in: Internal Report, DAMTP, Cambridge, UK (1990). [16] M. Casdagli, Physica D 35 (1989) 335. [17] D.S. Broomhead and D. Lowe, Complex Syst. 2 (1988) 321. [18] J.J. Healey, J. Fluid Mech. 255 (1993) 667. [19] O.E. Rossler, Phys. Lett. A 57 (1976) 397. [20] M. Gaster, Proc. R. Soc. A 430 (1990) 3.