Discrete Fourier transform computation for digital relaying

Discrete Fourier transform computation for digital relaying A G b m e z Expbsito and J A Rosendo Mac|as Departamento de Ingenieria El~ctrica, Universidad de Sevilla. Avda. Reina Mercedes s/n. 41 012 Sevilla, Spain

J L Ruis M a c i a s Sevillana de Electricidad, Avda. Borbolla 5, 41004 Sevilla, Spain

In this paper several ways of computing the discrete Fourier transform are presented. Starting from the conventional non-recursive version, three recursive algorithms are obtained. Even though these recursive filters are computationally efficient their steady-state response is not always accurate. So, an efficient non-recursive version is proposed which requires computing only the imaginary component of the transform. An example is included in order to compare the operation counts of the different schemes. Keywords: on-line computer systems, data gathering and analysis, software design

I. I n t r o d u c t i o n Currently, many protection devices based on digital technology are commercially available. It is hoped that, in the future, the use of microprocessor-based relays will increase gradually, as they offer several advantages such as reduced cost, better performance, flexibility, adaptability, self-checking capability and the possibility of integration in the digital environment of future substations 1'5. Computer relays rely on numerical methods to model the input waveforms and/or the differential equations of the system being protected ~. Some of the numerical techniques proposed for waveform modelling2'3'6'7 are •

• • • • •

II. Standard non-recursive DFT To begin with, the standard non-recursive form of the DFT will be reviewed. This form will be later used for comparison purposes. Let Xo, xl . . . . . xN_ 1 be the N real samples of the signal whose spectrum is required. The kth harmonic is given by 8

least-squares fitting, discrete Fourier transform (DFT), rectangular Fourier transform, Walsh transform, Haar transform, Kalman filtering.

N-I

~k=

The Kalman filter is a sophisticated tool that should be used in those situations where the covariance of the measurement noise is not constant (Reference 3, p 105). Received 30 December 1992; accepted 11 November 1993

Volume 1 6 Number 4 1 994

Curve fitting algorithms may be advantageous if certain parameters, such as DC offset time constants, are of interest. Otherwise, when only fundamental and harmonic phasors are required, transform-based methods are usually preferable in terms of computational cost and frequency response. Although rectangular-type transforms are implemented using only additions and subtractions, a certain number of multiplications must be carried out to obtain Fourier coefficients from the rectangular coefficients. In fact, based on the comparative analysis of Reference 7, it can be concluded that the DFT is the best choice for computing the fundamental and harmonic components of input waveforms. This paper is exclusively devoted to the DFT. Starting from the conventional, non-recursive version, several recursive algorithms are obtained. It is shown that the real and imaginary components of the recursive DFT lead to a certain type of second-order bandpass filters. These recursive expressions are finally used to obtain an improved non-recursive scheme which resorts to only two consecutive values of the imaginary component of the DFT. An example will be discussed along with each version of the DFT in order to compare the operation counts.

~ xie-~Ol

(1)

i=0

with 0-

2nk N

0142-0615/94/04/0229-05 © 1994 Butterworth-Heinemann Ltd

(2)

229

Discrete Fourier transform c o m p u t a t i o n . A. G r m e z E x p r s i t o et al

In computer relaying the samples are continuously updated as the data window moves. When a new sample is available the oldest one is discarded. Hence, at instant n, the window includes the N values x . _ N + ~ , x , _ u + z, . . . . x,, and equation (1) can be rewritten as ~--. =

xl e -j°"-"+N- '~

~

(3)

i=n-N+l

where, for simplicity of notation, the superscript k has been omitted (it is still implicit in the value of 0). If ~ . is expressed in rectangular coordinates, ~ . = C. + iS., its components are C, =

x i c o s ( ( / - n + N - 1)0)

~

(4)

i=n-N+l

S, = -

x i sin((/ - n + N -

~

1)0)

(5)

i=n-N+l

and its amplitude F,2 = C,2 + S,2

(6)

III. Improved non-recursive DFT Usually, only a few harmonics are of interest in computer relaying. Hence, F F T techniques find little application in their conventional form. However, noting that, because of the window shift and the symmetry of equations, some additions that must be carried out at a particular instant have already been computed at an earlier step, the operation count can be reduced by means of auxiliary storage. For simplicity, instead of devising general expressions, the former example will be used to illustrate this technique. Example

Defining the following intermediate variables w. = x. -- x._ 8

(13)

y. = x. + x._ 8

(14)

Vn = Yn -- Yn-4

(15)

equations (7) to (12) can be rewritten as C 1 = -w.-7

Exarnple

Throughout the paper it will be assumed that the amplitude of the first, second and fifth harmonics are of interest and that N = 16. This would be, for instance, the case if a current-restrained transformer differential relay were to be developed 6. For this particular case the harmonic components are

+ k2(wn_l - w . _ ~ )

(16)

+ k l ( w . - Wn_6) Jr k3(Wn_ 2 - Wn_4)

$1 = % - 3 + k z ( w . - 1 + w.-5) + k l ( w . - ~ + w.-4) + k3(w. + w.-6) C2 = - - Y . - 3

+ Y n - 7 + k2(yn -

= - - V n _ 3 Jr k 2 ( v n -

Y n - 2 --

(17)

Y.-4 +

Yn-6)

v._2)

(18)

C1 = X n - 1 5 -- X n - 7 Jr k 2 ( x n - 1 3 - X n - 9 -- X n - 5 + X n - 1 ) Jr k l ( X n - z ~ . -

$2 = Y.-1 - Y.-5 + k2(Y. + Y.-2 - Y . - , ~ - Y.-6)

x . _ 8 - - x n _ 6 Jr xn)

Jr k3(Xn_12 -- Xn_lO - - Xn_ 4 Jr X n _ 2 )

(7)

C~ = -w.-7 + kz(-w.-i S1 ~

--Xn_ll

+ k2(-x._13

- - Xn_ 9 Jr X . - 5 Jr X . - I )

Jr k l ( - - x n _ 1 2

- - X n _ l O Jr X n - 4 Jr X n - 2 )

+ k a ( - X . _ 1~ - x._ s + x~_ 6 + x.)

--

(10)

C 5 ~ Xn_15 -- Xn_ 7 Jr k z ( - - X n _ 1 3

85 ~

--Xn_ll

--

X n - l O -- X n - 4 Jr X . - 2 )

Jr

" " " Jr

Xn-1

e-J0(N-2)

-~- Xn - 1 e - jO(N- 1)

and subtracting the second of these equations from the first one multiplied by e -j°. This yields ~ - . e - j ° - o~._ 1 = x . - x . _ N

+ k 3 ( x n - x 2 + X n - I O -- Xn_ 4 -- X n _ 2 )

(12)

with k 1 = cos(22.5), k 2 = cos(45), k 3 = cos(67.5). Taking advantage of the common terms shared by the first and the fifth harmonic components, a total of 49 sums and 18 products are required. These include the 3 sums and 6 products implied by equation (6).

230

e-J°

~'~n-1 --- X n - N Jr X n - N + I e - J ° Jr " " ' Jr X n - 2 e - J ° ( N - 2 )

(11)

Jr k 2 ( X n - 1 3 Jr X n - 9 -- Xn_ 5 -- X n _ l )

- x . _ 8 + x._ 6 + x.)

"~- X n - N + 2

Jr X n e - jO(N- 1)

Jr X n - 3

+ k1(-x._l,~

The complex recursive form of the D F T can be easily obtained by writing equation (3) at two consecutive instants ~'~n ~- X n - N + I

Jr X n - 9 Jr X n - 5 - - X n _ l )

+ k 3 ( - - X n _ 1 4 + X n - 8 Jr X n - 6 -- Xn) Jr k l ( X . . - 1 2

The number of operations required now to compute (F.~)~, (F.2)2 and (FS.) z is 28 additions and 18 products, including the 3 additions of (13)-(15). This means almost a 50% reduction in the number of additions.

IV. Recursive forms of the DFT

Xn_12 Jr X n - I O

Jr Xn_ 8 -- Xn_ 6 -- Xn_ 4 Jr Xn_ 2 Jr Xn)

(21)

(9)

13 Jr X n - 9 -- X n - 5 Jr X n - 1

Jr k z ( - - X n _ 1 4

(20)

- w.-5)

+ k 3 ( - w n _ z - w . _ , , ) + k l ( w . + Wn_6)

(8)

Jr k 2 ( X n - 1 4 - - Xn_12 -- X . _ l O Jr X . - 8 Jr Xn_ 6 -- Xn_ 4 -- Xn_ 2 Jr Xn)

n + Wn_6) + k l ( W n _ ~ - - Wn_4)

S ~ = w.-3 + k ~ ( - w . - 1

C 2 = Xn_15 -- X n _ l l Jr X n - 7 -- Xn_ 3

--Xn

+ w._~)

Jr X n - 3 + k3(-w

82 ~

(19)

= V n _ 1 + k2(~) n + Vn_2)

or, equivalently, ~, = (ft.-

1 + x. -

x . _ ~ ¢ ) e j°

(22)

For practical purposes, equation (22) is decomposed into its real and imaginary components, C. = (C,_ 1 + Ax.) cos 0 - S._ 1 sin 0

(23)

Electrical P o w e r & E n e r g y S y s t e m s

Discrete Fourier transform computation: A. G6mez Exp6sito

S. = (C._ 1 + Ax.) sin 0 + S._ 1 cos 0

(24)

where Ax. = x. - x._u. Using matrix notation these equations can be written as

[C.] = ~ c o s 0 - s i n 0 1 F C , _ l l S.

Lsin0

+ Fcos01A x cos0 _J[_S._13 LsinOJ "

(25)

Equation (25) represents a first-order recursive linear system that can be easily transformed into two uncoupled second-order recursive equations. The application of the z-transform to (25) yields

[sC]= F c°s0 L sin 0 ÷

-sin01z-l[sC ] cos 0 _1

Fcos 01

(26)

and, rearranging,

[~]=(/

1

~cos0-sin01z_ -] - Lsin 0 cos0 _J 1./ x Fc°s~](1 - z-N)X Lsin

The following transfer functions are obtained from the above expression

C

--

X S X

(cos0-z-1)(1-z -u) 1-2cos0z

--

(27)

- l + z -=

sin 0(1 - z -u) 1 - 2 c o s O z - l + z -2

(28)

These transfer functions correspond to recursive bandpass filters made up of a comb filter (having N sampling delays) cascaded with a resonator. The combination comb-filter-resonator is a good example of a recursive design with finite impulse response (FIR) a. Sometimes, FIR filters are exclusively associated with non-recursive designs. In the discrete-time domain, equations (27) and (28) lead to

C. = cos 0(2C._ 1 + Ax.)

-

Cn_ 2 --

AXn_ 1

S, = 2 cos 0 S,_ x - S,_ 2 " J r sin OAx,

recursive version that gives directly the harmonic amplitude (equation (31)), a total of 14 additions, 6 products and 3 shifts are required. This reduction in the operation count is due to the fact that instead of using equation (6), the three harmonic amplitudes are directly obtained, without extra computations. Actually, the number of arithmetic operations given for these recursive filters is independent of the number of samples per period. Note, however, that those expressions are only valid for full-cycle DFT. In view of these results it could be concluded that recursive versions are preferable from a computational perspective. However, some problems are likely to appear in practical implementations making use of such recursive equations, particularly if integer arithmetic is adopted. On the one hand, the two complex-conjugate poles of equations (27) and (28) lie exactly on the unit circle. Hence, rounding the values of 2 cos 0 may lead to unstable algorithms. Fortunately, this problem can be solved with a careful design which takes into account two restrictions: the resonator poles must lie inside the unit circle and, at the same time, the comb filter zeros should cancel out these poles as perfectly as possible, so that the frequency response is not deteriorated. On the other hand, even though the stability problem is solved, the long-term behaviour of recursive filters is not always what could be expected for linear, time-invariant systems. For instance, Figure 1 shows the transient and steady-state response of the filter which computes the fundamental harmonic amplitude when the input is a fundamental frequency sine of constant amplitude. While the transient response is correct, a decaying term can be noted in the steady-state response. Ideally, the output should remain constant once the transient period has vanished. Other unexpected behaviours, such as limit cycles, residual outputs in absence of inputs, etc., have been also observed. A detailed discussion of these anomalous results are beyond the scope of this paper (the interested reader is referred for instance to Reference 9, Chapter 14). In summary, filter coefficient quantization (trunca-

(29) (30)

where Ax,_ 1 = x,_ t - x, _ ~u+ 1). Furthermore, a third recursive version can be obtained when only the harmonic amplitude is of interest. Substituting equations (23) and (24) for C, and S, in equation (6) yields F.~ = F.Z_l + Ax.(2C._ 1 + Ax.)

(

(31)

This equation allows the recursive computation of the squared harmonic amplitude based on its real component.

Example For the particular case under study, the coupled, first-order, recursive version (equations (23) and (24)) requires 13 additions and 16 products to compute the three harmonic amplitudes (the fact that cos 45 = sin 45 saves 2 products). On the other hand, using the decoupled, second-order equations ((29) and (30)) as many as 20 additions, 14 products and 3 shifts must be carried out. Finally, using equation (29) to compute C, and the

Volume 16 Number 4 1994

et al

/ / lo o

,/

!

lo 1

lo ~

lo ~

lo ~

SAMPLES

Figure 1. Response of the fundamental harmonic fundamental frequency sine (/V = 1 6 )

filter to a

231

Discrete Fourier transform computation. A. Gomez ExpOsito et al tion or round-off), imposed by finite wordlengths, gives rise to a nonlinear staircase function which forms part of the feedback loop. Hence, linear systems theory cannot completely predict filter responses in such situations. All these potential problems make the recursive DFT less reliable than the non-recursive versions. An acceptable combination could be interleaving the non-recursive DFT with several computations of the recursive algorithm so as to reset biasing errors in the filter output.

V. Non-recursive component

DFT based on a single

A compromise version of the D F T is presented in this section intended to be as reliable as the standard D F T but, at the same time, almost as fast as the recursive algorithms. The following ratio can be obtained from equations (27) and (28): C S

-

cos 0 - z-X

VI. H a l f - c y c l e DFT Throughout the paper it has been implicitly assumed that the sampling window spans a full cycle of the fundamental frequency. If the number of samples, m, is less than N, the counterpart of equation (22) is ~v. = ( ~ . _ 1 + x. e ~j,.0 _ x,,_,.) e j°

For an arbitrary value of rn it would be a cumbersome task to develop general expressions like those obtained in previous sections. Nonetheless, it is quite easy to repeat the former development for the practical case rn = N / 2 , i.e. for a half-cycle window. In such a case equation (36) reduces to ~n

C X S

F~ sin 2 0 = (C~ + S~) sin 2 0 2 = S~ + S,_ 1 - 2S~S~_ ~ cos 0

(34)

Except for a scaling factor, sin 2 0, equation (34) gives the harmonic amplitude with a minimum number of operations (S~_ 1 is available from the former step). Note that both equations (33) and (34) are non-recursive. Either recursive or non-recursive expressions can be adopted to compute S~ but, as was discussed earlier, non-recursive algorithms are more reliable in the long term. Equation (34) can be looked at from a completely different perspective. Assuming two samples of a sinusoidal signal are available S. ~ = F. sin q5

S. = F. sin(~b + 0) where 0 is the phase difference between the two values, equation (34) is also obtained by simply eliminating ~b. Resorting to equations (23) and (24) an alternative to equation (34) based on the real harmonic component, C., can also be obtained. F.2 sin 2 0 = C2. + (C._

1

"q-

Ax) 2

- 2C.(C._ ~ + A x . ) cos 0

eJ°

(37)

~

--

(z 1 - cos0)(1 + z -~/2) . . . . . . . . . . _ - 1-2cos0z-~ +z 2 - s i n 0(1 + z -N/2) 1-2cosOz-l+z

-2

(38) (39)

Hence, equations (33) and (34) are still applicable to the half-cycle window case. Depending on the expression used to compute S,, equation (34) will provide either the full-cycle or the half-cycle harmonic amplitude.

VII. Conclusions The important topic of D F T computation has been dealt with in this paper. The key points for a DFT algorithm to be chosen in computer relaying applications are speed and reliability. Even though recursive D F T filters are computationally efficient, their long-term responses may be incorrect owing to finite-wordlength effects. On the other hand, conventional, non-recursive D F T algorithms require a lot of computations which can be unacceptable depending on the number of samples, number of signals and computing power. After reviewing some recursive DFT formulae, a non-recursive version has been presented, both for the full-cycle and half-cycle cases, which requires computing only the imaginary component of the DFT. This results in a considerable saving in the number of operations, without penalizing reliability.

VIII. A c k n o w l e d g m e n t This work originated in a project financially supported by the Spanish OCIDE under PIE grant No. 132209.

IX. References (35)

However, it leads to a larger number of operations. Example

Using the improved non-recursive version of Section III to compute S.~, S.~, S.~, and equation (34) to compute the harmonic amplitudes, a total of 20 additions and 15 products are required, provided the scaling factor, sin ~ 0, is of no importance.

232

X

(33)

This equation provides a simple means of computing the real harmonic component based on two consecutive values of the imaginary component. For those cases in which only the harmonic amplitude is required, the following expression, based on equation (33), can be used

('~n - 1 -- Xn -- Xn-N/2)

are

which, in the discrete-time domain leads to, C, sin 0 = S, cos 0 - S,_ 1

=

Clearly, equations (23) and (24) remain valid if the definition of Ax. is changed to A x . = - - x . x._mz. Accordingly, the counterparts of equations (27) and (28)

(32)

sin 0

(36)

1 IEEE Tutorial Course 'Computer relaying' IEEE publication 79EHO148-7-PWR (1979) 2 IEEE Tutorial Course 'Microprocessor relays and protection systems' IEEE Publication 88 EHO269-1-P WR (1988) 3 Phadke, A and Thorp, J Computer relaying for power systems John Wiley (1988) 4

Slatem, R 'Modern trends in power system protection' IFA C Power Gener., Dist. and Protection Pretoria (1980) pp 81-96

Electrical Power & Energy Systems

Discrete Fourier transform computation: A. G6mez Exp6sito et al 5 Saha, M and Pfistner, R 'The state of the art of numerical relaying' Electr. Power Energy Syst. Vol 13 No 2 (1991) pp 91-99

relaying algorithms for the differential protection of three phase transformers' IEEE Trans. Power Syst. Vol 3 No 3 (1988) pp 1378-1384

6 Rahman, M and Jeyasurya, B 'A state-of-the-art review of transformer protection algorithms' IEEE Trans. Power Deliv. Vol 3 No 2 (1988) pp 534-544

Lynn, P and Fuerst, W Introductory digital signalprocessing with computer applications John Wiley (1989)

7 Habib, M and Matin, M 'A comparative analysis of digital

Volume 16 Number 4 1994

Phillips, C and Nagle, H Digital control system analysis and design Prentice-Hall (1984)

233