Analysis of noisy signals

PHYSICAi

Physica D 58 (1992) 284-287 North-Holland

Analysis of noisy signals A . R a b i n o v i t c h a a n d R. T h i e b e r g e r ab "Physics Department, Ben Gurion University, Beer Sheva, Israel bDepartment of Physics, N.R.C.N., P.O. Box 9001, Beer Sheva, Israel

Received 7 October 1991 Revised manuscript received 10 February 1992 Accepted 2 March 1992

A method is developed to obtain all possible pure signals which could have been the origins of a given noisy signal. The method also provides an explicit way to calculate power spectra of these signals.

1. Introduction

Autoregressive moving average ( A R M A ) models have been extensively used for predicting the time evolution of signals [1]. They constitute both the first approximation to use when the signal is not known to be chaotic, and a general m e t h o d of obtaining the pattern of the " n o n deterministic" part of a signal by Wold's decomposition [2]. Recently the A R M A method has been applied, in this sense, to the analysis of chaotic data [3, 4]. Such analysis should however be carried out with care. For chaotic and other essentially non-linear signals the non-deterministic part is only guaranteed to be " u n c o r r e l a t e d " , while the A R M A approach implies "independe n c e " [6]. Thus, e.g., an A R M A analysis of the chaotic logistic m a p x---~ b x ( 1 - x) for b < 4 yields some " p a t t e r n e d " part [4], while for b = 4 a completely " u n c o r r e l a t e d " generator (which can be interpreted as white noise) ensues. A further decomposition of the non-deterministic part can be based on general non-linear methods [5, 7] or on specific more direct methods [8] when chaotic attractors are known or assumed to exist. When a B o x - J e n k i n s [1] or similar analysis is

carried out, a specific A R M A filter is obtained as being responsible for "creating" the signal from a r a n d o m shock generator. Usually there appears an added white noise as a result of measuring instruments or as an intrinsic property of the process measured itself. The amount of noise included in the signal is generally unknown. This noise brings about alterations in the A R M A p a r a m e t e r s obtained [9] and hence also in the predicted time behaviour. Since the exact magnitude of such noise is not usually available, it is evident that a complete unravelling of the " p u r e " signal from the noisy one is impossible. It would seem however of advantage to havc even a partial possibility of "purifying" the signal, such as a knowledge of all possible pure signals which could have played the role of our signal's origin under the addition of different amounts of white noise. For the A R M A (2, 2) case a complete analysis was carried out previously [10]. By a non-linear transformation the loci of the parameters, for all possible pure signals, were shown to lie on a straight line. The Fourier transform and predicted time evolution could be calculated. H e r e we will discuss the general A R M A (m, l) case.

0167-2789/92/$05.00 © 1992- Elsevier Science Publishers B.V. All rights reserved

285

A. Rabinovitch, R. Thieberger / Analysis of noisy signals A

2. Method

W e consider a stationary and invertible signal of an A R M A (m, l) type, l <- m, given by

Z'

(1)

1 - q l B - q2 B2 . . . . . qtB l ~ = 1 -- - p ~ - p - ~ - : -- -Pm--ffm '

w h e r e t is the discrete time, 2, is the " c e n t e r e d " signal (i.e. Z, = Z , - ( Z , ) ) , Pl . . . . . P ~ , q~ . . . . . qt are the A R and M A p a r a m e t e r s , respectively, ~, is the r a n d o m shock g e n e r a t o r with zero m e a n and variance 02 and B is the shift o p e r a t o r . W e assume that this signal is c o m p o s e d of an original " p u r e " signal )(, with an a d d e d white noise 6, with variance 0-2:

"

- i

2, +/),.

=

i=0

- j

~ ,

L/]

--

Pm

-]

E U ~ 4 __

PiPi+i

i=0

qm 0

(2)

T h e p u r e signal )(,, has to be of the A R M A ( m , m ) type [11]. Since the addition of white noise does not c h a n g e the A R p a r a m e t e r s [1] the f o r m of X, is

m

~" qi,oqi+i,o /~/j __

t~.z

2t

A

for j = 0, 1 . . . . . m. H e r e q0 = P 0 = %,0 are defined as - 1 . These m+l equations for the m + 2 unk n o w n s ( qi,o, i = 1, 2 . . . . , m, 0-] and 0-2b) define the locus of all possible original processes. T h e s e e q u a t i o n s are non-linear and rather inconvenient to solve directly #1. Since they d e p e n d upon the a u t o c o r r e l a t i o n s they have a p p e a r e d r e p e a t e d l y in time series analysis [12]. We present here a t r a n s f o r m a t i o n which m a r k e d l y simplifies and facilitates the analysis. Define:

qiqi+j

j=0,1

i=0

.....

m-l,

qm and 2

j(,

=

1 - ql,o B

-

q2,0 B2

. . . . .

1 - p l B - pe B2 . . . . .

qm,o B ' ' pm Bm

,%

"

(3)

H e r e qi,0 are the p u r e signal M A p a r a m e t e r s and fi, is the g e n e r a t o r of the p u r e signal, having a v a r i a n c e 0-]. T h e values qi, P~, 0-2 are o b t a i n e d directly f r o m the signal, on using any A R M A p r o g r a m , and we wish to e s t i m a t e qi,o, 0-a and o-b which are u n k n o w n . L e t us d e n o t e ?, : (1 - p~B - p 2 B ~ . . . . .

pmBm)2~,

(4)

then, calculating E{Y,Y,+~} and using eqs. ( 1 ) (3) and the white noise c h a r a c t e r of fi,,/~, and ~,, we o b t a i n m-j 2

0-a E

qi,oqi+j.O

i=0

m-j

=0";2 E i=O

m j qiqi+y -- 0-2 E PiPi+y ' i

0

(5)

6_qm 0-c 2 • Pm O'b

(6)

By inserting eq. (6) into eq. (5) one obtains the locus as a straight line in m dimensions: u'=UJz+b(u~-U/z),

b=6/(6-1)

(7)

We thus see both the reduction of dimension (m vs. m + 1 in eq. (5)) and the linear character of the locus in u-space. T h e m a r k e d a d v a n t a g e of the t r a n s f o r m a t i o n is a p p a r e n t if the evaluation of the p o w e r s p e c t r u m ~2 is o n e ' s main interest. This point will n o w be discussed.

"~ Using MATHEMATICA, we had difficulties in solving the set of equations (5). For example, in the case m =4 discussed in fig. 1, no solution could be obtained on a 386 pc, because of the lack of sufficient memory after running for over 15 minutes• #2 Again, since the spectrum is based on the linear analysis, it can be considered as (an important) first approximation to the non-linear approach. In the latter, "polyspectra'" are used for handling higher comulants. See e.g. ref. [7].

286

A. Rabinovitch, R. Thieberger

Analysis of noisy signals

30 25 20 --

15

o

8.

10 I

0.05

011

0.i5

---0("i

0,25

0.3 "~....

0.35

014

0),5

0.5

f Fig. 1. Comparison of spectral densities, b = -0.005, -50000.

3. The power spectrum No direct subtraction of white noise from the p o w e r spectrum is possible. H o w e v e r , the power spectrum of the pure signal can be obtained directly f r o m the points in the u-space. Thus, this p o w e r spectrum is given by [1]: g(w) = constant ×

1 - ql,o B . . . . . ]-- p-~-

-

qm,oBml 2

--pm B

B = exp(-io~).

I '

(8)

Using eq. (6) we can transform eq. (8) to g(oJ) = constant × [u ° + 2u I cos(w) + - . .

+ 2u"

1 cos((m - 1)w) - 2 cos(mw)]

× [u ° + 2u z cos(w) + . . 71_ 2 U z

-1

cos((m - 1)w) - 2 cos(rnw)]-1 .

(9)

W h e r e the constant is determined by the normalization: 0.5

f 0

g(f) df=l,

f=w/2~r.

In fig. 1 we show the results for a demonstratior case. We chose A R M A (4, 4) with the followin~ parameters: ql = q3 = 0, q2 = - 1 . 7 7 , q4 = - 0 . 8 1 P~ = P3 = 0, P2 = - 1 . 5 7 , P4 = - 0 . 6 4 . We choose b = - 0 . 0 0 5 to represent a noisy signal and b = - 5 0 0 0 0 to represent a fairly cleaned up signal. We see quite clearly the "purifying" effect. It should be noted that the loss of invertibilit) or of stationarity can be detected from the powel spectrum directly, thus there is no need to transform back to the qi.0 for this purpose. This k, achieved by detecting a zero (node) at a point i~ the spectrum which indicates an intersection ol the zero of the q0 polynomial with the unit circle. Similarly, an appearance of negative parts in the " p o w e r s p e c t r u m " indicates a region of u-space which is not allowed (complex q/,o)"

(10)

References [1] G.E.P. Box, G.M. Jenkins, Time Series Analysb (Holden-Day, San Francisco, 1976). [2] H. Wold, A Study in the Analysis of Stationary Tim~ Series, 2nd ed. (Almquist and Wiksell, Stokholm, 1953)

A. Rabinovitch, R. Thieberger / Analysis of noisy signals [3] A. Rabinovitch and R. Thieberger, Physica D 28 (1987) 409. [4] A. Rabinovitch and R. Thieberger, J. Theor. Biol. 131 (1988) 509. [5] H. Tong, Non-Linear Time Series: a Dynamical System Approach (Clarendon, Oxford, 1990). [6] M.B. Priestley, Spectral Analysis and Time Series (Academic, New York, 1981). [7] M.B. Priestley, Non-linear and Non-stationary Time Series Analysis (Academic, New York, 1991). [8] M. Casdagli, Physica D 35 (1989) 335; G. Sugihara, B. Grenfell and R.M. May, Phil. Trans. R. Soc. Lond. B 330 (1990) 235;

[9] [10] [11] [12]

287

J.D. Scargle, Astrophys. J. 359 (1990) 469, and references therein. S.M. Kay, Modern Spectral Estimation (Prentice-Hall, Englewood Cliffs, NY, 1988). A. Rabinovitch and R. Thieberger, J. Time Series Anal. 13 (1992) 267. C.W.J. Granger and M.J. Morris, J. R. Star. Soc. A 139 1I (1976) 246. D.H. Anderson and E.C. Gartland JR., SIAM J. Sci. Stat. Comput. 6 (1985) 376.