A Robust Wavelet-Based Hybrid Method for the Simultaneous Measurement of Harmonic and Supraharmonic Distortion

This article introduces a new method for the simultaneous measurement of harmonic and supraharmonic distortion with improved robustness against amplitude and power frequency deviations. The proposed algorithm is based on wavelet analysis, and it is designed to analyze a ten-cycle measurement interval, something crucial when implementing methods to measure harmonic and supraharmonic content according to the International Electrotechnical Commission (IEC) 61000-4-7 standard. Discrete Fourier transform (DFT)-based methods, instead, must analyze a fixed 200-ms interval in order to avoid apparent shifts of the supraharmonic frequency components. The results of the experimental measurements presented in the article show that the DFT-based methodology suggested in the IEC 61000-4-7 standard produces the results affected by the value of the power frequency and amplitude, while the proposed method, instead, is insensitive to them. This feature also removes the need of a previous high-pass filtering stage, required for the IEC method. Moreover, the proposed method has the additional advantage of working on the same measurement interval to calculate both harmonics and supraharmonics, which reduces the complexity of data handling by avoiding multithreaded data acquisition (DAQ) operations.


I. INTRODUCTION
The evolution of smart grids, with distributed generation, enhanced communication features and modern power converters poses new challenges, such as the measurement of new Power Quality (PQ) phenomena [1]. Among them, the change in the harmonic distortion features is attracting the interest of the researchers and calls for new measurement methods and instruments able to measure the increasingly present distortion in the range from 2 kHz to 150 kHz, known as supraharmonics [2]. Moreover, the performance of new measurement methods must be carefully assessed under realistic situations, and the results must be comparable with the applicable IEC standards limits.
In principle, the measurement of supraharmonics could take advantage of the traditional characterization tools, but the reality is different as the higher frequency distortion is intrinsically different from the harmonic distortion. The harmonic emission is mainly made of waves whose frequency is an integer multiple of the fundamental frequency which, in power systems, is the power frequency (nominally 50 Hz or 60 Hz, although in this paper only 50 Hz systems are considered). The term supraharmonics, instead, refers to the emission beyond (supra) the harmonic range, but not necessarily multiple of the power frequency.
Currently, there is no standardized method for the measurement of grid disturbance levels in the supraharmonic range. Few suggestions are informatively described in the relevant standards [3], all of them based on the DFT, but no agreement has been reached so far. Moreover, these methods come with some drawbacks such as the unsuitability for field measurements or the need of employing a different measurement interval than for the harmonic range, which makes the DFT algorithms sensitive to fluctuations of the power frequency.
Moreover, the use of a different measurement interval for harmonics and supraharmonics analysis (10 cycles and 200 ms) requires flexible power quality analyzers, with two different processing threads, which increases their design complexity and, therefore, price. Moreover, the DFT shows important limitations when applied to fluctuating (non-stationary) signals, which is the normal situation in modern power systems. When non-stationary harmonics are analyzed with the DFT, severe energy leakage can be observed, producing large errors [4].
Extending the work presented by the authors in [5], this paper proposes a hybrid measurement method for harmonic and supraharmonics analysis based on the Wavelet Packet Decomposition (WPD). The method is compliant with the IEC 61000-4-7 standard for harmonics under stationary conditions and shows a better performance under non-stationary conditions. Moreover, it is comparable with the supraharmonic methods described in IEC 61000-4-7 and CISPR 16-1, but with the advantage of working with a 10 cycles measurement interval (synchronized measurement) without being affected by power frequency deviations. These features are accurately assessed through experimental measurements, proving that the method shows improved robustness against typical deviations from ideal conditions that occur in real smart grids. Finally, the proposed method does not require any pre-processing stage, while the IEC analysis should be preceded by a High Pass Filter (HPF) in order to achieve good accuracy under power frequency deviations. The rest of the paper is structured as follows: Section II presents the IEC standard measurement method for harmonics and reviews the existing methods for supraharmonic analysis. In Section III, the proposed method is described in detail, and it is then applied to real signals from the grid, in order to assess its performance under real conditions. Then, Section IV systematically studies the advantages of the proposed method when dealing with specific grid issues i.e., non-stationary conditions and deviations of the power frequency. Finally, the conclusions are presented in Section V.

A. Methods Existing in Standards
Harmonic distortion analysis, through DFT, was defined in the IEC 61000-4-7 standard [6] several years ago. The Fast Fourier Transform (FFT), a fast and efficient algorithm for performing the DFT, was the chosen tool for its superior performance under stationary conditions. However, this tool shows spectral leakage across the analysis bands when nonstationary harmonics are analyzed. In this regard, the IEC 61000-4-7 standard suggests to make an aggregation of bands to improve the accuracy under these working conditions, but leaves the door open for the adoption of any other kind of strategy compatible with the IEC requirements. In this case, the use of a wavelet-based measurement method can provide more accurate results compared to the grouping strategy proposed in the IEC 61000-4-7 standard, as it will be shown in Section III-B. The challenge is to maintain the compatibility with the IEC standard i.e., using a fixed rectangular sampling window of 10 power line cycles, frequency decomposition bands of 50 Hz centered at even and odd integer harmonics (for the harmonic groups), a measurement bandwidth up to at least 2000 Hz, and deviation results for stationary analysis below the limits specified in the IEC 61000-4-7 standard [6].
On the other hand, the field of supraharmonic distortion analysis still lacks a commonly accepted measurement method [2]. A lot of research has being carried out and the development of international standards is still ongoing [7]. Currently, the informative Annex C of the IEC 61000-4-30 standard on power quality measurement methods describes three different approaches. The first one is the only one suggested for the 2 kHz-9 kHz range and coexists with the other two options for the 9 kHz-150 kHz range. The first option (2 kHz-150 kHz) is the extension of the same measurement technique described in IEC 61000-4-7 i.e., DFT, up to 150 kHz. This method employs a 200 ms rectangular measurement interval. The obtained spectral components are then grouped into 200 Hz bands, centered in odd multiples of 100 Hz, via the Root Sum Square (RSS) [6]. The second option is to use the method of CISPR 16-1 [8], based on measurements performed in the frequency domain with a scanning receiver. Using a super-heterodyne analyzer, a narrowband filter is tuned on a 200 Hz frequency band and the CISPR peak, quasi-peak, and rms values are measured. The full bandwidth is covered by shifting the filter step-wise. Finally, the third option is a new approach proposed in IEC 61000-4-30, again based on DFT analysis, but with some differences to the IEC 61000-4-7 method. An analysis interval of 10 cycles is considered. Within this interval, 32 approximately equally spaced non-overlapping measurements are taken, collecting 512 samples per measurement, using a rectangular window. Since the specified sampling frequency is 1024 kHz, each interval has a duration of 0.5 ms. The 32 acquired intervals are then processed using DFT, resulting in 2 kHz-wide bands. For every 10 cycles interval, minimum, maximum and average values are recorded. This last method is less accurate (wider bands) and less complete (more than 90 % of each analysis interval is ignored) than the previous ones, but it has the advantage of requiring less data storage and being less expensive.
It must be noted that all IEC methods based on DFT have fixed frequency resolution, as the measurement intervals are not synchronized with the power frequency. This is an important difference with the IEC harmonics measurement method, which prescribes a 10 cycles, synchronized measurement interval. The reason is that the use of a synchronized interval is especially beneficial to the DFT for the harmonic analysis. By measuring exactly 10 cycles, indeed, the DFT resolution always matches the power frequency and its harmonic components. In fact, even if the power frequency (and consequently all its harmonics) deviates from its nominal value, the DFT resolution will change accordingly, continuing to match all the harmonic components. However, the supraharmonic distortion is independent of the power frequency, as it typically originates from devices whose operation is not correlated with the power frequency, like power electronic converters [9]. As a consequence, in case of a power frequency shift, the DFT resolution would shift accordingly, as said, but the supraharmonic components will not shift, as they are not linked to the power frequency. The resulting effect would be an apparent shift of the supraharmonic emission with respect to whole spectrum. In order to avoid this effect, the IEC 61000-4-7 standard suggests the use of a fixed 200 ms measurement interval, instead of the traditional 10 cycles interval [6], so that the frequency resolution is fixed and no apparent shift is observed. This solves the problem of the apparent shift, but creates another issue. The DFT is based on the assumption that the analyzed signal is infinitely periodic or, at least, strictly periodic inside the measurement interval. However, in case of a deviation of the power frequency from the nominal value of 50 Hz, a fixed 200 ms interval will not contain exactly 10 fundamental cycles anymore. The periodicity assumption of the DFT is therefore not valid, resulting in energy leakage across the frequency spectrum, especially of the fundamental component, which has the highest magnitude. As a possible mitigation strategy, the use of a previous filtering stage is suggested in the standard, although the exact filter specifications are not provided. In [10] a 3 rd order elliptic digital filter is proposed, while in [11] a 10 th order elliptic digital filter is proposed and implemented with a zero-phase filtering technique. Both HPFs have been used in this work when implementing the IEC measurement method, along with the unfiltered implementation, in order to provide a complete comparison. More details on the measurement methods can be found in [7], [12].

B. Other Existing Methods
Very few other measurement methods have been recently proposed in the literature for the supraharmonic analysis. In [13] a method based on analog filter banks was proposed with the purpose of reducing the sampling capability requirements of PQ analyzers, but at the same time increasing the requirements of the analog components. In [14], a small fraction of the supraharmonic range (up to only 10 kHz) was analyzed with a WPD approach. However, a frequency resolution of 1250 Hz and a measurement interval of (fixed) 80 ms were proposed, not comparable with the current standards. Multiple measurement vectors methods have also been suggested based on Baesyan Compressive Sensing and Orthogonal Matching Pursuit in [15], [16] as a strategy to improve the outputs of the DFT analysis and provide a sparse representation of the signal. The challenge, however, is represented by the lack of a priori knowledge about the sparsity of the signal which might affect the accuracy.

III. THE PROPOSED METHOD A. Description of the Method
The WPD is a wavelet decomposition, introduced by Coifman, Meyer and Wickerhauser [17], that generalizes the link between wavelets and Multiresolution Analysis (MRA). The WPD is particularly suitable for general harmonic analysis since an accurate selection of the decomposition structure can lead to uniform frequency bands of the desired size. For this reason, it has been proposed in the past as a tool for harmonics analysis, being especially effective under fluctuating conditions, where the basic requirements of the DFT (linearity and stationarity) do not hold [18], [19].
In this paper, a hybrid measurement method for the simultaneous assessment of both harmonic and supraharmonic distortion, based on WPD, is presented. The method is based on the work done by the authors in [5], extending the measurement capabilities to the harmonic range. By addressing both harmonic and supraharmonic analysis needs, the proposed method is not only suitable for the analysis of the whole frequency range (from 50 Hz to 150 kHz), but also provides a methodology compliant with the IEC standard for harmonic groups and comparable with the IEC and CISPR suggestions for supraharmonics.
Firstly, a rectangular window of 10 cycles of the power frequency is selected as measurement interval to be decomposed using the WPD. Then, similarly to the implementation proposed in [18], [20], at every decomposition step the signal is filtered at half of its bandwidth using a Quadrature Mirror Filter (QMF) pair, and then downsampled by a factor of 2, recursively. Therefore, at every step the number of nodes is doubled, and their bandwidth is halved, recursively. The original fullbandwidth signal represents the starting node (level 0). With an initial sampling frequency of 409 600 Hz (i.e. maximum theoretical bandwidth Bw = 204 800 Hz), the decomposition tree allows to obtain, for each node at level 10, a bandwidth of 200 Hz i.e., the same resolution as the IEC 61000-4-7 and CISPR 16-1 methods for the frequency range above 2 kHz. Taking advantage of the results of the supraharmonic analysis, a further decomposition of only 3 additional levels, and only for the frequenyc range 0 Hz to 3200 Hz, is performed, allowing to obtain the resolution required by the IEC 61000-4-7 standard for the harmonic range.
The resulting decomposition tree is schematically shown in Fig. 1. The green area identifies the supraharmonic calcula-   tion, while the red area shows the additional decomposition needed for the harmonic range. The characteristics of the decomposition structure are presented in Table I, which reports information about the number of nodes per level, samples per node, effective time resolution ∆t, and bandwidth Bw for every node at each level. It can be seen that, as an example, at level 10 there are 1024 nodes with a 200 Hz bandwidth each. Moreover, each node, after the iterative subsampling, contains 80 samples, which corresponds to a time resolution of ∆t = 2.5 ms. Level 13* in Table I describes the results obtained after the two-nodal-grouping, as explained in the following. The wavelet coefficients D at level j, node p, are obtained from the convolution of the input signal at the previous level with the filters coefficients and subsequent downsampling.
where p = 0, 1, . . . , 2 (j−1) − 1 identifies the nodes, j the level, and h(n) and g(n) are the low-pass and high-pass filter coefficients, respectively. The employed filters are a pair of two maximally-flat Infinite Impulse Response (IIR) Butterworth QMFs of order 29, selected following the procedure described in [18]. The filtering is performed employing the Zero-Phase Filtering (ZPF) technique described in [21], which offers the advantage of keeping phase information, although in this specific work it is not strictly necessary. For the supraharmonic analysis, the nodes of the WPD at level j = 10 contain the energy of the associated frequency intervals and can be used for measuring the rms amplitude as follows: where p is the node with a bandwidth of 200 Hz, N the total number of data, k the index counter, and j = 10 the level at which the supraharmonic assessment is performed.
On the other hand, as the nodes at the 13 th decomposition level have a limited bandwidth of just 25 Hz, the required 50 Hz bands were obtained by a two-nodal-grouping (level 13* in Table I) for a IEC compliant harmonic evaluation. According to [22], aggregated rms harmonic values can be obtained from individual node contents as follows: where p and q are the two nodes with a bandwidth of 25 Hz, N the total number of data, k the index counter, and j = 13 the level at which the rms calculation is performed. According to [18], in order to obtain the correct IEC grouped values, the grouping must satisfy: q = p + 1, p = 2, 4, 6, ... (the first node is discarded). These values can be compared directly with the rms grouped values (not subgrouped) obtained from the DFT analysis proposed in IEC standards. Beside being suitable for analyzing fluctuating signals, another valuable feature of the proposed WPD method is the use of a 10 cycles measurement interval, while the IEC method relies on a fixed 200 ms interval for supraharmonics which, as discussed, results in large energy leakage. For this work, the 10 cycles interval was obtained by detecting the zero crossing of the fundamental component, achieving a precision of 2 µs on the interval length. This feature removes the need of using different measurement intervals for low (harmonic) and high (supraharmonic) frequency distortion and has the potential to reduce the complexity of PQ instruments by only requiring one analysis thread, as shown in Fig. 2. Moreover, in Section III-B it will be shown that the proposed method does not require the initial filtering stage that is required by the IEC methods. The decomposition algorithm, as well as the methodology for rms calculations, were implemented in Matlab which eases the validation process. The sampling frequency has been kept the same for both methods, equal to 409 600 Hz, to meet the requirements of the proposed method and fulfill the Nyquist-Shannon theorem.
More information about the characteristics of the decomposition method as well as its general implementation and characterization can be found in the recent works carried out by the authors in [5], [20].

B. Assessment of The Proposed Method under Real Conditions
In this section, the performance of the proposed method is briefly illustrated, showing that the results match the performance of the DFT method presented in IEC 61000-4-7, with very similar results (when the IEC 61000-4-7 method has proper high pass filtering implemented) for both supraharmonic and harmonic evaluations. The frequency spectra were obtained after sampling real voltage waveforms, measured in the grid, of an Electric Vehicle (EV) charger and a Photovoltaic (PV) inverter, which are typical sources of distortion. The signal acquisition was performed by a National Instruments PXIe 6124 DAQ card, with a Pico TA044 70 MHz 700 V differential voltage probe at 1 MHz. The signal was resampled to meet the required sampling frequency of the proposed method, i.e. 409 600 Hz, and then analyzed as described in Section II.
Among the three options described in II-A for supraharmonic evaluation, this paper only considers the IEC 61000-4-7 method (here called IEC method for the sake of brevity) for comparison purposes because of its full signal coverage, frequency resolution and simplicity. The second method (CISPR method) was developed to test equipment emission in laboratory conditions and, as stated in the IEC 61000-4-30 standard, it is considered too complex and expensive for in-situ PQ measurements. The third option (IEC 61000-4-30 method), instead, only covers the 8 % of the signal and with worse frequency resolution, making it unsuitable for comparison. Fig. 3(a) shows the supraharmonic voltage emission of an EV charger connected to the grid whereas Fig. 3(b) includes the results regarding harmonic evaluations. The same comparison is presented in Fig. 4 for the harmonic and supraharmonic emissions of a PV inverter. Only for the supraharmonic comparison, and in order to better highlight the differences, the spectra are presented using the logarithmic unit dB µV , which is defined using 1 µV as reference level (it follows that 0 dB µV correspond to 1 µV and each decade is represented by 20 dB). It can be seen from Fig. 3(a) and 4(a) that the WPD performance is very similar to that of the IEC method with a previous filtering stage (and using grouped harmonic values). The IEC method without filtering, instead, has a worse performance, especially due to higher background levels, caused by energy leakage. The reason was found to be a deviation of the power frequency from its nominal value. It was indeed measured that the power frequency was 49.92 Hz in the PV case and 50.02 Hz in the case of the EV. This is not an uncommon situation in the grid, since small deviations are allowed. According to EN 50160 the power frequency must stay in the range between 49.5 Hz to 50.5 Hz during 99.5 % of a year for grids with synchronous connection to an interconnected system, and between 49 Hz to 51 Hz during 95 % of a week for grids with no synchronous connection to an interconnected system [23].
With regard to the harmonic analysis, Fig. 3(b) and 4(b) show the rms amplitude for the two considered methods (IEC and WPD). Although there are no limits specified in the IEC 61000-4-7 standard for this scenario (non-stationary conditions), slight differences are observed between the methods when looking at the first harmonic orders due to an excessive spectral leakage caused by the fluctuating nature of the waveforms. A more detailed characterization study is presented in the following sections through the evaluation of synthetic signals.

IV. PERFORMANCE OF THE PROPOSED METHOD UNDER NON-STATIONARY CONDITIONS
This section presents the analysis of the robustness of the proposed method under controlled amplitude fluctuations and power frequency deviations, and compares it to the IEC methods based on DFT for harmonic and supraharmonic analysis.

A. Metrics
In order to quantify the differences between two spectra for the harmonic evaluation, a specific index -called Root Mean Squared Error (RMSE) -is employed: where ε h is the error produced for order h, and N is the total number of bands, 40 in this case. This metric is normally used for statistical evaluations, however it can also be used to quantify the variation of errors observed on a certain test set, as it is the case. As it can be seen, the higher the RMSE, the lower the overall accuracy for the evaluated method. On the other hand, for the supraharmonic analysis, a different metric is employed. As it will be shown in section IV, the main effect of the power frequency deviation is observed in the background. Therefore, in order to accurately assess its content, giving the same importance to every bin of the spectrum, the Mean Spectra Difference (MSD) is defined as: where x i=1,...,N are the N bins of the first spectrum, while y i=1,...,N are the N bins of the second spectrum, taken as reference. Therefore, the obtained value expresses the mean difference between the two spectra per bin, in percentage. In this work, the reference spectrum was the one calculated at nominal power frequency 50 Hz.

B. Harmonic Performance under Non-Stationary Conditions
The sensitivity analysis of the robustness under changes in amplitude was carried out through the evaluation of signals with a known content of modelled harmonics which allowed the calculation of deviations. For this reason, one realistic signal, extracted from the National Physics Laboratory (NPL) PQ waveform library, was selected as the reference for this case study. Concretely, the selected waveform was the "IEC 61000-3-2 voltage limits for the power amplifier", a signal reflected in the Annex 1 of the IEC 61000-3-2 standard that is used for limiting the voltage supply distortion of power amplifiers on harmonic emission assessments. This signal is composed by a total of 38 harmonics with a Total Harmonic Distortion (THD) of 1.26 % which make it ideal for this benchmark exercise.
For this non-stationary study, a fluctuating pattern was applied to the selected waveform. A linear modulation, similar to a motor start, was applied to the last 115 ms of the 200 mslong signal. Concretely, the amplitude ranged from 1 p.u at instant 85 ms to an end value changed between 0.01 p.u and 2 p.u, in 0.05 p.u steps, which covers most of the situations that could be found in the grid (from voltage interruption levels, 1 %, to swell amplitudes, 200 %). Therefore, the nonstationary (fluctuating) character is obtained by changing the amplitude of the entire waveform, while the harmonic content is, in percentage, not changed.
After building up the waveform test bed and decomposing the signals with both IEC and WPD methods, the RMSE values were calculated. In this regard, Fig. 5 shows a plot of the results obtained with both methods. It can be seen that the performance of the IEC method is clearly affected by the variations of the amplitude within the analysis interval whereas the proposed method is insensitive to them. RMSE values close to 1.5 V at very low and high residual voltages (0.01 p.u and 2 p.u respectively) indicates that the FFT method suffers from very high inaccuracies when analysing important voltage fluctuations (this would correspond to 4 % if expressed with respect to the signal). When measuring waveforms with a medium voltage fluctuation level (1 ± 0.2 p.u), the FFT method reports RMSE values of about 0.23 V. On the other hand, the WPD method remains very stable under any level of voltage fluctuation, giving an average RMSE of 0.07 V and a maximum error of 0.11 V at the highest voltage fluctuation level (which would correspond to 0.2 % and 0.19 % of the signal, respectively). Only at values of residual voltage in the surroundings of 1 p.u (1±0.05 p.u), the FFT method performed slightly better than the WPD, as the voltage fluctuation can be considered negligible (similar to a non-fluctuating condition).
Currently, the IEC 61000-4-30 standard indicates that PQ measurement algorithms might produce an unreliable value in case of events with residual voltages from 0 p.u to 0.9 p.u (interruptions and dips) or larger than 1.1 p.u (swells). Given the results shown in Fig. 5, the WPD method performs accurately beyond the IEC limits so, if the proposed method is used, the flagging concept for harmonic estimation could be reconsidered for what it concerns the accuracy of the harmonic measurement. In any case, for a maximum residual voltage value of 0.9 p.u, the WPD method performs better when compared to the IEC method. In order to further support this, a deep study was carried out not only changing the amplitude but also the power frequency. Concretely, the residual voltage was varied linearly from 0.01 p.u to 2 p.u (with 0.05 p.u steps), whereas the power frequency ranged from 49 Hz to 51 Hz with 0.1 Hz steps. As a result, the plot in Fig. 6 shows the performance of both methods when dealing with non-stationary conditions as it might happen on real PQ campaigns. It can be noted that the IEC method reported RMSE values very close to the ones shown in Fig. 5 for any power frequency while the WPD remained stable, giving a very flat response over the entire test. For this reason, it can be concluded that the proposed method overcomes the current limitations of IEC method when measuring harmonics under fluctuating conditions.

C. Supraharmonic Performance under Frequency variations
The analysis is the result of a set of experimental measurements, taken at the point of connection of an EV charging system (different than those presented in Section III-B). The measurements were performed under controlled laboratory conditions, initially at nominal power frequency 50 Hz. In order to assess the robustness of the methods, the measurements were repeated with the same configuration but different values of power frequency f i , ranging from 49.0 Hz to 51.0 Hz, with 0.1 Hz steps. The results were then compared with the reference case of nominal power frequency and the deviations were quantified. The power supply, a pure sine wave at frequency f i , was generated using a Spitzenberger & Spies PAS 15000, with a voltage output accuracy < ±0.1 %, from DC to 15 kHz (−3 dB). For every value of the power frequency, a voltage waveform was acquired using a Dewe 2600 measurement system and its voltage module (Dewetron HSI-HV), sampling at 1 MS/s and setting the measurement range to ±400 V. The reported accuracy for this range is given by ± (0.016f ) ± 0.1 % of range, where f is the signal frequency in Hz. The measurements resulted therefore in 21 spectra for each measurement method.  On the other hand, the proposed method provides more stable results. High leakage is observed in the case of the DFT, as expected, since the analysis interval is not a multiple of the fundamental cycle. It can be noted that the WPD spectrum does not show any shift of the supraharmonic components when using a 10 cycles measurement interval. This confirms the discussion presented in Section II and the need of a HPF stage for the IEC method.
In order to quantitatively assess the variation between the results obtained at different power frequencies and the reference case, three approaches were employed: (i) by measuring the MSDs between the spectra obtained at frequency f i and the spectrum obtained at 50 Hz; (ii) by comparing the magnitude of the peaks located at the switching frequencies, for the different values of the power frequency; (iii) by measuring the MSDs between the spectra obtained at frequency f i and the spectrum obtained at 50 Hz, but excluding the peaks at the switching frequencies. This last approach allows to assess the variation of the leakage produced in the background.
The MSD between the spectra for every value of f i is shown in Fig. 8. The plot also includes the results obtained with the IEC method with the implementation of two different HPFs, according to [9] (identified by letter K) and [10] (identified by letter D). For the sake of clarity, the results of the IEC method without filtering are not plotted due to the large errors produced, one order of magnitude larger than the IEC method with the filtering stage. It can be observed that the results of the IEC method are largely affected by the power frequency. The larger the frequency deviation, the larger the difference between the spectra. It must be noted that the very high percentage values are due to the very low background values obtained in the reference case (50 Hz), which is at the denominator in Eq. 6. Moreover, it can be observed that the type of HPF produces little differences in the results. The proposed WPD method is the only one that does not show any correlation with the Another possibility to assess the variation is to only look at the peaks corresponding to the switching frequencies and their multiples. As shown in Fig. 7, peaks at 10 kHz and multiples were detected, suggesting that 10 kHz is the switching frequency of the EV charging system. The peak at 20 kHz, i.e. the second harmonic of the switching frequency, was selected for analysis and its magnitude was compared with the magnitude of the same peak under different values of the power frequency, in order to track its variation. It must be noted that the two adjacent bands are included in the calculation (total bandwidth of 600 Hz) in order to cover the whole energy of the emission band, as suggested in [9]. Fig. 9(a) shows the magnitude of these peaks at varying power frequency. In this case, the IEC method (not shown) still produces large, frequency dependent errors. The filtering stage, however, is able to remove this effect, producing results that are independent of the power frequency. The proposed WPD method is also able to provide results with no correlation to the power frequency, and the results are comparable with those of the IEC method with a previous HPF stage.
The last approach to assess the variation is complementary to the previous one. In this case the focus is on the variation of every frequency component except those related to the switching frequencies i.e., the background content. In order to do so, the MSD is calculated as in Eq. 6, but excluding the switching frequencies and their multiples from the calculation. Similarly to the procedure described before, the adjacent bands are excluded as well, as they contain the energy of the switching frequencies. The results are shown in Fig. 9(b) and they can be interpreted as the average background level. Any difference is due to energy leakage. It can be observed that the IEC method, even with the previous filtering stage, produces spectra whose background level is affected by leakage and largely dependent on the power frequency. The realization of a filtering stage before the analysis can attenuate this effect but cannot remove it. The proposed WPD method, instead, shows lower background levels for any power frequency variation.

V. CONCLUSIONS
The paper presented a novel WPD algorithm for the simultaneous measurement of both harmonic and supraharmonic contents. The proposed method analyzes a fixed 10 cycles interval synchronized with the power frequency, which has the potential of simplifying measurement equipment, and shows overall robustness against voltage fluctuations and power frequency deviations, which is an important drawback of the DFT-based methods proposed in IEC standards. The robustness has been compared with the IEC 61000-4-7 method through theoretical and experimental measurements. The results of the proposed method are constant and independent of the voltage amplitude modulation and the power frequency value. As expected, the IEC method performs worse, producing results that are largely affected by leakage and correlated with the value of the voltage amplitude and the power frequency. It has also been shown that the implementation of a previous filtering stage, as suggested by the IEC standard, improves the behavior of the method for supraharmonic estimations, especially for the peak values, but does not remove the dependence on the power frequency of the background values. This could represent an issue in case of low Signal-to-Noise Ratio (SNR). Regarding harmonic measurements, it was proved that the WPD proposed method is not affected by low residual voltages during any voltage event type. This characteristic makes the method suitable for a continuous monitoring of the harmonics regardless of the amplitude value. This could lead to modify the flagging limits specified in the IEC 61000-4-30 standard or even remove the harmonic mesurement from the list of algorithms affected by low residual voltages. Thus, the improved robustness of this implementation makes the proposed algorithm suitable for environments with high voltage and frequency volatility, such as weak grids or islanded operations.