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 and their associated Fourier transforms 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 both for the correlation and the 4th order cumulant, to investigate the state of asynchronous running of synchronous machines and the fault operation of inverter-fed induction motors. When the investigated signal is distorted by a coloured noise, more exact results can be achieved by applying cumulants

Leonowicz, Z

Lobos, T

Ruczewski, P

Szymanda, J

English

Repozytorium Eny Politechnika Wrocławska

 1  
Application of higher-order spectra for signal processing in electrical 
power engineering 
 
T. LOBOS, Z. LEONOWICZ, J. SZYMANDA, P. RUCZEWSKI
∗
 
 
During recent years higher order statistics (HOS) have found a wide applicability in many diverse fields, e.g.: 
biomedicine, 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 and their associated Fourier transforms 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 both for the correlation and the 4th order cumulant, to investigate the state of asynchronous running of 
synchronous machines and the fault operation of inverter-fed induction motors. 
When the investigated signal is distorted by a coloured noise, more exact results can be achieved by applying 
cumulants. 
 
Keywords: higher order statistics, spectral analysis, synchronous generators, induction motors, semiconductor 
inverters, fault detection. 
 
1.  Spectrum estimation methods 
 
The linear methods of spectrum estimation (Blackman-Tukey, direct methods based on FFT) suffer from a 
major problem of resolution (Naidu, 1996). 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” (Therrien, 1992). 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 techniques, based on the Pisarenko’s approach of separating the data 
into signal and noise subspaces, is the MUSIC method. 
 
2.  MUSIC (MUltiple SIgnal Classification) 
 
The MUSIC method (Therrien, 1992) involves projection of the signal vector: 
( )[ ]s i = −1 1e ej j N
T
i iω ω⋯             (1) 
onto the entire noise subspace.  
 
We consider a random sequence x made up of M independent signals in noise.  
x s= + =
=
∑ A A A ei
i
M
i i
j i
i
η φ
1
 ;                 (2) 
                                                           
∗
 Wroclaw University of Technology, pl. Grunwaldzki 13, 50-370 Wroclaw, Poland 
 2  
or: [ ]x = +S A A AM
T
1 2 ⋯ η             (3) 
where  [ ]S s s s= 1 2 M⋯                  
 (general model used in all “subspace methods”). If the noise (η ) is white, the correlation matrix is 
{ }R s s Ix
i
M
i iA A= +
=
∗∑E
1
0
2
i i
T σ             (4) 
or: R SP S Ix
T= +∗0 0
2σ                (5) 
where 
{ }
{ }
{ }
P0
1 1
2 2
0 0
0 0
0
0 0
=














∗
∗
∗
E
E
E
A A
A A
A AM 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 σ0
2
) corresponds to the signal subspace.  
We define the matrices of eigenvectors: 
[ ]E e e esignal M= 1 2 ⋯             (6) 
[ ]E e e enoise M M N= + +1 2 ⋯            (7) 
and also two matrices of eigenvalues: 
Λ signal
M
=












λ
λ
λ
1
2
0 0
0 0
0 0
⋯
⋯
⋮ ⋮ ⋱ ⋮
⋯
           (8) 
Λ noise
M
M
N
=
























+
+
λ
λ
λ
σ
σ
σ
1
2
0 0
0 0
0 0
0 0
0 0
0 0
⋯
⋯
⋮ ⋮ ⋱ ⋮
⋯
⋯
⋯
⋮ ⋮ ⋱ ⋮
…
=  
0
2
0
2
0
2
    
It is possible to write Rx as: 
R E E E Ex signal signal signal
T
noise noise noise
T= +∗ ∗Λ Λ        (9) 
Enoise can be used to form a projection matrix PX for the noise subspace  
P E E I Pnoise noise noise
T
signal= = −
∗
           (10) 
The squared magnitude of the projection of w (defined as in (1) ) onto the noise subspace is given by 
w P w w E E w∗ ∗ ∗=T noise
T
noise noise
T
          (11) 
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 
( ) [ ] [ ]⌢P e j T noise T noise noiseTω = =∗
− ∗ ∗ −
w P w w E E w
1 1
     (12) 
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 z e e z e N zi i i i N= + + + −− − −0 1 11 1…       (13) 
where ei[n] are components of the eigenvector ei. The equation (11) can be written using (13) as 
( ) ( )w e e w∗ ∗
= + = +
∗=∑ ∑T i i T
i M
N
i
j
i M
N
i
jE e E e
1 1
   ω ω          (14) 
The MUSIC pseudospectrum can therefore be expressed as 
 3  
( ) ( ) ( )⌢P e E z E zj i i
i M
N
z e j
ω
ω
=






∗ ∗
= +
−
=
∑ 1
1
1
/         (15) 
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 
( )⌢P e j T i
i M
N
ω =






∗
= +
−
∑ w e
2
1
1
             (16) 
 
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 ei j T iω = ∗w e                 (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 FFT’s 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 K N M M M N M= − − + −2 2 2lnρ     (18) 
( ) ( ) ( ) ( )MDL M K N M M M N M K= − − + −ln lnρ 12 2    (19) 
where: ( )
( )
( )
ρ
λ λ λ
λ λ λ
M
N M
M M N
N M
M M N
=
⋅ ⋅ ⋅
−
+ + +
+ +
−
+ +
1 2
1
1 2
1
…
…
        (20) 
(K-number of rows of the data matrix). The value of M is chosen to minimise (18) or (19). 
 
3.  Cumulant-based approach 
 
The second-, third-, and fourth-order cumulants of a zero-mean stationary random process x(t) are 
( ){ }C x tx1, = E = 0              (21) 
( ) ( ) ( ){ }C x t x tx2, τ τ= +E             (22) 
( ) ( ) ( ) ( ){ }C x t x t x tx3 1 2 1 2, ,τ τ τ τ= + +E         (23) 
( ) ( ) ( ) ( ) ( ){ }
( ) ( ) ( ) ( )
( ) ( ) ( )
C x t x t x t x t
C C C C
C C cum
x
x x x x
x x
4 1 2 3 1 2 3
2 1 2 2 3 2 2 2 3 1
2 3 2 1 2
,
, , , ,
, ,
, ,τ τ τ τ τ τ
τ τ τ τ τ τ
τ τ τ
= + + + +
− − − − +
− − =
E
      
      …
   (24) 
There are many signal processing applications where signals are either real or complex (sinusoidal) processes. 
This general class of signals can be described by 
( ) ( ) ( )x n a n si n i
i
p
=
=
∑ ω
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 (Mendel, 1991) 
( ) ( )[ ]x n j ni i i
i
p
= +
=
∑α ω φexp
1
           (26) 
where φi are independent identically distributed random variables uniformly distributed over [0, 2π], ωi≠ωj for i≠j 
and αi’s i  ωi’s are constants. 
The comparable model for real signals is  
 4  
( ) ( )y n ni i i
i
p
= +
=
∑α ω φcos
1
            (27) 
For random variables b=e
jψ
, 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 
( ) ( )[ ]C jx k k
k
p
4 1 2 3
4
1 2 3
1
, , ,τ τ τ α ω τ τ τ= − − + +
=
∑ exp     (28) 
For (27) in the case of real harmonics  
( ) ( )[
( )]
C x k k k
k
p
k
4
1
8
4
1 2 3 2 3 1
1
2 3 1
, = − − − + − − +
− −
=
∑α ω τ τ τ ω τ τ τ
ω τ τ τ
 cos cos
          + cos  
  (29) 
Additionally, ( )C x k k
k
p
2
1
2
2
1
, =
=
∑α ω τcos                 (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 (Mendel, 1991). 
Especially set τ τ τ τ
1 2 3
= = =  in (29) to see that 
( ) ( )C y k k
k
p
4
3
8
4
1
, τ α ω τ= −
=
∑ cos             (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 (Swami, 1991). 
 
4.  Experimental section 
 
4.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 doesn’t change during the 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 with 
frequencies that depend on the rotor slip.  
Assuming a steady-state asynchronous running of a synchronous machine (s=const) the stator winding current 
can be expressed (Kosty³a, 1994). 
( ) ( )[ ] ( )[ ]i I t I s t I s ts f s f b s b r s r= + + − + + − +cos cos cosω γ ω γ ω γ1 2 1  (32) 
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 
 5  
(toolbox hosa v.2.0). Coloured noise, generated by passing white Gaussian noise through a FIR filter [1,1], was 
added artificially (Fig. 3). 
  
0 200 400 600 800 1000
Sample No
-4
-2
0
2
4
A
m
p
lit
u
d
e
 
Fig. 1. Waveform of the stator winding current. 
 
 
40 42 44 46 48 50 52 54 56 58 60
Frequency [Hz]
M
a
g
n
it
u
d
e
MUSIC spectrum estimate
correlation-based
cumulant-based
 
Fig. 2. MUSIC spectrum estimates in noiseless case. 
 
 
10
6
 
10
-2
 
10
-1
 
10
0
 
10
1
 
10
2
 
10
3
 
10
4
 
10
5
 
 6  
40 42 44 46 48 50 52 54 56 58 60
Frequency [Hz]
M
a
g
n
it
u
d
e
MUSIC spectrum estimate
correlation-based
cumulant-based
 
Fig. 3. MUSIC spectrum estimates in coloured noise case (SNR=0 dB) . 
 
Table 1. Estimated signal frequencies. 
  Fig. 2  Fig. 3 
 corr. cum. corr. cum. 
f1  [Hz] 47.09 47.44 48.46 46.91 
f2  [Hz] 49.99 49.88 - 49.77 
 
 From the analysis it follows that the waveform of the stator current (Fig. 1) can be decomposed (resolved) 
into two sinusoidal components of frequencies 47 Hz and 50 Hz. Under noiseless or high-SNR conditions 
(Fig. 2) the correlation-based frequency estimates are more accurate (Table 1). Under poor SNR conditions only 
cumulant-based methods identify the correct number of signal components (Fig. 3). 
 
 
4.2  Fault operations of  inverter-fed induction motor drives 
Fault operations of inverter-fed induction motors have been also investigated. The investigated drive 
represents a typical configuration of industrial AC drives, consisting of three-phase asynchronous motor an a 
power converter composed by a single-phase half-controlled bridge rectifier and a voltage source inverter.   
 The waveforms of the inverter input current (Fig. 4) under various operation conditions, taken from 
experimental test  (Gentile, 1996), have been investigated using the correlation-based and the cumulant-based 
method. The main frequency at the inverter output was 40 Hz. Under normal conditions (Fig. 5) regular 
frequencies 240, 480, 720 and 960 Hz have been observed. Under fault conditions (switch failure) irregular 
frequencies 680 and 760 Hz have also been detected. The cumulant-based method shows a better accuracy and 
resolution than the correlation-based method. 
 
10
4
 
10
3
 
10
2
 
10
1
 
10
0
 
10
-1
 
10
-3
 
10
-2
 
 7  
0.00 0.01 0.02 0.03
[s]
0
1
2
3
[A
]
a
0.00 0.01 0.02 0.03
[s]
-2
0
2
4
[A
]
b
 
Fig. 4. Inverter input current. a) without fault  b) switch failure. 
 
0 200 400 600 800 1000
 Frequency [Hz] 
M
a
g
n
it
u
d
e
240
480 720
960
Spectrum estimates
power spectrum
4th order cumulant-based
 
 
0 200 400 600 800 1000
Frequency [Hz] 
M
a
g
n
it
u
d
e
40
80
120
680 720
760
Spectrum estimates
power spectrum
4th order cumulant-based
 
Fig. 5. Higher order spectrum estimates calculated by MUSIC method. a)without fault b)switch failure. 
10
3
 
a 
10
-8
 
10
-7
 
10
-6
 
10
-5
 
10
-4
 
10
-2
 
10
-3
 
10
0
 
10
-1
 
10
1
 
10
2
 
10
2
 
10
-7
 
10
-6
 
10
-5
 
10
-4
 
10
-3
 
10
-2
 
10
-1
 
10
0
 
10
1
 
b 
 8  
 
5.  Conclusions 
 
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 cumulant-based 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 cumulant-based methods indicate the correct 
number of signal components.  
When investigating the waveforms of the inverter input current of a inverter-fed AC drive, the cumulant-
based method shows a better accuracy and resolution than the correlation-based method. 
 
 
REFERENCES 
 
GENTILE G., ROTONDALE N., TURSINI M., FRANCESCHINI G., TASSONI C., 1996, Fault operations of inverter-fed 
induction motor. Proceedings of the 31
st
 Universities Power Engineering Conference, 2, 353-360. 
KOSTYLA P., LOBOS T., WACLAWEK Z. and PYTEL J., 1994, Identification of Asynchronous Running of 
Synchronous Machines using a Neural Network, XVII Seminar on Fundamentals of Electrotechnics and 
Circuit Theory, 2, 497-502. 
MENDEL J.M.,  1991, Tutorial on Higher-Order Statistics (Spectra) in Signal Processing and System Theory: 
Theoretical Results and Some Applications, Proceedings of the IEEE, 79, 278-305. 
NAIDU P. S., 1996,  Modern Spectrum Analysis of Time Series (Boca Raton: CRC Press), 241-336. 
SWAMI A., 1991, Cumulant-Based Approach to the Harmonic Retrieval and Related Problems, IEEE 
Transactions on Signal Processing, 39, 1099-1108. 
THERRIEN C.W., 1992, Discrete Random Signals and Statistical Signal Processing (Englewood Cliffs: Prentice-
Hall), 614-655. 
 


Application of higher-order spectra for signal processing in electrical power engineering

http://zet10.ipee.pwr.wroc.pl/getfile.py?docid=4&name=comp98&format=pdf&version=1

Application of higher-order spectra for signal processing in electrical power engineering

Abstract

Similar works

Full text

Available Versions

Repozytorium Eny Politechnika Wrocławska