Fast Fourier Transform

We have seen that N- point DFT is given by the following expression

N 1

X ( m) x n e

j 2 nm
N

n 0

j 2
N

th

Let a e , the N root of unit. Then the following relationships can be easily
derived.
1) 1 a a 2 a N 1 0
Proof: using Geometric progression series formula

x xr xr xr
2

We 1 a a a
2

N 1

N 1

x(1 r N )

(1 r )

1(1 a N ) (1 1)

0 because a N 1
(1 a )
(1 a )

2) From a m a N m 1 and a m (a m ) * 1 , we get

(a m ) * a N m
Hence a * a N 1 , (a 2 ) * a N 2 etc.
3) From the fact a 1 a 1 and a * a 1 , we have
2
* 2
N 2
etc.
a 1 a * a N 1 and a (a ) a

Now using the a-operator, DFT transformation m 0, , N 1 can be written as

follows.

X (0) x0 x1 x N 1
X (1) x 0 x1 a 1 x N 1 (a 1 ) N 1

X ( N 1) x 0 x1 a ( N 1) x N 1 (a 1 ) ( N 1)
Arranging it in a matrix format, we get the following

X (0)
X (1)

1
1

1
a 1

X ( N 1)

1 (a 1 ) N 1

x0
(a 1 ) N 1 x1
|

1 ( N 1) 2
(a )
x N 1
1

(A)

(or) stated more compactly

X P x
Where P (i, j ) ( a 1 ) (i 1) ( j 1)
The i th row of the matrix P indexes the i th frequency component while the j th
column of the P matrix indexes the j th sample. The matrix P enjoys a very
special property viz. its columns or rows are orthogonal to each other. If pi and p j
denote the i th and j th column of matrix P , then, it is easy to verify that

piH p j 0

(i j )

piH pi N
Where H- indicates Hermitian operator defined as

piH ( p iT ) * ( p i* ) T
i.e. each column is first transposed to a row and every element is then replaced by its
complex conjugate.
For a real number, the complex conjugate is identical to the original number. Hence
on real-valued vectors, Hermitian and transpose operators are one and the same.
However, for complex valued vectors the two differ.
It is now easy to verify that

P 1

1 H
P
N

1
1

N 1

H
and P

1 a ( N 1)

2 ( N 1)

( N 1) 2

Thus

X P x

(1)

(2)

1 H
P X
N

The invertibility of P , in essence captures the transform property of DFT.

In relaying, typically we are not interested in deriving all possible frequency
components. Our interest is primarily extracting the fundamental and sometimes 2 nd
and 5 th harmonic (differential protection). Since, matrix-vector product involves N 2
multiplications, and similar number of additions, we say that extracting all possible
frequency components by (1) would involve O( N 2 ) (read as order N 2 ) effort.
This effort is considered to be significantly high for real-time computing. However,
with some ingenuity, researchers have shown that the task can be achieved in
approximately

N
log 2 N computations. This fast approach to computing all possible
2

frequency transforms in discrete domain is called Fast Fourier Transform (FFT). For
example, with N 8 brute force implementation of (1) requires 64 complex
multiplications which can be reduced to 12 multiplications with FFT. As we would
not have much use for FFT in this course, we will not pursue this topic any further.
Rather, we now establish an equivalence between two very well known transforms
viz. multiple DFT (or) FFT and sequence component transformation. [for the
remaining text through out we will refer equation (1) as FFT]

Equivalence of FFT and N-phase sequence component transformation:

We will first review N-phase sequence transformation. Consider a N-phase
(balanced or unbalanced) system ( N 3 ). Let the phasors in phase domain e.g. (
Va , Vb ,Vc for 3-phase system) be represented by x 1 , x 2 ,........x N . Then, these

phasors can be expressed as a linear combination of N-set of balanced N-phase

systems as follows.
0-sequence component

x1 x 2 x N
X0
N

There are n- such phasors in the 0-sequence system

each of equal magnitude and angle

(+ve) 1-sequence component

x1 x 2 a x 3 a 2 x N a N 1
X1
N

(A1)

There are n- such phasors in 1-sequence system. If X 1,1 is taken as reference, then
X 1,2 is equal in magnitude to x 1 ( a ) but **** by angle

2
N 1
i.e. X 1,2 X 1,1a
N

(-ve) 2-sequence component

X2

x1 x 2 a 2 x 3 a 4 x N (a 2 ) N 1
N

(A2)

This system is obtained by relating lag in consecutive phases by operator a 2 (or a N 2

).
|
|
X N 1

x 1 x 2 (a 2 ) N 1 x 3 (a 4 ) N 1 x N a ( N 1)

(3)

Figure 1(a, b and c) visualize the system for a 3-phase system. Expressing these
equations in matrix format we obtain the following equations.

1
X0
1
X
1
a
1

1
| 1 a2
N

|
|

|
X N 1
1 a ( N 1)

x1

a N 1 x 2
a 2 ( N 1) |

( N 1) 2

a
xN

1
a2
a4
a 2 ( N 1)

Thus, writing compactly

X P 1 x

or

x P X

(B)

Where X X 0 , X 1 ,....., X N 1

x x1 , x 2 ,....., x N

Thus, from (1) and (B) we conclude that FFT and sequence transformation (defined
from sequence domain to the phase domain) involve the same transformation matrix
P. Hence, the two transforms are mathematically equivalent. In particular for N=3.
With a- phase as reference phasor, we see the following equivalence relationships.

Va 1 1
V 1 a2
b
Vc 1 a

1
a

V0
V
1

a 2 V2

and from DFT, we get

V ( 0) 1 1
V (1) 1 a 2

V (2) 1 a

1 Va
a Vb
a 2 Vc

This mathematical equivalence brings out an important concept viz. transformations

and decompositions done via orthogonal matrices can have multiple interpretations.
However, there is one important difference between the two transformations. The
samples x 0 ,......., x N 1 in signal processing are real numbers while corresponding
phasors in sequence analysis are complex numbers. Thus, while there is a redundancy
in information in DFT domain which leads to DFT symmetry property, there is no

such redundancy in sequence domain. Hence, in sequence domain we do not come

across such a property.1

Estimation of Frequency Lecture 37

So far while discussing the phasor estimation problem, we have assumed that
frequency of the power system remains fixed at its nominal value ( f 0 50 / 60 Hz ).
Hence, we have also fixed the sampling frequency to Nf 0 f s , where N-point DFT
is used in estimation. However, during disturbance and even in steady state to a
certain extent, the frequency varies. Thus, we expect the phasor estimation under
constant frequency assumption to be erroneous. Under such situations, how good is
our estimate of the phasor? We now plan to answer this question. As a by product of
the analysis, we will also develop a frequency estimation technique which can be
used in under frequency and rate of change of frequency relays. To simplify
presentation the analysis is done in multiple steps. First, we determine the DFT of
complex exponential signal at frequency f .

DFT of Complex Exponential

Let the signal be given by

x(t ) e

j 2 k f 0 t

Let N-sample be captured in P-cycles. P can be a positive integer or even a positive

real number. Then, sampling interval ( t ) is given by P Nf , sampling speed by
0
Nf 0

xn e

f s . The discrete sample at t nt is given by

nP
j 2 f 0
N f0

j 2 nP
N

(n 0)

Thus, mth DFT component corresponding to frequency mf 0 , is given by

N 1

X (m) xn e

j 2 n m
N

n0

The mathematical equivalence of two should put the reader at ease with both these transformations,
irrespective of which one he came across first.

N 1

j 2 n
( Pm)
N

---------- (3)

n0

It should now be fairly easy to compute the summation in (3)

2

Let q e j N ( P m ) , a constant once P and m are fixed.

------- (4)

N 1

n
N 1
Then X (m) q 1 q q

---------- (5)

n0

1 q N 1 e j 2 ( P m )

=
2
j
( P m)
1 q
1 e N

e j ( P m )
e

e j ( P m ) e j ( P m)

( P m )
N

( P m)
N

j ( P m )(1

1
)
N

j ( P m )(1

N e

1
)
N

1
)
N

( P m )
N

e j ( P m ) e j ( P m)
e

j ( P m )(1

( Pm)
N

2j

( P m )
N

2j

sin ( P m)

sin ( P m)
N

sinc( ( P m)) ( P m)

sinc( ( P m)) ( P m)
N
N

j ( P m )(1

1
)
N

sinc ( P m)

sinc( P m)
N

sin( x )
Note that MATLAB defines sin c( x)
x . So this definition is at variance
with our definition viz. sin c( x) sin( x ) x . While programming in MATLAB one
must carefully account for it.

Fig (1)
Fig (3)

Fig (2)

Fig (3)

Fig (1) shows the envelope of response for P 1 (or) f f 0 as a function of m

treated as a continuous variable. Since in DFT we can also choose m to be discrete
number, this plot is sampled at discrete points ( m 0,1,2,....., N 1 ). The plot
indicates that fundamental is extracted correctly as expected and at all harmonic
frequencies, the DFT gain i.e. X (m) is zero. Infact, fig(1) is identical to the
frequency response of full-cycle Fourier algorithm.
Fig (2) shows the envelope of response for 3 possible frequencies 49, 50, 51 Hz.
Where nominal frequency f 0 50 Hz . It is seen from the figure that magnitude
response is more or less identical.

Finally, Fig (3) shows a set of response when frequency deviation from nominal is
large enough i.e. f = 40, 50, 60 Hz. It can now be seen easily that if we sample the
envelope of 40Hz signal at 50 Hz ( m 1 ), a finite non-zero value different from unity
is returned. This effect is knows as DFT leakage. It is said that the energy in
frequency bin f 40 Hz has leaked into frequency bin f 50 Hz . Though, in DSP
literature, DFT leakage is considered to be undesirable from the relaying perspective,
it turns out to be a boon in disguise.
As seen in fig (3), the DFT magnitude leakage for 1 Hz frequency deviation around
nominal frequency say 50Hz, is practically negligible. For f = 51 Hz, it equals

51

sin c ( ( 50 1))
1 pu = -0.0065. This indicates that magnitude of phasor

sin c ( ( 50 1))
50

estimation in relaying applications by DFT with constant sampling frequency is

robust. We thus say that primarily magnitude leakage can be neglected. However, the
phase angle of DFT tells another story. Note that
For P

51
1
and N 12 X (m) ( P m)(1 )
50
N

Thus X (1)

180 0
1
180 11
(51 50)(1 )
3.3 0
50
12
50 12

Thus an error of 3.3 0 is introduced in phase computation.

DFT of real COSINE signal

In practice, we are interested in DFT of real sine (or) cosine signal and not the
complex exponential signal. However, deriving the response for real cosine (or) sine
signal from complex exponential is not very difficult. From the fact that

e j 2 f t e j 2 f t
and by following similar steps as in the case of
x (t ) cos 2 f t
2
complex exponential, we get x(n)

1 j
e
2

2 n P
N

2 n P
N

Thus

1 N 1 j 2N n ( P m ) j 2N n ( P m )
X (m) e
e

2 n 0

We have already deduced that

N 1

2 n
( P m)
N

Ne

j ( P m )(1

1
)
N

n 0

sinc ( P m)

sinc
( P m)
N

Similarly it can be shown that

N 1 j 2 n ( P m )
N

Ne

j ( P m ) ( P m )
N

n 0

sinc ( P m)
sinc ( P m)
N

Thus,

X c ( m)

N j ( P m )(1 N1 ) sinc[ ( P m)]

N j ( P m )(1 N1 ) sinc ( P m)
e
e
2
sinc[ ( P m) ] 2
sinc ( P m)
N
N

With P 1 , ( f f 0 1) , it can be seen that

sinc( ( P m)) sinc 2 0

while,

sinc( ( P m)) sinc (0) 1

and

N j ( P m )(1 N ) sinc[ ( P m)]

Pm
we
sinc(
) sinc (0) 1 ; Hence, X c ( m) 2 e
sinc[ ( P m) ]
N
N
1

conclude that even with real sine or cosine signal, any reasonable deviation of the
signal from the nominal value, leads to negligible magnitude leakage. However,
phase angle error for 1Hz deviation in frequency is approximately 3.3 which is not
negligible.
Figure shows the envelope of magnitude response of X c (m) for different frequencies
40, 50, 60 Hz.

f s T0

th
f 0 m P1 P 1 w

Estimation of Frequency
Our aim, now, is to develop a method of estimating power system frequency using the
recursive DFT approach to phasor estimation. We plan to show that in the moving
window approach; the phase angle estimated by recursive DFT approach rotates at a
speed proportional to the deviation from the nominal frequency. In turn, this deviation
can be assessed by measuring the rate of change of phase angle. In the previous
section, we have derived the DFT leakage for complex exponential and real cosine
signals. We will now generalize the DFT computation of complex exponential so that
it can handle moving window concept.
Generalized DFT of Complex Exponential
Following the methodology used while generalizing DFT, we can write generalized
DFT, we can write generalized DFT for wth window (window no. corresponds to the
sample no. of first sample in the window) as:
X w (m)

N w 1

n w

where x e
n

xn e

j 2 n
m
N

--------(9)

j 2 nP
N

Thus, xn in (4), we get

X (m)
w

N w 1

j 2 n
( Pm)
N

n w

j 2
( Pm) w
N

[1 q .... q N 1 ]

where q e

j 2
( Pm)
N
N 1

Summation

is the DFT corresponds to window no. 1 substituting from eqn. (9),

i 0

we obtain X w (m) e

j 2 n
( P m ) w
N

X c (m)

---------- (10)

If m = P, then it is clear that DFT X w (m) is stationary i.e. independent of window

number. In relaying applications where we are interested fundamental phasor

estimation, we choose f s Nf 0 . Hence, we can expect obtain stationary DFT in (10)

when m = P = 1.
With the knowledge, of the generalized DFT of the complex exponential, we can now
derive the generalized DFT expression for real cosine signal.
Generalized DFT of Real Cosine Signal
For a real cosine signal, with moving window concept eqn (7) can be generalized as
follows:
X w (m)

j 2 n
( P m )
1 N w1 j 2N n ( P m )
N
e

2 n w

--------- (11)

1 j 2N w ( P m ) N 1 j 2N n ( P m ) 1 j N2 w ( P m ) N 1 jN2 n ( P m )
e
e
e
e

2
2
n 0
n 0

1
N j 2N w ( P m ) j ( P m ) 1 N sin c P m
N j N2 w ( P m ) j ( P m ) 1 N sin c P m
e
e
e
e
2
2

P m
sin c
P m
sin c

(12)
Following the similar reasoning as outlined in the previous section of the DFT of real
cosine signal, we can neglect the second term in (12). This leads to following
approximation.

N j 2N w ( P m ) j ( P m ) 1 N sin c P m
X (m) e
e
2

sin c
P m
N

observe that P = m = 1, we obtain stationary DFT because e

----------- (13)
j 2 w
( P m )
N

= 1. However,

when frequency deviates from nominal frequency i.e P 1 , then X w (1) , the DFT
estimate for nominal frequency, will start rotating along with the window. From, (13)
j 2
( P 1)
X w1 (1)
N

e
we can derive that
w
X (1)

Since fundamental phasor estimated for window number w is given by

w 1
j 2
j2
( P 1)
w
X
N
X
X (1) , we deduce that

e
e
w
N 2
X
w

j 2 f f 0
N
f0

j 2 T0
N

w 1

e j 2 f t

----------- (14)

where t is the sampling time interval. Summarizing, if we set sampling frequency

for a sinusoidal signal to provide N-samples per cycle at the nominal frequency f 0 ,
(say 50 or 60hz), and keep the sampling frequency invariant of the frequency of
sinusoid, then the phasor estimated by moving window approach rotates at a speed
proportional to f . This rotation will be in anti-clockwise if

f >0 i.e. f > f 0 .

Conversely, it will be in the clockwise direction if f <0 i.e. f < f 0 .

If we monitor, this phase rotation, then from the proportionality relationship of (14),
we can estimate the frequency f. If denotes phase-angle, then from (14), we can
obtain the rate of change of as follows:
d w1 w

2f
dt
t
Discussion:
There are many advantages associated with the above algorithm. The method is not
based on zero crossing times and it is immune to noise and harmonics. Instead of
using single phase quartiles, one can estimate the +ve sequence component and derive
frequency from it. Such an approach will use all the three phase voltages and hence
will have better noise rejection properties. At harmonic of nominal frequency
X w (m) 0 and hence, the above approach will reject frequencies m f 0 completely.
Measurement of Frequency
f 0 accuracy measure frequency we integrate eqn. (15) until, a change equal to
approx. 0.5 radians is accumulate this change is T, then
f

0.5
2 T

If we assume that frequency computation will be further averaged over four

measurements to smoothen out noise, then time to compute deviation f is given by
T

0.32
sec
f

Thus, a frequency deviation of 1Hz is detected in 0.32sec while a frequency deviation

of 0.1Hz is detected in 3.2sec. This is an illustration of standard speed versus
accuracy conflict in relaying.
This technique can be recommended for development of under frequency and rate of
change of frequency relays.
mP