A Simple Algorithm for Power System

Frequency Estimation
Arghya Sarkar, Member, IEEE, and S. Sengupta, Senior Member, IEEE

Abstract--This paper presents development and II. FREQUENCY ESTIMATION ALGORITHM

implementation of a novel digital signal processing algorithm for Let, s(t) represents a continuously measured non-sinusoidal
on-line estimation of the fundamental frequency of non-
sinusoidal power system signal. The algorithm relies on stable signal (arbitrary voltage or current). In a typical situation, it
band-pass second degree digital integrator (BPSDDI) based has a general form of
rigorous mathematical deduction and recombines itself with s ( t ) = Smax 0 + ∑ Smax k sin ( 2π kft + α k ) + ξ ( t ) (1)
zero-crossing avoiding technique. The proposed algorithm k ≠0
provides high degree of immunity and insensitivity to harmonics where ξ(t) is zero mean random noise, Smax 0 is the dc
and noise, and fast response during time varying conditions.
Structural simplicity, wide range of application, controls over component of the signal, f is the fundamental frequency,
speed and accuracy, and parameter robustness are other salient Smax k are the peak values and α k are the phase angles, at kth
features of the method. The simulation results confirm the
validity and accurate performance of the proposed approach harmonic. In practice, all the parameters in (1) may undergo
under different operating conditions. variations with respect to time.
Pre-filtering of non-sinusoidal signal defined in (1) by a
Index Terms-- Algorithm, filter, frequency estimation, suitable band-pass filter (BPF) can eliminates dc offset and
harmonics, second degree digital integrator. higher order harmonics keeping fundamental component
almost unaltered and leads the following observation model
I. INTRODUCTION sF ( t ) = Smax1 sin ( 2π ft + α F ) + ξ F ( t ) (2)

F REQUENCY information takes a great concern in power where αF(t) and ξF(t) are the new phase angle and noise of
system operation, control and protection. While this is sF(t), respectively, due to band-pass filtering operation.
straightforward under sinusoidal environments, the task of If all the parameters remain constant within the observation
frequency estimation becomes challenging and interesting in window, the second degree integration of (2) takes the
the presence of harmonic distortion or noise. A variety of following form
techniques and algorithms have been proposed in different t ⎡t ⎤
literature to estimate frequency, for example modified zero sFI ( t ) = ∫ ⎢ ∫ sF ( t ) dt ⎥ dt
0 ⎣0 ⎦
crossing technique [1], discrete Fourier transform (DFT) [2],
extended Kalman filtering [3], least mean square (LMS) Smax1 sin (α F ) Smax1 t cos (α F )
= + (3)
algorithm [4], statistical methods [5] and orthogonal finite 4π 2 f 2 2π f
impulse response (FIR) filters based recursive algorithm [6].
Smax1 sin ( 2π ft + α F )
For real-time use, most of the aforementioned methods have − + ξ FI ( t )
trade-off between accuracy and speed [2]. 4π 2 f 2
In this paper, a novel digital signal processing algorithm where ξ FI ( t ) is the noise signal.
has been proposed to estimate the fundamental frequency of
the distorted power system signal. Theoretical aspects of the The second term of (3) is a ramp signal and makes the
proposed algorithm have been discussed, the measurement output of the second degree integrator totally unstable. An
accuracy and response time are evaluated, and the simulation appropriate choice of BPF can eliminate dc part (first term)
results are presented. and ramp component of (3) and gives the following stable
sinusoidal component with noise signal ξ FIBPF ( t )
Smax1 sin ( 2π ft + α F )
Arghya Sarkar is with the MCKV Institute of Engineering, 243 G. T.
sFIBPF ( t ) = − + ξ FIBPF ( t ) (4)
Road (N), Liluah, Howrah-711204,(e-mail- [email protected]) 4π 2 f 2
S. Sengupta is with the Department of Applied Physics, University of If, the distorted signal contains a finite number of
Calcutta, 92, A. P. C. Road, Kolkata - 700009, ( e-mail - significant harmonics and uniformly sampled at sampling
[email protected]).
frequency greater than the Nyquist rate, then (2) and (4) can
be discretized at any arbitrary sample instant, n, as
978-1-4244-7781-4/10/$26.00 ©2010 IEEE
 sF [ n ] = Smax1 sin ( 2π fn + α F ) + ξ F [ n ] readings sFIBPFT [ n] are taken and would thus tend to be
(5)
= sFT [ n ] + ξ F [ n ] worst where sFIBPFT [ n] is small, i.e., in the region of a zero
Smax1 sin ( 2π fn + α F ) crossing. Hence, according to following theorem [8]
sFIBPF [ n] = − + ξ FIBPF [ n ] If
4π 2 f 2 (6)
a c e
= sFIBPFT [ n ] + ξ FIBPF [ n ] = = = ...... > 0 (12)
b d f
where sFT [ n ] = Smax1 sin ( 2π fn + α F ) (7) Then
Smax1 sin ( 2π fn + α F ) a a + c + e + ......
and sFIBPFT [ n ] = − (8) = (13)
b b + d + f + ......
4π 2 f 2
sFT [ n ] and sFIBPFT [ n] are the true samples of sF [ n] and r̂ [ n ] can be modified as
A −1
sFIBPF [ n] , respectively, contain theoretical sinusoidal values.
∑ s FT [ n ]
ξ F [ n ] and ξ FIBPF [ n ] are referred to the values ξ F ( t ) and rˆ [ n ] =
j =0 (14)
A −1
ξ FIBPF ( t ) , respectively, at discrete time index n, with ∑ s FIBPFT [ n ]
j =0
additional noise due to quantization.
From (14) it has been observed that the problem of the
If, the true samples sFT [ n ] and sFIBPFT [ n] have been
proximity to zero crossings has been completely eliminated
considered, the fundamental frequency can easily be obtained due to the sum of modulus of consecutive samples of
from (8) and (9), as sFIBPFT [ n] .
1
f [ n] = r [ n] (9) Since, the theoretical values of sFT [ n ] and sFIBPFT [ n]
2π
s FT [ n ] are unknown, the obtained samples sF [ n] and sFIBPF [ n] ,
where r [ n] = (10)
− sFIBPFT [ n ] are used in their place to obtain the final expression of the
fundamental frequency. Hence, the fundamental frequency can
Since, the true samples sFT [ n ] and sFIBPFT [ n] are be estimated from (9) and (14) as
subjected to additive noise ξ F [ n ] and ξ FIBPF [ n ] , (10) A −1

becomes an estimator, r̂ [ n ] , of r [ n ] . If, the variance of 1

∑ s FT [ n ]
f [n ] =
j =0 (15)
2π A −1
sFT [ n ] and sFIBPFT [ n] are 2 2
σ FIBPFT
σ FT and , respectively,
∑ s FIBPFT [ n ]
and σ FT , FIBPFT is the corresponding covariance, then j =0
Although, (15) has been derived considering all the
approximate expression for the variance of r̂ [ n ] , Var ( rˆ [ n ]) ,
instantaneous parameters (amplitude, frequency, and phase) of
can be obtained from the standard approximations via Taylor’s the input signal constant, simulation results (Section IV-F)
theorem [7] as show that the developed algorithm also provides excellent
2 2 accuracy during time varying conditions.
⎛ ∂r [ n] ⎞ 2 ⎛ ∂r [ n] ⎞
Var ( rˆ [ n]) = ⎜ ⎟⎟ σ FT +⎜ 2
⎟ σ FIBPFT
⎝ ∂sFT [ n] ⎠
⎜ ⎜ ∂ ( sFIBPFT [ n]) ⎟ III. DESIGN CONSIDERATION
⎝ ⎠
For a practical implementation of the fundamental
⎛ ∂r [ n] ⎞ ⎛ ∂r [ n] ⎞
frequency estimation algorithm as developed in the previous
+2 ⎜ ⎟ ⎜ ⎟σ FT , FIBPFT
⎜ ∂sFT [ n] ⎟ ⎜ ∂ ( s
⎝ ⎠ ⎝ FIBPFT [ n]) ⎟⎠ section, it is necessary to (A) design stable BPSDDI with high
accuracy, and (B) choose suitable value of A.
2
⎛
=⎜
1 ⎞ 2
σ FT +
2
( sFT [ n]) σ 2 A. Design of BPSDDI
⎟
⎜ sFIBPFT [ n] ⎟ 4 FIBPFT (11)
⎝ ⎠ ( sFIBPFT [ n]) The basic requirement is to design a suitable second degree
digital integrator so that there is a π phase difference between
⎛ 1 ⎞ sFT [ n] sF [ n ] and sFIBPF [ n ] . At the same time both the signals
−2 ⎜ σ
⎜ sFIBPFT [ n] ⎟⎟ 2 FT , FIBPFT
⎝ ⎠ ( sFIBPFT [ n]) should not have any dc offset. The proposed design procedure
starts from conventional Simpson rule based first degree
2
2
σ FT + ( r [ n]) σ FIBPFT
2
+ 2r [ n]σ FT , FIBPFT digital integrator (FDDI), then extend it to second degree case.
= The detailed design procedure has been explained below.
2
( sFIBPFT [n]) The transfer function of a Simpson FDDI is given by [9]
( )
Clearly, Var rˆ [ n ] depends on the region where the
 H FDDI ( z ) =
(
Ts 1 + 4 z −1 + z −2 ) (16)
BPSDDI, the distorted power system signals s [ n ] can directly
be fed to BPSDDI to get sFIBPF [ n] . At the same time, to
(
3 1 − z −2 ) maintain a phase difference, –π, between sF [ n ] and
where Ts is the sampling interval.
The simplest way to determine the transfer function of a sFIBPF [ n] , s [ n ] is filtered by the same 4 order Chebyshev th

double integrator is to cascade two FDDI. Hence, using (16) I BPF as characterized by (18), to get sF [ n] . These design
the transfer function of the second degree digital integrator
considerations and (19) lead the construction of following
(SDDI) can be obtained as:
block diagram (Fig. 1) for fundamental frequency estimation

H SDDI ( z ) = Ts2
(1 + 8z −1
+ 18 z −2 + 8 z −3 + z −4 ) (17)
algorithm. However, if distorted power system signals
contain a d.c component, a high-pass pre-filter is required to
(
9 1− z )
−2 2
diminish its effect.
The system described by (17) provides excellent accuracy B. Choice of A
within specified frequency range, but it is unstable, since two
Choice of A plays an important role in the accuracy and the
poles lie on the unit circle. In order to get a stable system, it is
computational complexity of the proposed algorithm. Larger
cascaded with an appropriate band-pass filter so that all poles
value of A will produce more accurate and smooth results but
of HSDDI are cancelled by the zeros of BPF. The selected BPF
at the cost of increased computational load. For variance
should also have ability of filtering the fundamental
reduction, A was found to be the best choice in terms of
component of power system signal with steep roll-off at
computational burden and smoothing criteria when it is almost
transition bands; and reduced delay and computational
equal to half of number of samples per fundamental cycle.
complexity. To that effect, a fourth-order Chebyshev I digital
Hence, it has been defined as
BPF with cut-off frequencies 48 and 52 Hz has been chosen
(as presented in [10]) whose transfer function is of the form A [ n ] = round ( f s 2 f [ n − 1]) to the nearest integer (20)

( )
2
GBPF 1 − z −2
H BPF ( z ) = (18) IV. PERFORMANCE ANALYSIS
a ( 0 ) + a (1) z −1 + a ( 2 ) z −2 + a ( 3) z −3 +a ( 4 ) z −4 A set of simulation test has been performed in MATLAB
where GBPF is the filter gain and a(0), a(1),….,a(4) are the environment to estimate the validity and performance of the
denominator coefficients. Since, the numerator of (18) proposed algorithm under different operating conditions.

( )
2
contains 1 − z −2 term, pole-zero cancellation has been A. Static Sinusoidal Test
In this test, sinusoidal signals with amplitude 1 p.u. and
occurred. Hence, from (17) and (18) the transfer function of frequencies in the range from 30 to 70 Hz in steps of 5 Hz
the stable BPSDDI can be achieved as: have been provided as inputs to the algorithms. The absolute
H BPSDDI ( z ) = HSDDI ( z ) H BPF ( z ) error in frequency estimates has been shown in Fig. 2 which

=
(
T 2GBPF 1 + 8z−1 + 18z−2 + 8z −3 + z−4 ) (19) reveals that a very small (~μHz) error is present within the
measured frequency. This result also confirms the capability
(
9 a ( 0) + a (1) z−1 + a ( 2) z−2 + a ( 3) z−3 +a ( 4) z−4 ) of the proposed algorithm to measure the fundamental
The system presented in (19) is stable because all poles of frequency over a wide range.
-5
its transfer function lie within unit circle. 1.2
x 10
Absolute Error in Frequency Estimation (Hz)

0.8

0.6

0.4

0.2

0
30 35 40 45 50 55 60 65 70
Input Signal Frequency (Hz)

Fig. 2. Absolute steady-state errors in frequency estimation under static

sinusoidal conditions.
Fig. 1. Block diagram of the fundamental frequency estimation algorithm.
Since, BPF is an integral part of the designed stable
B. Static Harmonic Test 0.025

To study the effect of harmonics on the performance of the

proposed frequency estimator, an fundamental frequency of 50 0.02

Hz, a third harmonic component in the range from 0% to 60%,

and a fifth harmonic component equal to a half of the third
0.015

|Error| (Hz)
component have been utilized. The absolute maximum errors 0.01
for the proposed technique has been shown in Fig.3, from
which it can be concluded that high accuracy can be achieved 0.005

utilizing the proposed scheme.

C. Noise Test
0
0 10 20 30 40 50 60
3 rd Harmonic (%)

Sinusoidal 50-Hz signals with the superimposed white

zero-mean Gaussian noise has been utilized as input test Fig. 3. Absolute maximum steady-state errors in terms of the third harmonic.
signals. A range from a highly noisy signal (SNR=20 dB) to a
low noisy signal (SNR=60 dB) is covered. The absolute peak 0.06

of oscillating steady state error has been depicted in Fig.4,

which exhibits that the error rapidly drops from SNR=20 dB 0.05

to SNR=60 dB and proposed algorithm provides high 0.04

accuracy at the presence of noise.

|Error| (Hz)
0.03
D. Dynamic Response during Step Variation of Frequency
A sinusoid with fixed amplitude of 1 pu has been utilized 0.02

for this test of which frequency abruptly drops from 51 to 48 0.01

Hz at 2s. Corresponding to the original system frequency, Fig.
5 shows the frequency estimation during step frequency 0
20 30 40 50 60

change. From figure it has been observed that the proposed SNR (dB)

scheme provides fast convergence with transient time equal to

Fig. 4. Absolute maximum steady-state errors in terms of the SNR.
60 ms (3 cycle).
E. Dynamic Response During Step Variation of Amplitude 51.5

As the response of a sudden increase in the amplitude of 50 51 Proposed

Hz sinusoids at 2s from 1 to 1.6 pu, tracking of frequency has Actual

50.5
been depicted in Fig.6, which exhibits that the proposed
approach again provides fast convergence (transient time = 80 50

ms). 49.5

F. Tracking of Instantaneous Frequency 49

Since, the basic assumption of the proposed algorithm is to 48.5

keep all the instantaneous parameters of input signal constant, 48

the performance of the presented technique during time- 47.5

varying condition is of great interest. Real non-stationary 1.98 2 2.02 2.04 2.06 2.08 2.1 2.12 2.14 2.16 2.18 2.2

signals can be characterized by slowly varying amplitude and

a polynomial phase as [11] Fig. 5. Frequency estimation during step frequency change.
s ( t ) = B ( t ) cos ( Φ ( t ) ) (21)
55
IB
where B ( t ) = ∑ bit i Proposed
is a slowly varying amplitude (IB is the 54 Actual

i =0 53

highest degree of the amplitude polynomial), called

Frequency (Hz)

52
IΦ
instantaneous amplitude and Φ (t ) = ∑φi t i is the 51

i =0 50

instantaneous phase ( I Φ is the highest degree of the phase 49

polynomial). 48
1.98 2 2.02 2.04 2.06 2.08 2.1 2.12 2.14 2.16 2.18 2.2
The instantaneous frequency is then given by [11] Time (s)

Fig. 6: Frequency estimation during step change in amplitude.

 I VI. REFERENCES
1 dΦ 1 Φ
f (t ) =
2π dt
=
2π ∑
iφi t i −1 (22) [1] M. M. Begovic, P. M. Djuric, S. Dunlap, and A. G. Phadke, "Frequency
tracking in power networks in the presence of harmonics," IEEE Trans.
i =0
Power Delivery, vol. 8, no. 2, pp. 480–486, Apr. 1993.
In this test, both amplitude and frequency of the input [2] J. Z. Yang, and C. W. Liu, "A precise calculation of power system
signal undergoes a polynomial variation with following frequency and phasor," IEEE Trans. Power Delivery, vol. 15, no. 2, pp.
polynomial coefficients 494-499, Apr. 2000.
[3] I. Kamwa and R. Grondin, "Fast adaptive schemes for tracking voltage
b1 = 0.5, b0 = 1, φ3 = −2π × 5, φ2 = 2π ×10, phasor and local frequency in power transmission and distribution
φ1 = 2π × 50, φ0 = (π 6 ) systems," IEEE Trans. Power Delivery, vol. 7, no. 2, pp. 789–795, Apr.
1992.
Fig. 7 corresponds to the tracking of the polynomial phase [4] M. S. Sachdev, and M. M. Giray, "A least error squares technique for
determining power system frequency," IEEE Trans. on Power
signal by the proposed algorithm. Close tracking of the Apparatus and Systems, vol. PAS-104, no. 2, pp. 437–443, Feb. 1985.
instantaneous frequency has been observed utilizing the [5] R. M. Adelson, "Frequency estimation from few measurements," Digital
proposed algorithm. Signal Processing, vol. 7, no. 1, pp. 47-54., Jan. 1997.
[6] M. D. Kusljevic, "A simple recursive algorithm for frequency
55 estimation," IEEE Trans. Instrum. Meas., vol. 53, no. 2, pp. 335-340.,
Apr. 2004.
54 Proposed [7] K. M. Wolter, Introduction to variance estimation, Springer, 2007.
Actual [8] Xu. Q. Q, Suonan. J. L, Ge. Y. Z. "Real-time measurement of mean
53 frequency in two-machine system during power swings," IEEE Trans.
Power Delivery, vol. 19, no. 3, pp.1018-1023, July 2004.
52 [9] M. A. Al-Alaoui, "Novel IIR differentiator from the simpson integration
rule," IEEE Trans. Circuits and Systems-I: Fundamental Theory and
Frequency (Hz)

51
Applications, vol. 41, no. 2, pp. 186-187, Feb. 1994.
50 [10] A. Sarkar and S. Sengupta, "Band-pass second degree digital integrator
based power system frequency estimation under non-sinusoidal
49 conditions," Accepted for publication in IEEE Trans. Instrum. Meas.
[11] B. Boashash, "Estimating and interpreting the instantaneous frequency
48 of a signal-Part 2: Algorithms and applications, " Proce. of IEEE, vol.
80, no. 4, pp. 540–568, Apr. 1992.
47

46 VII. BIOGRAPHIES
45
0.2 0.4 0.6 0.8 1 1.2 Arghya Sarkar (M’06) was born in West
Time (s) Bengal, India, on December 25, 1974. He
received the B.Sc. (Hons.) degree in physics and
the B.Tech. M.Tech. and Ph.D. degrees in
Fig. 7. Tracking of the polynomial phase signal. electrical engineering from the University of
Calcutta, Calcutta, India. He is currently an
V. CONCLUSIONS Associate Professor with the MCKV Institute of
Engineering, Howrah, India. His research
A novel digital signal processing algorithm for estimation interests are concerned with the application of
of the fundamental frequency of distorted power system digital methods to electrical power quality
t
signals have been presented and its performance is evaluated
Samarjit Sengupta (M’04, SM’10) received the
by means of simulation studies. High accuracy and B.Sc. (Hons.) degree in physics and the B.Tech.,
insensitivity to harmonics and noise, and fast response during M.Tech., and Ph.D. degrees in electrical
step parameter changes or instantaneous frequency tracking engineering from the University of Calcutta,
Calcutta, India. He is currently a Professor of
have been observed. Structural simplicity of the proposed electrical engineering with the Department of
estimator makes it suitable for digital implementation in both Applied Physics, University of Calcutta. His main
software environment, e.g., a DSP, and a digital hardware research interests include power quality
instrumentation, power system stability, and
environment, e.g., FPGA or ASIC. power system protection