**1. Introduction**

The fast evolution of electrical and electronic technologies is generating an exponential growth in the use of non-linear loads, such as switch-mode power supplies, electronic lighting ballasts [1], and other harmonic sources whose penetration can be expected to increase substantially in the future [2]. Harmonic distortions coming from new-generation sources and loads are generally larger and less regular than those due to traditional sources and loads, making power and energy more difficult to measure [3,4]. International standards have been developed to define acceptable limits and suitable measurements techniques for power systems' harmonics, e.g., [5,6]. Fourier spectral analysis is the most widely used approach for the evaluation of the harmonic content in power systems, due to its excellent accuracy under stationary conditions. This choice is reflected in the IEC 61000-4-7 [5] standard, which defines harmonics measurements and instrumentation requirements. This standard proposes the use of the Fast Fourier Transform (FFT), a fast and efficient algorithm for performing the Discrete Fourier Transform (DFT). Although the DFT is valid under general conditions, it is well known that limitations arise when it is applied to non-stationary signals, which is the normal situation in modern power systems. When fluctuating harmonics are analyzed with a DFT, the energy contained in the harmonics is distributed across the spectrum, affecting the accuracy [7]. For this reason, the IEC 61000-4-7 standard proposes a grouping strategy in order to improve the accuracy in non-stationary conditions and, at the same time, explicitly mentions the possibility to employ more advanced processing tools. During the last few years, alternatives to DFT have been proposed in the literature. A common workaround is to use the Short Time Fourier Transform (STFT), which assumes the signal to be piecewise stationary, but other methods have been suggested as well [8]. Among them, there is Wavelet Transform (WT), typically in the form of Discrete Wavelet Transform (DWT). The main advantage of DWT is that, while DFT requires a stationary, perfectly periodic sinusoidal signal to work properly, this is not a requirement for DWT. This signal processing tool has a wide range of applications, including power systems, where it proved to be extremely useful in signal denoising, short time predictions, fault detection, and energy managemen<sup>t</sup> [9–14], as well as for Power Quality [15–19]. From WT, Wavelet Packet Decomposition (WPD) was derived, generalizing the link between wavelets and Multiresolution Analysis (MRA). Thanks to the uniform frequency bands that can be obtained with WPD, its application in Power Quality has been proposed for power calculation [20,21] and harmonic measurements [22–24]. However, the WPD harmonic analyses proposed in the literature do not look complete from the point of view of compliance with international standards and capabilities under stationary conditions, as will be discussed in Section 2.

To complete the current scenario, this paper presents a WPD method for a full harmonic analysis, in contrast with the valuable, but partial work developed in recent years, e.g., [22–25]. The proposed method takes advantage of the superior performance of WPD with fluctuating harmonics with respect to Fourier analysis. On the other hand, it is for the first time compliant with the IEC 61000-4-7 standard. This result, never achieved before, is possible thanks to the accurate selection of the filtering implementation and to the uniform frequency bandwidth division of WPD, as studied by the authors in [26]. To assess the method's implementation and validate its superior performance within the IEC standard requirements, several voltage signals with known harmonic content, both in stationary and fluctuating conditions, were tested. The accurate assessment of individual harmonic values was the main focus of the proposed method, and therefore, they were the objective of the validation. Other indices, like for instance Total Harmonic Distortion (THD), were not considered since they are the result of a further processing stage that can be performed in the same way as for FFT, as explained in Section 2.2. In Section 4.4, a comparison of the computational effort is provided, and finally, in Section 5, the method is tested with real voltage and current signals measured in a high harmonic distortion environment.

### **2. The Proposed WPD Method**

### *2.1. Wavelet Transforms in Power Systems*

WT is a mathematical transform that provides a representation of a function by a series of orthonormal functions, generated by scaling and translating a wavelet. It is particularly powerful in the measurement of the time-frequency variations of spectral components, but not suitable for power harmonics analysis because of the non-uniform frequency bands obtained. WPD is an extension of WT where the frequency axis is divided into equally wide intervals, instead of the traditional division in dyadic intervals, whose sizes have an exponential growth [27]. The advantage is that, with an accurate selection of the sampling frequency and the decomposition structure, it can lead to uniform frequency bands of the desired size and can therefore be employed for harmonics analysis. Previous works on WPD based algorithms for harmonics analysis can be found in the literature, but up to this moment, no IEC compliant approach, with complete harmonic analysis, has been proposed. Hamid et al. presented in [21] a WPD method for Root Mean Squared (RMS) and power measurements, from which harmonic distortion can be obtained. However, with this approach, only odd harmonics can be evaluated, and the wide frequency bands that can be obtained are not compliant with IEC standards. Eren et al. proposed a decomposition scheme in [22], but again, only bands centered in odd harmonics were taken into account. Moreover, results were shown only up to the seventh order. A more refined strategy was proposed by Diego and Barros in [24], where the 50 Hz wide frequency bands were compatible with IEC 61000-4-7 [5], but the decomposition was performed only up to the 15th order, far from compliance with the IEC standards.

Additionally, a major drawback of WPD based algorithms is their difficulty to compete with FFT in terms of computational speed. The FFT algorithm is the result of many years of developments and advances of Fourier analysis, while newly proposed methods (such as WPD) are not ye<sup>t</sup> optimized for high performances. However, this gap is reducing, thanks to the constant increase in available computational power and the possibility of hardware implementation. In [23], it was shown that WPD can be efficiently implemented and be employed online for harmonics analysis. However, although the work showed the feasibility of online measurements, the same limitations on the results were found, i.e., the evaluation of only odd harmonics, limited total bandwidth (only up to 15th order), and too wide frequency bands (even harmonics contributions were included in the odd ones).

WPD algorithms continue to be interesting for power system researchers, as suggested by recent works [20,23], but in order to move forward and make WPD a feasible and better alternative to DFT, a complete algorithm for harmonic assessment must be designed. Modern power systems, with the increasing presence of non-linear loads, require a method compatible with the IEC standards, able to measure not only odd harmonics, but also even harmonics, and up to the highest possible order.

### *2.2. Characteristics of the Proposed Method*

The past limitations of WPD arose from the difficulty of adjusting the algorithm to the constraints given in the IEC 61000-4-7 standard, which is not trivial. Among these constraints, there are: a fixed rectangular sampling window of 10 cycles of the fundamental (200 ms at 50 Hz), frequency bands of 50 Hz centered in integer multiples of the fundamental (both even and odd), a bandwidth up to the 50th harmonic order (2500 Hz), and mainly, the fulfillment of the accuracy requirements for the stationary case. Therefore, the implementation of an IEC compliant algorithm is not just an extension of previous studies, since it requires overcoming the previous limitations and keeping the accuracy within the prescribed limits.

In this paper, it is shown that this can be achieved by an accurate study of the WPD implementation. The novel algorithm hereafter presented is the continuation of the work performed by the authors in [26], where a detailed study was conducted to define a methodology for the selection of the most suitable filter and convolution technique for WPD applied to power systems' harmonics. The methodology proposed in [26] is here employed to develop and characterize a complete measurement method. In order to obtain the IEC frequency bands, up to at least the 50th order, a sampling rate of 6400 Hz is necessary, and seven levels of decomposition must be implemented. These necessary features set high requirements on the frequency selectivity of the decomposition filters. The work performed in [26] showed that a maximally-flat Butterworth Infinite Impulse Response (IIR) filter of order 29 was the most suitable choice, and for this reason, it was employed in this paper as the mother wavelet. Its superior behavior in frequency selectivity, number of coefficients, and spectral leakage under realistic conditions, compared with traditional mother wavelet functions (e.g., db20, db10, Vaidyanathan, etc.) provides the high level of accuracy required to comply with the IEC standards in both stationary and non-stationary conditions. It must be noted that, although IIR filters have a non-linear phase, it does not represent a problem because the purpose of the method is to extract the RMS content of each frequency band, and the convolution technique presented in [26] allows preserving this information even with non-linear phase filtering.

In the following, the filtering procedure at each node of the decomposition tree is explained: the input signal *f* is recursively convolved by a pair of Quadrature Mirror Filters (QMFs), high pass and low pass filters (*h*(*n*) and *g*(*n*), respectively), and then, the filtered signals are downsampled by a factor of two, taking into account that the employed filters are half band filters (see Equations (1) and (2)). Thus, the high pass and low pass QMFs used were obtained from the corresponding IIR Butterworth filter of order 29, as discussed before. This WPD scheme was designed and implemented using MATLAB, allowing the authors to speed up the validation and testing procedure under realistic

input conditions, i.e., with the evaluation of real voltage waveforms. The following equations describe the basic filtering stage of the WPD method, in which the wavelet coefficients at level *l*, node *i*, are obtained from the convolution of the input signal at the previous level with the filters' coefficients and subsequent downsampling.

$$D\_1^{2i}(k) = \sum\_{n} g(n) D\_{1-1}^i(2k - n) \tag{1}$$

$$D\_1^{2i+1}(k) = \stackrel{\cdots}{\sum\_n} h(n) D\_{1-1}^i(2k - n) \tag{2}$$

where *i* = 0, 1, ... , 2(*l* − 1) − 1, and *h*(*n*) and *g*(*n*) are the high pass and low pass filter coefficients, respectively. In this regard, Figure 1 shows the decomposition tree scheme employed.

**Figure 1.** Wavelet Packet Decomposition (WPD) filtering and decimation scheme.

Table 1 summarizes the characteristics of the decomposition structure at each node, including information about the number of nodes per level, samples per node, effective time resolution, and bandwidth for every node at each level.


**Table 1.** Decomposition tree summary.

It is possible to observe that nodes at the seventh decomposition level have a bandwidth of 25 Hz. The required 50 Hz bands were obtained by a two-nodal-grouping (7\*, last level in Table 1). Aggregated RMS harmonic values are derived from individual node content as follows [15]:

$$\chi\_{\text{rms}}(j, p+q) = \sqrt{\frac{\sum\_{k} \left(d\_{j,k}^{p}\right)^{2} + \sum\_{k} \left(d\_{j,k}^{q}\right)^{2}}{N}} \tag{3}$$

where *p* and *q* are the two grouped nodes, *N* the total number of data, *k* the index counter, and *j* the level at which the RMS calculation is performed. In order to obtain the correct IEC bands, the grouping must satisfy: *q* = *p* + 1, *p* = 2, 4, 6, ... (the first node is discarded).

Thanks to this two-nodal-grouping, each band corresponded exactly to an IEC harmonic group, and therefore, any further processing stage, e.g., THD calculation or time aggregation according to IEC 61000-4-30, is perfectly possible. The algorithm was indeed designed to represent an alternative to FFT for power quality analyzers, assuming sampled values as input and producing compliant individual harmonic values as output. As a result, the method did not have any additional requirements than those of FFT.

### **3. Compliance with IEC Accuracy Requirements**

This section presents a test of accuracy, showing that the proposed method fulfills the IEC precision requirements for the stationary case and overperforms the IEC standard method for the non-stationary case.

IEC 61000-4-7 prescribes, for Class I accuracy, a maximum error of 5% in the evaluation of harmonic components, in case of single frequency and steady-state signals. In order to test the compliance of the proposed method, fifty steady-state signals *<sup>u</sup>*i(*t*), of single harmonic frequency ranging from harmonic order one to 50, were generated with a MATLAB script and subsequently analyzed using both the IEC method and the proposed WPD method. The input signals have the following expression:

$$u\_i(t) = 100\sqrt{2}\sin\left(i\,2\pi f\_0 t\right) \tag{4}$$

where *i* is the *i*th test containing the *i*th pure tone, *f*0 is the power frequency (50 Hz), and *t* is the time.

The top graph in Figure 2 shows the errors obtained with the proposed method for each i*th* harmonic order, with respect to the expected values, which in this case was always 100. It can be seen that the errors were well below the limit for Class I instruments for every harmonic order. The slightly higher errors for order 31st and 33rd were due to an intrinsic characteristic of WPD, which affected the central decomposition bands, as explained in [26]. This, however, did not compromise the accuracy, and compliance with IEC 61000-4-7 was proven.

The same test was repeated for non-stationary single frequency signals. To obtain non-stationary signals, the fluctuation pattern proposed in IEC 61000-4-7 Annex C.3 was employed, i.e., a reduction of the RMS of the signal from 100% to 20% at *t* = 85 ms in a 200 ms window. Since FFT is not suitable for measuring fluctuating harmonics, the standard proposes a grouping strategy, where the FFT bins are summed into harmonic groups to reduce the uncertainty. However, no accuracy limits are prescribed for this scenario. The bottom graph in Figure 2 shows the errors obtained in the case of fluctuating signals using the conventional FFT, the grouped IEC strategy, and the proposed WPD method. It can be seen that the conventional FFT produced errors of the order of 20%, while the IEC grouping strategy allowed a reduction of the error below the 5% limit. In this case, the proposed WPD method further reduced the error, for every order but 31 and 33, never exceeding the IEC stationary-state limits. The proposed method was therefore equally valid as the grouped FFT for single frequency signals, both stationary and fluctuating, and significantly superior to standard FFT, as expected.

**Figure 2. Top**: errors of the proposed method with steady-state single frequency signals. **Bottom**: comparison of the error produced by the conventional FFT, the grouped IEC, and the WPD methods in case of single frequency signals, under fluctuating conditions (as explained in the text). The grey lines show the limits provided by IEC 61000-4-7 for stationary conditions.

### **4. Validation of the Method**

The validation of the method required a more complex harmonic content than single frequency signals. Moreover, in order to be able to assess the accuracy of the method and compare it with the accuracy of Fourier analysis, the true harmonic content of the signals needed to be known a priori, i.e., synthetic signals must be employed. For this reason, realistic signals were generated, starting from the National Physics Laboratory (NPL) Power Quality Waveform Library [28], a public repository of waveforms with known harmonic content, representative of voltage signals found in the grid or used for calibration of equipment. This approach allowed obtaining waveforms whose true harmonic content was known and, at the same time, was realistic harmonic content.
