APPLICATION OF HIGHER-ORDER SPECTRA FOR SIGNAL PROCESSING IN ELECTRICAL POWER ENGINEERING

Zbigniew Leonowicz, Tadeusz obos, Jarosaw Szymada, Piotr Ruczewski

Technical University of Wroclaw pl. Grunwaldzki 13, PL-50-370 Wroclaw, Poland Tel: (+4871)3203448; Fax: (+4871)229725 E-mail: [email protected] ABSTRACT During recent years higher order statistics (HOS) have found a wide applicability in many diverse fields, e.g.: biomedicine, seismic data processing, harmonic retrieval and adaptive filtering. In power spectrum estimation, the signal under consideration is processed in such a way, that the distribution of power among its frequency is estimated and phase relations between the frequency components are suppressed. Higher order statistics known as cumulants, and their associated Fourier transforms, known as polyspectra, reveal not only amplitude information about a signal, but also phase information. If a non-Gaussian signal is received along with additive Gaussian noise, a transformation to higher order cumulant domain eliminates the noise. These are some methods for estimation of signal components, based on HOS. In the paper we apply the MUSIC method (Multiple Signal Classification), both for the correlation and the 4th order cumulant, to investigate the state of asynchronous running of synchronous machines. When the investigated signal is distorted by a coloured noise the exact results can be achieved by applying cumulants. SPECTRUM ESTIMATION METHODS The linear methods of spectrum estimation (BlackmanTukey, direct methods based on FFT) suffer from a major problem of resolution [4]. Because of some invalid assumptions (zero data or repetitive data outside the duration of observation) made in these methods, the estimated spectrum will be a smeared version of the true spectrum. Closely spaced spectral peaks get merged into a broader peak, leading to false estimates of spectral power. Even the best window becomes ineffective when the dynamic range of the spectrum is high. Methods using a data adaptive window (ML, Capon), invented to alleviate this problem, deliver high resolution at moderate SNR only. If we are able to extrapolate the available data outside the duration of the data (it means to fit the data with a finite rational model), we can completely eliminate the use of the window itself. Once the model parameters are estimated, the spectrum may be derived from the model parameters or from the extrapolated data. The most important time series models are: AR(autoregressive), MA (moving average) and ARMA. The most recent methods of spectrum estimation are based on the linear algebraic concepts of subspaces and so have been called subspace methods [3]. Its resolution is theoretically independent of the SNR. The model of the signal in this case is a sum of random sinusoids in the background of noise of a known covariance function. Pisarenko first observed that the zeros of the z-transform of the eigenvector corresponding to the minimum eigenvalue of the covariance matrix lie on the unit circle, and their angular positions correspond to the frequencies of the sinusoids. In a later development it was shown that the eigenvectors may be divided into two groups, namely, the eigenvectors spanning the signal space and eigenvectors spanning the orthogonal noise space. The eigenvectors spanning the noise space are the ones whose eigenvalues are the smallest and equal to the noise power. One of the most important technique, based on the Pisarenkos approach of separating the data into signal and noise subspaces, is the MUSIC method. MUSIC (MULTIPLE SIGNAL CLASSIFICATION) The MUSIC method [3] involves projection of the signal vector:

s i = 1 e j

L ej

( N 1) i

(1)

onto the entire noise subspace. We consider a random sequence x made up of M independent signals in noise.

x = Ai s i + ;
i =1

Ai = Ai e j

(2) (3)

or: x = S A1 where S = [ s1

A2

L AM

s2 L s M ]

(general model used in all subspace methods). If the noise is white, the correlation matrix is

R x = Ai Ai s i s iT + 02 I
i =1

(4)
(5)

or: R x = SP0 S T + 02 I

 A1 A1 0 L 0 0 L 0 A2 A2 where P0 = M M O 0 0 0 L AM A M N-M smallest eigenvalues of the correlation matrix (matrix dimension N>M+1) correspond to the noise subspace and M largest (all greater than 2 0 ) corresponds to the signal subspace.

The denominator polynomial (15) has M double roots lying on the unit circle, which correspond to the signal frequencies. If the frequencies are to be found by plotting the pseudospectrum, then the spectrum can be written in the form 1 2 ) j N T ( ) P e = w ei (16) i = M +1 Let the components of the vector ei be regarded as a sequence ei[0], ei[1],..., ei[N-1]. Then the Fourier transform of this sequence is E i ( e j ) = w T e i (17) The L values of (17) corresponding to the L desired frequency points can be found by zero-padding the sequence ei[n] to a length L and using an L-point FFT. A total of N such FFTs is needed to evaluate (16). In the most common situation where the number of signals is not known at the outset, it must be estimated from the data. When the SNR is high, it is easy to determine M from the inspection of the eigenvalues. In other cases the AIC (Akaike Information Criterion) or the MDL (Minimum Description Length) are useful. AIC( M ) = 2 K ( N M ) ln ( M ) + 2 M ( 2 N M ) (18)
1 MDL( M ) = K ( N M ) ln( M ) + 2 M ( 2 N M ) lnK (19)

We define the matrices of eigenvectors: E signal = [e 1 e 2 L e M ]

(6)

E noise = e M +1

e M +2

L eN

(7)

and also two matrices of eigenvalues: 1 0 0 2 signal = M M 0 0

L 0 L 0 O M L M 0

(8)

noise

M +1 0 = M 0

M +2
M 0

L 0 02 L 0 =0 O M M L N 0

2 0

M 0

L 0 L 0 O M K 02 (9)

It is possible to write Rx as: T T R x = E signal signal E + E noise noise E signal noise

Enoise can be used to form a projection matrix PX for the noise subspace T Pnoise = E noise E = I Psignal (10) noise The squared magnitude of the projection of w (defined as in (1) ) onto the noise subspace is given by T w T Pnoise w = w T E noise E w (11) noise Since each of the elements of the signal vector is orthogonal to the noise subspace, the quantity (11) goes to zero for values of the frequency where w=si. The MUSIC pseudospectrum is defined as 1 1 ) T P( e j ) = w T Pnoise w = w T E noise E w (12) noise

(20) 1 + M + 2 +K+ N N M M +1 (K-number of rows of the data matrix). The value of M is chosen to minimise (18) or (19).

where: ( M ) =

KN M +1 M + 2

1 NM

CUMULANT-BASED APPROACH The second-, third-, and fourth-order cumulants of a zeromean stationary random process x(t) are C1, x = E x( t ) = 0 (21)

] [

{ }

C 2 , x ( ) = E { x( t ) x(t + )}

(22)

and it exhibits sharp peaks at the signal frequencies where w=si. An alternative variation of the method of pseudospectrum estimation can be developed as follows. Define the eigenfilter ( ) E i ( z ) = ei [0] + ei [1]z 1 +K+ei [ N 1]z N 1 (13) where ei[n] are components of the eigenvector ei. The equation (15) can be written using (17) as
i = M +1

C 3, x 1 , 2 = E x ( t ) x t + 1 x t + 2

C 4 , x 1 , 2 , 3 = E x(t ) x t + 1 x t + 2 x t + 3 C2 , x C2 , x
1 2,x 2 3 2,x 2 2,x 3 1

) { ( )C ( ( )C (
3 2, x

)(

)}

(23)

( )( )( ) C ( )C ( ) = cum(K)
2

)} + )+

(24)

T ei e w= i

i = M +1

E (e ) E (e )
j i

(14)

There are many signal processing applications where signals are either real or complex (sinusoidal) processes. This general class of signals can be described by

The MUSIC pseudospectrum can therefore be expressed as 1 ) N P( e j ) = E i ( z ) E i (1 / z ) (15) i = M +1

z = e j

x(n) = ai (n) sn ( i )
i =1

(25)

where sn(.) are signal waveshapes, i are constants and ai(n)s are random variables.

In case of complex harmonic processes (25) can be expressed as[2]

x( n) = i exp j i n + i
i =1

[(

)]

(26)

where i are independent identically distributed random variables uniformly distributed over [0, 2], ij for ij and is i is are constants. The comparable model for real signals is

short time of the fault and the rotor gains speed to store the excess energy. If the fault persists long enough, the rotor angle will increase continuously and the synchronism will be lost. After switching off the fault the current waveforms contain three main components which frequencies depend on the rotor slip. Assuming a steady-state asynchronous running of a synchronous machine (s=const) the stator winding current can be expressed [5]

y( n) = i cos i n + i
i =1

(27)

is = I f cos s t + f + I b cos (1 2 s) s t + b + + I r cos (1 s) s t + r

(32)

For random variables b=ej, where is uniformly distributed over [0, 2], was shown the following. 1) All third-order cumulants of complex harmonic b are always zero. 2) Of the three different ways to define a fourth-order cumulant of a complex harmonic only one always yields a nonzero value , i.e., cum(b, b, b, b)=0, cum(b*, b, b, b)=0 but cum(b*, b*, b, b)= 1. For (26) in the case of complex harmonics the fourth-order cumulant is given by

where If, f - amplitude and phase angle of the forward field component, Ib, b - amplitude and phase angle of the backward field component, Ir, r - amplitude and phase angle of the rotor field component, s - rotor slip, s angular velocity of the rotating stator field. At the beginning of an asynchronous running, differences between frequencies are small and their identification is rather difficult. EXPERIMENTAL SETUP A small synchronous generator running at 47 Hz was switched on the 50 Hz power supply. The stator winding current samples (fs=1000 Hz)(Fig. 1) were analysed with the help of the MATLAB (v.4.2c.1) function harmest (toolbox hosa v.2.0). Coloured noise, generated by passing white Gaussian noise through a FIR filter [1,1], was added artificially (Fig. 3).
4

C 4 , x 1 , 2 , 3 = k4 exp j k 1 + 2 + 3
k =1

[ (

)]

(28)

For (27) in the case of real harmonics

1 C4 , x = 8 k4 cos k 1 2 3 + cos k 2 3 1 + k =1 p

+ cos k 2 3 1
Additionally, C 2 , x =
1 2

)]
p 2 k

(29) cos k

k =1

( )

(30)

If x(n) is a sum of p real-valued sinusoids, then the diagonal slice of the fourth-order cumulant retains all of the pertinent information about the number of harmonics, their amplitudes and frequencies [2]. Especially set 1 = 2 = 3 = in (29) to see that
k =1

Amplitude

3 C4 , y ( ) = 8 k4 cos( k )

(31)

which is nearly identical with the autocorrelation of the signal. It is known that signal parameters can be estimated from the output correlation in the white noise case, they can also be estimated from the fourth-order cumulant, which is useful in the additive coloured noise case. It means that existing high-resolution methods, such as MUSIC, can be applied by replacing correlation quantities with fourth-order cumulant [1]. ASYNCHRONOUS RUNNING OF SYNCHRONOUS GENERATORS When a fault occurs at the terminals of a synchronous generator the power output of the machine is reduced but the input power from the turbine doesnt change during the

-2

-4 0 200 400 600 800 1000

Sample No

Fig. 1. Current waveforms.

1E+6

Table 1. Estimated signal frequencies.

MUSIC spectrum estimate correlation-based cumulant-based

1E+5

1E+4

f1 [Hz] f2 [Hz]

Fig. 2 corr. cum. 47.09 47.44 49.99 49.88 CONCLUSIONS

Fig. 3 corr. cum. 48.46 46.91 - 49.77

1E+3

Magnitude

1E+2

1E+1

1E+0

1E-1

1E-2 40 42 44 46 48 50 52 54 56 58 60

Frequency [Hz]

Fig. 2. MUSIC spectrum estimates in noiseless case.

In the paper it was shown that high-resolution spectrum estimation methods, such as MUSIC, can be effectively used to identify the loss of synchronism of synchronous machines. Under high-SNR conditions, correlation- and cumulant-based results are similar, although cumulantbased estimates show more variability and require approximately two times longer computation times that those correlation-based. Under poor SNR conditions, when the additive Gaussian noise is coloured, only cumulantbased methods indicate the correct number of signal components. REFERENCES

1E+4

MUSIC spectrum estimate

1E+3

correlation-based cumulant-based

1E+2

1E+1

1E+0

1E-1

1E-2

1E-3 40 42 44 46 48 50 52 54 56 58 60

[1] A. Swami, Cumulant-Based Approach to the Harmonic Retrieval and Related Problems, IEEE Trans. on Signal Processing, vol. 39, pp. 1099-1108, May 1991. [2] J.M. Mendel, Tutorial on Higher-Order Statistics (Spectra) in Signal Processing and System Theory: Theoretical Results and Some Applications, Proceedings of the IEEE, vol. 79, pp. 278-305, March 1991. [3] C.W. Therrien, Discrete Random Signals and Statistical Signal Processing, Englewood Cliffs, New Jersey, Prentice-Hall, 1992, pp. 614-655. [4] P.S. Naidu, Modern Spectrum Analysis of Time Series, Boca Raton, Florida, CRC Press, 1996, pp. 241-336. [5] P. Kostya, T. obos, Z. Wacawek and J. Pytel, Identification of Asynchronous Running of Synchronous Machines Using A Neural Network, XVII Seminar on Fundamentals of Electrotechnics and Circuit Theory, vol. 2, pp. 497-502, May 1994. ACKNOWLEDGMENT

Magnitude

Frequency [Hz]

Fig. 3. MUSIC spectrum estimates in coloured noise case (SNR=0 dB).

This work was supported by the polish State Committee for Scientific Research - KBN (Grant No. 8T10A07210)