3RD INTERNATIONAL CONFERENCE ON MODERN POWER SYSTEMS MPS 2010, 18-21 MAY 2010, CLUJ-NAPOCA, ROMANIA

Fundamental Frequency Estimation Using

Wavelet Denoising Techniques
D. Gheorghe, M. Chindriş, A. Cziker and R. B. Vasiliu

1
interharmonics, which strongly deteriorate the quality of the
Abstract—Harmonic and interharmonic analysis in power power supply voltage. Periodicity intervals in the presence of
systems are usually based on the nominal frequency of 50Hz or interharmonics can be very long. Parameter estimation of the
60Hz, which is an approximation of the fundamental frequency.
components is very important for control and protection tasks.
Many spectrum analyzers consider this frequency constant and
therefore, their measurement is not very accurate. This paper The design of harmonics filters relies on the measurement of
presents a new fundamental frequency estimation method based distortions in both current and voltage waveforms [1].
on an improved zero crossing algorithm. The electric signal is Digital control and protection of power systems require the
first filtered using the wavelet transform and the frequency is estimation of supply frequency and its variation in real-time.
then computed using a method that counts the number of zero Variations in system frequency from its normal value indicate
crossings in time using an adaptive window that detects the false the occurrence of a corrective action for its restoration. A
zero crossings. The paper continues with a virtual instrument large number of numerical methods is available for frequency
that computes the frequency from a sampled signal. Several
estimation from the digitized samples of the system voltage.
computer simulations test results are presented in the paper to
highlight the usefulness of this approach in estimating near Discrete Fourier transforms , Least error squares technique,
nominal power system frequencies. Kalman filtering , Recursive Newton-type algorithm , adaptive
notch filters etc. are known signal processing techniques used
Index Terms—Fundamental Frequency, Harmonics and for frequency measurements of power system signals [2].
Interharmonics, Virtual Instrument, Wavelet Denoising, Zero The real-time performance of a fundamental frequency
Crossings estimation algorithm depends not only on its computational
efficiency but also on its ability to obtain accurate estimates
I. INTRODUCTION from short signal segments.

S PECTRUM estimation of discretely sampled processes is The algorithm proposed in this paper is based on a Discrete
usually based on procedures employing the Fast Fourier Wavelet Transform filter that attenuates the high frequency
harmonics which would create false zero crossings near the
Transform (FFT). This approach is computationally efficient
real one, an adaptive window of search and an algorithm that
and produces reasonable results for a large class of signal
tracks the fundamental frequency and approximates it each
processes. However, there are several performance limitations
period.
of the FFT.
The most prominent limitation is that of frequency resolution, II. WAVELET FILTER DESIGN
i.e. the ability to distinguish the spectral responses of two or
more signals. A second limitation is caused by data A. Discrete Wavelet Transform
windowing, which manifests as leakage in the spectral The Wavelet Series is just a sampled version of Continuous
domain. These performance limitations are particularly Wavelet Transform and its computation may consume
troublesome when analyzing short data records, which occur significant amount of time and resources, depending on the
frequently in practice, because many measured processes are resolution required. The Discrete Wavelet Transform (DWT),
brief. which is based on sub-band coding is found to yield a fast
Modern frequency converters generate a wide spectrum of computation of the Wavelet Transform. It is easy to
harmonic components. Large converter systems and arc implement and reduces the computation time and resources
furnaces can also generate non-characteristic harmonics and required.
Filters are one of the most widely used signal processing
All authors are with the Electrical Power Systems Department, Technical functions. Wavelets can be realized by iteration of filters with
University of Cluj-Napoca, 15 C. Daicoviciu St., RO 400020, Cluj – Napoca. rescaling. The resolution of the signal, which is a measure of
D. Gheorghe and R.B. Vasiliu are PhD students (e-mails: the amount of detail information in the signal, is determined
[email protected], [email protected]), A.
Cziker is an associate professor (e-mail:[email protected]) and
by the filtering operations, and the scale is determined by
M. Chindriş is a full professor (e-mail: [email protected]). upsampling and downsampling (subsampling) operations[3].
The firs author D. Gheorghe received the License Degree in Electrical The DWT is computed by successive lowpass and highpass
Engineering from the Technical University of Cluj-Napoca, Cluj-Napoca, filtering of the discrete time-domain signal as shown in figure
Romania in 2009. Since 2009 he is working towards his Ph.D. at the Electrical
1. This is called the Mallat algorithm or Mallat-tree
Power Systems Department of Technical University of Cluj-Napoca. His
research project is concentrated on studying the power quality parameters. decomposition. Its significance is in the manner it connects

148
 3RD INTERNATIONAL CONFERENCE ON MODERN POWER SYSTEMS MPS 2010, 18-21 MAY 2010, CLUJ-NAPOCA, ROMANIA

the continuous-time mutiresolution to discrete-time filters. In B. Analytical Approach

the figure, the signal is denoted by the sequence x[n], where n Wavelets are similar to continuous wavelets, but the scale
is an integer. The low pass filter is denoted by G while the a and the location parameter b are measured in discrete
0
high pass filter is denoted by H . At each level, the high pass intervals. The Wavelet used must be Orthogonal translation on
0
filter produces detail information, d[n], while the low pass itself and in dilation, the discrete wavelet is defined as:
filter associated with scaling function produces coarse
⎛ t − nb a m ⎞
approximations, a[n]. 1 ⎜ 0 0 ⎟
ψ (t ) = ψ⎜ ⎟ (1)
m, n m m
a ⎜ a ⎟
0 ⎝ 0 ⎠
Discrete Wavelets involve a scaling function or “father
wavelets” this function must be orthogonal on translation and
orthogonal to other wavelets in its family. The scaling
function is defined as:
m
−
φ (t ) = 2 2 φ ⎛⎜ 2 − m t − n ⎞⎟ (2)
m, n ⎝ ⎠
Fig. 1. Three-level wavelet decomposition tree.
Using the inner product, the projection of a function onto
the wavelets is found. The same rule applies for the scaling
At each decomposition level, the half band filters produce
function. This is generated using the following formulas:
signals spanning only half the frequency band. This doubles
the frequency resolution as the uncertainty in frequency is ∞
reduced by half. In accordance with Nyquist’s rule if the T = ∫ x (t )ψ (t ) dt (3)
m, n n, m
original signal has a highest frequency of ω, which requires a −∞
sampling frequency of 2ω radians, then it now has a highest The discrete wavelet transform.
frequency of ω/2 radians. It can now be sampled at a ∞
frequency of ω radians thus discarding half the samples with S = ∫ x (t )φ (t ) dt (4)
no loss of information. This decimation by 2 halves the time m, n n, m
−∞
resolution as the entire signal is now represented by only half Discrete transform of scaling functions.
the number of samples. Thus, while the half band low pass
filtering removes half of the frequencies and thus halves the The wavelets and scaling functions are orthogonal and form
resolution, the decimation by 2 doubles the scale. a suitable basis. The function can now be reconstructed as a
With this approach, the time resolution becomes arbitrarily linear combination of wavelet and scaling function and their
good at high frequencies, while the frequency resolution corresponding coefficients:
becomes arbitrarily good at low frequencies. ∞ m0 ∞
The filtering and decimation process is continued until the x(t ) = ∑S φ
m ,n n , m (t ) + ∑ ∑T ψ n ,m (t )
m ,n (5)
desired level is reached. The maximum number of levels n = −∞ m = −∞ n = −∞

depends on the length of the signal. The DWT of the original The n indices represent the translation through time will the
signal is then obtained by concatenating all the coefficients, m indices represent the Dilation. The actual wavelets behave
a[n] and d[n], starting from the last level of decomposition. like a band pass filter and the scaling behaves like a low pass
filter.
∞
d (t ) = ∑ Tm,nψ n ,m (t ) (6)
m −∞
∞
a (t ) = ∑ S m,nφn ,m (t ) (7)
m −∞

Intuitively as m increases the wavelet is compressed and

fluctuates more rapidly, like a sine wave with increasing
frequency. The signal can be reconstructed in time with
Fig. 2. Three-level wavelet reconstruction tree. different levels of detail, hence the reconstructed function
called the m detail. The higher values of m dictate higher
Figure 2 shows the reconstruction of the original signal
levels of detail and much of the high frequency components.
from the wavelet coefficients. Basically, the reconstruction is
the reverse process of decomposition. The approximation and C. Case Study
detail coefficients at every level are upsampled by two, passed The following signal contains a fundamental component of
through the low pass and high pass synthesis filters and then 50hz and two harmonics, a low frequency harmonic (rank 3)
added. This process is continued through the same number of and a high frequency one (rank 25).
levels as in the decomposition process to obtain the original ⎛ π⎞
f (t ) = 230 sin(ωt ) + 23 sin ⎜ 3ωt + ⎟ + 17 sin(25ωt ) (8)
signal. ⎝ 2⎠

149
 3RD INTERNATIONAL CONFERENCE ON MODERN POWER SYSTEMS MPS 2010, 18-21 MAY 2010, CLUJ-NAPOCA, ROMANIA

This signal will be sampled at 512 points and it will be

300
stored in vector F:
0.02
n = 9 i = 0...2 n − 1 ti = n i Fi = f (ti ) (9)

Magnitude [V]
200
2
100
200 Voltage Waveform
Magnitude [V]

0
− 100
0 200 400

Time [s]
− 200
Coefficient 2 + 250V
0 0.005 0.01 0.015 0.02 Coefficient 4 +100V
Coefficient 5
Time [s]

Fig. 3. Voltage Waveform Fig. 6. Discrete Wavelet Transform coefficients (2, 4 and 5)

260 The signal is now decomposed in 10 levels of coefficients,

230
200
each representing a frequency band. The first two elements, 0
Magnitude [V]

170 and 1, are called approximation coefficients. The remaining

140 elements are the detail coefficients. The transform DWT
110 contains 8 levels of detail. The last 256 entries represent
80 information at the smallest scale, the preceding 128 entries
50
represent a scale twice as large, and so on.
20
− 10 The parameter scale in the wavelet analysis is similar to the
− 40 scale used in maps. As in the case of maps, high scales
−1 2 5 8 11 14 17 20 23 26 29
correspond to a non-detailed global view (of the signal), and
low scales correspond to a detailed view. Similarly, in terms
Rank
of frequency, low frequencies (high scales) correspond to a
Fig. 4. Magnitude spectrum of the original voltage waveform global information of a signal (that usually spans the entire
signal), whereas high frequencies (low scales) correspond to a
The discrete wavelet transform of function f(t): detailed information of a hidden pattern in the signal (that
x + 2n usually lasts a relatively short time).
DWT ( f , n, x) = ∫ f (t )W (t − x, n)dt (10) The large scale coefficients will contain the fundamental
frequency and the low frequency harmonics and the high
x −2n
The mother wavelet chosen is Daubechies 14: frequency harmonics will be stored in the low scale
coefficients. Therefore, if the low scale coefficients are
truncated, the remaining coefficients will contain only low
frequencies. The signal is reconstructed by taking the Inverse
Discrete Wavelet Transform of the remaining coefficients.

Original Waveform
200 Filtered Waveform
Magnitude [V]

Fig. 5. Daubechies 14 wavelet and scaling function

− 200

Where S is the scaling function and W is the wavelet 0 0.005 0.01 0.015 0.02

function. Time [s]

The transform coefficients are computed :
p = 0...n − 2 A = DWT (11)
p, i p +1 ⎡ i ⎤ Fig. 7. Original voltage waveform and filtered voltage waveform
2 +⎢ ⎥
n − p −1 ⎥
⎣⎢ 2 ⎦ Magnitude spectrum of the reconstructed signal taken with
the Discrete Fourier Transform:

150
 3RD INTERNATIONAL CONFERENCE ON MODERN POWER SYSTEMS MPS 2010, 18-21 MAY 2010, CLUJ-NAPOCA, ROMANIA

260 The rectangular search window is dimensioned after the

230 first zero crossing when the width and the translation are
200
Magnitude [V]

defined. This window will be then shifted along the whole

170
140 signal.
110 Multiple zero crossing can occur on a single period duration
80 leading to major measuring errors. Usually these problems are
50 caused by harmonics and the most difficult ones to detect are
20
the ones very close to the fundamental zero crossing point.
− 10
− 40 Luckily, if the “false” zero crossing are very close to the
−1 2 5 8 11 14 17 20 23 26 29 “real” crossing points, the harmonic frequency that produces
them is also very high. Since the harmonic frequency is high,
Rank it will not pass through the low pass wavelet filter. If the false
Fig. 8. Magnitude spectrum of the filtered waveform zero crossings are further than the real ones, they are produced
by low rank harmonics with high amplitudes. These effects
III. FREQUENCY ESTIMATION METHOD are solved by adjusting a rectangular window of search around
the expected fundamental frequency value. Furthermore, the
The zero crossing detection algorithm is based on the sign
window of search will considerably reduce the number of
difference between two consecutive samples. When a group of
iterations per period since the samples outside the window are
ti and ti+1 samples are found, a linear interpolation is
not being read.
performed in order to find an approximate coordinate in time
when the signal value is 0.
IV. VIRTUAL INSTRUMENT
Zero crossing condition:
sign(ti ) + (sign(ti+1 ) − 2 ) = 0 A virtual instrument has been implemented with the help of
2
(12)
LabView. The current or voltage waveform can be either
Linear interpolation method : inputted as an array of samples or as an analytical expression.
The instrument requires the width of the search window, the
total analysis time (if the waveform is an analytical
expression) and the sampling frequency.
The instrument returns a figure of the frequency
approximation in time, the frequency being approximated
every cycle, the frequency deviation from the real frequency
(analytical case), the total integer cycles time, the number of
detected false zero crossings and the total analysis duration.
The input waveform and the wavelet filtered waveform are
also displayed.
Fig. 9. Linear interpolation The analytical test waveform contains a 50hz fundamental
frequency, a low frequency interharmonic (rank 3.3) and a
f (ti ) ⋅ (ti +1 − t i ) high frequency interharmonic (rank 25.7), each component
ti 0 = ti − (13)
f (ti +1 ) − f (t i ) with a phase shift of 0.36 rad, 0.7 rad and 1.6 rad.
After the first zero crossing, the axes are translated into that
location in time so the first non integer cycle is truncated. The
frequency estimation is determined by dividing the number of
zero crosses counted and the total duration of the integer
cycles.
number of crosses
f = (14)
total time

Fig. 11. Virtual instrument front panel. The wavelet filter is activated.
Fig. 10. First zero crossing and the adaptive search window

151
 3RD INTERNATIONAL CONFERENCE ON MODERN POWER SYSTEMS MPS 2010, 18-21 MAY 2010, CLUJ-NAPOCA, ROMANIA

The signal has been sampled over 10 seconds, which

corresponds to 500 cycles of 20ms each, with a sampling
frequency of 6400 samples/second.
The total numbers of samples was 64000 but the program
performed only 5163 iterations. This reduced calculation
complexity is given by the rectangular search window. The
window slides along the sampled signal and the virtual
instrument reads only the samples that are inside the window.

TABLE I
FREQUENCY ESTIMATION RESULTS

Input Estimated Input Estimated

Error Error
Frequency Frequency Frequency Frequency
[mHz] [mHz]
[Hz] [Hz] [Hz] [Hz]
42.5 42.497 2.731 50.5 50.496 3.245
43 42.997 2.763 51 50.996 3.277
43.5 43.497 2.795 51.5 51.496 3.309
44 43.997 2.827 52 51.996 3.342
44.5 44.497 2.860 52.5 52.496 3.374
45 44.997 2.892 53 52.996 3.406 Fig. 13. Virtual instrument front panel. Wavelet filter disabled.
45.5 45.497 2.924 53.5 53.496 3.438
46 45.497 2.956 54 53.996 3.470 The second test has been performed with the same input
46.5 46.497 2.988 54.5 54.496 3.502
dates but without the wavelet filter in order to point out the
47 46.996 3.020 55 54.996 3.534
47.5 47.496 3.052 55.5 55.496 3.566 difference in error.
48 47.996 3.084 56 55.996 3.599 The measurement error without the wavelet filter for a
48.5 48.496 3.117 56.5 56.496 3.631 50Hz input signal is -95.34mHz
49 48.996 3.149 57 56.996 3.663
49.5 49.496 3.181 57.5 57.496 3.695 V. CONCLUSIONS
50 49.996 3.213 58 57.996 3.727
In the last decades power quality estimation has become a
very important issue in power systems so the necessity of
accurate parameters measurement has grown in importance
too.
The algorithm presented and the results obtained are
satisfactory. The new signal processing techniques, like the
discrete wavelet transform offered important accuracy
improvements and solved most of the false zero crossings
instances.
The paper describes briefly the wavelet denoising process
which can be useful in many other applications in electrical
engineering.
The virtual instrument developed performed well in high
harmonic polluted conditions. The total analysis duration is
Fig. 12. Measurement error on the 42.5Hz – 57.5Hz domain less than the total waveform duration so the instrument can be
used online.
The figure above presents a set of measurements on the
interval 42.5Hz to 58Hz, comparing the analytically input VI. REFERENCES
frequency and the estimated frequency. The frequency
deviation is calculated as the difference between the estimated [1] C. Mattavelli, “Analysis of Interharmonics in DC Arc Furnace
frequency and the real frequency. As seen from the table, the Installations”, 8th Int. Conf. on Harmonics and Quality of Power”,
maximum error for the analyzed domain is 3.72mHz. Athens, Greece, vol.2, pp.1092-1099, 1998.
According to the 61000-4-30 standard, the instrument [2] P. K. Dash, "Frequency Estimation of Distorted Power System
complies with class A equipment where the maximum Signals Using Extended Complex Kalaman Filter", IEEE Transactions
on Power Delivery, Vol. 14, No. 3, July 1999
allowable deviation for a 10 seconds analysis is ±10mHz [5].
[3] C. Taswell, "The What, How, and Why of Wavelet Shrinkage
The total analysis duration for a 10 seconds waveform is Denoising," Computing In Science And Engineering, vol. 2, no. 3,
about 5 seconds so the instrument is able to perform a real May/June 2000, p. 12-19.
time frequency measurement. [4] Golovanov, Carmen ş.a. „Metode moderne de măsurare în
electroenergetică.” Bucharest: Editura Tehnică, 2001.
[5] IEC 61000-4-30 Ed.2: Electromagnetic compatibility (EMC) – Part
4-30: Testing and measurement techniques –Power quality measurement
methods , April 2007

152