Digital Power System Protection
Digital Power System Protection
Digital Power System Protection
SystemProtection
S.R. Bhide
Associate Professor of Electrical Engineering
Visvesvaraya National Institute of Technology
Nagpur
Delhi 110092
2014
Digital power system protection
S.R. Bhide
2014 by PHI Learning Private Limited, Delhi. All rights reserved. No part of this book may be
reproduced in any form, by mimeograph or any other means, without permission in writing from the
publisher.
ISBN-978-81-203-4979-7
The export rights of this book are vested solely with the publisher.
Published by Asoke K. Ghosh, PHI Learning Private Limited, Rimjhim House, 111, Patparganj
Industrial Estate, Delhi-110092 and Printed by Star Print-O-Bind, F-31, Okhla Industrial Area
Phase I, New Delhi-110020.
To
the memomy of My Parents
Pushpa and Ramkrishna Bhide
Contents
Prefacexi
Appendix 247251
References 253256
Index 257259
Preface
On the one hand the electronic technology is changing at a breath taking speed. The relentless
march of Moors law is making it not only possible to put computing power into every
imaginable existing gadget but is also giving birth to entirely new gadgets and gizmos. Todays
note-book puts yesterdays supercomputer to shame. The internet and the digital computer
have become ubiquitous. Interestingly, advances in the VLSI technology have also made the
digital computer largely invisible by embedding it deeply into the innards of the system that
it is controlling. You have digital computers embedded in every imaginable object.
Thus, the current trend is to replace the mechanical parts with electronic hardware and the
electronic hardware, with software. It is no wonder, therefore, that the field of power system
protection has also been caught up in the throes of the digital revolution. The result is that, the
protective relay has been completely transformed into an intelligent electronic device (IED).
On the other hand, equally dramatic changes have taken place in the general power
system scene. The electric utility business has been deregulated in large parts of the world.
Environmental concerns are giving rise to the proliferation of non-conventional energy sources
like wind, solar, bio-mass etc. There is more variety in the sources of power generation today
than a decade back, leading to challenges and opportunities for the protection engineer.
All these factors make the present times very exciting and fascinating for the practicing
protection engineers. However, for the new comer, to put it mildly, it becomes confusing and
chaotic. This is because an appreciation of three major disciplines is required before one can
make sense of all the digital brouhaha. Firstly, one must never lose sight of the underlying
electrical principles. The various entities that are being protected; the generators, bus-bars,
transmission lines and transformers have their own idiosyncrasies, i.e., unique set of responses
and behaviours that set them apart from every other entity. These must first be thoroughly
understood. Secondly, since the implementation of all the ideas is through the digital signal
processing technology, one needs thorough appreciation of the art and science of DSP at the
level of algorithms. Thirdly, we need a machine; the DSP microprocessor to practically
xi
xiiPreface
implement the algorithms. Thus, appreciation of this physical layer of digital technology is
equally important.
It is not possible to do justice to all the three areas mentioned above in a text like this.
We have mainly attempted to introduce the reader to the middle-ware, i.e., various digital
algorithms which are at the heart of the digital protective relays.
The text is aimed at the final year undergraduate and postgraduate students of Electrical
Engineering who have already undergone courses on Power System Protection, Signals and
Systems and Microprocessors. A course on Numerical Methods will be a definite advantage.
The practicing engineers, who routinely install, operate and maintain these systems, will
also find the text useful in deepening their insight into the ideas that make the digital relays
intelligent.
The author feels indebted to all the researchers, inventors and pioneers who have enriched
the state of the art by their contributions and made this field extremely fascinating.
The author wishes to thank his teacher and guide Prof. Y.G. Paithankar for introducing this
fascinating field to him in the early 1990s. His love for the subject, enthusiasm for learning,
his simplicity and pure humanism has been very inspiring. The author will always cherish the
memories of his association with Dr. Paithankar.
The author would like to gratefully acknowledge the affection and encouragement received
from all his colleagues at the Electrical Engineering Department, VNIT, Nagpur.
Last but not the least, the numerous students who have taken the courses on DSP
Applications to Power Systems and Digital Protection, over the years, deserve special thanks
for being catalysts in enriching authors understanding of the subject and enhancing his patience.
S.R. Bhide
1
Evolution of Power System Protection and
the Emergence of Digital Relaying
1
2Digital Power System Protection
The technology used for making ammeters, at that time was based on force experienced
by a current carrying coil placed in a magnetic field. We call this the electromechanical
technology which made measurements through the balancing of deflecting torque against a
restoring or restraining torque. The restraining torque could be developed by a spring, by
gravity or by another current etc.
It is easy to see that electromechanical measuring instruments evolved into their new avatar
of dedicated and very sensitive and precision mechanisms called protective relays. Somewhere
along the evolutionary path, the indicating needles were dropped. It is interesting to note that
modern numerical trelays do store the instantaneous values of voltage and current and can
display them on request of the user. Thus, the measuring instrument from which the relay
has evolved is still there. This reinforces the idea that modern relays can be looked up on as
measuring instruments endowed with numerical computational and logical decision making
capabilities.
electronic devices were not as reliable as their present day counterparts. Further, user experience
and confidence in the new technology was lacking. It was only during the early 1960s that
solid-state analog technology was accepted by the industry.
Even though digital computers were around right from the early 1940s, they were very
bulky, in fact monstrous by todays standards, consumed astounding amounts of power and
lacked speed for real time applications demanded by a protective relay. Not only that, they
were so expensive that only big corporations and government defence establishments could own
them. In those days, it was unthinkable to use a digital computer in a relay. This is borne out
by the fact that early research papers on use of digital computer for protective relaying go to
great lengths to justify their use and almost sound apologetic about using digital computers for
such a mundane task as protective relaying!
Application layer
Presentation layer
Session layer
Transport layer
Network layer
Data link layer
Physical layer
In this book, we will mainly focus on the top 3 layers with special emphasis on algorithms
for extraction of information useful for relaying purposes.
2
Digital Signal Processing Basics and
Architecture of Numerical Relay
It can be seen from Figure 2.4 that for a 3-bit ADC with an input voltage range of 0 to V
volts, all voltages in the range 0 to (V/8) are assigned the same code 000. All voltages in the
range from (V/8) to (2V/8) are assigned the same code 001, and so on. Thus, the resolution
of the ADC is (V/8). As the input voltage changes from 0 to slightly more than (V/8) volts,
Digital Signal Processing Basics and Architecture of Numerical Relay9
the code changes from 000 to 001. Thus, the LSB corresponds to a change of (V/8) volts. The
quantum of (V/8) volts is therefore called 1 LSB, and the resolution of the ADC is said to
be 1 LSB. Alternatively, ADC resolution is expressed simply as number of bits in the output.
Let us indicate the total excursion in input voltage that the ADC can handle by Vspan where
Vspan can be written as:
Vspan = Vin, maxmim Vin, minimum
Magnitude of 1 LSB for an N-bit bipolar ADC with an input span of Vspan is given by:
Vspan (Vm - ( -Vm )) 2Vm
1 LSBbipolar = N
= N
=
2 2 2N
Note that for a bipolar ADC, Vin , minimum = Vm while for a unipolar ADC, Vin, minimum = zero.
Let us assume Vin, maximum for both ADCs as +Vm. Then, the span of a bipolar ADC will
be 2Vm while that of unipolar ADC will be Vm.
If we apply a ramp to the input of the ADC, the output of the ADC will be a quantized
digital output. If we feed this to a DAC, it will generate a staircase type waveform as shown
in Figure 2.5. We define the quantization error of the ADC as:
Quantization error = output of ADC (as seen through DAC) analog input
Note that even an ideal ADC having no other imperfections cannot escape from quantization
error. We can only make the quantization error smaller and smaller by increasing the number
of bits in the ADC, but we can never completely eliminate it.
We can plot the error signal waveform as shown in Figure 2.5. The error in the output
of the ADC can be considered as a noise signal attributable to the quantization phenomenon.
Hence, it is customary to quantify the error due to quantization in terms of signal to noise ratio
defined as ratio of signal voltage to noise voltage usually expressed in decibels.
Figure 2.5 shows that the actual transfer curve of the ADC, always hovers on the lower side
of the ideal transfer curve. Statistically it will be more appropriate if the actual transfer curve
symmetrically encompasses the ideal transfer curve. This can be easily achieved by adding a
LSB offset to all the analog levels at which the code makes transition as shown in Figure
2.6. This results in the error curve also being symmetric about the x-axis.
Figure 2.6 Transfer characteristics of a 3-bit ADC with half LSB offset and quantization error.
We are now in a position to express the quantization error in terms of signal to noise ratio.
Let us assume that a sinusoidal wave with a peak value of Vm is applied to the ADC. The RMS
value of the signal is thus Vm / 2 . Note that the signal excursion is from Vm to +Vm, which
is a change of 2Vm requiring an ADC with a span of 2Vm. Thus, the ADC is representing a
voltage range of 2Vm using N-bits.
The RMS value of the saw-tooth wave, generated because of quantization phenomenon,
which is considered as noise can be written as:
Digital Signal Processing Basics and Architecture of Numerical Relay11
1 T/2 2
Vnoise = ve dt
T T/2
where ve, the error signal, having the saw-tooth waveform, is given by:
Ve = (VHLSB/(T/2))t, where T is the time period of the saw-tooth error waveform and VHLSB
is the magnitude of voltage corresponding to LSB given by:
1 1 2Vm Vm
VHLSB = LSB = = N
2 2 2N 2
As noted earlier, the value of 2Vm in the numerator of above equation is due to the fact that
a sine wave with a peak of |Vm| needs a bipolar ADC which can accommodate input voltage
excursion from Vm to +Vm which is actually an excursion of 2Vm.
Hence, we can write:
2
1 T/2 -VHSLB
Vnoise =
T T/2 (T/2)
t dt
VHSLB 1 T/2 2
Vnoise =
T/2
(t ) dt
T T/2
T/2
V 1 t3 V
Vnoise = HSLB = 2 HLSB
T/2 T 3 -T/2 12
Hence,
Vm Vm
Vnoise, RMS = 2 N
= ( N -1)
2 12 2 12
SNR = 2 N 1 6
Let us express it in decibels:
Vsignal,RMS
SNR dB = 20 log10
Vnoise,RMS
12Digital Power System Protection
SNR dB = 20 log10 (2 N -1 6 )
The reason for calling it a flash converter is that the conversion from analog to digital
value in this topology is very fast since it involves only the switching time of the comparators
which is very small. In fact all other methods of analog to digital conversion are slower than
the flash method. Since all the comparators work in parallel, it is also known as parallel
comparator type ADC.
14Digital Power System Protection
It can be seen from Figure 2.7 that the comparators have thresholds of V/22, 2V/22, and
3V/22, i.e., V/4, 2V/4 and 3V/4. Further, the output digital code generated by the quantizer is
not a natural binary code. Thus, the ADC needs a code converter or the so-called encoder.
The encoder is a simple combinatorial circuit. The encoder can be easily designed as shown
in Figure 2.7 using K-map method.
2.4Sampling
An ADC cannot handle a continuously changing analog signal. We must first sample the analog
signal and then feed one sampled analog value to the ADC at a time. Figure 2.8(a) shows the
sampling of an analog signal at 16 samples per cycle.
In the above discussion we had assumed an ideal ADC. However, a real-life ADC has
many imperfections. One such important imperfection is that all ADCs require a finite amount
of time to convert the input analog voltage into its digital value.
this will solve the problem of finite conversion time of ADC. To illustrate how things get
more complicated, consider Figure 2.8(b) where a unipolar ADC is assumed and three analog
voltage waveforms with different frequencies but same peak value are shown. Using the rule
that signal should not change by more than LSB during conversion we see that for a given
number of ADC bits and peak value of input signal, the frequency of the signal that can be
handled without error goes down as we employ a better ADC. A better ADC is the one with
a smaller conversion time and greater number of bits. We must therefore hold the input to the
ADC constant for a time equal to conversion time of ADC.
Figure 2.8(b) shows three voltages of different frequencies LF, MF and HF having the
same peak value. It can be seen the LF < MF < HF. In a time duration equal to conversion
time Tc, each voltage changes by a different quantum. The quantum being highest for highest
frequency signal since all signals are assumed to have same amplitude. The figure also shows
the voltage corresponding to LSB. For error free operation of ADC, the input voltage should
not change by more than this value of LSB during conversion time Tc. It can be seen that
the condition that input voltage should not change by more than LSB is satisfied only by
the lowest frequency signal LF.
Now if we employ a better ADC, i.e., an ADC with higher number of bits N and smaller
conversion time Tc, the maximum frequency that can be handled goes down still further than
LF. This can be seen from Figure 2.8(c).
Thus, we cannot improve upon the maximum frequency that can be handled by using ADCs
with faster conversion time and more number of bits. This can be easily appreciated by finding
out the maximum frequency that can be handled by an ADC with excellent specifications.
Consider a 16-bit ADC with a conversion time of 10ms. Let us find the maximum frequency
that can be handled by an ADC with so good specifications.
16Digital Power System Protection
(dv/dt)ADC,max i.e., the maximum rate of change of voltage that the ADC can accommodate is
given by (dv/dt)ADC,max = (Dv/Dt)ADC,maxwhere and Dv = (Vm/2N) and Dt = Tc
dv V
= Nm
dt ADC,max 2 Tc
Equating (dv/dt)signal,max and (dv/dt)ADC,max we get:
Vm
2p f maxVm =
2 N Tc
Hence,
1
f max = ( N +1)
2 p Tc
Digital Signal Processing Basics and Architecture of Numerical Relay17
This is too low a frequency to be of any use in power system protection work! And the
frequency will become still lower if we use a better ADC! Thus we have to think of a device
in between the analog input signal and the ADC which will be able to cope up with the
high (dv/dt) rates involved. Fortunately, a simple solution to this problem exists. If we use a
sampler followed by a hold circuit, i.e., a sample and hold circuit (S/H) then we can go
for sampling much higher frequencies. This becomes possible because now it is the aperture
time of the hold circuit, Ta, that has to meet the (dv/dt) requirement rather than the ADC.
Fortunately, the aperture times of S/H circuits are many orders of magnitude smaller than the
ADC conversion times. Thus, fmax with sample and hold can be obtained by substituting Ta in
place of Tc in the expression for fmax derived earlier.
1
f max = ( N +1)
without hold circuit
2 p Tc
1
f max = ( N +1)
with hold circuit with aperture time Ta
2 p Ta
Let us rework the maximum frequency by assuming a realistic aperture time of 250 ps.
1
f max = (16 +1)
with hold circuit
2 p (250)(10 -12 )
Thus, we get improvement in maximum frequency from 0.2428 Hz to 9.714 kHz because
of simple expedient of using a sample and hold circuit ahead of the ADC. We make suitable
modification in the DSP signal chain as shown in Figure 2.9(a). Figure 2.9(b) shows an analog
signal , the sampling impulses, the sampled signal and sampled-and-held versions of an analog
signal. Note that and fsignal = (1/Tsignal) and fsampling = (1/Tsampling)
Figure 2.9(a) Modified DSP signal chain with sample and hold circuit included.
18Digital Power System Protection
The sampling theorem is pictorially shown in Figure 2.10. It can be seen from Figure 2.10
that as dictated by the sampling theorem, there is a band of frequencies from zero to 2fsignal,max
which are forbidden for the sampling frequency.
This naturally leads to the next question: What happens if we violate the sampling theorem
and sample at a forbidden frequency, i.e., at a frequency less than that stipulated by Shannon?
We encounter the phenomenon of aliasing in such an eventuality, wherein a higher
frequency signal appears as a lower frequency signal. When aliasing takes place, we cannot
recover correct information from the signal. Thus, aliasing is to be avoided at all costs.
Figure 2.11(a) shows the effect of violating the sampling theorem. It can be seen that
the sampling is at rate less than twice the signal frequency since there are no samples during
alternate halves of the input sine wave. The wave that will be recreated from the samples will
be a low frequency aliased wave different from the input.
As expected, there is aliasing at sampling frequencies of 49 Hz, 51 Hz and even at 100 Hz.
Sampling at 49 Hz and 51 Hz give an aliased signal of 1 Hz with phase angles of zero and
180 respectively. At 100 Hz, since the sampling is synchronised with the zero crossing of
the waveform, we get a DC with zero magnitude as the aliased wave. However, at sampling
frequency of 500 Hz we are able to correctly recreate a 50 Hz wave, since there are 50 complete
oscillations in 1 second.
In DSP based systems, sampling plays very important role as what we perceive
through the sampling process strongly depends upon the relation between the signal
bandwidth and the sampling frequency. Figure 2.11(c) shows perceived frequency Vs
the actual input frequency for a fixed sampling frequency. It can be seen that only in the
region where Shannons sampling theorem is not violated, do we correctly perceive the
input frequency. Outside this range we actually either have no perception or are working
in the illusory realm because of aliasing.
Digital Signal Processing Basics and Architecture of Numerical Relay21
Fig. 2.11(c) Perceived frequency Vs input frequency for a fixed sampling frequency.
Figure 2.12(b) shows the characteristics of an ideal anti-aliasing filter that will entirely
filter out all frequency components above fsampling/2 . It can be seen that such an ideal anti-
aliasing filter will have to have its cut-off frequency at fsampling/2 and exhibit an infinite roll-off.
However, all practical filters have finite roll-off. Hence in case of a practical anti-aliasing filter
we will have to push-back the cut-off frequency so that the magnitude of the signal falls below
LSB in the stop band as shown in Figure 2.12(c).
1
A( factual ) =
2n
f
1 + actual
fcut-off
If we represent the ratio [(fsampling/2)/fcut-off] by a, then we can write the above equation as:
Vm 1 2Vm
1 + (a ) 2n 2 2N
Solving the above equation gives order of the low-pass filter as:
loge (4 N - 1)
n
2 loge (a )
loge (412 - 1)
n (3.61 rounded to 4)
2 loge (10)
Hence, 4th order low-pass anti-aliasing Butterworth filter should be used for the given system.
24Digital Power System Protection
Review Questions
1.
What do you mean by resolution of an ADC?
2.
What do you mean by 1 LSB, LSB w.r.t. ADCs?
3.
In a 10-bit ADC with input range of 10 V to 10 V, what is the value of 1 LSB and
LSB?
4.
In an ADC with input range of +/5 V, a resolution (1 LSB) of 10 mV is required.
What is the minimum number of bits it must have?
5.
Find the signal to noise ratio (in dB) due to quantization in an DSP system which has
ADC of 16-bits and sampling is done at the Nyquist rate. What will be the error if
sampling is done at 10 times the Nyquist rate?
6.
Find the SNR (in dB) due to quantization in a DSP system with the following
specifications:
Number of bits in the ADC = 12, signal bandwidth = 1000 Hz and sampling frequency
= 5000 Hz.
7.
How does one decide the minimum sampling frequency in a DSP system?
8.
A power system signal is given by:
v(t ) = 100 sin(2p 50t + 45) + 50 sin(2p 100t + 25) + sin(2p 150t + 15)
Suggest minimum sampling frequency, in hertz to avoid aliasing errors.
9. What, if any, are the disadvantages of sampling at very high frequencies , much above
the Nyquist sampling rate?
10.
Give at least two examples where aliasing is encountered in daily life.
(Hint: Movies, Table fan)
11.
The frequency spectra of signal is as shown in Figure Q.11. If the signal is ideally
sampled at intervals of 1 ms, then state with justification whether this will lead to
aliasing. Sketch the frequency spectra of the sampled signal.
Figure Q.11
28Digital Power System Protection
12.
Find the maximum frequency that can be sampled without using hold circuit for a DSP
system with the following specifications:
Conversion time of ADC = 5 s and number of bits in the ADC = 16
13.
Find the order of low-pass anti-aliasing filter of Butterworth type for a DSP system
with the following specifications:-
ADC resolution = 12-bits , sampling frequency = 6000 Hz and cut-off frequency = 300 Hz
14.
Provide a frequency domain explanation of the aliasing phenomenon.
3
Algorithms Based on Undistorted
Single Frequency Sine Wave
3.1.2 Derivation of the Sample and Derivative (Mann & Morrison) Algorithm
In this algorithm, we model the signal in a very simplistic manner. It is assumed that the signal
can be represented as a pure sinusoid whose amplitude as well as frequency is constant during
the period under consideration as shown in Figure 3.1.
Thus, voltage samples collected by the relay are assumed to belong to a signals of the
type:
v(tk ) = Vm sin(w tk + q v ) (3.1)
where we assume Vm, w and qv to be constant. Further we assume that the frequency w is
known. Thus in the above equations we have two unknowns, namely the peak value of the
signal and its phase angle. Hence, we need one more equation to solve for the unknowns. Let
us generate another equation by taking derivative of the first equation giving:
v (tk ) = wVm cos(w tk + q v ) (3.2)
from which we get:
v (tk )
= Vm cos(w tk + q v ) (3.3)
w
Combining Eqs. (3.1) and (3.3), we get:
2
v (tk ) (3.4)
Vm2 = (v(tk ))2 +
w
Hence, we can find the amplitude of voltage signal as:
v (tk )
2
Vm = (v(tk ))2 + (3.5)
w
Similarly, we can write for the current signal as:
i(tk ) = I m sin(w tk + q i ) (3.6)
Algorithms Based on Undistorted Single Frequency Sine Wave31
i (tk )
2
I m = (i(tk )) +
2
w
v(tk ) (3.9)
tan(w tk + q v ) =
v (tk ) /w
v(tk )
w tk + q v = tan -1 (3.10)
v (tk )/w
v (t k )
q v = tan -1 ' - w tk (3.11)
v ( t k )
w
Similarly, we can get:
i(tk )
q i = tan -1 - w tk (3.12)
i (tk )
w
The derivatives can be computed numerically using the central difference formula as:
(vk +1 - vk -1 )
vk = (3.13)
2 Dt
and (ik +1 - ik -1 ) (3.14)
ik =
2t
Thus, computation of one derivative requires three samples. Now, we cannot base
the trip decision on a single estimate (which involves only three samples), we have to
continuously keep estimating the peak and the phase angle and hence the phasor, in order to
be sure about what way the phasor is really going. We should, thus, base our trip decision
on a sufficiently large number of estimates rather than a single estimate. This will make
the decision more robust and provide some measure of noise immunity as discussed in the
following section.
32Digital Power System Protection
Let us summarise the results of the Mann and Morrison algorithm as shown in Table 3.1.
v (tk )
2
( v (t k )) 2
+ w
Vm =
(vk +1 - vk -1 )
where vk =
2 Dt
v(tk )
qv = tan -1 - w tk radians
v (tk ) / w
i (tk )
2
(i(tk )) +
2
w
Im =
(ik +1 - ik -1 )
where ik =
2 Dt
i(tk )
qi = tan -1 - w tk radians
i (tk ) / w
The threshold value to be used for the counter will be based on severity of noise at a
particular location. The general logic of such an implementation is shown in Figure 3.2.
Digital Power System Protection
reactance as seen from the relay location.
Figure 3.4 Spread sheet for simulation of Mann and Morrison algorithm.
samples of current, slides through these samples and computes the apparent resistance and
current which are generated for a cycle. A window of three samples of voltage and three
Algorithms Based on Undistorted Single Frequency Sine Wave35
Vm_est(k)=sqrt(v(k)^2 + (v_dash(k)/w)^2);
Im_est(k)=sqrt(i(k)^2 + (i_dash(k)/w)^2);
Phi_V_est_rad(k) = (atan2(((v(k))/((v_dash(k))/w)),1))-(w*k*dT);
Phi_V_est_deg(k)= (Phi_V_est_rad(k))*(180/pi);
Phi_I_est_rad(k) = (atan2(((i(k))/((i_dash(k))/w)),1))-(w*k*dT);
Phi_I_est_deg(k)= (Phi_I_est_rad(k))*(180/pi);
R_est(k) = Zmag_est(k)*cos(Zang_est_rad(k));
X_est(k) = Zmag_est(k)*sin(Zang_est_rad(k));
L_est(k) = X_est(k)/w;
end
fprintf(==============================================================
=======================================================\n)
fprintf(Samp v i v_dash Vm_est i_dash Im_est Phi_V_est
Phi_I_est R_est L_est\n)
fprintf(True Value = %5.2f %5.4f %5.2f %5.4f
%5.2f% 5.4f\n,Vm,Im,Phi_v,Phi_v-Phi_i,R,L )
fprintf(==============================================================
=======================================================\n)
for n = 1:no_of_samples-1
if(n==1)
fprintf(%2i | %f %f\n,n,v(n),i(n))
continue
end
fprintf(%2i | %f %f| %f| %f| %f| %f| %f| %f| %f| %f\n,n,v(n),i(n),v_
dash(n),Vm_est(n),i_dash(n),Im_est(n),Phi_V_est_deg(n),Phi_I_est_
deg(n),R_est(k),L_est(k) )
end
% Vm_est
% Im_est
%
% Phi_V_est_rad
% Phi_V_est_deg
% Phi_I_est_deg
Algorithms Based on Undistorted Single Frequency Sine Wave37
Let the three signal samples, starting from an arbitrary index k be:
vk = A sin(2p f (k Dt ) + q ) (3.17)
vk +1 = A sin(2p f (k Dt + Dt ) + q ) (3.18)
vk + 2 = A sin(2p f (k Dt + 2 Dt ) + q ) (3.19)
X - Y = 2 p f (k Dt + Dt ) - 2 p f Dt = 2 p f (k Dt ) + q (3.22)
Hence vk + vk+2 can be written as:
vk + vk + 2 = A [sin( X - Y ) + sin( X + Y )] (3.23)
We can use the following trigonometric identity:
[sin( X - Y ) + sin( X + Y )] = 2 sin( X ) cos(Y ) (3.24)
and we can write:
vk + vk + 2 = A [2 sin( X ) cos(Y )] (3.25)
vk + vk + 2 = A [2 sin(2 p f (k Dt + Dt ) + q ) cos(2 p f Dt )] (3.26)
However,
A sin(2 p f (k Dt + Dt ) + q ) = vk +1 (3.27)
Hence, we can write
vk + vk + 2 = [2 vk +1 cos(2 p f Dt )] (3.28)
vk + vk + 2
cos(2 p f Dt ) = (3.29)
2 vk +1
vk + vk + 2 (3.30)
2 p f Dt = cos -1
2 vk +1
1 v + vk + 2
f= cos -1 k (3.31)
2 p Dt 2 vk +1
Thus, since Dt, the sampling time period which is 1/fsampling is known, the unknown
frequency can be estimated with the help of three consecutive samples.
1 - cos(2(2p f (k Dt + Dt ) + q )) (3.34)
(sin(2p f (k Dt + Dt ) + q )2 =
2
further, if we write
(2p f (k Dt + 2 Dt ) + q ) = X (3.35)
40Digital Power System Protection
= A2 1 - cos 2p f (2 Dt ) (3.44)
2
Using the trigonometric identity cos (2 X ) = 2(cos ( X )) - 1
2
vk
= sin(2p f (k Dt ) + q ) (3.52)
A
v
2p f (k D t ) + q = sin -1 k (3.53)
A
v
q = sin -1 k - 2p f (k Dt ) (3.54)
A
In the above expression vk, A, k and Dt all are known. Hence, phase angle q can be estimated.
Signal v = Vm sin(w t + q )
vk = A sin(2p f (k Dt ) + q )
Samples vk +1 = A sin(2p f (k Dt + Dt ) + q )
vk + 2 = A sin(2p f (k Dt + 2 Dt ) + q )
1 v + vk + 2
Frequency f f (k ) = cos -1 k
2 p Dt 2 vk +1
v
Phase q q (k ) = sin -1 k - 2p f (k Dt )
A
42
Digital Power System Protection
for n=2:numel(x)-1
fprintf(%d %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n,n,w_est(n),f_
est(n),A_est(n), theta_est_deg(n),err_A(n),err_w(n),err_theta(n) )
end
fprintf(\n\n)
for p=1:3
fprintf(Problem No. %d \n,p)
fprintf(The samples of a single frequency signal are as shown below.\n)
fprintf(Make at least three estimates of the amplitude,phase and frequency\n)
fprintf( 0 1 2 3 4 5 6 \n)
for n=p+1:p+7
fprintf( %f ,x(n) )
end
fprintf(\n)
fprintf(===================================================== \n)
fprintf(\n\n)
end
===============================================
The true values are:-
A=104.000000 fsig=50.195001 theta=39.000000
===============================================
======The estimated values and errors are ======================================================
n Omega freq Ampltd PhAng Err_A% Err_w% Err_Ph_Ang%
=======================================================================================
2 315.384 50.195 104.000 39.000 0.000 0.000 0.000
3 315.384 50.195 104.000 39.000 0.000 0.000 0.000
4 315.384 50.195 104.000 39.000 0.000 0.000 0.000
5 315.384 50.195 104.000 39.000 0.000 0.000 0.000
6 315.384 50.195 104.000 39.000 0.000 0.000 0.000
7 315.384 50.195 104.000 39.000 0.000 0.000 0.000
8 315.384 50.195 104.000 39.000 0.000 0.000 0.000
where we assume Vm, w and qv to be constant. Further we assume that the frequency w is
known. Thus in the above equations, we have two unknowns, namely the peak value of the
signal and its phase angle. Hence, we need one more equation to solve for the unknowns. Let
us generate another equation by taking derivative of the first equation giving:
v (t ) 2 v (t ) 2
Vm = (3.58)
+
k k
2
w w
i (t ) 2 i (t ) 2
I m = k + 2k (3.59)
w w
Further, on the lines of Mann and Morrison algorithm, it can be easily shown that:
v (tk )/w 2
q v = tan -1 - - w tk (3.60)
v (tk )/w
i '' (t )/w 2
q i = tan -1 - ' k - w tk
(3.61)
i (tk )/w
v (t ) 2 v (t ) 2
Vm = k
+
k
2
w w
Vm (vk +1 - vk -1 )
where, vk =
2 Dt
v(k + 1) - 2 v(k ) + v(k - 1)
v (k ) @
Dt2
v (tk )/w 2
qv q v = tan -1 - - w tk radians
v (tk )/w
i ' (t ) 2 i '' (t ) 2
I m = k + 2k
w w
Im (ik +1 - ik -1 )
where, ik =
2 Dt
i(k + 1) - 2i(k ) + i(k - 1)
i (k ) @
Dt 2
i (tk )/w 2
qi q i = tan -1 - - w tk radians
i (tk ) /w
Algorithm
Figure 3.8 Spreadsheet for implementation of first and second derivative algorithm.
Algorithms Based on Undistorted Single Frequency Sine Wave47
f_samp = 2*OverSamplingFactor*f ;
X = R*XbyR_ratio ;
L = X/(2*pi*f) ;
Z = sqrt( R*R + X* X) ;
line_ang_deg = (180/pi)*atan2(X,R) ;
w = 2* pi* f ;
T = 1/f ;
dT = 1/f_samp ;
Vm = 100 ;
Im = Vm / Z ;
Qv_deg = 30 ;
Qv = Qv_deg*(pi/180) ;
Qi = Qv - atan2(X,R) ;
Qvdeg = Qv * (180/pi) ;
Qi_deg = Qi*(180/pi) ;
k=0;
sss=50;
%---------------------- Generation of Samples ----------------------------
for tim = 0:dT: T
k =k+1 ;
t(k) = tim ;
v(k) = Vm * sin( 2* pi* f * t(k) + Qv) ;
i(k) = Im * sin( 2* pi* f * t(k) + Qi) ;
end
% ---------------------------------------------------------------------
samples =k ;
for n = 2:samples-1
Algorithms Based on Undistorted Single Frequency Sine Wave49
%---------- Estimation of voltage phasor -------------------------------
v_d(n) = ( v(n+1)-v(n-1)) / (2*dT) ;% First derivative
v_dd(n) = (v(n+1)-2*v(n)+ v(n-1))/(dT*dT) ;% 2nd Derivative
Vm_est(n) = sqrt( (v_d(n)/w)^2 + ( v_dd(n)/(w*w) )^2 ) ;
num = (v_dd(n))/((w*w)) ;
den = v_d(n)/w ;
y = num/den;
theta_v_est(n) = (atan2(- y,1))-((w*(n-1)*dT));
theta_v_deg_est(n)=theta_v_est(n)*(180/pi);
err_Vpc(n)=(((Vm_est(n))-(Vm))/(Vm))*100;
err_theta_v_deg(n)=(((theta_v_deg_est(n))-(Qv_deg))/(Qv_deg))*100;
%----------------------------------------------------------------------
%---------- Estimation of Current phasor -------------------------------
i_d(n) = ( i(n+1)-i(n-1)) / (2*dT) ;% First derivative
i_dd(n) = (i(n+1)-2*i(n)+ i(n-1))/(dT*dT) ;% 2nd Derivative
Im_est(n) = sqrt( (i_d(n)/w)^2 + ( i_dd(n)/(w*w) )^2 ) ;
numi = (i_dd(n))/((w*w)) ;
deni = i_d(n)/w ;
y = numi/deni;
theta_i_est(n) = (atan2(- y,1))-((w*(n-1)*dT));
theta_i_deg_est(n)=theta_i_est(n)*(180/pi);
err_ipc(n)=(((Im_est(n))-(Im))/(Im))*100;
err_theta_i_deg(n)=(((theta_i_deg_est(n))-(Qi_deg))/(Qi_deg))*100;
%----------------------------------------------------------------------
R_est(n)=(Vm_est(n)/Im_est(n))*cos(theta_v_est(n)-theta_i_est(n));
X_est(n)=(Vm_est(n)/Im_est(n))*sin(theta_v_est(n)-theta_i_est(n));
err_R_pc(n)=(((R_est(n))-(R))/(R))*100;
err_X_pc(n)=(((X_est(n))-(X))/(X))*100;
end
% ------------------------------------------------------------------
fprintf(------------------First and Second Derivative Algo.-------------
----------\n)
fprintf(--------------------------------------------------------------
------------\n)
fprintf(Vm Vm_est(n) Error in Vm Ph Ang Estimate_of_ph Error in Ph\n)
fprintf(--------------------------------------------------------------
------------\n)
for n=2:numel(Vm_est)
fprintf(%f %f %f %f %f %f\n,Vm,Vm_est(n),err_Vpc(n),Qvdeg, theta_v_
deg_est(n),err_theta_v_deg(n))
end
%----------------------------------------------------------------------
fprintf(--------------------------------------------------------------
------------\n)
fprintf(Im Im_est(n) Error in Im Ph Ang Estimate_of_ph Error in Ph\n)
fprintf(--------------------------------------------------------------
------------\n)
for n=2:numel(Im_est)
50Digital Power System Protection
(v(k + 1))2 - 2 v(k )v(k + 1) cos wDt + (v(k ))2 cos2 (wDt )
(v(k ))2 + = Vm2 (3.71)
sin (wDt )
2
(v(k ))2 sin 2 (wDt ) + (v(k ))2 cos2 (wDt ) + (v(k + 1))2 - 2 v(k )v(k + 1) cos(wDt )
Vm2 = (3.72)
sin 2 (wDt )
Simplifying further, we get:
Note that all the terms on the RHS of Eq. (3.74) are known.
Similarly,
%==============================================
mag_z_set = 30 ; % Relay Reach Setting
theta_n = pi/2 ; % Relay MTA setting
%=====================================================
%============= Numerical Problem Setting ===================
f = 50 ; % System frequency
OverSamplingFactor = 20 ; % fsamp = 2 * OSF * fsig
km = 100 ; % kilometers
R = 0.037 ; % ohms per km for 345kV line
XbyR_ratio = 9.9189 ; % for 345 kV line
%=====================================================
f_samp = 2*OverSamplingFactor*f ;
X = R*XbyR_ratio ;
L = X/(2*pi*f) ;
Z = sqrt( R*R + X* X) ;
line_ang_deg = (180/pi)*atan2(X,R) ;
w = 2* pi* f ;
T = 1/f ;
dT = 1/f_samp ;
Vm = 100 ;
Im = Vm / Z ;
Qv_deg =0 ;% Theta_V = 0
Qv = Qv_deg*(pi/180) ;
Qi = Qv - atan2(X,R) ;
Qi = 1 * Qi ; %In the algorithm ph. angle of
current is taken positive
Qvdeg = Qv * (180/pi) ;
Qi_deg = Qi*(180/pi) ;
k =0;
sss =50;
%========== Generation of Samples ================
for tim = 0:dT: T
k =k+1 ;
t(k) = tim ;
v(k) = Vm * sin( 2* pi* f * t(k) ) ;
end
Algorithms Based on Undistorted Single Frequency Sine Wave55
% ===============================================================
samples =k ;
for n = 1:samples-1
%---------- Estimation of voltage phasor -------------------------------
Vm_est(n)=sqrt((((v(n))^2)+((v(n+1))^2)-(2*v(n)*v(n+1)*cos(w*dT)))/
((sin(w*dT))^2));
Im_est(n)=sqrt((((i(n))^2)+((i(n+1))^2)-(2*i(n)*i(n+1)*cos(w*dT)))/
((sin(w*dT))^2));
%---------------------------------------------------------------
%---------- Estimation of Current phasor -------------------------------
num=(i(n)*v(n))+(i(n+1)*v(n+1))-(((i(n)*v(n+1))+(i(n+1)*v(n)))*cos(w*
dT));
den=(Im_est(n)*Vm_est(n)*((sin(w*dT))^2));
Ph_i_est(n)=acos(num/den);
Ph_i_est_deg(n)=(Ph_i_est(n))*(180/pi);
%---------------------------------------------------------------
%-------- Error in estimation of voltage and current phasors
err_Vm(n)=(((Vm_est(n))-(Vm))/(Vm))*100;
err_Im(n)=(((Im_est(n))-(Im))/(Im))*100;
err_Ph_i_deg(n)=(((Ph_i_est_deg(n))-(Qi_deg))/(Qi_deg))*100;
R_est(n)=(Vm_est(n)/Im_est(n))*cos(Ph_i_est(n));
X_est(n)=(Vm_est(n)/Im_est(n))*sin(Ph_i_est(n));
err_R_pc(n)=(((R_est(n))-(R))/(R))*100;
err_X_pc(n)=(((X_est(n))-(X))/(X))*100;
end
% ---------------------------------------------------------------
fprintf(------------------Two Sample Technique -----------------------
\n)
fprintf(--------------------------------------------------------------
------------\n)
fprintf(Vm Vm_est(n) Error in Vm \n)
fprintf(--------------------------------------------------------------
------------\n)
for n=1:numel(Vm_est)
fprintf(%f %f %f \n,Vm,Vm_est(n),err_Vm(n))
end
%---------------------------------------------------------------
fprintf(--------------------------------------------------------------
------------\n)
fprintf(Im Im_est(n) Error in Im Ph Ang Estimate_of_ph Error in Ph\n)
fprintf(--------------------------------------------------------------
------------\n)
for n=2:numel(Im_est)
fprintf(%f %f %f %f %f %f\n,Im,Im_est(n),err_Im(n),Qi_deg,Ph_i_est_
deg(n),err_Ph_i_deg(n))
end
fprintf(--------------------------------------------------------------
------------\n)
56Digital Power System Protection
Vm Vm_est(n) Error in Vm
Review Questions
1. The voltage and current samples at a relay location are shown below. The system
frequency is known to be 50 Hz and the sampling frequency is 10 kHz. Using Mann
and Morrison algorithm find the apparent resistance and reactance as seen from the
relay location.
Sample set no. 1
v i
(i) 55.000000 8.887405
(ii) 57.965138 8.682345
(iii) 60.873070 8.468716
(iv) 63.720929 8.246729
Sample set no. 2
(i) 66.505903 8.016604
(ii) 69.225243 7.778567
(iii) 71.876266 7.532854
(iv) 74.456357 7.279707
Sample set no. 3
(i) 76.962967 7.019376
(ii) 79.393625 6.752117
(iii) 81.745931 6.478195
(iv) 84.017563 6.197880
Sample set no. 4
(i) 86.206280 5.911448
(ii) 88.309922 5.619183
(iii) 90.326413 5.321371
(iv) 92.253762 5.018309
Sample set no. 5
(i) 94.090069 4.710293
(ii) 95.833519 4.397630
(iii) 97.482394 4.080626
(iv) 99.035065 3.759595
58Digital Power System Protection
2.
What are the assumptions in the Mann and Morrison algorithm?
3.
Derive the Mann and Morrison algorithm.
4.
What are the drawbacks of using numerical derivative?
5.
Assuming that amplitude of signal is known, can you find its frequency by using an
algorithm similar to Mann and Morrison algorithm.
An ac signal v = Vm sin(2 pft + f) has three unknown Vm, f and f. Thus, we need
6.
three equations to solve for the three unknown. Using this idea can you derive an
algorithm?
7.
Three-sample algorithm assumes a sine waveform. Derive three-sample algorithm when
cosine waveform is taken as reference.
8.
What could be the applications of three-sample algorithm in general and in the field
of power system protective relaying?
9.
Discuss the importance of measurement of power system frequency.
10.
What are the assumptions underlying three-sample algorithm?
11.
What is the minimum number of samples required for making one estimate each of
amplitude, phase and frequency, assuming the signal to be undistorted sine wave?
12.
What happens to the estimation using three-sample algorithm if one of the samples
happens to be zero or very low in magnitude?
13.
Compare and contrast three-sample algorithm with Mann and Morrison algorithm.
14.
Can three-sample algorithm be applied under non-stationary conditions?
15. Numerical problem on three-sample technique:
(i) The samples of a single frequency signal are as shown below:
Make at least three estimates of the amplitude,phase and frequency
0 1 2 3 4 5 6
56.642460 26.917181 5.442939 37.270267 65.449321 87.221739 100.456286
=================================================================
The sampling frequency is 1003.90002 Hz.
(ii) The samples of a single frequency signal are as shown below:
Make at least three estimates of the amplitude,phase and frequency
0 1 2 3 4 5 6
5.442939 -37.270267 65.449321 87.221739 100.456286 103.857472 97.092364
=================================================================
The sampling frequency is 1003.90002 Hz.
Algorithms Based on Undistorted Single Frequency Sine Wave59
16.
Derive the first derivative and second derivative algorithm.
17.
What are the assumptions underlying the first and second derivative algorithm?
18.
What could be the motivation for such algorithms involving higher derivatives?
19.
Derive the two-sample algorithm. What are the advantages of using this algorithm?
20.
What are the assumptions underlying the two-sample algorithm?
21.
Develop a spreadsheet for simulating the first and second derivative and two-sample
algorithm.
22.
Write MATLAB programs for simulating the first and second derivative algorithm and
two-sample algorithm.
4
Algorithms Based on Solution of
Differential Equation
It may be noted that we have ignored the shunt capacitance of the line and it is assumed
that fault resistance is zero. Note that v t) and i t) are the voltage and current available at the
relay location through CT and PT which are not shown in Figure 4.1. Hence, both the CT and
the PT ratio may be considered to be equal to 1. Figure 4.2 shows the sampling instants and
sampled values of voltage and current.
If we apply KVL around the loop formed because of fault, at instant tk , we can write:-
di(tk )
v(tk ) = R i(tk ) + L
dt
It can be easily seen that in the above differential equation, the knowns are v(tk) and i(tk)
, and the unknowns are R and L, i.e., the resistance and the inductance of the transmission
line from relay location upto the fault point. The rate of change of current term, even though
60
Algorithms Based on Solution of Differential Equation61
not measured directly, can be computed from the current measurements. In order to solve for
the two unknowns, we need two equations. The second equation can be generated by applying
KVL at another instant of time say (k + 1)th instant as shown below:
di(tk )
v(tk ) = R i(tk ) + L (4.1)
dt
di(tk +1 )
v(tk +1 ) = R i(tk +1 ) + L (4.2)
dt
Let us represent di(t)/dt by i(t) and rewrite Eqs. (4.1) and (4.2) as:
v(k ) = R i(k ) + Li (k ) (4.3)
v(k + 1) = R i(k + 1) + Li (k + 1) (4.4)
Note that when the numerical values of the derivative of the current are inserted in the
above equations, they are no longer differential equations but become simultaneous algebraic
equations with constant coefficients.
We can write the above two equations in matrix form as:
v ( k ) i( k ) i ( k ) R
= (4.5)
v ( k + 1) i(k + 1) i (k + 1) L
Which gives us:
-1
R i(k ) i ' ( k ) v( k )
= (4.6)
L i(k + 1) i ' (k + 1) v(k + 1)
R 1 i (k + 1) -i (k ) v(k )
L = i(k )i (k + 1) - i (k )i(k + 1) -i(k + 1) i(k ) v(k + 1) (4.7)
62Digital Power System Protection
Finally we can express R and L as a function of various samples of voltage and current
and computed derivatives of current
i (k + 1)v(k ) - i (k )v(k + 1)
R = (4.8)
i(k )i (k + 1) - i (k )i(k + 1)
i(k )v(k + 1) - i(k + 1)v(k )
L= (4.9)
i(k )i (k + 1) - i (k )i(k + 1)
Note that the numerical derivatives can be calculated using central difference formula as:
v(k + 1) - v(k - 1)
v (k ) = (4.10)
2 Dt
i(k + 1) - i(k - 1) (4.11)
i ( k ) =
2 Dt
where Dt is the sampling time period given by 1/fsampling.
The block-diagram showing the input to the algorithm and the block required for
implementation of specific distance relay characteristics is shown in Figure 4.3.
Figure 4.3Block-diagram showing inputs to differential equation algorithm and block for
implementation of specific distance relay.
It can be seen from Figure 4.3 that after estimating the apparent R and L as seen from
the relay location, we will have to check whether the apparent impedance lies within the trip
region of the relay characteristics being implemented. Thus, in addition to the differential
equation algorithm, we will have to write another algorithm for implementing specific distance
relay characteristics.
of 100 km, 345 kV line. We will use the data from Table 4.1 which happens to be for a 60
Hz power system.
Table 4.1 Data of a 60 Hz power system
230 kV 345 kV 500 kV 765 kV 1100 kV
R /km 0.050 0.037 0.028 0.012 0.005
XL 0.488 0.367 0.326 0.329 0.292
XL/R 9.76 9.9189 11.6428 27.4166 58.4
1
tan (xL/R) deg 84.15 84.2430 85.09 87.91 89.019
S s/km 3.371 4.518 5.20 4.978 5.544
Zc 380 285 250 257 230
4.518 10 6
Hence, C= = 0.01198436721 micro-farad (mF)
(2)(p )(60)
Consider 100 km length of line, with total capacitance lumped at one end.
1 1
Zshunt = XC = = = j2210.486 W
jw C j ((2)(p )60(1.2 10 6 ))
For p model, of shunt capacitance is connected at each end. Hence, XC becomes double,
i.e., j 4420.97 W.
Figure 4.4 Entire shunt capacitance modelled at sending end with fault on.
64Digital Power System Protection
Table 4.2
Consider of the 100 km length of the 345 kV line modelled by:
Series lumped model Total capacitance lumped at sending p model. with capacitance
end lumped at each end
XL = (100) (0.367) = j 36.7 W Zseen = (3.7 + j36.7)||(j 2210.486) Zseen = (3.7 + j36.7)||(j 4420.97)
R = (100) (0.037) = 3.7 W (3.7 + j 36.7) ( - j 2210.486) (3.7 + j 36.7) ( - j 4420.97)
= =
Zseen = (3.7 + j 36.7) W (3.7 + j 36.7)( - j 2210.486) (3.7 + j 36.7)( - j 4420.97)
Thus, we see that the error introduced because of ignoring shunt capacitance is less than
1%, which is not at all serious. Hence, our representation of the transmission line as a series
R-L circuit is justified.
4.3 E
Equation Algorithm
f_samp = 2*OverSamplingFactor*f ;
X = R*XbyR_ratio ;
L = X/(2*pi*f) ;
Z = sqrt( R*R + X* X) ;
line_ang_deg = (180/pi)*atan2(X,R) ;
w = 2* pi* f ;
T = 1/f ;
dT = 1/f_samp ;
Vm = 400*sqrt(2) ;
Im = Vm / Z ;
Qv_deg = -15 ;
Qv = Qv_deg*(pi/180) ;
Qi = Qv - atan2(X,R) ;
Qvdeg = Qv * (180/pi) ;
Qideg = Qi*(180/pi) ;
k =0 ;
sss = 50 ;
for tim = 0:dT: T
k = k+1 ;
t(k) = tim ;
v(k) = Vm * sin( 2* pi* f * t(k) + Qv) ;
i(k) = Im * sin( 2* pi* f * t(k) + Qi) ;
end
samples = k ;
m = 0 ;
for n = 2:samples-1
m = m+1 ;
Algorithms Based on Solution of Differential Equation67
fprintf(================================================= \n)
disp(Numerical Problem)
fprintf( A Mho Relay with following settings is to be implemented \n) ;
fprintf( Z Setting = %f Angle Z setting = %f\n , mag_z_set , (180/pi)* theta_n)
fprintf( Use differential algorithm to make two estimates apparent R and L
\n)
fprintf( Given that the system freq is %f \n,f )
fprintf( and given that the sampling freq is %f \n,f_samp )
%fprintf(\n Sampling Interval =%f Sampling Freq = %f\n , dT , 1/dT)
%---------------------------------------------------------------------
--
fprintf(The Voltage and Current Samples are :- \n);
fprintf(Voltage Samples )
fprintf(\n==================================================\n);
nn =8 ;
for kk = 1 : nn
fprintf( %f | ,v(kk)) ;
end
fprintf(\n===================================================\n);
fprintf(\nCurrent Samples)
fprintf(\n===================================================\n);
for kk = 1 :nn
fprintf(%f | ,i(kk)) ;
end
fprintf(\n===================================================\n);
for k=2:num_calc
den(k) = (i(k)*i_dash(k+1) )- ( i(k+1)*i_dash(k) ) ;
RR(k) = ((i_dash(k+1)*v(k))-(i_dash(k)*v(k+1)))/(den(k)) ;
LL(k) = ((i(k)*v(k+1))-(i(k+1)*v(k)))/(den(k)) ;
err_R_pc(k) = (((RR(k))-(R))/(R))*100 ;
err_L_pc(k) = (((LL(k))-(L))/(L))*100 ;
end
R,RR,err_R_pc
L,LL,err_L_pc
RR = 00.03760.03760.03760.03760.03760.03760.0376
err_R_pc = 01.55811.55811.55811.55811.55811.55811.5581
L = 0.0012
LL = 00.0012 0.0012 0.00120.00120.00120.00120.0012
err_L_pc = 1.0e003 *
00.16450.16450.16450.16450.16450.16450.1645
4.5
Solution of Differential Equation Algorithm Using Numerical
Integration
We use the same model of the faulted transmission line, as before, as shown in Figure 4.7.
We have already seen that the differential equation relating the voltage and current at the relay
location v(t ) = R i(t ) + L di(t )/dt gets converted into an algebraic equation when we substitute
the numerical value of the derivative term. However, a derivative has the problem that it tends
to amplify noise because of its high-pass nature. This is very easy to see. Let the signal be
v = Vm,signal sin(wsignalt) + Vm,noise sin(wnoiset) (4.12)
Note that wnoise >> wsignal and Vm,noise << Vm,signal
Algorithms Based on Solution of Differential Equation69
Thus, the noise amplitude may be smaller than the signal amplitude to start with, but since
the noise frequency is many orders of magnitude higher than the signal frequency, after taking
derivative it will dominate over the signal as can be seen from the expression for derivative,
v= (wVm) cos(wsignalt) + (wnoiseVm,noise) cos(wnoiset) (4.13)
Look at the amplitude of the noise. It is amplified by wnoise. This is problematic because
real life signals are always replete with noise. Thus, in the process of differentiation the noise
will get amplified. This may happen so much so that many times the noise will mask out the
signal of interest. Thus, wherever possible we should try to avoid taking derivative. This leads
to implementation of the differential equation algorithm using numerical integration. We start
with the basic equation relating the relay voltage and current:
di
v = Ri + L (4.14)
dt
Integrating both sides between two instants of time t1 and t2 and then between t3 and t4,
we get:
t2 t2 t2
di(t ) (4.15)
v(t )dt = R i(t ) dt + L dt
dt
t1 t1 t1
t4 t4 t4
di(t ) (4.16)
v(t )dt = R i(t ) dt + L dt
t3 t3 t3 dt
t4 t4 t4
(4.18)
v(t )dt = R i(t ) dt + L di(t ) dt
t3 t3 t3
Note that the last term in each equation does not involve finding the numerical integration
since it reduces simply to:
70Digital Power System Protection
t2
t4
Note, further, that the other integrals can be evaluated numerically. Thus, when we
substitute the numerical values of the integrals, we get algebraic equations. To simplify the
derivation further, let us use the following representation for the numeric values of the integrals:
t2 t2 t2
v(t )dt = E i(t )dt = A di(t )dt = B (4.21)
t1 t1 t1
t4 t4 t4
v(t )dt = F i(t ) dt = C di(t )dt = D (4.22)
t3 t3 t3
Thus, we can write the two differential equations as the following two simultaneous
algebraic equations with constant coefficients:
E = R A + L B (4.23)
F = R C + L D (4.24)
Let us represent the system of equations in matrix form as shown below:
E A B R
F = C D L (4.25)
Thus, we can find the unknowns R and L as follows:
-1
R A B E (4.26)
L = C D F
R 1 D -B E (4.27)
L = ( AD - BC ) -C A F
Solving the above, we get:
( DE - BF )
R= (4.28)
( AD - BC )
( AF - CE )
and L= (4.29)
( AD - BC )
Note that A, C, E and F are the values returned by the numerical integrals defined
above.
Algorithms Based on Solution of Differential Equation71
f ( x ) dx = 2 [ y0 + y1 + y1 + y2 + y2 + + yn -1 + yn -1 + yn ]
h
(4.38)
x0
xn
f ( x ) dx = 2 [ y0 + 2[ y1 + y2 + + yn -1 ] + yn ]
h
(4.39)
x0
Figure 4.9 Simpsons rule for numerical integration.
Algorithms Based on Solution of Differential Equation73
x n
h
f ( x )dx = [first sample + 4(sum of odd indexed samples) + 2(sum of even indexed samples)]
x0 3
(4.42)
% [A ] = [B C ] [R ]
% [D ] = [E F ] [L ]
% [R ] = [ F/d -C/d ] [A ]
% [L ] = [ -E/d B/d] [D ]
% d = 1/(BF-EC)
% R = AF-CD /BF-CE
% L = BD-AE/ BF-CE
% ----------------------------------------------------------
clf ;
clear ,
clc
% ----------------------------------------------------------
f = 50 ;
OverSamplingFactor = 500 ;
Vm = 110 ;
R =1 ;
XbyR = 20 ;
% ----------------------------------------------------------
fsampling = 2 * f * OverSamplingFactor ;
delta_t = 1/fsampling ;
L = XbyR /(2*pi*f) ;
X = 2*pi*f*L ;
phi = atan2(X,R) ;
74Digital Power System Protection
phi_deg = phi*180/pi ;
theta = phi ; % Switching angle
tau = L/R ;
X = 2*pi*f*L ;
Z = R+j*X ;
Zmag = abs(Z) ;
Im = Vm/Zmag ;
phi = atan2(X,R) ;
% v = Vm sin(2*pi*f*t)
% i = Imsin(2*pi*f*t + theta - phi) - Im exp(-t / tau) sin(theta -phi);
k =0 ;
%delta_t = 1e-5 ;
T = 1/50 ;
h = delta_t ; % for Simpsons Rule
% ------- Generation of samples of v and i ----------------------
% -----------------------------------------------------------------
for t=0:delta_t:T
k = k +1 ;
v(k) = Vm * sin(2 * pi * f * t ) ;
i_ss(k) = Im*sin((2*pi*f*t)-phi) ;
i_tran(k) = Im*(exp(-t/tau))*(sin(theta-phi)) ;
%i(k) = Im*sin((2*pi*f*t)-phi)- Im*(exp(-t/tau))*(sin(theta-phi));
i(k) = i_ss(k) - i_tran(k) ;
subplot(3,1,3)
plot(i_tran)
axis([1 length(v) -Im Im])
xlabel(Sample no.)
ylabel(DC offset in current)
Algorithms Based on Solution of Differential Equation75
% --------------------------------------------------------------------
% ------- Creation of two arrays of v() and i() for integration -----
% --------------------------------------------------------------------
samp = 64 ;
for n = 1:samp
v1(n) = v(n) ;
v2(n) = v(n+samp) ;
i1(n) = i(n) ;
i2(n) = i(n+samp) ;
end
v1_intg = my_simpson(v1,h) ;
v2_intg = my_simpson(v2,h) ;
i1_intg = my_simpson(i1,h) ;
i2_intg = my_simpson(i2,h) ;
A = v1_intg; B = i1_intg;C = i1(samp)-i1(1) ;
D = v2_intg; E = i2_intg; F =i2(samp)-i2(1) ;
% ----------------------------------------------------------------
den = B*F - C*E ;
RR = ( A*F - C*D) / (den ) ; % Computed value of L
LL = (B*D - A*E) / (den ) ; % Computed value of R
XX = 2 * pi * f * LL ;
Error_R = ((RR-R)/(R))*100 ;
Error_L = ((LL - L)/(L))*100 ;
% ------------------------------------------------------------------
fprintf(Sampling frequency = %f \n, fsampling);
sprintf(Samples per integration interval %d,samp)
sprintf(Setting R = %f Setting L = %f Setting X=%f,R,L,X)
sprintf(Computed R = %f Computed L = %f Computed X=%f,RR,LL,XX)
sprintf(Error in R = %f %% Error in L = %f %%,Error_R,Error_L)
di di di
dVaa = ( Ra dx )ia + ( La dx ) a + ( Lab dx ) b + ( Lac dx ) c (4.43)
dt dt dt
78Digital Power System Protection
Figure 4.12 Model of three-phase line for developing DEA for three-phase lines.
Assuming that there is a single line to ground bolted fault at length x from the sending
end, the voltage across length x can then be written as:
di di di
Vaa = ( Ra x )ia + ( La x ) a + ( Lab x ) b + ( Lac x ) c (4.44)
dt dt dt
di di di
Vaa = x ( Ra )ia + x ( La ) a + x ( Lab ) b + x ( Lac ) c (4.45)
dt dt dt
d L L
Vaa = x ( Ra )ia + x ( La ) ia + ab ib + ac ic (4.46)
dt La La
Thus, voltage of phase a at the relay location not only depends upon the current in phase
a but also upon currents in the other phases as well as the self-inductance of phase a and
mutual inductances between the phases a and b, a and c.
Note that vaa = va for line to ground fault on phase a. Hence, we can write:
diq L L
va = x ( Ra )i p + x ( La ) where ip = ia and iq = ia + ab ib + ac ic (4.47)
dt La La
Similarly, for b-g fault and c-g fault we can write:
diq L L
vb = x ( Rb )i p + x ( Lb ) where ip = ib and iq = ab ia + ib + bc ic (4.48)
dt Lb Lb
Algorithms Based on Solution of Differential Equation79
diq L L
vc = x ( Rc )i p + x ( Lc ) where ip = ic and iq = ca ia + cb ib + ic (4.49)
dt Lc Lc
In each of these equations, there are two unknowns (xRa) and (xLa) and thus we again have
di(t )
equation of the form v = R i(t ) + L albeit the current terms i(t) and di(t)/dt are different.
dt
Thus, we can follow the same methodologies that we followed for the single phase case.
Figure 4.13 shows the input and output of the phase a ground fault measuring unit. Similar
measuring units in block-diagram form for b to g and c to g faults will have to be set
up to monitor and provide protection against b-g and c-g faults.
Note that we have exclusively used all the phase quantities in the above measurement.
We can instead make use of sequence (symmetrical component) quantities. We will take this
up in a later section.
di d
va = ( Rx )ia + ( LS x ) a + ( LM x ) (3i0 - ia ) (4.53)
dt dt
di di di
va = ( Rx )ia + ( LS x ) a + 3 xLM 0 - xLM a (4.54)
dt dt dt
di di
va = (Rx )ia + x ( LS - LM ) a + 3 xLM 0 (4.55)
dt dt
We have positive and zero sequence inductance of the line L1 and L0 given by:
L1 = LS LM (4.56)
L0 = LS + 2 LM (4.57)
L0 L1 = 3 LM (4.58)
di di
va = ( Rx )ia + x (L1 ) a + x ( L0 - L1 ) 0 (4.59)
dt dt
d L - L1
va = ( Rx )ia + x ( L1 ) ia + 0 i0 (4.60)
dt L1
d
va = x R ia + xL1 (ia + ki0 ) (4.61)
dt
L - L1
where, k = 0 (4.62)
L1
Hence for a three-phase line, for a-g fault we get an equation which is exactly of the form
di(t ) ; albeit the terms in the i(t) and di(t ) are different but known currents.
v(t ) = Ri(t ) + L
dt dt
Hence we can use the differential equation algorithm developed for single phase line.
For catering to b-g and c-g faults we can similarly prove that:
d
vb = xRib + x L1 (ib + ki0 ) (4.63)
dt
d
vc = x R ic + x L1 (ic + k i0 ) (4.64)
dt
Algorithms Based on Solution of Differential Equation81
Figure 4.14Input and output of ag ground fault Figure 4.15Input and output of bg ground fault
measuring unit. measuring unit.
We can write:
dia di di
va = x ( Ra )ia + x ( La ) + x Lab b + x Lac c (4.65)
dt dt dt
dia di di
vb = x ( Rb )ib + x ( Lab ) + x Lb b + x Lbc c (4.66)
dt dt dt
di di
va - vb = x ( Ra )ia - x ( Rb )ib + x ( La ) a - x ( Lab ) a
dt dt
di di di di
+ x ( Lab ) b - x ( Lb ) b + x ( Lac ) c - x ( Lbc ) c (4.67)
dt dt dt dt
R di di di
va - vb = x Ra ia - b ib + x ( La - Lab ) a + x ( Lab - Lb ) b + x ( Lac - Lbc ) c (4.68)
Ra dt dt dt
For a fully transposed line it can be assumed that Lac Lbc, hence the last term in the
above expression vanishes, giving:
R di L - Lb dib
Va - Vb = x Ra ia - b ib + x ( La - Lab ) a + ab (4.69)
Ra dt La - Lab dt
R di L - Lab dib
Va - Vb = x Ra ia - b ib + x ( La - Lab ) a - b (4.70)
Ra dt La - Lab dt
R d L - Lab
Va - Vb = x Ra ia - b ib + x ( La - Lab ) ia - b ib (4.71)
Ra dt La - Lab
diq
Va - Vb = x Ra (i p ) + x ( La - Lab ) (4.72)
dt
Rb Lb Lab
where i p = ia - ib and iq = ia - (ib ) (4.73)
Ra La - Lab
L - Lab R
Note that for a symmetric line b and b will be equal to 1 and the above equation
reduces to: La - Lab Ra
d (ia - ib )
Va - Vb = x Ra (ia - ib ) + x ( La - Lab ) (4.74)
dt
Algorithms Based on Solution of Differential Equation83
Figure 4.18 Inputs and outputs of a-b phase fault measuring unit.
4.7.4 P
hase Fault Protection of Three-Phase Lines Using Sequence
Quantities
We have seen that during a-b phase fault, the voltage (Va Vb) can be expressed as:
d (ia - ib )
Va - Vb = x Ra (ia - ib ) + x ( La - Lab ) (4.75)
dt
diq
Va - Vb = x Ra (i p ) + x ( La - Lab ) (4.76)
dt
Rb L - Lab
where i p = ia - ib and iq = ia - b
L - L (ib ) (4.77)
Ra a ab
L - Lab R
Note that for a symmetric line, b and b will be equal to 1 and the above
equation reduces to: La - Lab Ra
d (ia - ib )
Va - Vb = x Ra (ia - ib ) + x ( La - Lab ) (4.78)
dt
We have: Ra = R, La = LS, the self-inductance and Lab = LM, the mutual inductance
d (ia - ib ) (4.79)
Va - Vb = x R(ia - ib ) + x ( LS - LM )
dt
d (ia - ib )
Va - Vb = x R(ia - ib ) + x ( L1 ) (4.80)
dt
84Digital Power System Protection
Figure 4.19 Inputs and outputs of measuring unit for a-b phase fault.
4.8 E
limination of Specific Harmonics by Integration between
Selected Limits
Again consider the basic differential equation between voltage and current at the relay
location
diq
v = R ip + L (4.81)
dt
Referring to Figure 4.20, if we choose the first limits of integration as 0 to a and second
limit as (p/3) to ((p/3) + a) we get the following two simultaneous equations:
a /w a /w a /w
v dt = R i p dt + L diq (4.82)
0 0 0
p p p
+a w +a w +a w
3 3 3
v dt = R i p dt + L diq (4.83)
p / 3w p / 3w p / 3w
Algorithms Based on Solution of Differential Equation85
p p
+a w
+a w 3
a /w 3 a /w
E= vdt + vdt; A = i p dt + i p dt (4.86)
p / 3w 0 p / 3w
0
p
3 + p
B = iq - iq [0 ] + iq - iq (4.87)
w w 3w
p
+b w p
+b w
b /w 3 b /w 3
F= vdt + vdt; C = i p dt + i p dt (4.88)
0 p / 3w 0 p / 3w
p
b 3 +b p
D = iq - iq [0 ] + iq - iq (4.89)
w w 3w
Review Questions
We are interested in finding the value of x which minimises the sum of the squares of
errors. Hence, this value will be that for which, the following two conditions are met:
d (SSE)
(i) = 0 and
dx
d 2 (SSE)
(ii) = positive
dx 2
Let us find the value of x for which the first derivative of SSE w.r.t. x is zero, i.e.,
d d2
(( x - 2)2 + ( x - 3)2 ) = 0 and 2
(( x - 2)2 + ( x - 3)2 ) > 0
dx dx
2( x - 2) + 2( x - 3) = 0
10 5
2 x - 4 + 2 x - 6 = 0 fi 4 x - 10 = 0 fi 4 x = 10 fi x = fi x = = 2.5
4 2
Algorithms Based on Least Squared Error (LSQ)89
d2 d
Further, (( x - 2)2 + ( x - 3)2 ) = (4 x - 10) = 4 which is > 0
dx dx
5
Hence, x = = 2.5 is the value at which the sum of the squares of errors is indeed minimum.
2
5 3+2
Now, it is easy to see that x = = = 2.5
2 2
Thus, the estimated value with minimum error that we arrived at by minimising the sum
of the squares of the errors is the same as the mean or the average value. Hence, there is
greater significance to the average value as a value which minimises the sum of squares of
errors, which is the closest that we can approach the truth.
1 2
[ x ] = [1 1]
2 3
1 2
[ x ] = [1 1]
2 3
1 1 2
[ x] =
2 2 3
2 3
[ x] = +
2 2
[ x ] = 2 = 2.5
5
Representing the problem in symbolic fashion, we started with the matrix equation:
[A]21 [x]11 = [c]21
T
Pre-multiplying both sides by [A]
[A]T[A]21[x]11 = [A]T [c]21
Pre-multiplying both sides by {[A]T[A]21}1
{[A]T[A]21}1 {[A]T[A]21}[x]11 = {[A]T[A]21}1[A]T[c]21
[I]11[x]11 = {[A]T[A]21}1[A]T[c]21
[x]11 = [A+][c]21
where [A+] = {[A]T[A]21}1[A]T is called the pseudo-inverse of the rectangular matrix A whose
number of rows is greater than number of columns. Pseudo-inverse has properties very similar
to inverse. Some of the properties of the pseudo-inverse are given below:
(i) (P P+) P = P
(ii) (P+ P) P+ = P+
(iii) (P+ P) = (P+ P)T
(iv) P P+ = (P P+)T
Pseudo-inverse can also be determined for under-determined system as:
[A]T {[A]21[A]T}1
[n T ][n ] = ( Z T - x T H T )([ Z ] - [ H ][ x ])
(SSE)11 = ( Z T [ Z ] - [ x ]T [ H T ][ Z ] - Z T [ H ][ x ] + [ x ]T [ H T ][ H ][ x ])
Note: [ x ]T [ H T ] [ Z ] = Z T [ H ][ x ]
SSE = ( Z T [ Z ] - 2[ x ]T [ H T ][ Z ] + [ x ]T [ H T ][ H ][ x ])
As already seen d(SSE)/d x = 0 and d2(SSE)/d( x)2 > 0 for minimising the sum of squares
of errors. Thus, differentiating the SSE and equating it to zero gives:
d (SSE)
= ( - 2[ H T ][ Z ] + 2[ H T ][ H ][ x ]) = 0
d ( x )
([ H T ][ H ][ x ]) = [ H T ] [ Z ]
[ x ] = [ H ]+ [ Z ]
+ -1
where [ H ] = [[ H ][ H ]] is the pseudo-inverse of H.
T
d 2 (SSE)
then 2
= [ H T ][ H ]
dx
T
Now, square matrix of the type [H ][H] is always positive definite. Hence both, the first
order and second order conditions for minima are met.
Hence, we can say that the pseudo-inverse method minimises in the LSE or LSQ sense.
Thus, we are justified in calling the algorithm by Sachdev as an LSQ/LSE algorithm.
x2
e x =1 + x + +
2!
I0t I0t 2
writing (I0et/t ) using only three terms of the expansion as: I 0 e - t / t @ I 0 -
+ 2.
t 2t
Thus, the model of the fault current with which we will be trying to fit the samples of the
fault current becomes:
I0t I0t 2
+ 2 + Im1 sin(w0t + q1) + Im3 sin(3w0t + q3)
y(t) = f(t) = I 0 -
t 2t
Further we will expand the sin(w0t + q0) and sin(3w0t + q3) terms to get:
I0t I0t 2
i(t ) = I 0 - + 2 + Im1 sin(w0t) cos(q1) + Im1 cos(w0t) sin(q1)
t 2t
+ Im3 sin(3w0t) cos(q3) + Im3 cos(3w0t) sin(q3)
By Rearranging, we get:
I0t I0t 2
i(t ) = I 0 - + 2 + Im1 cos(q1) sin(w0t) + Im1 sin(q1) cos(w0t)
t 2t
+ Im3 cos(q3) sin(3w0t) + Im3 sin(q3) cos(3w0t)
Again by rearranging, we get:
i(t) = I0 + Im1 cos(q1) sin(w0t) + Im1 sin(q1) cos(w0t)
I0t I0t 2
+ Im3 cos(q3) sin(3w0t) + Im3 sin(q3) cos(3w0t) - + 2
t 2t
Let us denote i(t) since it is a sample value by s(t).
s(t) = I0 + Im1 cos(q1) sin(w0t) + Im1 sin(q1) cos(w0t)
I 0t I 0t 2
+ Im3 cos(q3) sin(3w0t) + Im3 sin(q3) cos(3w0t) - + 2
t 2t
In the above equation we have seven unknowns. Hence, we need to generate a total of
seven equations to be able to solve for the seven unknowns. This is shown below with unknown
quantities in each equation shown in bold.
94Digital Power System Protection
In each equation there are seven knowns, which can be pre-computed. For example, in the
first equation these are:
a11 = 1, a12 = sin(w0t1), a13 = cos(w0t1), a14 = sin(3w0t1), a15 = cos(3w0t1), a16 = t1 and a17 = t12
The seven unknowns are:
x1 = I0, x2 = Im1 cos(q1), x3 = Im1 sin(q1), x4 = Im3 cos(q3),
I I
x5 = Im3 sin(q3), x6 = 0 , x7 = 02
t 2t
The amplitudes and phase angles of the fundamental and the third harmonic can be found
as shown below:
IIm3,theta3,theta3_calc)
plot(s)
%---------------------------------------------------------------
%---------------------- End of Program ----------------------
%--------------------------------------------------------------%
Review Questions
For the sake of completeness, we list the frequencies of various visible colours (in terra-
hertz) and audible sounds (in Hz) as shown in Tables 6.1 and 6.2 respectively, which by the
way, shows that we are almost deaf as well as blind since we can perceive only a small fraction
of sound and light spectra.
99
100Digital Power System Protection
Deaf we may be, but it can be seen from Table 6.2 that all of us, whether we are musically
trained or not, have the ability to correctly decode the ratio of two frequencies. The ratio of
frequency of adjacent (pure) notes on the piano is exactly 21/7, so that the eighth note (first
note in the next octave) is double in frequency as compared to the first note! In fact that is the
origin of the word octave (Latin for eight). This is a tribute to our ears, which over millions
of years, evolved as specialised instruments for doing real-time frequency analysis of very high
complexity. In fact our ears have evolved dedicated hardware to do frequency analysis before
sending the signal to the brain.
Thus, sound, which is information in time domain, encoded in time varying pressure wave
in air, is actually converted into a frequency domain signal by the mechanism of the ear.
What we finally perceive is not the instantaneous variations of pressure but a tone, i.e., a
frequency with a certain magnitude which we perceive as weak or strong. This is exactly the
job of an instrument called spectrum analyser. Thus, we have a spectrum analyser for sound
in the audible frequency range of 20 Hz to 20 kHz.
It is interesting to note that it is also possible to trick our ears into believing what we
want it to believe if we program the sounds correctly. We can create the illusion of a stereo
sound or surround sound by manipulating the time of arrival of sound signals. Similarly, we
can substantially reduce perception of noise in a recording by using some more tricks like
DolbyTM techniques.
Similarly our eyes perceive the time domain light signals as colours, i.e., frequencies with
their respective strengths (lightness or darkness of the colour) before passing it to the brain.
We cannot tell the instantaneous value of the light signal. Thus, our eyes constitute a spectrum
analyser working in the visible frequency range.
This is a fascinating topic but we are constrained to leave it here and concentrate on
electrical engineering!
Today our thinking is so much dominated by digital computers that it is difficult to imagine
Fourier analysis being done by any other method than by using digital computers. However,
Discrete Fourier Transform101
using computers for Fourier analysis is of comparatively recent origin. For tens of years,
scientists have done Fourier analysis using mechanical, acoustic or optical methods. Prior to
1965 such devices were in wide use for Fourier analysis. However, with the publication of
the FFT algorithm by Cooley and Tuckey in 1965, all other methods have been superseded by
digital computer oriented methods.
Now, how to find the coefficients C0, C1, C2 etc.? This is where the property of orthogonality
becomes extremely useful. If we multiply both sides of above equation by f0 and integrate
between a and b, we get:
b b b b
f (t )f0 (t )dt = C0 f0 (t )dt + C1 f0 (t )f1 (t )dt + + Cn f0 (t )fn (t )dt +
2
a a a a
a a
Thus, we get:
b
f (t )f0 (t )dt = k0 C0
a
b
1
C0 =
k0
f (t )f0 (t ) dt
a
In general,
1 b
Cn = f (t )fn (t ) dt
kn a
Thus, multiplying by fn(t) and integrating between a and b has the effect of filtering
out all the other terms except the n the term.
Note that the above form is phase angle free. The phase angle is not explicitly stated but
is encoded in the ratio of bn and an. Another form which is sometimes used is a cosine wave
with phase angle, explicitly stated. To develop that form consider:
Cn cos(n w 0 t + q n )
i.e., Cn cos(n w 0 t + q n ) = Cn cos(n w 0 t ) cos(q n ) - Cn sin(n w 0 t )sin(q n )
b
Let C0 = a0 , then Cn = (an2 ) + (bn2 ) and q n = - tan -1 n (Note: - ve sign)
an
Then, Cn cos(n w 0 t + q n ) = an cos(n w 0 t ) + bn sin(n w 0 t )
The coefficients of the Fourier series can be found by evaluating the following integrals:
1T
0 T f (t )dt
a = is the average value of the function over a cycle
0
2T
an = f (t )cos(n w 0 t )dt
T0
2T
bn = f (t )sin(n w 0 t )dt
T0
The above simple expressions for the coefficient become possible because of the orthogonal
nature of the following set of functions, which form the basis functions of Fourier series:
{1, cos(w0t), cos(2w0t), cos(3w0t),, cos(nw0t), sin(w0t), sin(2w0t), sin(3w0t),, sin(nw0t)}
The orthogonal relationship between the basis functions is captured in Table 6.3.
104Digital Power System Protection
T
T
sin(nw 0t ) zero zero
2
=p zero zero
0
T
T
cos(m w 0t ) zero zero zero
2
=p zero
0
T
T
cos(nw 0t ) zero zero zero zero
2
=p
0
1 1
cos(nw 0 t ) = (e jnw 0 t + e - jnw 0 t ) and sin(nw 0 t ) = (e jnw 0 t - e - jnw 0 t )
2 2j
Hence, we can write f(t) as:
where Cn is the coefficient in the cosine-phase angle form of Fourier series. f(t) can now be
expressed in compact form
f (t ) = [Cn~ e jnw t ]
0
n =-
1T
Cn~ = f (t ) [cos(nw 0 t ) - jsin(nw 0 t )] dt
T0
1T - jnw t
Cn~ = f (t )e 0 dt
T0
Fourier analysis is a search for the values of and bn coefficients which fit the given
curve with the help of sines and cosines. We should try to find the values of an and bn which
are optimal in the LSQ sense. Assuming N number of terms in the series, we can express
the error as:
N
error = f (t ) - a0 + (an cos(n w 0 t ) + bn sin(n w 0 t ))
n =1
Sum of squared errors can be written as:
T 2
N
SSE = f (t ) - a0 + (an cos(n w 0 t ) + bn sin(n w 0 t )) dt
0 n =1
For minimising the SSE, we impose the conditions that:
T T T N
0 = 2 f (t )dt - 2 a0 dt - 2 (an cos(n w 0 t ) + bn sin(n w 0 t ))dt
0 0 0 n =1
Note that each of the term within the summation sign integrates to zero. Hence, we are
left with:
T T
0 = 2 f (t )dt - 2 a0 dt
0 0
T T
0 = f (t )dt - a0 dt
0 0
T T
a0 dt = f (t )dt
0 0
T
a0 T = f (t )dt
0
1T
a0 = f (t )dt
T0
.... which is the same expression derived earlier. Hence, a0 is best in
the LSQ sense.
T N N T
- bn sin(n w 0 t )cos(n w 0 t ) dt - bn cos(m w 0 t )cos(n w 0 t ) dt
0 n =1 0 m =1, m n
T
T
N N
n b sin( n w 0 t )cos( n w 0
t ) dt = bn cos(m w 0 t )cos(n w 0 t ) dt = 0
0 n =1 0 m =1, m n
T
T
and (cos(n w 0 t )) dt =
2
; thus we are left with:
0 2
Discrete Fourier Transform107
T
T
0 = f (t )cos(nw 0 t )dt - an
0 2
2T
an = f (t )cos(nw 0 t )dt
T0
which is the same expression as obtained earlier.
1T 2T 2T
a0 = f (t )dt
T0
an = f (t )cos(n w 0 t )dt
T0
bn = f (t )sin(n w 0 t )dt
T0
a0
f (t ) = + (an cos(n w 0 t ) + bn sin(n w 0 t ))
2 n =1
2T 2T
an = f (t )cos(nw 0 t )dt
T0
bn = f (t )sin(n w 0 t )dt
T0
n = 1 to
3. Trigonometric Form of Fourier Series with Cosine Wave and Phase Angle
f (t ) = C0 + Cn cos(n w 0 t + q n )
n=0
b
q n = - tan -1 n
C0 = a0 Cn = (an2 ) + (bn2 ) an
(Note: ve sign)
108Digital Power System Protection
an - jbn an + jbn
a0 = C0 = Cn~ = C-~n
2 2
1T - jnw t
Cn~ = f (t )e 0 dt
T0
Fourier analysis is useful in a very large number of fields. In fact the double helix form
of DNA was discovered in 1962 through X-ray diffraction techniques and Fourier analysis.
NASA improves the clarity and detail of picture of celestial objects taken in space by means
of Fourier analysis. Fourier analysis has found application in plasma physics, semiconductor
physics, microwave acoustics, seismography, oceanography, radar mapping, medical imaging,
chemical analysis and biology to name a few areas.
First let us find the Fourier series representation for the pulse and plot its amplitude and
phase spectrum. The complex Fourier coefficients Cn~ are given by:
T
Cn~ = f (t ) e - jnw 0 t dt
0
a/2
V0
Cn~ = e - jnw 0 t dt
T -a / 2
a/2
V e - jnw 0 t V0 e - jnw 0 a / 2 - e jnw 0 a / 2
Cn~ = 0 =
T - jnw 0 - a / 2 T ( - jnw 0 ) 1
1 V0 2 e jnw 0 a / 2 - e - jnw 0 a / 2
Cn~ =
T n 2 p f0 1 2j
110Digital Power System Protection
V0 e jnw 0 a / 2 - e - jnw 0 a / 2 w a
Cn~ = ; multiplying and dividing by n 0 , we get:
np 2j 2
V0 n w 0 a sin(nw 0 a/2)
Cn~ =
np 1 2 (nw 0 a/2)
V0 a sin(nw 0 a/2)
Cn~ =
T (nw 0 a/2)
V0 a sin( x ) nw 0 a
Cn~ = where x =
T ( x) 2
Discrete Fourier Transform111
V0 a
Cn~ = sinc( x )
T
sin( x ) sin( x )
where sinc( x ) is defined as . It can be shown that the value of is 1 when x = 0.
( x) ( x)
The c in sinc() comes from the Latin cardinalis. The full name of the sinc() function
is sinus-cardinalis. The plot of sinc() function using MATLAB ezplot() function is shown in
Figure 6.4. In Figure 6.5 plot of Cn~ and magnitude and phase spectrum of sinc() function are shown.
It can be shown that as T , the discrete frequency spectrum remains the same in shape
sin(nw 0 a/2)
but assumes the envelope of the spectrum of the recurring pulse. The envelope of
sin(w a/2) (nw 0 a/2)
becomes . Note that in the first expression we have w0, a discrete frequency
(w a/2)
whereas in the second expression we have w the continuous frequency variable.
The mathematical treatment of the non-recurring pulse can thus be derived from the
recurring pulse by letting the time period T approach in the limit. The changes that take place
in the representation as we transit from recurring to non-recurring pulse are listed in Table 6.4.
112Digital Power System Protection
Table 6.4
Recurring pulse Single pulse
Time period Time period T is finite Time period T
Frequency Fundamental frequency is (1/T) Fundamental frequency ffunda is zero
Nature of Spectrum is discrete (spaced at ffunda Spectrum is continuous (spaced at zero Hz
frequency from each other) from each other)
spectrum
Summation/ Fourier sum Fourier integral
Integration 1
f (t ) = Cn~e jw n t
0
f (t ) = jw t
F ( jw )e dt
n =- 2p -
1T - jnw t
Coefficients Cn~ = f (t )e 0 dt F ( jw ) = f (t )e - jw t dt
T0 -
Last row of Table 6.4 shows that the difference between the frequency spectrum of the
recurring and non-recurring versions of a pulse is not in the shape. Both versions have the
same shape. The difference is that the recurring pulse has discrete (discontinuous) frequency
spectrum while the non-recurring pulse has continuous frequency spectrum. Thus, we can make
the following observations:
1. The shape of the continuous amplitude and phase spectra for a non-recurring f(t) is
identical with the envelopes of the amplitude and phase line spectra for the same pulse
recurring.
2. All frequencies are present in the continuous amplitude spectrum since F(jw) is defined
for all w from to +.
3. Lightning, which is a very short duration pulse in time, makes its presence felt
irrespective of the frequency of the radio station to which we may be tuned. (We do
not hear the lightening crashes in our mobile phones because, we are not listening to
the raw audio signal, the phone processes the audio signal to filter out noise using
DSP techniques.)
We draw an important inference from the above discussion that in case of non-recurring
pulse, we have to use Fourier integral instead of the Fourier summation. Fourier integral is also
popularly known as Fourier transform. Note that the Fourier transform that has been described
is a continuous transform. As can be expected, this continuous transform will have to give way
to discrete transform. But this will be taken up in the next section.
n =-
Discrete Fourier Transform113
1T 1T
Cn~ =
T0
f (t )cos( nw 0 t ) dt - j f (t )sin(nw 0 t ) dt
T0
1T
Cn~ = f (t )[cos(nw 0 t ) - jsin(nw 0 t )] dt
T0
1T - jnw t
Cn~ = f (t )e 0 dt
T0
Table 6.5
Fourier sum (Recurring signal) Fourier integral (Windowed signal)
1
f (t ) = Cn~e jw n t
0
f (t ) = F ( jw )e
jw t
dt
n =- 2p -
We develop the DFT by noting the expression for F(jw) for windowed signal:
F ( jw ) = f (t )e - jw t dt
-
Note that w refers to angular frequency of the signal in the continuous domain, i.e., wsignal.
wsignal = 2pfsignal, where fsignal is the signal frequency in the continuous domain. Note that time
t can take only discrete values t = nDt; where Dt = 1/fsampling in the discrete domain and fsignal
can take only discrete values corresponding to mffunda. Thus changing over to the discrete
domain where X(m) represents the Fourier coefficient at the mth harmonic frequency, we can
write the discrete version of the Fourier integral as:
N -1
- j 2p fsignal nDt
X ( m) = x (n)e
n=0
- j 2p m ffunda n
N -1
x (n)e
fsampling
X ( m) =
n=0
n=0
114Digital Power System Protection
We summarise the transition from CFT to windowed DFT as shown in Table 6.6.
j 2p mn
f (t ) = Cn~e jw n t
0
1 jw t 1 N -1
n =-
f (t ) =
2p
F ( jw )e dt x ( n) = X ( m)e N
- N m =0
- j 2p mn
1T - jnw t
N -1
Cn~ = f (t )e 0 dt F ( jw ) = f (t )e - jw t dt X ( m) = x ( n) e N
T0 - n=0
It is worth noting that Fourier analysis assumes that we know the frequency of the signal
apriori because otherwise we would not have been able to decide the sampling frequency.
The primary purpose of Fourier analysis is not frequency measurement. Its purpose is to
estimate magnitudes of various frequency components. The frequencies of these components
automatically are tied up with the sampling frequency and therefore known beforehand (apriori).
Thus we can write:
1 1 fsampling
Tfunda = N Tsampling =N fsampling = Nffunda ffunda =
ffunda fsampling N
This is illustrated in Figure 6.9.
n=0
n=0
Discrete Fourier Transform117
N -1 - j 2p ( N ) n j 2p m n
X ( N - m) = x ( n) e N e N
n=0
N -1 + j 2p m n
X ( N - m) = x ( n) e - j 2 p n e N
n=0
n=0
Thus, the only difference between X(m) and X(N m) is the sign of the exponents, which
is reversed. Hence, X(N m) terms are complex conjugate of X(m) terms. Thus, the magnitudes
of X(N m) and X(m) will be identical but their phase angles will be negative of each other.
Hence, we can write:
X(N m) = {X(m)}* or {X(N m)}* = X(m) We are assuming here that the signal x(n) is real.
Table 6.7
For N = 8; assuming that the signal x(n) is real.
m=0 m=1 m=2 m=3 m=4 m=5 m=6 m=7
Nm=8 Nm=7 Nm=6 Nm=5 N m =4 Nm=3 Nm=2 Nm=1
X(0) = X(1) = X(2) = X(3) = X(4) = X(5) = X(6) = X(7) =
X(8)* X(7)* X(6)* X(5)* X(4)* X(3)* X(2)* X(1)*
Note: X(8) does not exist and X(0) is always real. Since X(4)* = X(4), X(4) is also real.
Remaining Fourier coefficients are complex.
Sample no. 0 1 2 3 4 5
Value 253.000 275.2891 161.2628 77.0000 74.0526 65.7109
Sample no. 6 7 8 9 10 11
Value 33.0000 31.9212 35.9474 11.0000 29.2628 89.0788
Note: The DFT formula does not include the multipliers of 1/N for DC term
and 2/N for other terms. This is to speed-up the DFT calculation by avoiding
an extra division. The additional division is time consuming and hence avoided.
The result is not wrong because we can always factor in the 1/N and 2/N terms.
Sometimes these may also get cancelled out during division. For example, if
we are finding impedance then 2/N factor in voltage and 2/N factor in current
will get cancelled.
120Digital Power System Protection
n =11 p p p p p
- j n - j 0 - j 1 - j 2 - j 3
6 6 6 6 6
X (1) = x (n)e = x (0)e + x (1)e + x (2)e + x (3)e
n=0
p p p p p p
- j 4 - j 5 - j 6 - j 7 - j 8 - j 9
6 6 6 6 6 6
+ x (4)e + x (5)e + x (6)e + x (7)e + x (8)e + x (9)e
p p
- j 10 - j 11
6 6
+ x (10)e + x (11)e
On the Casio fx-991ES and similar calculators, we can easily evaluate the above
equation by noting that re jq = rq . Hence, we can write the expression for X(1) as
follows:
n =11 p
- j ( )n -1.p -2.p -3.p
X (1) = x ( n) e 6 = x (0) + x (1)
6 + x (2)
6 + x (3)
6
n=0
n =11 p
- j n
6
X (1) = x (n)e = x (0) + x (1)- 30 + x (2)- 60 + x (3)- 90 + x (4)- 120
n=0
6. If all is well you will see the answer on the LCD as:
461.999789 396.0000809i
2 2
C1 = (461.999789 396.0000809 i) = (a1 - jb1 )
N 12
Hence,
a1 = (461.999789/6) = 76.99996483 77 (Answer)
and b1 = (396.000809/6) = 66.00001348 66 (Answer)
(iii) To find a2 and b2 of 2nd harmonic for m = 2 and n = 0 to 11 in the general expression
for DFT:
n =11 n =11 n =11 p
- j n
3
X (2) = x (n)e- j 2p 2 n /12 = x (n)e- jp n / 3 = x ( n) e
n=0 n=0 n=0
n =11 p p p p p
- j n - j 0 - j 1 - j 2 - j 3
3 3 3 3 3
X (2) = x ( n) e = x (0)e + x (1)e + x (2)e + x (3)e
n=0
p p p p p
- j 4 - j 5 - j 6 - j 7 - j 8
3 3 3 3 3
+ x (4)e + x (5)e + x (6)e + x (7)e + x (8)e
p p p
- j 9 - j 10 - j 11
3 3 3
+ x (9)e + x (10)e + x (11)e
n =11 p
- j n -1.p -2.p -3.p
x (n) e
3
X (2) = = x (0) + x (1) + x (2) + x (3)
n =0 3 3 3
Hence,
330
a2 = = 55 (Answer)
6
264.00005011
b2 = = 44 (Answer)
6
Let us tabulate the results:
DC a1 b1 a2 b2 a3 b3
88 77 66 55 44 33 22
max_error = 5.8524e-013
a0_calc = 88.000000
a1_calc = 76.999996
b1_calc = 66.000013
a2_calc = 55.000000
b2_calc = 44.000008
a3_calc=33.000000
b3_calc=22.000000
n = N + k -1 - j 2 p (n-k )m
X k ( m) = x (n)e N
n=k
126Digital Power System Protection
where k is the window number and m is the order of the DFT; m = 0 DC, m = 1
fundamental etc.
Let us run through a few values of the indices to verify that the basic DFT expression
remains same.
Let N = 8
window k n = k to (N + k1) (nk)
0 07 07
1 18 07
2 29 07
3 3 10 07
4 4 11 07
n=k
n = k +1
If we start the above summation from n = k, instead of n = k + 1, we will get an extra term:
- j 2p ( -1) m j 2p m
x ( k )e N = x ( k )e N
Thus, if we modify the lower limit of DFT for (k + 1)th window to (n = k) then we will
j 2p m
get the above extra term, hence we can keep the sum unaltered by subtracting x (k )e N
from
Discrete Fourier Transform127
the sum with lower limit (n = k) instead of (n = k + 1). Hence, we can write it as:
n = N + k - j 2p ( n - k -1) m j 2p m
X k +1 (m) = x (n)e - N
N
x ( k ) e
n = k
If we take out the last term of the above summation corresponding to n = (N + k), we will
have to reduce the upper index of the summation to (N + k 1) and the term corresponding
to n = (N + k) will be:
- j 2p ( N + k - k -1) m - j 2p ( N -1) m
x ( N + k )e N = x ( N + k )e N
- j 2p ( N -1) m - j 2 p ( N ) m j 2p m
x ( N + k )e N = x ( N + k )e N e N
j 2p m
- j 2p m
= x( N + k ) e e N
j 2p m
= x(N + k ) e N since e - j 2p m is equal to 1 for all values of m.
Hence, we can write:
n = N + k -1 - j 2 p ( n - k -1) m j 2p m j 2p m
N - x ( k )e N
X k +1 (m) = x (n)e N
+ x ( N + k ) e
n = k
n = N + k -1 - j 2 p ( n - k ) m j 2p m j 2p m j 2p m
X k +1 (m) = x (n)e N e N
+ x( N + k ) e N - x ( k )e N
n = k
j 2p m n = N + k -1 - j 2p ( n - k ) m
X k +1 (m) =e N x ( n) e N
+ x ( N + k ) - x (k )
n = k
j 2p m
X k +1 (m) = e N
{X k (m) + x( N + k ) - x(k )}
j 2p m
DFT(m) k +1 = e N
{DFT(m)k + [Newest sample - Oldest sample ]}
Consider the sliding DFT computation of an undistorted sine wave shown in Figure 6.14.
In the first window, the samples from x(0) to x(7) will be used for computing DFT(1)1. Next
we slide the widow by one sample discarding x(0) and including x(8). Let us find how the
DFT(1)2 is related with DFT(1)1.
j 2p 1
DFT(1)2 = e 8
{DFT(1)1 + [ x(8) - x(0)]}
128Digital Power System Protection
However, in this case, since x(8) = x(0), magnitude of DFT remains unaffected but there
is an apparent phase shift of p/4 radians or 45 as the window slides along even if the signal
is stationary. Thus:
jp
DFT(1)2 = e4 {DFT(1)1 }
This is an unwanted artefact and provides motivation for modifying the DFT further as
described in the next section.
{X k (m) + x( N + k ) - x(k )}
X k +1 ( m ) = e N
- j 2p ( k +1)
Multiplying both sides of equation for DFT by e N , we get:
- j 2p ( k +1) - j 2p ( k +1) j 2p
X k +1 (1)e N =e N e N {X k (1) + x( N + k ) - x(k )}
- j 2p ( k )
(1) = e
X k +1
N
{X k (1) + x( N + k ) - x(k )}
- j 2p ( k ) - j 2p ( k )
(1) = X (1)e N + [ x ( N + k ) - x (k )]e N
X k +1 k
Using the definition for the modified DFT, it can be seen that:
- j 2p ( k )
X k (1)e N (1)
=X
k
clear , clc
fsig = 50 ; T=1/fsig;
N = 12 ;
fs = N * fsig ;
dt = 1/fs ;
% -----------------------------------------------------------
%---------- Lets generate samples of signl ----------------
Im1 = 200 ; Im2 = 1000;
phi1 = 45*(pi/180) ; phi2 = 0 ;
n = 0 ;
for t=0:dt:(T-dt)
n = n+1;
130Digital Power System Protection
x(n)=Im1*sin(2*pi*fsig*t+phi1);
end
fprintf(Next 12 Samples\n)
disp(x(13:24))
% Next 12 samples are from 50 Hz signal with peak value
of 1000
%--------------------------------------------------------
%---------------- DFT calculation ----------------------------
% Initialization of coeffs. to zero before start of
calculations
for m = 0:1 % m is harmoic number 0=dc , 1 = fundamental
X(m+1)= 0 ;
end
X_hat_mag =(2/N)*abs(X_hat)
% ----- Plotting of Sliding DFT -----
subplot(2,1,2)
stem(X_hat_mag)
xmin1=0;xmax1=13;ymin1=0;ymax1=1001;
axis([xmin1 xmax1 ymin1 ymax1])
grid on
xlabel(Sliding Window Number --->)
ylabel(Sliding DFT Magnitude)
title(Sliding DFT)
% ------------- End of program ------------------------
Review Questions
1. What are the conditions to be fulfilled by a signal for being represented by Fourier
series?
2. What are the Dirichlet conditions?
3. State Fourier series in trigonometric form without phase angle.
132Digital Power System Protection
4. State Fourier series in trigonometric form as a cosine wave with phase angle.
5. State Fourier series in trigonometric form as a sine wave with phase angle.
6. State the relationship between coefficients of Fourier series in Q.4 with Fourier series
in Q.3.
7. State the relationship between coefficients of Fourier series in Q.5 with Fourier series
in Q.3.
8. Derive expression for Fourier series in complex exponential form.
9. State the relationship between coefficients of Fourier series in Q.8 with Fourier series
in Q.3.
10. Prove that the an and bn coefficients of the Fourier series are optimal in LSQ sense.
11. Compare the expressions for CFT and DFT and point out similarities between them.
12. If the DFT has N frequency terms, what is the highest harmonic frequency in terms
of multiples of fundamental?
13. In a protective relay hardware, 8 samples at a sampling frequency of 1000 Hz are
collected in the DFT window. What is the frequency of the fundamental, 2nd, 3rd,
4th harmonic etc.?
14. All samples inside the DFT window have same magnitude. What frequency content
is indicated?
15. Only one sample is inside the DFT window, at the centre is non-zero, all other samples
are zero. What frequency content is indicated?
16. If there are N samples in time domain, how many samples will be there in the
frequency domain?
17. What is the effect of increasing number of samples in a given window on the precision
of DFT?
18. What is special about X(0)?
19. What is special about X(N/2); assuming N is even.
20. Explain the statement: all N DFT coefficients are not unique.
21. If x(n) are all real, is it guaranteed that X(n) will all be real?
22. Why conventional Fourier series theory fails when we try to apply it to non-recurring
waveform?
7
FFT and Goertzel Algorithm
The expression for DFT is normally written in the form shown below:
n = N -1
X (m) mm == 0N -1 = x (n)e - j 2p mn/N
n=0
Let us represent ej2p/N = WN (Note that negative sign is included in WN) and rewrite expression
for DFT as shown below:
n = N -1
X (m) mm == 0N -1 = x (n)WNmn
n=0
n=7
X (m) mm == 07 = x (n)W8mn
n=0
This is actually a matrix equation since we have a N 1 column vector on LHS and
product of a (N N) matrix with a (N 1) column vector on the RHS. Hence, we can write
the equation as:
133
134Digital Power System Protection
Table 7.2 shows the substantial saving in computational effort offered by FFT assuming
N = 213 = 8192.
The FFT algorithm was invented by J.W. Cooley and J.W. Tuckey in 1965. In those days,
digital computer capabilities were rather limited in terms of processing speed and memory
storage etc. Hence, FFT algorithm had a dramatic impact upon the scientific community.
In fact, it started a revolution in signal processing. Many applications where one would not
think of applying FFT (because of the high cost) could now benefit from FFT. Eventually, all
other methods of doing Fourier analysis were superseded by digital computer based FFT. It is
interesting to note that like most other mathematical inventions even FFT can be traced back
to Carl Fredrich Gauss, the 18th century mathematical genius.
Consider an 8-point DFT. We have taken out the relevant portion of the matrix equation
for extraction of fundamental. Another way to look at it is as a dot product or inner product
of two vectors. One vector consists of the 8 weights and the other vector is composed of
the 8 samples in time domain. Weights for computation of fundamental can be visualised
as shown in Figure 7.1.
Figure 7.1 Visualisation of weights for computation of fundamental in an 8-point DFT (contd.).
136Digital Power System Protection
To develop the decimation-in-time FFT algorithm, we start with the expression for DFT:
n = N -1 n = N -1
X (m) mm == 0N -1 = x (n)e - j 2p mn/N = x (n)WNmn
n=0 n=0
Decomposing the sum as that of even indexed and odd indexed time domain samples:
N N
n= -1 n= -1
2 2
X (m) mm == 0N -1 = x (2n)e - j 2p m (2 n )/N
+ x (2n + 1)e - j 2p m (2 n +1) /N
n=0 n=0
N N
n= -1 n= -1
2 2
X (m) mm == 0N -1 = x (2n)e - j 2p m (2 n )/N
+e - j (2p m )/N
x (2n + 1)e - j 2p m (2 n)/N
n=0 n=0
Consider the first part of the sum. By comparing it with standard DFT expression, we can
conclude that, the expression represents an N/2 point DFT of the even indexed time domain
samples.
N
n= -1 (2 n)
2 - j 2p m N
x (2 n)e N This is the
2
point DFT of even indexed samples
n=0
Consider the summation in the second part of the expression for DFT. Again it can
be seen that it is another N/2 point DFT, but this time, of the odd indexed time domain
samples.
138Digital Power System Protection
N
n= -1
2 N
x (2 n + 1)e - j 2p m (2 n)/N This is the
2
point DFT of odd indexed samples
n=0
Recall the definition of WN; WN = ej2p/N. Hence the remaining term, i.e., ej(2pm)/N can be
written as WNm. Thus, what we have done is the decomposition of the N-point DFT as:
[mth term of N-point DFT] = [N/2 point DFT(even indexed samples)]+ WNm [N/2
point. DFT(odd indexed samples)]
Thus, N-point DFT is not simply a sum of two N/2 point DFTs, we have to use a weighting
factor of WNm before combining the two DFTs. Using more compact notation we can write:
N
- j 2p mn /
- j 2p m (2 n )/N 2
e =e = WNmn
/2
N N N
X m + ; m = 0 to - 1 fi X to X ( N - 1)
2 2 2
N N N
m= -1 n= -1 N N n = -1 N
N 2 2 m + n m + 2 m + n
X m +
2 m =0
= x (2n)WN/2 2 + WN 2 x (2 n + 1)WN/2 2
n=0 n=0
N
m+ n
Consider WN/2 2
N N - j 2p nN/2
m + n n
WN/2 2
= WN( mn ) 2
/2 WN/2 = WN( mn )
/2 e
N/2
- j 2p nN/2
e N/ 2 = e - j 2p n = 1 (always)
Hence,
N
m + n
WN/2 2
= WN( mn
/2
)
FFT and Goertzel Algorithm139
N
m +
Similarly, considerWN 2
N N
m +
WN 2
= WN( m )WN 2
N
Consider WN 2
N - j 2p N/2
WN 2 = e N
= e - jp = -1
N
m +
Hence, WN 2
= -WN( m )
i.e., weights for the second half of DFT terms are negative of those for the first half. Thus,
we can write:
N N N
m= -1 n= -1 n= -1
N 2 2 2
X m +
2 m =0
= x (2n)WN( m/2) n - WN( m ) x (2n + 1)WN( m/2) n
n=0 n=0
N N
n= -1 n= -1
2 2
X (m) m =(NN-1) = /2 - WN
x (2n)WNmn ( m)
x (2 n + 1)WNmn
/2
m = n=0 n=0
2
N
n= -1
2
If we let /2 = Am = DFT of all even indexed samples
x (2n)WNmn
n=0
N
n= -1
2
and x (2n + 1)WNmn
/2 = Bm = DFT of all odd indexed samples
n=0
140Digital Power System Protection
N
m = -1
then we can write X ( m) m = 0 2 = A(m) + WNm B(m)
We take a look at the whole process in a step by step manner in Figures 7.5 and 7.6.
We can continue to follow the same mathematical logic to break down the 4-point DFT
into a combination of two 2-point DFTs and each 2-point DFT as a combination of two 1-point
DFTs. However, instead, we illustrate graphically, in Figures 7.7 and 7.8 how a 4-point DFT
is synthesised from two 2-point DFTs. Figure 7.9 shows the synthesis of a 2-point DFT using
two 1-point DFTs. Figure 7.10 depicts the concept of a butterfly. Finally Figure 7.11 combines
all these steps into a full-fledged flow-chart of an 8-point FFT with decimation-in-time.
FFT and Goertzel Algorithm141
It is to be noted that the original sequence of the indices which was in natural progression,
i.e.,
01234567
has become scrambled. It is now, 0 4 2 6 1 5 3 7.
Table 7.3 shows the binary representation of the indices.
Thus, it is seen that the sequence of indices for the FFT can be easily generated by
bit reversal of the naturally ordered indices. To take advantage of this property, DSP
FFT and Goertzel Algorithm145
microprocessors include facilities for bit reverse addressing. Note that bit-reversal is not the
same as bit complementation.
Decomposing the sum as that of first-half and second-half of the time domain samples:
N N
n= -1 n= -1 N
2 2 N - j 2p m 2 + n /N
X (m) mm == 0N -1 = x ( n) e - j 2p m ( n )/N
+ x + n e
n=0 n=0
2
e - j 2p m ( n) / N = WNmn
N
- j 2 p m + n /N
2
e = e - j 2 p mN/2 N e - j 2p mn/N
146Digital Power System Protection
N
- j 2p m + n / N
2
e = e - j 2p m/2 e - j 2p mn/N ; in the first term the constant 2 in numerator is not
cancelled with that in the denominator so that the term remains in standard form WN = ej2p/N.
Hence,
N
- j 2p m + n / N
2
e = W2mWNmn
N N
n= -1 n= -1
2 2 N
X (m) mm == 0N -1 = x (n)WNmn + W2m x + n WNmn
2
n=0 n=0
In order to introduce further simplification by way of recursion, consider the even indexed
terms of the DFT, i.e., X(2m):
N N
n= -1 n= -1
2 2 N
X (2m) mm == 0N / 2 -1 = x (n)WN(2 m ) n + W2(2 m ) x + n WN(2 m ) n
n=0 n=0
2
N N
n= -1 n= -1
2 2 N
X (2m) mm == 0N / 2 -1 = x (n)WNmn/ 2 + W2(2 m ) x + n WNmn/ 2
n=0 n=0
2
What is W2(2 m ) ?
- j 2p 2 m
W2(2 m ) = e 2 = e - j 2p m = 1
Hence, we get:
N N
n= -1 n= -1
2 2 N
X (2m) mm == 0N / 2 -1 = x (n)WNmn/ 2 + x + n WNmn/ 2
n=0 n=0
2
We can club the summation as:
N
n= -1
2 N mn
X (2m) mm == 0N/2 -1 = x (n) + x 2 + n WN/2
n=0
The above is an N/2 point DFT of sum of signals from the two halves, i.e., for an 8-point
DFT
[x(0) + x(4)], [x(1) + x(5)], [x(2) + x(6)] and [x(3) + x(7)]
But remember that the above expression gives only even indexed terms of the DFT. We
need to find expression for remaining half of the DFT terms, i.e., the odd indexed terms of
the DFT.
148Digital Power System Protection
What is WN(
2 m +1) n
?
WN(2 m +1) n = WN2 mnWNn = WNmn/ 2WNn
What is W2(2 m +1) ?
- j 2p 2 m - j 2p 1
W2(2 m +1) = W22 mW21 = e 2 e 2 = e - j 2p m e - jp = 1 ( -1) = -1
Hence, we get:
N N
n= -1 n= -1
2 2 N
X (2m + 1) mm == 0N / 2 -1 = x (n)WNmn/ 2WNn - x + n WNmn/ 2WNn
n=0 n=0
2
n=0 2
The above is a weighted N/2 point DFT of difference of signals from the two halves,
weighted by factor WNn , i.e., for an 8-point DFT:
[x(0) x(4)], [x(1) x(5)], [x(2) x(6)] and [x(3) x(7)]
Thus, on the input (time domain) side we have to club the first half of samples with the
second half of samples and on the output side (frequency domain) we get a sequence of even
indexed DFT terms followed by odd indexed DFT terms. If we compare this with decimation-
in-time FFT then by analogy we can call this as decimation-in-frequency FFT. Hence, the term
called decimation-in-frequency FFT.
The process of breaking down the bigger DFT into DFT of sum and difference of smaller
sequences is continued till there is only one sample left in the sequence. The idea can be
FFT and Goertzel Algorithm149
grasped intuitively by Figures 7.13 to 7.15. Finally Figure 7.16 shows the full-fledged signal
flow graph for the FFT with decimation-in-frequency.
The process of breaking down the DFT can be continued further till we get to one point
DFT which is the signal itself. This is shown for an 8-point FFT in Figures 7.14 and 7.15.
Figure 7.14 Breaking down a 4-point DFT into two number of 2-point DFTs.
150Digital Power System Protection
Figure 7.15 Further breaking down of 2-point DFT into 1-point DFTs.
clear, clc
x=[1 2 3 4 5 6 7 8];
XX = myFFTDIT(x)
YY = fft(x), % For checking results against MATLABs fft()
function
% the user created function myFFTDIT() is include here for ease in
documentation
% function [X] = myFFTDIT(x)
% Implementation of decimation-in-time FFT
% N = numel(x);
% if (N == 1 )
% X = x ;
% else
% xeven = x(1:2:N-1);
% xodd = x(2:2:N);
% WNm=exp( (-j * 2 * pi / N)*(0:(N/2)-1 ) );
FFT and Goertzel Algorithm151
%
% Am = myFFTDIT(xeven);
% Bm = myFFTDIT(xodd);
% X1 = Am + Bm .* WNm;
% X2 = Am - Bm .* WNm;
% X = [X1 X2]
%
%
% end
XX =
Columns 1 through 6
36.0000 4.0000 + 9.6569i 4.0000 + 4.0000i 4.0000 + 1.6569i 4.0000 4.0000 1.6569i
Columns 7 through 8
YY =
Columns 1 through 6
36.0000 4.0000 + 9.6569i 4.0000 + 4.0000i 4.0000 + 1.6569i 4.0000 4.0000 1.6569i
Columns 7 through 8
clear, clc
N = 8 ;
x =[ 1 2 3 4 5 6 7 8] ; % for testing
x
XX = myFFT(x);% FFT output in XX() is scrambled.
%-------------------- XXX() collects unscrambled FFT
XXX(1)=XX(1);% 0
XXX(5)=XX(2);% 4
XXX(3)=XX(3);% 2 Note:
XXX(7)=XX(4);% 6 Output of FFT
DIF is
XXX(2)=XX(5);% 1 scrambled.
XXX(6)=XX(6);% 5
XXX(4)=XX(7);% 3
XXX(8)=XX(8);% 7
XXX
YY = fft(x) % Compare the results with MATLABs fft() function
% X1 = myFFT(sum);
% X2 = myFFT(diffW);
%
% X=[X1 X2];
%
%
% end
36.0000 4.0000 + 9.6569i 4.0000 + 4.0000i 4.0000 + 1.6569i 4.0000 4.0000 1.6569i
Columns 7 through 8
YY = Columns 1 through 6
36.0000 4.0000 + 9.6569i 4.0000 + 4.0000i 4.0000 + 1.6569i 4.0000 4.0000 1.6569i
Columns 7 through 8
n =3 - j 2p m n
X (m) mm ==30 = x (n)e 4
n=0
- j 2p - j 2p mn
Let e 4 = W4 hence, e 4 = W4mn
Hence, we can write the 4-point DFT as:
n =3
X (m) mm ==30 = x (n)W4mn
n=0
Let us say, we are interested in finding only the fundamental component; X(1), then
m = 1 and we can write:
n =3
X (1) = x (n)W4n
n=0
expanding the summation:
X (1) = x (0)W40 + x (1)W41 + x (2)W42 + x (3)W43
since W40 = 1; we can write:
X (1) = x (0) + x (1)W41 + x (2)W42 + x (3)W43
X (1) = x (0) +W41 ( x (1) + x (2)W41 + x (3)W42 )
X (1) = x (0) +W41 ( x (1) + W41 ( x (2) + x (3)W41 ))
X (1) = x (0) +W41[ x (1) + W41 ( x (2) +W41 x (3))]
The algorithm can be implemented as shown in Figure 7.18.
Figure 7.18 Flow chart of Goertzel algorithm for finding fundamental component of DFT.
FFT and Goertzel Algorithm157
Figure 7.19 Flow chart of Goertzel algorithm for finding fundamental component of DFT.
We can also appreciate the Goertzel algorithm with the help of the signal flow graph
shown in Figure 7.20.
We can run through the signal flow graph by following entries as in Table 7.4.
Step Time x(n) Output of delay block Output of adder block Input to delay block
2 1 x(2) W41 x (3) x(2) + W41 x (3) W41 ( x (2) + W41 x (3))
3 2 x(1) W41 ( x (2) + W41 x (3)) x (1) + W41 ( x (2) + W41 x (3)) W41 ( x (1) + W41 ( x (2)
+W41 x (3)))
4 3 x(0) W41 ( x (1) + W41 ( x (2) x (0) + W41 ( x (1)
+W41 x (3))) +W41 ( x (2) + W41 x (3)))
(4-point fundamental DFT
complete)
X1=0;
x=[0 1 2 3 4 5 6 7] ;
for k=N:-1:1
X1= X1 * WN1 + x(k)
end
X1
%------- Lets check our result with MATLAB fft()
XX= fft(x);
XX1=XX(2) % First term of FFT is index 1 and not zero and is
the DC term, hence XX(2)
X1 = 7
FFT and Goertzel Algorithm159
X1 = 10.9497 4.9497i
X1 = 9.2426 11.2426i
X1 = 2.5858 14.4853i
X1 = 5.4142 12.0711i
X1 = 10.3640 4.7071i
X1 = 9.6569 + 4.0000i
X1 = 4.0000 + 9.6569i
Review Questions
161
162Digital Power System Protection
is an integral multiple of the signal frequency and the window length is adjusted to contain
integral number of cycles of the signal then the various frequency components get collected
in the correct bin. However, if the sampling frequency is not an integral multiple of the signal
frequency or if the window length is such that it does not contain integer number of cycles
of the signal, all the bins collect something, even if the signal has only single frequency
component. Thus, the frequency components are said to leak into those bins where they are
not supposed to show up.
8. Sampled data system has vision only at integral multiples of fsampling/N, i.e., at analysis
frequencies of mfsampling/N; where m is from 0 to N/2. Thus, we have vision only
at discrete points in the frequency domain, up to a maximum frequency of fsampling/2
at m = N/2. At all other points in the frequency domain, we remain blind.
9. Frequency response of DFT is zero at integral multiples of fundamental. And these
are the points at which the DSP system has vision.
10. Frequency response of DFT is non-zero at non-integral multiples of fundamental, but
at these points DSP does not have vision.
It is to be noted that we are performing the mathematical operation of multiplication
between signal and window in the time domain. Multiplication in time domain has effect of
convolution in frequency domain. Hence, spectrum of time windowed signal can be found by
convolution of spectra of signal and window function.
Spectrum of the rectangular function is the well-known sinc() function. Hence, the spectrum
of the rectangular windowed signal is obtained by convolution of sinc() function with spectrum
of the signal.
Figure 8.2 The sampled and windowed signal in time and frequency domains.
164Digital Power System Protection
Figure 8.2 shows that in spite of sampling and windowing, the spectrum of the signal
has remained intact! Careful observation of Figures 8.2 and 8.3 will reveal that this feat was
possible only because of the following:
1. Selection of the window size to match one complete cycle of the signal.
2. Selecting sampling frequency as an integral multiple of signal frequency.
The convolution operation between the window spectrum and signal spectrum can be
verified with the help of Figure 8.4.
It can now be easy to appreciate that in real life these conditions can seldom, if ever, be
satisfied. This raises the natural question: What will happen if:
1. a fraction of cycle of signal appears in the window?
2. sampling frequency is not an integral multiple of the signal frequency?
The answer to both above questions is that in such a circumstance there will be spectral
leakage. We take up spectral leakage in the next section.
Windowing and Spectral Leakage165
Figure 8.4 Discrete convolution between window spectra and signal spectra.
distortions have appeared at the window edges! This is bound to cause additional frequency
components in the frequency spectrum of the windowed signal. These additional spurious
frequency components are referred to as spectral leakage in the signal processing literature.
to zero at the edges of the window. If we look at the spectrum of a window function, we
will notice that in general there is an inverse relationship between width of main lobe and
height of side lobe, i.e., side lobe ripple. Ideally we would like the width of main lobe to be
zero and side lobe ripple to be zero as well, which makes it an impulse at 0 Hz. However,
such a window shape (impulse) in frequency domain corresponds to infinite window length
in time domain, obviously an impractical proposition. Thus, we have to live with ill effects
of windowing, and the best that we can do is to choose a window shape which will minimise
spectral leakage. Table 8.1 lists a few popularly used windows.
N -1 2p n
1. Hanning window = 0.5 - 0.5cos
n = 0 w ( n)
N
N -1 2p n
2. Hamming window = 0.54 - 0.46 cos
n = 0 w ( n)
N
n
N /2
=
n = 0 w ( n)
N/2
3. Triangular window
N -1 n
=2-
n=
N w ( n)
+1 N/2
2
Many other windows are reported in the literature, which are named after their inventors.
Some of these are: Bartlett, Blackmann, Chebyshev, Gaussian, Kaiser, Nuttal, Parzen and
Tuckey.
Figure 8.7 shows three windows in time and frequency domain. It can be observed that
in the frequency domain, the peak of rectangular window is shown normalised to 1. This is
known as processing gain or coherent gain of the window. This is in fact the DC signal gain
of the window. For all other windows, the gain is reduced because the window smoothly goes
to zero near the boundaries. The processing gain of the triangular window and the Hamming
window are 0.5 and 0.5328 respectively. This is a known bias on spectral amplitudes of the
windowed signal.
Review Questions
It can be proved that the computational process that goes on inside the box corresponds to
some kind of filtering operation like high-pass, low-pass, band-pass or band-stop filtering and
differentiating, integrating, convolution or co-relation. In fact every act of processing on time
domain samples creates some or the other effect in the frequency domain, be it intentional or
unintentional. We can also implement other types of digital filters which cannot be classified
into one of the basic four types listed above, such as adaptive filters used for noise cancelling
and beam steering. For example we can devise an adaptive digital filter for picking out an
EEG signal buried deep in 50 Hz line noise. We can design a filter for predicting tomorrows
share price index or the electrical load on the power system during the next hour and so on.
Keeping in mind that an antenna is a spatial filter, we can design it in such a way as to have a
null in the direction of an interfering signal, without, of course physically moving the antenna!
(Adaptive beam steering).
Thus, any mathematical operation is an act of filtering. Even if you do not do anything to
samples of a signal and just pass them from input to output, it is an all-pass filtering operation.
However, there are some drawbacks of digital filters. These are listed in Table 9.2.
The IIR filters can be naturally implemented using recursion because of the inherent
feedback. Hence, they are sometimes referred to in the literature as recursive filters. There
is no apparent recursion in the FIR filters and are naturally implemented using non-recursive
procedures. Hence, they are referred to in the literature as non-recursive filter. However, it can
be shown that both FIR as well as IIR filters can be expressed recursively or non-recursively.
Hence, it is better to classify them as FIR and IIR instead of recursive and no-recursive.
FIR filters are the simplest of all. Hence, we first take up the explanation of FIR filters.
In FIR filters as shown in Table 9.3, the present sample of output is produced by linearly
combining the present input with a finite number of past input samples. FIR filters do not use
any of the past output samples, they use past input samples only.
Figure 9.2 shows a block diagram representation of FIR and IIR filters. The block diagram
also suggests the hardware that will be required in case the filter is realised in hardware.
Figure 9.2 Block diagrams of FIR and IIR filters suggesting hardware implementation.
Introduction to Digital Filtering177
Thus, the block-diagram serves a dual purpose of a flowchart as well as hardware diagram of
the digital filter! It can be seen that multiplier, adder, delay (MAD) and memory are the basic
elements from which any digital filter of the FIR or IIR type can be constructed.
The outputs of the FIR and IIR filters, as can be seen from the block diagram, will be
given by:
FIR y(n) = b(0) x(n) + b(1) x(n 1) + b(2) x(n 2) + + b(k) x(n k)
IIR y(n) = b(0) x(n) + b(1) x(n 1) + b(2) x(n 2) + + b(Q) x(n Q)
a(1) y(n 1) a(2)y(n 2) a(P) y(n P)
We will answer all the above questions as we progress. Let us write the equation for y(n)
as follows:
y(n) = h(0) x(n) + h(1) x(n 1) + h(2) x(n 2) + h(3) x(n 3) + h(4) x(n 4)
Which can be written in the compact form as:
k =4
y ( n) = h( k ) x ( n - k )
k =0
Does the above equation look familiar? It surely does, as we can easily see that it is the
operation of convolution. It is a convolution between two sequences h( ) and x( ). This is a
very important result because it gives us a procedure for implementing digital filter and provides
an insight into its working.
Input Output
Time x(n) x(n 1) x(n 2) x(n 3) x(n 4) y(n)
0 1 0 0 0 0 1/5 = h(0)
of impulse
Ts 0 1 0 0 0 1/5 = h(1)
Duration
response
impulse lasts for a finite duration after which it becomes zero. Hence, it is named as finite
(time) impulse response or FIR filter. We can also see that the impulse causes the coefficients
of the filter to be outputted one-by-one. Since the coefficients are all that can be known about
a filter, the impulse response reveals everything about the filter. The impulse response is
sketched in Figure 9.5.
1 fsignal
j 2 p fsignal n j 2p n
j 2p fsignal nTsampling fsampling fsampling j 2p ( f pu ) n jw pu n
x ( n) = e =e =e =e =e
Dropping the subscript pu from wpu but keeping in mind that w is per unit radian frequency,
we can write the expression for x(n) as:
x (n) = e jw n
where, w is the pu radian frequency of the signal.
Since output is the average of present sample and previous samples, so we can write:
1 jw n 1 jw ( n -1)
y ( n) = e + e
2 2
1 jw n 1 jw n 1 - jw
y ( n) = e + e e
2 2 2
1
y(n) = (1 + e - jw )e jw n , but e jw n = x (n)
2
Hence, the frequency versus magnitude response is a cosine curve and phase versus
frequency is linear, with negative slope, in the fourth quadrant as shown in Figure 9.6. Thus,
phase shift is always negative as the output of the filter always lags the input because of the
delay inherent in the filter operation.
Figure 9.6 Frequency response (magnitude and phase response) of 2-point running average filter.
For the 2-point running average filter, we can find the output by running through the
sequence and finding the average of every two consecutive points.
(1 + 1) (1 + 1) (1 + 1)
y(n) = * *
2 2 2
y(n) = [* 1 1 1 *]
x (n) = e jp n = [1 -1 1 -1 1]
Again, calculating the 2-point running average of the sequence, we get another sequence:
182Digital Power System Protection
(1 - 1) (1 - 1) (1 - 1)
y(n) = * *
2 2 2
Hence, we get:
y(n) = [* 0 0 0 *]
Thus, the signal at the highest (per unit) frequency is totally attenuated to zero.
Substituting z = ejw
1
H ( jw ) = (1 + e - jw )
2
1
H ( jw ) = (e + jw /2 e - jw /2 + e - jw /2 e - jw /2 )
2
1
H ( jw ) = (e + jw /2 +e - jw /2 )e - jw /2
2
e + jw /2 +e - jw /2 - jw /2
H ( jw ) = e
2
w
H ( jw ) = cos e - jw /2
2
w
H (w ) = cos is the magnitude response and H ( jw ) = e - jw /2 is the phase angle response.
2
These responses are exactly same as those derived without using Z-transform.
Introduction to Digital Filtering183
Let x (n) = e jw n
Hence, we can write y(n) as:
1 1 1
y(n) = e jw n + e jw ( n -1) + e jw ( n - 2)
3 3 3
1 1 1
y(n) = e jw n + e jw n e - jw + e jw n e - jw 2
3 3 3
1
y(n) = [1 + e - jw + e - j 2w ] e jw n
3
but e jw n = x (n); hence we can write:
1
y(n) = [1 + e - jw + e - j 2w ]x (n)
3
y ( n) 1
= H ( jw ) = 1 + e - jw + e - j 2w
x (n) 3
Recall that for a 2-point running average filter the cut-off frequency was p which can be
written as 2p/2. One can expect that the cut-off frequency for a 3-point running average filter
will be 2p/3. Hence, substituting w = 2p/3, we get
1
H ( j) = 3 + 4 cos(w ) + 2 cos(2w )
3
2p 1 2p 2p
H j = 3 + 4 cos + 2 cos(2)
3 3 3 3
2p 1 -1 -1
H j = 3+4 +2
3 3 2 2
2p 1
H j = 3 - 2 -1 = 0
3 3
To find H(jw):
(sin(w ) + sin(2w ))
H ( jw ) = - tan -1
(1 + cos(2w ) + cos(2w ))
(sin(w )
H (jw ) = - tan -1
cos(w )
H (jw ) = -w
Hence,
2p 2p
H j =-
3 3
Introduction to Digital Filtering185
The frequency response of the 3-point running average filter is shown in Figure 9.7. Note
that at DC, i.e., w = 0;
1 1 1
H ( jw ) = 3 + 4 cos(w ) + 2 cos(2w ) = 3 + 4 .1 + 2 .1 = 9 =1
3 3 3
Substituting z = ejw
1
H (w ) = (1 + e - jw + e - j 2w )
3
jw
Taking out common factor e , we get:
1
H (w ) = (e jw + 1 + e - jw )e - jw
3
1 e jw + e - jw - jw
H (w ) = 1 + 2 e
3 2
1
H (w ) = (1 + 2cos(w ))e - jw = Magnitude angle
3
1
Thus, magnitude of output = (1 + 2cos(w )) and phase angle = w.
3
However, (1 + 2cos(w)) cannot be taken as magnitude since as cos(w) becomes 1, the
1
term becomes negative. Magnitude of (1 + 2cos(w )) can be written as:
3
1
H (w ) = (1 + 2cos(w ))2
3
1
H (w ) = (1 + 4 cos(w ) + 4(cos(w ))2
3
1
H (w ) = (1 + 4 cos(w ) + (2 cos(2w ) + 2))
3
H (w ) =
1
((3 + 4 cos(w ) + (2 cos(2w )))
3
and H (w ) = -w
Note that we have used the following:
1 + 2cos(w ) = (1 + 2cos(w ))2
Let us find the frequency at which |H(w)| is zero by solving the following:
1 + 2 cos(w ) = 0
2 cos(w ) = -1
1
cos(w ) = -
2
1
w = cos -1 -
2
2p
w= = 120 ( deg.)
3
This is the same as derived earlier without using Z-transform.
Introduction to Digital Filtering187
9.7 F
requency Response of 4-Point Running Average Filter Using
Z-transform
Output of the 4-point running average filter will be given by average of present sample taken
with the last three samples:
1 1 1 1
y ( n) = x (n) + x (n - 1) + x (n - 2) + x (n - 3)
4 4 4 4
1 1 1 1
Y (z) = X ( z ) + z -1 X ( z ) + z -2 X ( z ) + z -3 X ( z )
4 4 4 4
1
Y (z) = (1 + z -1 + z -2 + z -3 ) X ( z )
4
Y (z) 1
H (z) = = (1 + z -1 + z -2 + z -3 )
X (z) 4
Y (z) 1
H ( jw ) = = (1 + e - jw + e - j 2w + e - j 3w )
X (z) 4
Y (z) 1
H ( jw ) = = (1 + cos(w ) - j sin(w ) + cos(2w ) - j sin(2w ) + cos(3w ) - j sin(3w ))
X (z) 4
Y (z) 1
H ( jw ) = = [{1 + cos(w ) + cos(2w ) + cos(3w )} - j{sin(w ) + sin(2w ) +sin(3w )}]
X (z) 4
1
H ( jw ) = (1 + cos(w ) + cos(2w ) + cos(3w ))2 + (sin(w ) + sin(2w ) +sin(3w ))2
4
As the expressions for magnitude and phase response of the 4-point running average digital
filter is somewhat involved and not easy to sketch manually, let us write a small MATLAB
script to study the magnitude and phase response as shown in Figure 9.8. The program is
shown below:
188Digital Power System Protection
% --------------------------------------------------------------
% ----- Freq Response of 4-Point Running Average Digital Filter
% --------------------------------------------------------------
clear , clc,clf
k=0;
step=pi/100;
for w = 0 : step : pi
k = k+1;
freq(k) = w;
x = 1 + cos(w)+ cos(2*w)+ cos(3*w);
y = sin(w) + sin(2*w) + sin(3*w);
H(k) =sqrt(x*x + y*y);
Ang_H(k) = atan2(-y, x);
end
subplot(2,1,1)
plot(freq,H)
title(Magnitude response of 4-point running average digital filter)
xlabel(Frequency , radians , w,max,pu=pi)
subplot(2,1,2)
plot(freq,Ang_H)
title(Phase response of 4-point running average digital filter)
xlabel(Frequency , radians , w, max , pu = pi)
In the next section, we write a MATLAB program for plotting the frequency response of
an N-point running average filter.
Introduction to Digital Filtering189
clear,clc,clf
%--------------------------------------------------------
%------ Enter the number of points in the R.A. filter ---
N = 10;
%--------------------------------------------------------
x= 1
y=0
w=1;
for n=1:N-1
s=num2str(n);
ss=strcat(+cos(,s,,*,w,))
sss= strcat(+sin(,s,,*,w,))
x = [ x strcat(ss)]
y = [ y strcat(sss)]
end
x
y
step = pi/100;
k=0;
for w=0:step:pi
k =k+1;
ang(k)=w;
190Digital Power System Protection
xx=eval(x);
yy=eval(y);
magH(k)=(sqrt(xx*xx+yy*yy))/(N);
phH(k)=atan2(-yy,xx)
end
ang_deg=(180/pi).*ang;
xmin=0;xmax=180;ymin=0;ymax=1;
xmin1=0;xmax1=180;ymin1=-pi;ymax1=pi;
subplot(2,1,1)
plot(ang_deg,magH)
axis([xmin xmax ymin ymax])
title(Magnitude Response of Running Average Filter)
xlabel( Freuency in pu deg ( max=pi rad))
subplot(2,1,2)
plot(ang_deg,phH)
axis([xmin1 xmax1 ymin1 ymax1])
title(Phase Response of Running Average Filter)
xlabel( Freuency in pu deg (max =pi rad) )
x;
y;
mag_of_H =[|H(jw)| =sqrt ( x ) ^2 + ( y )
^2];
ph_of_H =[ ang(H(jw) =atan2 (,-,( y ,),/,x ,)]
formula1=num2str(mag_of_H)
formula2=num2str(ph_of_H)
fprintf( |H(jw)| = %s \n ,formula1 )
fprintf( Angle(H(jw)) = %s \n,formula2 )
text(1,10,formula1)
text(1,0.8,formula2)
text(1,11,No. of points=)
text(40,11,num2str(N))
% subplot(2,1,2)
% text(1,
Figure 9.9 Frequency response of 10-point running average filter.
The response of the IIR filter to the impulse will be of the form:
[from 0 0 0 0 y1 y2 y3 y4 y5 to ]
where y1, y2, are all finite values. The basic cause for continuation of the impulse response
for infinite time is the feedback that exists in the IIR filter. This feedback causes the signal to
get recycled again and again even after the input has vanished. Note also that the basic reason
for instability of IIR filters also stems from the feedback. FIR filters do not have feedback.
Hence, impulse response of FIR filter lasts for a finite time and there is no possibility of the
FIR filter becoming unstable.
We repeat the time domain expression for output of the IIR filter below, for convenience.
Q P
y ( n) = bk x (n - k ) - ak y(n - k )
k =0 k =1
The block-diagram of Figure 9.10 can be used to directly implement an IIR filter in either
hardware or software form. Hence, this method of implementation or the topology of the IIR
filter is called direct form IIR filter.
Taking the Z-transform of the above equation, gives us:
Q P
Y (z) = bk z - k X (z) - ak z - k Y (z)
k =0 k =1
Rearranging further gives:
P Q
Y ( z ) + Y ( z ) ak z - k = X ( z ) bk z - k
k =1 k =0
P Q
Y ( z ) 1 + ak z - k = X ( z ) bk z - k
k =1 k =0
Introduction to Digital Filtering193
Thus, we get the expression for H(z), the transfer function of the filter in frequency domain:
Q
Y (z)
bk z - k
= H ( z ) = k = 0P
X (z) -k
1 + ak z
k =1
b0 + b1 z -1 + b2 z -2 + b3 z -3 + + bQ z -Q
H (z) =
(1 + a1 z -1 + a2 z -2 + a3 z -3 + + aP z - P )
Polynomial of order Q in z
H (z) =
Polynomial of order P in z
Factored polynomial of order Q in z
H (z) =
Factored polynomial of order P in z
( z - q1 )( z - q2 )( z - q3 ) ( z - qQ )
H (z) =
( z - p1 )( z - p2 )( z - p3 ) ( z - pP )
Note that we can express the transfer function in various ways. For example we can write:
H(z) = H1(z).H2(z)Hk(z)
The above form can be implemented using a cascade of lower order transfer functions.
Alternatively we can express H(z) as:
k l
H ( z ) = 1 + H I ( z ) + H II ( z )
i =1 i =1
where HI(z) are k first order transfer functions and HII(z) are l second order transfer functions.
This type of representation leads to parallel form of IIR filter realisation. A question that needs
to be answered at this stage is: What is the motivation for alternate forms of implementation
when we can implement the filter easily in the direct form? The answer to this question is
not intuitively obvious. The answer lies in the error analysis of various topologies and their
sensitivity to quantization of coefficient values. It can be shown that direct form is not always
the best, from the above consideration and cascade and parallel form prove to be better under
certain situations. Hence, large numbers of topologies have been developed in the literature
on digital filtering.
Review Questions
1.
FIR stands for
(a) Finite magnitude impulse response
(b) Finite duration impulse response
(c) Finite duration, finite magnitude impulse response
(d) Finite magnitude, infinite duration impulse response.
194Digital Power System Protection
2.
FIR filter makes use of
(a) Samples of past input only
(b) Samples of past output only
(c) Samples of past input as well as samples of past output
(d) Samples of past input along with present sample of input.
3.
IIR filter makes use of
(a) Samples of past input only
(b) Samples of past output only
(c) Present sample, samples of past input as well as samples of past output
(d) Samples of past input along with present sample of input.
4.
What are the basic hardware elements required for implementation of FIR or IIR filters?
5.
Discuss why FIR has finite-time and IIR has infinite-time impulse response.
6.
Why IIR filters are referred to as recursive filters?
7.
Why FIR filters are called non-recursive?
8. Show how an FIR filter can be implemented recursively.
9. FIR filters are referred to as filters with only zeros and a Qth order pole at the origin.
Justify/Explain.
10. IIR filters have both poles as well as zeros. Explain.
11. FIR filters are free from possibility of becoming unstable. Explain.
12. Explain why IIR filters need to be carefully designed in order to avoid instability.
13.
Explain intuitively how running average can have low-pass filtering effect.
14.
By corollary, which mathematical operation will cause high-pass filtering?
15.
What will be the effect of increasing the number of points in the running average?
16.
Analyse a 4-point running average filter and find is cut-off frequency and phase shift
at cut-off frequency.
17.
Perform the analysis of 2-point and 3-point running average filters using Z-transform.
18.
State whether the running average filter is of FIR type or IIR type.
19.
Explain the concept of per unit frequency.
20.
What is the maximum per unit radian frequency allowable in digital signal processing?
21.
What is the maximum per unit hertz frequency allowable in digital signal processing?
22.
What happens if a (pu) frequency greater than p appears at the input of a digital
processing system?
23.
Write the equation in time domain for output of a 2-point running average filter.
24.
Write the equation in time domain for output of a 3-point running average filter.
Introduction to Digital Filtering195
25.
Write the equation in time domain for output of a n-point running average filter.
26.
Convert above equations into Z-domain equations.
27.
Convert above equations into transfer functions of respective filters.
28.
Write a MATLAB script to plot the magnitude and frequency response of an n-point
running average filter.
10
Digital Filter Design
10.2 R
elation between Linear Phase and Symmetry of Filter
Coefficients for FIR Filter [52]
Before we delve into design procedure for FIR filters, the following points need to be kept
in mind:
196
Digital Filter Design197
Figure 10.3 Impulse response symmetry patterns for type I, II, III and IV FIR filters.
Digital Filter Design199
Note: The differences between DFT and IDFT formulae are as follows:
1. 1/N factor is included in IDFT formula. No such factor is in DFT formula.
2. Sign of exponent is positive in IDFT and negative in DFT formula.
N -1 j 2p ik jp mi
1
h( k ) =
N
Ar ( fi )e N
i =0
N -1 j 2p i ( k 0.5 m )
1
h( k ) =
N
Ar ( fi )e N
i =0
N -1
1 2p i(k - 0.5 m) 2p i(k - 0.5 m)
h( k ) = Ar ( fi ) cos + j sin
N i =0 N N
We are designing a filter with all h(k) values real. Hence, the j sin() terms cancel out, giving:
N -1
1 2 p i (k - 0.5m)
h( k ) =
N
Ar ( fi ) cos N
i =0
Using the symmetry property of DFT, i.e., the contributions of the ith and (Ni)th terms are
the same. Thus, we can sum half of the terms and double the result. Further recalling that b(k)
= h(k) for FIR filters and noting that m the order of the filter is given by:
m = N 1; where N is the number of time/frequency domain samples
N = m + 1; where m is the order of the filter
We get the final expression for filter coefficients as:
m
floor
k = m -1 Ar (0) 2 2
2 p i ( k - 0.5m )
k = 0 [ bk ] = +
(m + 1) (m + 1)
Ar ( fi ) cos
m +1
i =1
Note that for type-I filter, the order of the filter m is even, in which case floor (m/2) = m/2.
But for type-II filter, m is odd, hence m/2 will not be an integer but the floor() function of
MATLAB will generate an integer, thus allowing the summation.
10.3.1 N
umerical Problem on Design of FIR Filter Using Frequency Sampling
Method
Let us illustrate the design of FIR filter using the frequency sampling method with the help
of a numerical problem.
Example 10.1 Design a low-pass FIR filter to mimic an ideal low-pass filter with a cut-off
frequency of fsampling/4. Assume order of the filter to be 20.
Digital Filter Design201
Solution Let us design a type-I filter. We have order of the filter m = 20. We have m = N
1; where N = number of frequency domain samples. Hence N = m + 1 = 20 + 1 = 21. Thus,
number of frequency domain samples = 21. Let us first sketch the desired frequency response
of the filter as shown in Figure 10.4.
Figure 10.5 FIR filter block diagram.
m
floor
2
A (0) 2 2 p i (k - 0.5m)
[b(k )] = r +
(m + 1) (m + 1)
Ar ( fi ) cos
m +1
i =1
From the sketch of the desired frequency response we can see that:
i = 5 i = 20
Ar ( fi ) = 1 and Ar ( fi ) =0
i = 0 i=6
Further, we have m = 20, hence 0.5m = 10. Thus, we can write the above expression as:
1 2 10 2 p i(k - 10)
b( k ) = + Ar ( fi )cos
21 21 i =1 21
Let us write a small MATLAB script to compute the filter coefficients and verify that the
coefficients indeed implement the desired filter by finding the FFT of the coeffiients:
clear , clc
%--------------------------------------------
%---------Order of filter -----------------
m=20;
%-------------------------------------------
for k = 0 : m
fpu(k+1)= k/m;
b(k+1)=(1/(m+1));
for ii =1:5
b(k+1)=b(k+1) + (2/(m+1))* cos( 2 *pi*ii*(k-10)/(m+1)) ;
end
Digital Filter Design203
sT 2
1- + -
2 2!
If we neglect higher order terms, we can write:
1 + (sT/2) 2 + sT
z@ =
1 - (sT/2) 2 - sT
Hence, we can write the following transformation as:
2 ( z - 1)
s=
T ( z + 1)
The relation between analog and digital frequencies can be derived as follows:
Noting that: s = jW and z = ejw
We have the bilinear transform given by:
2 ( z - 1) 2 (1 - z -1 )
s= =
T ( z + 1) T (1 + z -1 )
Hence, we can write:
2 (1 - e - jw )
s = jW =
T (1 + e - jw )
which can be simplified as follows:
2 e - jw /2 (e jw /2 - e - jw /2 )
jW =
T e - jw /2 (e jw /2 + e - jw /2 )
2 (e jw /2 - e - jw /2 )
jW =
T (e jw /2 + e - jw /2 )
2 (e jw /2 - e - jw /2 ) 1
W= jw /2
T 2j (e + e - jw /2 ) / 2
w
sin
2 2
W=
T w
cos
2
2 w
W= tan
T 2
WT
w = 2 tan -1
2
2 w
W max = tan max
T 2
As digital frequency w reaches its maximum value, i.e., wmax = p, it gets translated as:
2 p
W max = tan = infinity
T 2
Thus, the entire range of analog frequencies from 0 to gets compressed into digital
frequencies from 0 to p; giving rise to the warping effect as shown in Figure 10.7.
206Digital Power System Protection
2 ( z - 1) 2 (1 - z -1 )
s= =
T ( z + 1) T (1 + z -1 )
Digital Filter Design207
Table 10.2
2 ( z - 1) 2 (1 - z -1 )
Bilinear transform s= =
T ( z + 1) T (1 + z -1 )
WT
W
Analog frequency to digital frequency
w w = 2 tan -1
2
2 w
w
Digital frequency to analog frequency
W W= tan
T 2
Pre-warping 2 W
Pre-warping W tan
T 2
2 0.25p 2 p 2 0.828
Wc ( pre-warped ) = tan = tan = 0.414 =
T 2 T 8 T T
Thus, the transfer function of the prototype filter which was H a (s) = 1(1 + s/Wc ) becomes:
1 1
H a (s) = =
s sT
1+ 1+
0.828 0.828
T
208Digital Power System Protection
2 z 1
Applying bilinear transformation, i.e., substituting, s =
T z +1
1
H (z) = ; which can be simplified as shown in the following steps:
2 ( z - 1) T
1+
T ( z + 1) 0.828
1
H (z) =
2( z - 1)
1+
0.828( z + 1)
1
=
0.828( z + 1) + 2( z - 1)
0.828( z + 1)
0.828( z + 1)
=
0.828z + 0.828 + 2 z - 2
0.828( z + 1)
=
2.828z - 1.172
0.828( z + 1)
=
2.828( z - 1.172 / 2.828)
0.2927( z + 1)
=
( z - 0.414)
0.2927(1 + z -1 ) z
=
(1 - 0.414 z -1 ) z
0.2927(1 + z -1 )
=
(1 - 0.414 z -1 )
Y ( z ) 0.2927(1 + z -1 )
H (z) = =
X ( z ) (1 - 0.414 z -1 )
Thus, we can transform the above z domain equation into a time domain difference
equation as shown below:
y(n) - 0.414 y(n - 1) = 0.2927 x (n) + 0.2927 x (n - 1)
y(n) = 0.414 y(n - 1) + 0.2927 x (n) + 0.2927 x (n - 1)
The above time domain equation can be directly implemented with the help of multiplier,
adder and delay blocks, as shown in Figure 10.8.
j 2p f pu
z = re
jw pu
z = re
Normally, the subscript pu is dropped from wpu. But we should keep in mind that w is per
unit radian frequency whose maximum value is p. Thus,
z = r ejw
Remembering that w is an angle, say q which varies from 0 to p.
z = r ejq
z = r e jq (angle from 0 to p)
The frequency response is found by evaluating the transfer function H(z) on the unit circle,
i.e., |z| = 1 or r = 1. Hence, z = ejq. Thus, we see that W which is a linear distance (in
the vertical direction) on the s plane, becomes an angle q on the z plane as shown in
Figure 10.9. Thus, as frequency tends to infinity in the s domain, the angle tends to p in
the z domain as shown in Table 10.3.
Mapping of z plane from s plane from Figure 10.9 is shown in Table 10.4.
Table 10.4 Mapping of z plane from s plane
s plane z plane
Left half Inside unit circle
+/ j w axis Circumference of unit circle
Zero frequency z=1
Infinite frequency z = +/ p or 0.5 fs
Observation of the desired frequency response shown in Figure 10.10 leads to the following
conclusions:
1. The transfer function must contain a pole.
2. The pole should be at frequency F0 q0.
3. The pole should be located inside unit circle on the z plane, i.e., r < 1 for stability.
Now, the mapping is as follows:
W = 0 to infinity analoge frequency (Hz)
w = 0 to p digital frequency (pu)
f = 0 to fs /2 signal frequency (Hz)
Hence, we can write:
F0
q0 = p
f s /2
2p F0
q0 =
fs
The z domain transfer function must have a pole at z = rejq 0, where r < 1
212Digital Power System Protection
Now, since we want real coefficients, poles must always appear in conjugate pairs. Hence,
we can now directly write the expression for the z domain transfer function as:
Numerator
H (z) =
( z - re jq0 )( z - re - jq0 )
b0 ( z 2 - 1)
H (z) =
z 2 - z r (e jq0 + e - jq0 ) + r 2
b0 ( z 2 - 1)
H (z) =
z 2 - 2 z r (cos(q 0 )) + r 2
If we design the gain to be unity at w = w0, then it will be automatically less than 1 at all
other frequencies. The gain has already been designed to be zero at w = 0 and w = p.
Imposing this constraint on H(z), we can write:
b0 e2 jq0 - 1
H ( z ) z = e jq = 1 =
e2 jq0 - 2e jq0 r cos (q 0 ) + r 2
which gives b0 as:
e2 jq0 - 2re jq0 cos(q 0 ) + r 2
b0 =
e2 jq0 - 1
The value of r should be nearly equal to 1 and on the lower side. The empirical formula
for r is given by:
DF p
r =1-
fs
where DF is change in frequency from F0 at which the gain is 3 dB down.
Digital Filter Design213
b0 ( z 2 - 1)
H (z) =
z 2 - 2 z r (cos(q 0 )) + r 2
b0 e2 jq0 - 1
H ( z ) z = e jq = 1 =
e2 jq0 - 2e jq0 r cos(q 0 ) + r 2
Example 10.3 Design a digital resonator for a frequency of 200 Hz. Assume sampling
frequency to be 1200 Hz and DF = 6 Hz.
Solution We have:
2p F0 200 2p p
q0 = = 2p = =
fs 1200 6 3
We have
DF p 6p
r =1- =1- = 0.9843
fs 1200
0.0156( z 2 - 1)
=
z 2 - 2 z (0.9843)(cos(p / 3)) + (0.9843)2
214Digital Power System Protection
0.0156( z 2 - 1)
H (z) =
z 2 - 2 z (0.9843)(0.5) + 0.9688
0.0156( z 2 - 1)
=
z 2 - (0.9843) z +0.9688
0.0156(1 - z -2 )
=
1 - (0.9843) z -1 + 0.9688z -2
Y (z) 0.0156(1 - z -2 )
= =
X ( z ) 1 - (0.9843) z -1 + 0.9688z -2
Y ( z ) - (0.9843) z -1Y ( z ) + 0.9688z -2Y ( z ) = 0.0156 X ( z ) - 0.0156 z -2 X ( z )
The above difference equation leads to the direct form implementation as shown in
Figure 10.11.
Review Questions
1.
Sketch the ideal low-pass, high-pass, band-pass and band-stop filter magnitude
responses.
2.
Why phase response of digital filters is always negative?
3.
What are the advantages of having linear phase response?
4.
What constraints get imposed on the FIR filter impulse response (filter coefficients)
when it is required to have linear phase response?
5.
Describe the frequency sampling method of FIR filter design.
6.
Compare FIR and IIR filters.
7.
Which type of filter (FIR/IIR) is more suitable for numerical protective relaying
application?
8.
Design a low-pass FIR filter with a cut-off frequency of fsampling/6. Assume order of
the filter to be equal to 32.
9.
In IIR filter, state whether the magnitude of the impulse response is infinite or its
duration is infinite.
10.
What is the basic reason for infinite impulse response?
11.
Write expression for direct time domain implementation of the IIR filter.
12.
Write the frequency domain transfer function, H(z) of the IIR filter.
13.
What are the various alternate ways in which H(z) can be expressed?
14.
What is the motivation for implementing the filter in other than direct form?
15.
Derive bilinear transformation.
16.
What is the range of digital frequencies in Hz and in (pu) radians?
17.
What do you mean by warping of frequency?
18.
How frequency warping can be pre-compensated by pre-warping?
19.
State transfer functions of low-pass Butterworth, Chebyshev and elliptical filters in the
s domain.
20.
Numerical problems:
(i) Design an IIR digital filter based on first order low-pass Butterworth filter prototype
with 3 dB cut-off frequency wc = 0.1p ; by using bilinear transformation.
(ii) Design a digital resonator for a frequency of 100 Hz. Assume sampling frequency
to be 1500 Hz and DF = 4 Hz.
11
Synchrophasors
11.1Introduction [1,4,7,9,12,18,32,35,36,37]
Synchrophasors are marching inexorably into the field of power systems and seem to
have caught the imagination of the protection engineers. They have become the mainstay
of the so called smart-grid and penetrated considerably into the existing power system
transmission network. They have made feasible the concept of wide area protection, control
and monitoring. It is said that If SCADA system is akin to X-Ray then synchrophasor is
the MRI! One single factor which has made the synchrophasor revolution possible is the
affordability GPS technology.
Synchrophasors combine the concept of phasor with that of synchronisation in time. The
phasors were conceptualised by Charles Proteus Steinmetz in 1893. In 1988 synchrophasor
measurement units were invented by Arun G. Phadke and James S. Thorp at Virginia Tech
University in the U.S.A. Thus, Steinmetzs technique of phasor calculation has evolved
into the calculation of real time phasors that are synchronised to an absolute time reference
distributed by the global positioning system (GPS). Early prototypes of the PMU were
built at Virginia Tech. Macrodyne Corporation (U.S.A.) built the first commercial phasor
measurement unit (PMU) in 1992 almost a hundred years after Steinmetzs conceptualisation
of phasor. The concept of synchrophasor was standardised by IEEE in 1995 as IEEE Std
13341995 which was reaffirmed in 2001. A working group was established in 2001 to
update the standard. The new synchrophasor standard was issued in 2005 as IEEE Std
C37.1182005. This was modified in 2011 and broken up into two standards C37.118.1
which specified the measurement part of the PMU and C37.118.2 which specified the
communication aspects. The new standard now includes specifications of PMU output under
transient conditions, about which its previous version was silent. In the following sections
we will first visit the concept of phasor and then develop the idea of synchrophasor as
defined by the IEEE.
216
Synchrophasors217
11.2The Phasor
An electrical signal represented by the equation:
vx = 2 VRMS cos (2p f t + f ) where w = 2 p f
or vx = 2 VRMS cos (w t + f )
vx = Re( 2 VRMS e j (w t + f ) )
vx = Re{ 2(VRMS e jf )e j w t }
vx
= Re 2
VRMS e jf e j w t
Instantaneous value Phasor
Thus, while defining the phasor, we drop the term ejwt in the above expression and consider
only the RMS value and phase angle to express the phasorV as:
V = VRMS e jf
with the tacit understanding that frequency w is constant. In phasor calculations we replace
the signal vx by the phasorV . Thus, we have made the reversible phasor transformation denoted
by P( ).
P (vx ) V
Yields
The reason, for popularity of phasors, is that in the steady state, it is very convenient to
deal with the signal represented in the phasor form rather than in instantaneous form. We
should note, therefore, that phasor is inherently a steady state concept. Conventional phasor
representation is valid for single frequency undistorted sinusoids of constant amplitude, i.e.,
for stationary signals. In a power system, immediately after occurrence of faults or switching
events or during a power swing, signals are not stationary. During these times, it is difficult to
apply classical concept of phasor. However, the recent modification to the IEEE synchrophasor
standard for synchrophasor C37.118.1, provides definition of the synchrophasor under non-
stationary conditions, i.e., when both its amplitude and phase angle are being modulated by a
low frequency signal.
Figure 11.1 shows, the snap-shot of a length Vm rotating at w rad/s. As the length Vm rotates,
its instant to instant projection along the horizontal axis is the instantaneous value of the phasor.
If the length Vm happened to be aligned with the horizonal axis at t = 0 then we associate
218Digital Power System Protection
a phase angle of 0 with it. Thus, expressing the phasor as V1 = (Vm / 2 )0 . Suppose at time
t = 0 the rotating length Vm happened to be at an angle f with the horizontal axis, then we
associate a phase angle of f with it, thus, expressing the phasor as V2 = (Vm / 2)f . Hence,
the instantaneous value, the peak value, the RMS value and the phase angle are all linked
with the same entity, i.e., the sinusoidal waveform, the phasor captures all this information in
a compact form.
Figure 11.2 shows the reference cosine wave, given by vx = Vm cos(w t ), i.e., a cosine wave
with a phase angle of zero, with which any given 50 Hz wave will be compared for its
synchrophasor representation. As defined in the IEEE Std C37.118-2011, the time count starts
from t = 0 instant, which is defined as 00 hr : 00 min : 00 sec on January 1st, 1970. Precise
time-keeping has been maintained from this instant with the help of highly precise 1 pulse per
Synchrophasors219
second signal. At the start of every second (which can be traced back to this t = 0 instant),
a signal known as the 1 PPS signal is broadcasted from the geographical positioning system
satellites, and is available throughout the power system over widely separated geographical
locations, indeed throughout the globe. A small fist sized GPS antenna connected to a receiver
is used to pick up this signal and decode it. Figure 11.2 also shows the phasor diagram for the
reference wave, which will be expressed mathematically as:
V V
Vref = m 0 = m e j 0
2 2
Figure 11.3 shows a signal with amplitude of Vm and phase angle of f leading (i.e.,
positive phase angle). It can be seen that the phase angle of the measured waveform is given
by measuring the occurrence of its peak with respect the rising edge of the 1 PPS signal.
Figure 11.3 Measured signal with peak value Vm and phase angle f.
220Digital Power System Protection
Note that {w(t)t} represents the phase which continuously changes with time while f o is the
fixed phase offset.
v(t ) = Vm (t ) cos(f (t ) + fo )
We have
d
w= (f (t ))
dt
d (f (t ))
2p f =
dt
d (f (t )) = 2p f dt
d (f (t )) = 2p f dt
Phase change with time = 2p f dt
Thus, when the input frequency is varying, we can write the input signal as:
(
v(t ) = Vm (t ) cos 2 p f dt + fo )
v(t ) = Vm (t ) cos (2 p ( g(t ) + f )dt + f )
0 o
As can be seen from Figure 11.4, the maximum error is 1% when the phase angle error is zero
and maximum error in angle is tan -1 (.01OP/OP) = tan -1 (.01) = 0.572938698 @ 0.573.
d 2 j (t ) d
ROCOF(t ) = = [ D f (t )]
dt 2 2 p dt
Frequency in phasor measurements may be reported as the actual frequency f(t) or the
deviation of frequency from nominal, Df(t). In steady state conditions, Df(t) can be represented
as a scalar number.
The SOC counter started filling up at 00:00 on January 1st, 1970 and is getting incremented
every second. One may wonder how long it will take for the counter to become full and hence
over-run in the next second. It can be calculated as:
Seconds sec min hr day
= = (60)(60)(24)(365) = 31536000
Year min hr day year
Thus, it will be 136 years from January 1st, 1970, i.e., 2106 A.D., by the time SOC counter
over-runs. In all probability long before this happens, the technology will have changed and
hopefully will not cause something like the infamous Y2K problem during roll-over from 1999
to 2000.
The exact number of seconds, including the fractional part, since January 1st, 1970 can be
computed from the time stamp as:
FRACSEC
Exact number of seconds = SOC +
224
The PMU has to communicate the measurements with other devices such as phasor data
concentrators. Some protocol has to be used for such communications. The protocol allows
for necessary identifying information, such as the PMU IDCODE and status to be embedded
in data frame for proper interpretation of the measured data. Four message types are defined,
i.e., data, configuration, header, and command. The first three message types are transmitted
from the PMU and the last, i.e., command is received by the PMU. Configuration is a machine
readable message describing the data that the PMU sends and providing calibration factors.
Header is human-readable descriptive information sent from the PMU but provided by the
user. Commands are machine-readable codes sent to the PMU for control and configuration.
Information can be stored in any convenient form in the PMU but when transmitted it shall be
formatted as frames. Figure 11.6 shows the configuration of a data frame sent out by a PMU.
System frequency 50 Hz 60 Hz
Reporting rates in frames per second 10 25 10 12 15 20 30
Figure 11.7 Block diagram of a synchrophasor enabled relay or phasor measurement unit.
Synchrophasors227
approach, we keep the sampling frequency locked to a fictitious cosine wave whose positive
peak synchronised with 00 hr: 00 min : 00 sec on January 1st, 1970 as shown in the figure.
11.8.1 A
Demonstration of Error Caused by Off-nominal Frequency
Operation
Let the nominal frequency be 50 Hz and the sampling be done at four times the nominal
1 1
frequency, i.e., at 200 Hz. Hence, the sampling time period = Dt = = = (5)10 -3 s.
fsampling 200
Now let us consider two phasors, both of value 1 0 but one having frequency of
50 Hz and another of 51 Hz and sample both of them at 200 Hz for one cycle thus drawing
four samples of each. From these four samples let us evaluate the DFT and compare the
computed and actual values, thus finding errors, if any. The samples of the two signals
are shown in Table 11.2.
228Digital Power System Protection
V50 = v50 (0) ( -90(0)) + v50 (1) ( -90(1)) + v50 (2) ( -90(2)) + v50 (3) ( -90(3))
V50 = 1 ( -90(0)) + 0 ( -90(1)) - 1 ( -90(2)) + 0 ( -90(3))
V50 = 1 ((0)) - 1 ( -180) = 2 0
2 2
Hence, actual (peak) value of V50 = V50 = 20 = 10 .
N 4
Thus, there is no error, either in magnitude or in phase angle when the 50 Hz phasor is sampled
at 4 50 Hz frequency. Let us do similar calculations on the samples of the 51 Hz signal as
shown below:
3 - j 2 p 1n 3 - jp n 3
v51 (n)e 4 = v51 (n)e 2 = v51 (n) ( -90(n))
n=0 n=0 n=0
V51 = v51 (0) ( -90(0)) + v51 (1) ( -90(1)) + v51 (2) ( -90(2)) + v51 (3) ( -90(3))
V51 = 1 ( -90(0)) - 0.03141 ( -90(1)) - 0.9980 ( -90(2)) + 0.0941 ( -90(3))
V51 = 2.001965 ( 3.5946)
2 2
Hence, actual (peak) value of V51= V51 = 2.001965(3.5946) = 1.00098273.5946 .
N 4
Thus, there is error, both in magnitude and in phase angle when the 51 Hz phasor is sampled
at 4 50 Hz frequency. It can be seen that the error in magnitude is only 0.098% which
is less than 0.1% while the error in phase angle is 3.5946. Thus, phase angle error in
synchrophasors when off-nominal frequency signals are encountered is more serious than
the error in magnitude.
where fnom is the nominal frequency (50 Hz or 60 Hz). It is assumed that N is even. Let the
sample belong to an off-nominal frequency f . Hence, from the definition of a phasor we can
express the sample as:
x[ nDt ] = 2 Real [ X e j 2p f nDt ]
where X is the true phasor value of the samples.
j 2p f n
x[ nD t ] = 2 Real X e N fnom
The one cycle phasor estimate can then be written using the DFT equation for an off-
nominal frequency signal as:
2 N -1
j 2p f n - j 2p n
X=
N
2 Real X e Nfnom e N
n=0
Now for any phasor P , the Real ( P) can be written as
P + P*
Real (P) =
2
Hence, we can write the above equation during kth window as:
k + N -1 - j 2p n 1 k + N -1 - j 2 p f n - j 2p n
j 2p f n
1
X= X e N f nom e N + X * e N fnom e N
N n=k N n=k
1 k + N -1 1 k + N -1 - j 2 p f n - j 2 p n
j 2 p f n j 2p n
-
X = X e N f nom N + X * e N fnom N
N n=k N n=k
j 2p n f k + N -1 - j 2p n f +1
1 k + N -1 N -1
f nom + 1 X * e N fnom
X = X e
N n=k N n=k
n n
j 2p f k + N -1 - j 2p f
1 k + N -1 -1
N f nom f +1
X = X e + 1 X *
e
N nom
N n=k N n=k
230Digital Power System Protection
k q k q
j 2p f N -1
j 2p f - j 2 p f +1 N -1
- j 2p f
1 -1
N f nom
-1
N f nom + 1 X * e N fnom
+1
N f nom
X = X e
e e
N q =0 N q =0
1 2p f
N sin f - 1
2 N nom
N -2 p f
k
- j 2 p f +1 sin 2 N f + 1 -1 2 p f
j ( N -1) 2 N f +1
Second term = X * e
N f nom nom nom
e
-1 2 p f
N sin f + 1
2 N nom
Noting that
f f - f nom w - w nom f w + w nom
- 1 = = and + 1 =
f nom f nom w nom f nom w nom
Synchrophasors231
w - w nom
j 2 p k w -w nom sin p p w -w nom
N w w nom e
j ( N -1)
w nom
First term = X e nom N
p w - w nom
N sin
N w nom
Similarly,
w + w nom
- j 2 p k w +w nom sin p - j ( N -1) p w +w nom
w w nom e w
nom nom
Second term = X e
N N
p w + w nom
N sin
N w nom
Hence, expression for the off-nominal frequency DFT in terms of the true DFT X can be
written as:
j 2p k w -w nom - j 2p k w +w nom
N w nom N w nom
X = AXe + BX * e
w - w nom
sin p w -w nom
w nom
p
j ( N -1)
where A= e
N w
nom
p w - w nom
N sin
N w nom
w + w nom
sin p w +w nom
w nom - j ( N -1) N
p
and B = e
w
nom
p w + w nom
N sin
N w nom
The variation of magnitudes of A and B as fnom changes from 45 Hz to 55 Hz, i.e., a
change of 5 Hz as shown in Figures 11.10 and 11.11.
Further it can be shown that the locus of all phasor estimates, for a given phasor amplitude
and variable phase angle which is a circle at the nominal frequency happens to be an ellipse
at off-nominal frequency.
v50;
voff_nom1;
V50spectrum = fft(v50);
V50_complex(p) = (2/N)*V50spectrum(2);
V50mag(p) = abs(V50_complex(p));
voff_nom1spectrum = fft(voff_nom1);
voff_nom1_complex(p) = (2/N)*voff_nom1spectrum(2);
voff_nom1mag(p) = abs(voff_nom1_complex(p));
err_mag(p) = voff_nom1mag(p)-100;
voff_nom1angle(p) = (angle(voff_nom1_complex(p)))*(180/pi);
voff_nom2spectrum = fft(voff_nom2);
voff_nom2_complex(p) = (2/N)*voff_nom2spectrum(2);
voff_nom2mag(p) = abs(voff_nom2_complex(p));
err_mag(p) = voff_nom2mag(p)-100;
voff_nom2angle(p) = (angle(voff_nom2_complex(p)))*(180/pi);
end
figure(1)
polar(ph_angle,V50mag)
title(Variation magnitude w.r.t. phase for 50Hz signal sampled at
4*50 Hz)
figure(2)
polar(ph_angle,voff_nom1mag)
title(Variation magnitude w.r.t. phase for 45Hz signal sampled at
4*50 Hz)
figure(3)
polar(ph_angle,voff_nom2mag)
title(Variation magnitude w.r.t. phase for 55Hz signal sampled at
4*50 Hz)
Figure 11.12 Variation of magnitude w.r.t. phase for 50 Hz signal sampled at 4 50 Hz.
234Digital Power System Protection
Figure 11.13 Variation of magnitude w.r.t. phase for 45 Hz signal sampled at 4 50 Hz.
Figure 11.14 Variation of magnitude w.r.t. phase for 55 Hz signal sampled at 4 50 Hz.
Thus, the calculated state is at best an approximation to an averaged system state. Hence, the
estimates that are produced by the state estimation program are called static state estimates.
Synchronised phasor measurements of positive sequence bus voltage and currents most
naturally lead to state estimates. In the utopian view of the power system, every location is
equipped with synchrophasor measurements and information from all is available at a central
location. Thus, we have is a real-time state measurement rather than state estimation. However,
it is not necessary to have synchrophasor at every location to perform state estimation. It
can be shown that we can take measurements from a few key locations and perform state
estimation. This is possible because we can divide the network into observable locations where
synchrophasor are located and unobservable locations where synchrophasor are not located
but estimates can be indirectly calculated. From estimated values at observable locations, we
can make close estimates at unobservable locations.
2. Power system control: Prior to synchrophasor measurements, all control in power systems
used local measurements and a mathematical model of the larger power system was used to
arrive at control decisions. It was recognised that such controllers are rarely optimum and
could produce completely unacceptable response to system phenomenon when the models are
inaccurate. Synchrophasor measurements can improve such control actions. Typical applications
to problems where the control objectives were global in nature, i.e., for example an HVDC
controller may be called upon to damp electromechanical oscillations between two widely
separated areas of the power system. The high speed synchronised phasor measurements
provided by PMU offers a very attractive alternative to the problem of power system controllers.
Remote measurements could be brought to the controlled device at high speed, and the remote
inputs used as feedback signals in the controller. Phasor data are being used increasingly by
individual utilities to manage real-time grid operations.
3.Wide area power system protection: An additional protection area where phasor
measurements play a role is that of adaptive security and dependability. The existing protection
system has redundant primary protection coupled with multiple backup schemes. The resulting
system is highly dependable in that virtually every fault is ultimately cleared. The trade-off
is false trips which are tolerated. As the system has evolved, however, the fact that false
trips during large disturbances exacerbate the disturbance and allow it to cascade has been
recognised. The solution is to adaptively alter the security dependability balance during times
of stress. Remote phasor measurements can be used to determine when this happens. The
mechanism is to alter the relaying logic which normally lets any relay which sees fault to trip
the breaker to logic that demands a vote such as two out of three. This controlled dependability-
security system is an example of adaptive relaying. Other adaptive relaying applications
include distance relay zone adjustments, transformer protections, reclosing operations, out-of-
step relays, etc. Another popular innovation in modern protection systems is the increasing use
of remedial actions schemes (RAS) now renamed as system integrity protection (SIPS). These
schemes use wide area measurements to identify system conditions under which extraordinary
actions must be taken to avoid major system outages. Since the synchrophasor measurement
provides very detailed and precise time synchronised measurements of the state of the electrical
grid, the analysis of past electricity blackouts has indicated that monitoring synchrophasors
has the potential to detect system instability much earlier than current supervisory control and
236Digital Power System Protection
data acquisition (SCADA) systems. This will prevent large-scale blackouts by giving system
operators more time to respond to disruptive events enabling them to take corrective action
before system stability is affected. Synchrophasor projects involving the installation of hundreds
of phasor measurement units (PMUs) are currently taking place throughout the globe.
4.Forensic analysis: The dictionary meaning of the word forensic is relating to the
application of science to decide questions arising from crime or litigation. In case of major
system events like blackouts, the data provided by PMUs becomes the basis on which the
disturbance can be analysed. Since, PMUs collect and store high volumes of high-speed time
synchronised data about the conditions across an interconnection, those data can be compiled
quickly and analysed to determine the sequence of events which caused the disturbance. Phasor
data were essential in the recent investigations of blackouts in Europe and the USA.
5. Smart grid: At the transmission and generation level, synchrophasor systems are the single
most effective technology element to realise and implement smart grid, because synchrophasor
systems collect, distribute and analyse critical data and convert it in real time into information
and insight that improve grid automation and operation.
Review Questions
The current in the above circuit is given by writing the KVL around the loop as:
di(t )
v(t ) = R i(t ) + L
dt
237
238Digital Power System Protection
iDC = - I m e - t / t sin(a - f )
Hence, the RMS value of the current can be written as:
I RMS = I AC
2
+ I DC
2
Thus, it can be seen that the RMS value of current with DC offset is more than the RMS
value of the pure ac current. This can cause fast acting over current as well as distance relays
to over-reach. Further, presence of decaying DC offset is troublesome for the DFT algorithm.
The classical DFT algorithm which is based on Fourier analysis gives correct results when a
constant DC is present in the signal but gives erroneous results when decaying DC offset is
present. This is because while a constant DC (in time domain) shows up at only zero frequency
(in frequency domain), a decaying dc offset in the time domain has a wider frequency spectrum
and thus corrupts all the DFT bins. This can be seen from the Fourier transform of eat, denoted
by FT( ), as shown below:
1
FT(e -a t ) =
a + jw
Numerous researchers have worked on this problem and there is a large body of literature
on this topic. We discuss three approaches is the next sections. Interested reader is urged to
follow the literature on this subject.
d I 0 e-t / t
vmimic = KR I 0 e - t / t + K L
dt
1
vmimic = KR I 0 e - t / t - KL I 0 e - t / t
t
240Digital Power System Protection
Noting that t = L
R R -t / t
vmimic = KR I 0 e - t / t - KL I 0 e
L
vmimic = KR I 0 e - t / t - KR I 0 e - t / t = 0
Thus, the mimic impedance completely rejects the decaying dc offset. However, one must
remember that this complete rejection was made possible because the mimic circuit parameters
had been exactly matched with the time constant of the decaying DC current. This may not
always be possible and may lead to incomplete rejection of the DC offset.
{
V (z) = K R I (z) +
L I (z)
R Dt } L
[1 - Z -1 ] ; Noting that = t
R
{
V (z) = K R I (z) +
t I ( z)
Dt }
[1 - Z -1 ] ; representing
t
Dt
= t1
V ( z ) = K R I ( z ){1 + t 1[1 - Z -1 ]}
V (z)
= K R {1 + t 1[1 - Z -1 ]}
I (z)
V (z)
= H ( z ) = K R {(1 + t 1 ) - t 1 Z -1 ]}
I (z)
which is written as:
V (z)
= H ( z ) = K m {(1 + t 1 ) - t 1 Z -1 ]}
I (z)
where Km = K R can be thought of as a gain term.
The gain Km is set in such a way that at rated frequency, the filter gain will be one.
Substituting Z = e j w T = e j 2 p f T = e j 2 p 50 T ; we can write:
Gain at 50 Hz = K m {(1 + t 1 ) - t 1 e - j 2 p 50 T ]}
where T is the time period corresponding to the sampling frequency which can be expressed
as T = (1/fsampling). Hence,
- j 2 p 50 /fsampling
K m {(1 + t 1 ) - t 1 e }=1
2 p 50 2 p 50
K m (1 + t 1 ) - t 1 cos + j t 1 sin = 1
fsampling fsampling
2 2
2 p 50 2 p 50
Km (1 + t 1 ) - t 1 cos + t 1 sin =1
fsampling fsampling
242Digital Power System Protection
1
Km =
2 2
2 p 50 2 p 50
(1 + t 1 ) - t 1 cos + t 1 sin
fsampling fsampling
It can be easily appreciated that the mimic filter is a high-pass filter will therefore amplify
any noise contained in the CT secondary current.
Thus, giving
N /2
PS1 = A0 r 2 k -1
k =1
which can be written in the form highlighting the underlying geometric series as:
N /2
PS1 = A0 r -1 r 2 k
k =1
N /2
PS1 = A0 r -1 (r 2 ) k
k =1
r N r2 - r2
PS1 = A0 r -1
r2 - 1
rNr - r
PS1 = A0
r2 - 1
r (r N - 1)
PS1 = A0
r2 - 1
Using similar reasoning, the second term in the above equation turns out to be zero, giving:
N
-1
2
PS2 = A0 r 2 k
k =0
(r 2 ) N/2 - (r 2 )0
PS2 = A0
r2 - 1
rN - 1
PS2 = A0
r2 - 1
From equations for PS1 and PS2, we can write:
PS1
r=
PS2
r2 - 1
and A0 = PS2
(r N - 1)
Having the values of r and A0, the samples can be modified as:
Sknew = Sk A0rk; k = 0 to N 1
This, new set of samples, is free from decaying DC offset and can be used for further
processing.
244Digital Power System Protection
where I0 is the peak value of the dc offset, t is the time constant of the decaying dc offset,
k is the harmonic order, Ik is the peak value of the kth harmonic and qk is the phase angle of
the kth harmonic component and p is the maximum harmonic order. If the above expression
is integrated for one time period T of the fundamental, then second term evaluates to zero,
and we are left with integral of the decaying dc offset term as shown below:
t t
p
i(t )dt = I
0 e -t / t
+ ( I k sin(k w1 t + q k )) dt
t -T t -T k =1
t
-t / t
= I0 e dt = - I 0t [e - t / t ]tt -T
t -T
= - I 0t (e - t / t - e - (t -T ) / t ) = - I 0t (e - t / t - e - t / t eT / t )
= - I 0t e - t / t (1 - eT / t )
Let us call the above result of the integral as Z(t). Hence,
Z (t ) = - I 0t e - t / t (1 - eT / t )
Now, the expression for Z(t + Dt), which is the integral of the dc component Z(t) extended
by a small interval of time Dt can be written as:
Z (t + Dt ) = - I 0t e - (t + Dt ) / t (1 - eT / t )
Z (t + Dt ) = - I 0t e - t
/t
(1 - eT / t) e -Dt / t
Z (t )
Z (t + Dt ) = Z (t )e Dt/t
Simplifying further,
Dt Z (t + Dt )
- = ln
t Z (t )
Thus, the time constant of the decaying dc offset t is found out as:
Dt
t =-
Z (t + Dt )
ln
Z (t )
Removal of DC Offset245
Thus, once I0 and t are found out, we can completely characterise the decaying dc offset
and compensate the raw signals to extract decaying offset free part of the signal as explained
earlier.
Review Questions
1.
Clearly state the difference between (constant) dc offset and decaying DC offset.
2.
Sketch the frequency spectra of a decaying dc offset, thus explain why it corrupts all
the DFT bins.
3.
DFT cannot filter out the effect of decaying dc offset while it does so for constant
dc offset. Explain.
4.
Describe how decaying dc offsets are created in ac power system.
5.
What is the problem caused by existence of decaying dc offsets?
6.
Explain whether decaying DC offset is present in output of CT or output of PT or
both.
7.
Out of current and voltage signals which one is more likely to exhibit decaying dc
offset?
8.
Explain the effect of decaying DC offset on OC relays.
9.
Explain the effect of decaying dc offset on distance relays.
10.
Explain the mimic impedance method of filtering decaying dc offset.
11.
Explain the digital implementation of the mimic impedance for filtering out decaying
dc offset.
12.
Explain why the mimic impedance acts as a high-pass filter. What are the implications
of this?
13.
Explain the method of partial sums to filter out decaying dc offset.
14.
Explain how dc offset can be characterised by integration over the fundamental time
period, and this can be used to filter out the decaying dc offset.
Appendix
Range and time is broadcast on three L-band frequencies. Time can be derived from
the coarse/acquisition (C/A) code transmitted at 1575 MHz. All satellites broadcast at the
same frequencies. Signals are encoded using code division multiple access (CDMA) allowing
messages from individual satellites to be distinguished from each other based on unique
encodings for each satellite (that the receiver must be aware of). Two distinct types of CDMA
247
248Appendix
encodings are used, i.e., the coarse/acquisition (C/A) code, which is accessible by the general
public, and the precise (P(Y)) code, which is encrypted so that only the U.S. military can
access it. The signal can be received by a simple omnidirectional antenna. The spread-spectrum
technique that is used, makes the signal resistant to interference. At the heart of the GPS is
a ground based cesium clock ensemble that itself is referenced to the universal coordinated
time (UTC), which is the world time standard. Each satellite provides a correction to the UTC
time that the receiver automatically applies to the outputs. With continuous adjustments, the
time accuracy is limited only by short term signal reception whose basic accuracy is 0.2 s.
Baseline accuracy can be improved by advanced decoding and processing techniques. For
general applications a 0.5 s is a safer figure to depend on. The short term frequency accuracy
of the received signal is of the order of 10E-11.
One might wonder how ordinary public throughout the world can access the U.S. Department
of Defense satellite but access issues were largely resolved in the 1990s with agreements to
add new civilian frequencies and culminating in termination of the selective availability signal
degradation of the C/A code in the year 2000. The U.S. Department of Defence and Department
of Transportation have committed to make GPS available to civilian users at all-time except
in national emergency.
A.2 Comtrade
COMTRADE (Common format for Transient Data Exchange for power systems) is a file
formatfor storing oscillography and status data related to transient power system disturbances.
COMTRADE files are typically generated by Intelligent Electronic Devices (IEDs), such
as a protective relay, in electrical substations. Analysis tools can download the COMTRADE
file and calculate useful information related to the power system. COMTRADE files from
multiple substations can be used collectively to perform forensic analysis of large scale power
disturbance events (e.g., blackouts) to determine the root cause of the disturbance, help to
improve system protection and guide future mitigation strategies.
The COMTRADE file format has been standardised by the Power System Relaying
Committee (PSRC) of the IEEE Power & Energy Society as C37.111. The most widely used
version of the COMTRADE standard is C37.111-1999. This version specified a file format that
consisted of multiple file types designated by the assigned file extensions of *.CFG, *.INF,
*.HDR, and *.DAT. The *.DAT file contains the digitised sample data in an ASCII text format.
The *.CFG file contains configuration data as in the *.DAT file including information such as
signal names, start time of the samples, number of samples, min/max values, and more. Only
the *.CFG and *.DAT files were mandatory. Although the values of the digitised samples in the
*.DAT file are viewable without the *.CFG file, the value of the data is greatly diminished as
it would be very difficult to fully reconstruct the meaning of the data without the *.CFG file.
As applications for using COMTRADE grew the limitations of the file format began to
cause problems. For instance, the North American Synchro Phasor Initiative (NASPI) sponsored
by the North American Electric Reliability Corporation (NERC) wanted to use the COMTRADE
file format to exchange synchrophasor data. The C37.1111999 version of the COMTRADE
format only supported ASCII encoded data, supported limited analog and digital status data
Appendix249
types, and did not convey any power system configuration information. This would have
made large scale exchange of synchrophasordata more difficult. Beginning in 2011, the IEEE
PSRC began work on a new version of the COMTRADE standard that has been published as
C37.111-2013. This version of the standard supports a variety of new features including XML
encoding configuration and data, a single file type (*.CFF) that combines all the C37.111-1999
types into a single file, additional data types for 32-bit binary and floating point values and a
partial binary encoding to reduce file size.
The IEEE Power Quality Subcommittee of the IEEE Power & Energy Society has also
developed a file format called the Power Quality Data Interchange Format (PQDIF) that is
similar to COMTRADE in structure but is used primarily to convey power quality data instead
of transient disturbance data.
There are numerous COMTRADE viewers available from many commercial companies as
well as free viewers like Comcalc Pro, TOP, GTPPLOT and others.
81 Frequency Relay
83 Automatic Selective Control or Transfer Relay
84 Operating Mechanism
85 Pilot Communications, Carrier or Pilot-Wire Relay
86 Lockout Relay
87 Differential Protective Relay
89 Line Switch
90 Regulating Device
91 Voltage Directional Relay
92 Voltage and Power Directional Relay
94 Tripping or Trip-Free Relay
B Bus
F Field
G Ground or generator
N Neutral
T Transformer
References
253
254References
[12] Guo, Yong, Mladen Kezunovic, and Deshu Chen, Simplified algorithms for removal
of the effect of exponentially decaying DC-offset on the Fourier algorithm, Power
Delivery, IEEE Transactions on 18.3, 2003: 711717.
[13] Hays, H.H., Schaums Outline of Theory and Problems of Digital Signal Processing,
Tata McGraw-Hill, New Delhi, 2004.
[14] Holbert, Keith E., G.I. Heydt, and Hui Ni, Use of satellite technologies for power
system measurements, command, and control, Proceedings of the IEEE 93.5, 2005:
947955.
[15] Hope, G.S. and V.S. Umamaheswaran, Sampling for computer protection of
transmission lines, Power Apparatus and Systems, IEEE Transactions on 5, 1974:
15221534.
[16] Horowitz, S H., A.G. Phadke, and J.S. Thorp, Adaptive transmission system
relaying, Power Delivery, IEEE Transactions on 3.4, 1988: 14361445.
[17] Horton, John W, The use of Walsh functions for high speed digital relaying, IEEE
Publication 75CH1034-8, 1975: 19.
[18] IEEE Standard for Synchrophasor Measurements for Power Systems, IEEE Std
C37.118.1TN 2011 and IEEE Std C37.118.2TN 2011.
[19] Ifeachor, Emmanuel C. and B.W. Jervis, Digital Signal ProcessingA Practical
Approach, Pearson, 2007.
[20] Jeyasurya, B. and W.J. Smolinski, Identification of a best algorihthm for digital
distance protection of transmission lines, Power Apparatus and Systems, IEEE
Transactions on 10, 1983: 33583369.
[21] Johns, A.T. and M.A. Martin, Fundamental digital approach to the distance protection
of EHV transmission lines, Proceedings of the Institution of Electrical Engineers,
Vol. 125, No. 5, IET Digital Library, 1978.
[22] Lyons, Richard G., Understanding Digital Signal Processing, 3rd ed., Prentice Hall,
2010.
[23] Madhav Rao, T.S., Power System Protection: Static Relays with Microprocessor
Application, 2nd ed., Tata McGraw-Hill, New Delhi, 2001.
[24] Mahmood, Mahir K., Janan E. Allos, and Majid AH Abdul-Karim, Microprocessor
implementation of a fast and simultaneous amplitude and frequency detector for
sinusoidal signals, Instrumentation and Measurement, IEEE Transactions on 34.3,
1985: 413417.
[25] Makino, J. and Y. Miki, Study of operating principles and digital filters for protective
relays with digital computer, IEEE PES Winter Power Meeting, 1975.
[26] Mann, B.J. and I.F. Morrison, Digital calculation of impedance for transmission line
protection, Power Apparatus and Systems, IEEE Transactions on 1, 1971: 270279.
[27] Mann, B.J. and I.F. Morrison, Relaying a three-phase transmission line with a digital
computer, Power Apparatus and Systems, IEEE Transactions on 2, 1971: 742750.
[28] Martin, K.E., et al., IEEE Standard for synchrophasors for power systems,Power
Delivery, IEEE Transactions on 13.1, 1998: 7377.
References255
[29] Mason, C. Russell, The Art and Science of Protective Relaying, Wiley, 1956.
[30] McInnes, A.D. and I.F. Morrison, Real time calculation of resistance and reactance
for transmission line protection by digital computer, IEEE Trans, 1971: 1623.
[31] McLaren, Peter Grant and M.A. Redfern, Fourier-series techniques applied to distance
protection, Proceedings of the Institution of Electrical Engineers, Vol. 122, No. 11,
IET Digital Library, 1975.
[32] North American Synchrophasor Initiative: Synchrophasor System Benefits Factsheet:
http://www.naspi.org.
[33] Paithankar, Y.G. and S.R. Bhide, Fundamentals of Power System Protection, PHI
Learning, New Delhi, 2010.
[34] Paithankar, Y.G., Transmission Network Protection, Theory and Practise, CRC Press,
1997.
[35] Phadke, A.G. and J.S. Thorp, Synchronised Phasor Measurements and Their
Applications, Springer, 2008.
[36] Phadke, A.G. and Bogdan Kasztenny, Synchronised phasor and frequency measurement
under transient conditions, Power Delivery, IEEE Transactions on24.1, 2009: 8995.
[37] Phadke, A.G. and J.S. Thorp, History and applications of phasor measurements,Power
Systems Conference and Exposition, 2006, PSCE06, 2006, IEEE PES, IEEE, 2006.
[38] Phadke, A.G., et al., Synchronised sampling and phasor measurements for relaying
and control,Power Delivery, IEEE Transactions on 9.1, 1994: 442452.
[39] Phadke, A.G., J.S. Thorp, and M.G. Adamiak, A new measurement technique
for tracking voltage phasors, local system frequency, and rate of change of
frequency,Power Apparatus and Systems, IEEE Transactions on5, 1983: 10251038.
[40] Phadke, A.G., T. Hlibka, and M. Ibrahim, A digital computer system for EHV
substations: analysis and field tests,Power Apparatus and Systems, IEEE Transactions
on 95.1, 1976: 291301.
[41] Phadke, Arun G. and James S. Thorp, Computer Relaying for Power Systems, 2nd
ed., John Wiley & Sons, 2009.
[42] Pratap, Rudra, Getting Started with MATLAB 7, Oxford University Press (Indian
Edition), 2006.
[43] Rahman, M.A. and B. Jeyasurya, A state-of-the-art review of transformer protection
algorithms, Power Delivery, IEEE Transactions on 3.2, 1988: 534544.
[44] Ram, Badri, Power System Protection and Switchgear, Tata McGraw-Hill, 2011.
[45] Ramamoorty, M., Application of digital computers to power system protection,
Journal of the Institution of Engineers, India, 52, 1972: 235238.
[46] Ranjbar, A.M. and B.J. Cory, Algorithms for distance protection, IEEE Conf. on
Developments in Power System Protection, Institution of Electrical Engineers, London,
1975.
[47] Ranjbar, A.M. and B.J. Cory, An improved method for the digital protection of
high voltage transmission lines,Power Apparatus and Systems, IEEE Transactions
on 94.2, 1975: 544550.
256References
[48] Rockefeller, G.D. and E.A. Udren, High-speed distance relaying using a digital
computer II-test results,Power Apparatus and Systems, IEEE Transactions on 3,
1972: 12441258.
[49] Rockefeller, George D., Fault protection with a digital computer,Power Apparatus
and Systems, IEEE Transactions on 4, 1969: 438464.
[50] Sachdev, M.S. and M.A. Baribeau, A new algorithm for digital impedance
relays, Power Apparatus and Systems, IEEE Transactions on 6, 1979: 22322240.
[51] Scarborough, James Blaine, Numerical Mathematical Analysis, Oxford and IBH
Publishing, 1955.
[52] Schilling, Robert Joseph and Sandra L. Harris, Introduction to Digital Signal
Processing, 2nd ed., Cengage Learning, 2011.
[53] Segui, Trinidad, et al., Fundamental basis for distance relaying with parametrical
estimation, Power Delivery, IEEE Transactions on 16.1, 2001: 99104.
[54] Sidhu, T.S., et al., Discrete-Fourier-transform-based technique for removal of decaying
DC offset from phasor estimates, IEEE Proceedings-Generation, Transmission and
Distribution 150.6, 2003: 745752.
[55] Singh, Lankeshwar Prakash, Digital Protection: Protective Relaying from
Electromechanical to Microprocessor, New Age International, 1997.
[56] Smolinski, W.J., An algorithm for digital impedance calculation using a single PI
section transmission line model, Power Apparatus and Systems, IEEE Transactions
on 5, 1979: 15461551.
[57] Soman, S.A., Power System Protection Course, NPTEL Website.
[58] Strang, Gilbert, Linear Algebra and Its Applications, 4th ed., Cengage Learning, 2006.
[59] Thorp, J.S. and A.G. Phadke, A microprocessor based three-phase transformer
differential relay,Power Apparatus and Systems, IEEE Transactions on 2, 1982:
426432.
[60] Thorp, J.S., et al., Limits to impedance relaying, Power Apparatus and Systems,
IEEE Transactions on 1, 1979: 246260.
[61] Van Valkenburg, M.E., Network Analysis, Prentice Hall Inc., 1974.
[62] Warrington, AR van C. (Ed.), Protective Relays Their Theory and Practise: Volume
I, 1962.
[63] Warrington, AR van C. (Ed.), Protective Relays Their Theory and Practise: Volume
II, 1969.
[64] Wilson, Robert E., Uses of precise time and frequency in power systems, Proceedings
of the IEEE 79.7, 1991: 10091018.
[65] Yu, Sun-Li and Jyh-Cherng Gu, Removal of decaying DC in current and voltage
signals using a modified Fourier filter algorithm,Power Delivery, IEEE Transactions
on 16.3, 2001: 372379.
[66] Ziegler, Gerhard, Numerical Distance Protection, Siemens, Publicis, Germany, 2006.
Index