CZT-Based Harmonic Analysis in Smart Grid Using Low-Cost Electronic Measurement Boards

Abstract

This paper validates the use of a harmonic analysis algorithm on a microcontroller to perform measurements of non-stationary signals in the context of smart grids. The increasing presence of electronic devices such as inverters of distributed generators (DG), power converters of charging stations for electric vehicles, etc. can drain non-stationary currents during their operation. A classical fast Fourier transform (FFT) algorithm may not have sufficient spectral resolution for the evaluation of harmonics and inter-harmonics. Thus, in this paper, the implementation of a chirp-Z transform (CZT) algorithm is suggested, which has a spectral resolution independent from the observation window. The CZT is implemented on a low-cost commercial microcontroller, and the absolute error is evaluated with respect to the same algorithm implemented in the LabVIEW environment. The results of the tests show that the CZT implementation on a low-cost microcontroller allows for accurate measurement results, demonstrating the feasibility of reliable harmonic analysis measurements even in non-stationary conditions on smart grids.

1. Introduction

Today’s power grid is fully automated, easy-to-access, reliable, and efficient services, thanks to the introduction of smart grids and their continuous development. Therefore, the evolution of appropriate measurement systems supported by a reliable communication infrastructure is the basis for the development of smart grids to increase remote service monitoring, such as power quality analysis, anti-islanding detection, arc fault detection, and so on. In this context, the number of devices connected to the grid has increased dramatically, and the proliferation of distributed generation (DG) from renewable sources has encouraged users to become active by installing power plants that can independently power their utilities [1,2,3,4,5,6,7]. In this context, power systems have also become more complex, exhibiting time-varying harmonic behaviour due to the increasing non-linearity of electronic loads. For example, variable speed motor drives cause increased harmonic pollution, by generating and feeding harmonics and inter-harmonics into distribution. The integration of DGs also has negative effects on the total harmonic distortion at the point of connection (POC) [8,9,10,11]. Furthermore, it is crucial to consider the rising number of electric vehicles (EV) and what this means in terms of harmonic generation during charging phases, which are sometimes concentrated in the same time slots. B. Grasel et al. [12] investigated EV chargers and their harmonic patterns, evaluating their potential impact on the grid.

Power grids have not been designed to operate with unpredictable renewable energy sources, and this is causing several issues with power systems that can no longer be ignored. For this reason, it is necessary to think of an advanced metering architecture solution that allows the distributor to know in real-time the status of the network and to take actions to avoid instability and unwanted voltage or frequency variations. The observability of the network can be obtained only by increasing both local and distributed measurement capabilities [13]. In this framework, the authors are developing a research project named SInERT (“Solutions innovantes pour l’intégration des énergies renouvelables sur le réseau électrique tunisien”), funded by the European Union in the framework of the Italy-Tunisia cross-border cooperation program 2014–2020. The project aims at developing innovative solutions for integrating distributed generator and storage systems into distribution networks. To achieve this aim, a distributed measurement architecture is proposed [14]. It allows monitoring different quantities: active and reactive power flows, voltage and frequency (also for islanding detection), harmonics, and other power quality phenomena. A demonstrator of the proposed solutions is under construction in the Elgazala Technopark (Ariana, Tunisia). A photovoltaic power plant with storage will be installed using innovative power converters. They will be remotely controlled by a supervisory control and data acquisition (SCADA) system, which will also collect monitoring quantities both from power converters and from distributed measurement systems installed at different points of the low-voltage (LV) network. A picture of the pilot site architecture is shown in Figure 1, where the presence of the metering devices is highlighted. They will be able to perform different measurement, including those related to advanced power quality and harmonic analysis. In the framework of a distributed knowledge of power quality issues, these measurements systems should have advanced functionalities with the aim of estimating harmonic pollution, detecting disturbing sources, etc. [15,16,17,18]. Thus, one of the goals of the project is to investigate the feasibility of new algorithms for harmonic analysis.

Performing power quality measurements means having the possibility to measure the spectral quantities that characterize the acquired signal. Consequently, several mathematical methods have been studied and developed to perform harmonic analysis and monitor the harmonic variation in power systems. In very brief, for these applications, they can be essentially divided into time-domain and frequency-domain methods [19]. Some of them are preferably proposed for power system measurements and related applications. DFT (Discrete Fourier Transform) or FFT (Fast Fourier Transform) methods are the most commonly used for power quality analysis [20,21], and they are specifically addressed in the Standard IEC 61000-4-7, which is devoted to harmonics and inter-harmonics measurements and instrumentation for power systems and equipment connected thereto. Apart from this, other methods have been provided, such as filter banks or digital filtering techniques [22,23], wavelet-based methods [24,25], compressive sensing [26], STFT (Short Time Fourier Transform) [27], and sliding-window recursive-discrete Fourier transform (SWR-DFT) [28,29]. Each of these methods has its own advantages and disadvantages in terms of accuracy, computational effort, and so on. For example, the SWR-DFT uses the Fourier transform and a sliding window to compute time-varying harmonics. Consecutive decompositions are carried out until all harmonics have been fully extracted in this iterative technique. This approach has a low computational cost, no phase delay, and a short transient time. On the other hand, all DFT-based methods can suffer from spectral leakage problems, poor resolution with short observation windows, and they are not suitable when dealing with nonstationary signals. In more detail, this class of methods can be used to break down a complex signal into its harmonic components until it exhibits elements of stationarity. In fact, harmonic analysis algorithms that rely on Fourier harmonic decomposition require the signal to remain stationary throughout the duration of the observation. Moreover, the frequency resolution in algorithms such as DFT, or derivatives, depends on the observation window: $f_{R_{FFT}}=1/T_w$ where Tw is the length of the temporal observation window. For example, to have a resolution of 1 Hz, an observation window of 1 s is required, but due to the variability of the loads and the unpredictability with which the signal can become non-stationary, the signal is unlikely to remain stationary for long periods. The signal’s non-stationarity makes any frequency domain transformation impossible. To face this issue, it is possible to use a different harmonic analysis technique, which can allow obtaining a resolution independent from the observation window; in this way, the observation window can be properly reduced to better cope with nonstationary conditions without losing information on the signal. In this approach, a temporal interval in which the signal remains stationary can be identified, and a frequency domain transformation can be performed. A suitable algorithm that can work with short observation windows is a particular type of Z-transform, i.e., the Chirp-Z Transform (CZT) [30,31,32]. In fact, the main advantage of the CZT algorithm is related to the possibility of reaching a better resolution (if compared with the DFT/FFT algorithms) in a specified frequency range, even with small observation windows.

In this framework, the study presented in this paper is aimed at investigating the feasibility of using the CZT algorithm for harmonic analysis, from the viewpoint of both its accuracy in harmonic measurement, its capability of operating with nonstationary signals, and the possibility of implementing the algorithm on smart metering platforms. To test the performance in a real-world case study, the algorithm has been implemented on the STMicroelectronics NUCLEO—STM32H723ZG development board [33]. This board was chosen because it is similar to that inside the electronic meters normally used in distribution power systems, such as those addressed in the SInERT project or even used in most smart metering and household applications around the world. For example, the power meters shown in Figure 1 represent the measurement nodes whose metrological features could be improved by the algorithm developed in this paper.

The paper is organized as follows: in Section 2, the CZT algorithm is presented and analysed; in Section 3, the CZT algorithm on board implementation is explained; and the experimental results are finally reported in Section 4.

2. Case Study Algorithm: CZT

Harmonic analysis aims to perform a signal transformation from the time domain to the frequency domain in order to obtain the spectrum of the input signal and evaluate the possible presence of harmonics. In power quality applications, this type of analysis is particularly important because it allows for the evaluation of significant parameters such as total harmonic distortion (THD), unwanted harmonic presence, harmonic pollution, etc. As mentioned in Section 1, signals can be stationary or not. Depending on the stationarity of the underlying signal, one spectral analysis algorithm could be preferred to another for its performance.

In fact, for FFT or DFT, the spectral resolution is inversely proportional to the observation window:

$$f_{R_{FFT}}=\frac{1}{T_w}$$

(1)

where $f_{R_{FFT}}$ is the spectral resolution or the bin size, and $T_w$ is the observation window.

To improve the spectral resolution, a longer observation window is needed. For example, to obtain a spectral resolution of 5 Hz, $T_w$ must be 200 ms; if a better spectral resolution is needed, for example, equal to 1 Hz, the observation window must be enlarged to 1 s. If such spectral resolution is required, it is necessary to assume that the signal remains stationary for the entire duration of the observation window. On the other hand, if the signal is nonstationary, the DFT or FFT algorithm cannot be used unless a shorter observation window is used, in which case the assumption of signal stationarity can still be considered valid; however, the use of a DF or FFT algorithm with a short observation window may not result in a suitable frequency resolution. In this case, other spectral analysis algorithms could give better performances.

In this paper, the feasibility of the chirp Z-transform (CZT) algorithm is investigated because it does not entail the relationship of Equation (1) between the spectral resolution and the observation time. In the following, the basics of the CZT algorithm are summarized, putting into evidence its features in terms of spectral resolution improvement with respect to the DFT/FFT algorithm.

2.1. CZT Fundamentals

The Z-chirp transform, also known as the chirp Z-transform (CZT), is a transform that allows for information to be transformed from the time to frequency domains. The CZT algorithm offers a higher resolution than FFT, with flexibility in choosing the frequency bandwidth. Generally, in the field of digital signal analysis, the most commonly used algorithm is the FFT to identify tones, frequencies, harmonics, and study the frequency behaviour of a given signal. To understand the main difference between CZT and FFT, let’s consider the FFT and analyse the differences on equal terms. In the case of FFT, the frequency range goes from 0 Hz to fs/2 where fs is the sampling frequency with equally spaced samples; this involves a constant distance between two successive samples (Figure 2).

So, if a signal that has two harmonics close to each other is considered, the spectral resolution would be fundamental for their identification. If the resolution is not sufficient, a single lobe will be displayed, as shown in Figure 3a. The figure shows an input signal consisting of two ideal Dirac deltas at frequencies between two “bin sizes”, whose frequency signal is in the band of interest, so the transform approximates the two deltas in a single lobe due to poor resolution. In the FFT algorithm, the frequency resolution depends on the observation window; the greater the observation window, the better the frequency resolution, as their ratio is inversely proportional (1). If the resolution is not adequate, it may happen that the signal is not displayed correctly, and the transform shown in Figure 3a can be obtained. Indeed, the two close Dirac deltas appear inside the main lobe, so it is impossible to discriminate between them. Now if a CZT is applied to the same signal, the two Dirac deltas could be more efficiently captured, and the two lobes can be seen as shown in Figure 3b. Of course, not even with the CZT can we see two ideal deltas, but the increase in resolution permits us to identify two lobes instead of one.

2.2. CZT Algorithm

To understand how it is possible to increase the resolution using the CZT, the goniometric circle must be used. In more details, starting from a sequence x[n] of N discrete elements, the Fourier transform of the input signal can be written as $X\left(e^{j\omega }\right)$. In Figure 4, a representation of M elements among the $\left(e^{j\omega }\right)$ is shown in the goniometric circle of unit radius.

The Fourier transform equally distributes the points on the goniometric circle. The distance between one point and the next represents the spectral resolution. Starting from A(1;0), it is possible to consider:

$${\omega }_k={\omega }_0+k\Delta \omega , k=0, 1, \dots , M-1$$

(2)

where ${\omega }_0$ is the starting frequency and $\Delta \omega$ is the increment. For DFT and FFT, ${\omega }_0=0$, M = N, $\Delta \omega =2\pi /N$. The index k represents the k-th element of the output sequence. Starting from the following expression of DFT:

$$X\left(e^{j{\omega }_k}\right)={\displaystyle \sum }_{n=0}^{N-1}x\left[n\right]e^{-j{\omega }_{k}n}, k=0, 1, \dots , M-1$$

(3)

and substituting (2) into (3), CZT is obtained in (4).

$$X\left(e^{j{\omega }_k}\right)={\displaystyle \sum }_{n=0}^{N-1}x\left[n\right]e^{-j{\omega }_{0}n}{e}^{-j\Delta \omega kn}, k=0, 1, \dots , M-1.$$

(4)

If W is defined as the following:

$$W=e^{-j\Delta \omega }$$

(5)

then, Equation (4) can be written as (6),

$$X\left(e^{j{\omega }_k}\right)={\displaystyle \sum }_{n=0}^{N-1}x\left[n\right]e^{-j{\omega }_{0}n}{W}^{nk}.$$

(6)

Considering $nk=\frac{1}{2}\left[n^2+k^2-{\left(k-n\right)}^2\right]$, (6) can be written as indicated in (7),

$$X\left(e^{j{\omega }_k}\right)={\displaystyle \sum }_{n=0}^{N-1}x\left[n\right]e^{-j{\omega }_{0}n}{W}^{n^2/2}{W}^{k^2/2}{W}^{-{\left(n-k\right)}^2/2}.$$

(7)

If we denote $g\left[n\right]=x\left[n\right]e^{-j{\omega }_{0}n}{W}^{n^2/2}$, (7) changes into (8).

$$X\left(e^{j{\omega }_k}\right)=W^{k^2/2}\left({\displaystyle \sum }_{n=0}^{N-1}g\left[n\right]W^{-{\left(n-k\right)}^2/2}\right),k=0, 1, \dots , M-1$$

(8)

Equation (8) shows that $X\left(e^{j{\omega }_k}\right)$ elements correspond to the convolution between g[n] samples and the terms $W^{-{\left(n-k\right)}^2/2}$ followed by the product with $W^{k^2/2}$. Starting from these considerations, (8) can be represented as shown in Figure 5.

Considering (4), the frequency band of interest can be selected by choosing $\Delta \omega =2\pi \Delta f$ value, so frequency resolution can be written as follows:

$$f_{R_{CZT}}=\frac{\Delta f}{N}.$$

(9)

2.3. Example of Comparison between CZT and FFT

To appreciate the differences between FFT and CZT algorithms, a simple numerical example can be provided. It refers to the test case conditions of the experimental results reported in the following, i.e., the processing of signals acquired with a number of samples equal to 800, a sampling frequency $f_s$ = 10 kHz, and a frequency band of interest between 0 Hz and 500 Hz. In the case of the FFT, the resolution is given by (10); for the CZT, the resolution is given by (11), instead.

$$f_{R_{FFT}}=\frac{1}{T_w}=\frac{f_s}{N}=\frac{10000}{800}=12.5 \mathrm{Hz}$$

(10)

$$f_{R_{CZT}}=\frac{\Delta f}{N}=\frac{500}{800}=0.625 \mathrm{Hz}$$

(11)

A sampling frequency of 10 kHz and 800 points correspond to an observation window of 80 ms (four cycles of the fundamental power system frequency of 50 Hz). The spectral resolution that we obtain with the FFT is 12.5 Hz. The CZT resolution does not depend on the observation window, allowing for a higher resolution of 0.625 Hz with the same number of samples.

This leads to a different capability of detecting near spectral components between the two algorithms, as shown in Figure 6. In the figure, the FFT and CZT spectra are compared for a 50 Hz signal with near spectral components (sampling conditions of the previous example, i.e., $f_s$ = 10 kHz, number of samples equal to 800, frequency band of interest between 0 Hz and 500 Hz). From the figure, it can be observed that the frequency components are distinguishable in the CZT plot while they are not in the FFT plot.

3. Architecture and Hardware Implementation

The CZT algorithm has been implemented on a microcontroller development board, the NUCLEO-STM32H723ZG. This kind of platform has been chosen, instead of others, such as FPGA, as it is similar to those typically used in electronic smart meters used in power distribution systems. In fact, FPGA boards are typically more expensive than microcontroller boards, and they are not mostly used in common commercial devices for smart metering applications, including power quality monitoring and analysis. Even if FPGA increase the computational capabilities and be very useful for many other applications, in the aim of this paper was to investigate the feasibility of implementing power quality metrics (with proper accuracy and resolution features even with nonstationary signals) on low cost existing devices, to enable the large-scale use of such devices in different smart metering applications, such as those related to distributed generation and/or storage systems, EV charging stations, and so on.

3.1. CZT Architecture

The signal is windowed to reduce scallop loss and leakage errors. Precisely, a Hanning window is chosen and thus expressed as follows:

$$w\left(n\right)=0.5\left[1-cos\left(\frac{2\pi n}{N}\right)\right],0\le n\le N$$

(12)

The retained sampling frequency $f_s$ is 10 kHz. The latter value is compatible with commercial devices frequencies. An observation window $T_w$ of 80 ms was chosen. Considering these parameters, the number of points $N$ for the CZT algorithm is 800. The frequency interval $\Delta f$ was set to 500 Hz. Hence, it is possible to summarize the parameters as follows:

$$\Delta f=500 \mathrm{Hz}$$

(13)

$$f_s=10 \mathrm{kHz}$$

(14)

$$T_w=80 \mathrm{ms}$$

(15)

$$N=f_s·T_w=800$$

(16)

Furthermore, a number of output samples equal to the number of input samples was chosen, i.e., k = n.

Equation (6) can be written as (17):

$$X\left(k\right)={\displaystyle \sum }_{n=0}^{N-1}x\left[n\right]A^{-n}{W}^{nk}$$

(17)

where $A=e^{-j{\omega }_{0}n}$. Considering module and phase, (17) can be written as indicated in (18).

$$X\left(k\right)={\displaystyle \sum }_{n=0}^{N-1}x\left[n\right]·\left(A_0^{-n}{e}^{-j{\theta }_{0}n}\right)·{\left(W_0^{-n}{e}^{-j\Delta \omega n}\right)}^k.$$

(18)

If the starting coordinate point (1; 0) of the goniometric circle is followed by the path along the circumference, the radius will be constant and equal to 1. Furthermore, we set the terms $A_0=1$, ${\theta }_0=0$, $W_0=1$, and we obtain (19).

$$X\left(k\right)={\displaystyle \sum }_{n=0}^{N-1}x\left[n\right]e^{-j\Delta \omega nk}$$

(19)

It is possible to write $\Delta \omega$ as (20),

$$\Delta \omega =\frac{2\pi \Delta f}{N{f}_s}$$

(20)

So, (19) becomes (21),

$$X\left(k\right)={\displaystyle \sum }_{n=0}^{N-1}x\left[n\right]e^{-j\frac{2\pi \Delta f}{N{f}_s}nk}$$

(21)

Consider $n·k=\frac{1}{2}\left[n^2+k^2-{\left(k-n\right)}^2\right],$ (21) can be written as (22).

$$X\left(k\right)={\displaystyle \sum }_{n=0}^{N-1}x\left[n\right]e^{-j\frac{2\pi \Delta f}{N{f}_s}·\frac{1}{2}\left[n^2+k^2-{\left(k-n\right)}^2\right]}.$$

(22)

Because $\frac{\pi \Delta f}{N{f}_s}$ is a constant, it can be defined as a parameter $C=\frac{\pi \Delta f}{N{f}_s}$, then (22) becomes (23).

$$X\left(k\right)=e^{-jC{k}^2}{\displaystyle \sum }_{n=0}^{N-1}\left\{x\left[n\right]e^{-jC{n}^2}\right\}·e^{jC{\left(k-n\right)}^2}.$$

(23)

So, the CZT implementation can be summarized in four “stages”, which can be drawn as in the block diagram shown in Figure 7. A computational cost estimation can be done by observing Figure 7. Stage 1 is composed of Hanning register multiplication, the input signal, and a complex register. This stage consists of N complex operations, meaning four multiplications and two sums. Stage 2 consists of k × n multiplications, so the total number depends on k. In our algorithm, we chose k = n, so $N^2$ operations will be executed. Stage 3 is composed of N sums, while Stage 4 is similar to Stage 1 in terms of the number of operations, that is N complex operations will be performed. This means that Stage 2 operates with the largest number of operations—that is $N^2$. Thus, the final computational cost can be compared to the DFT one.

On the basis of the aforementioned definition relating to the CZT algorithm, the system input-output can be summarized as shown in Figure 8. Particularly, the system inputs are: the system clock (CLK), the system reset (which clears both outputs and internal registers), and two information inputs that can be used either for two real signals or for one complex signal. Finally, the system outputs are the CZT real and imaginary parts.

3.2. Implementation of CZT on Nucleo STM32H723ZG

The STM32H7 Nucleo-144 development kit is an affordable and easy-to-use board designed for rapid development and evaluation of STM32H7 microcontrollers. It features a powerful Arm Cortex-M7 processor, which offers exceptional performance while maintaining energy efficiency and user friendliness. The Cortex-M7 processor has a range of features, including a six-stage dual-issue pipeline, dynamic branch prediction, Harvard architecture with L1 caches, a 64-bit AXI interface, and two 32-bit DTCM interfaces. It also supports various memory interfaces, including separate instruction and data buses, a tightly coupled memory (TCM) interface, an AXI bus interface, and a dedicated low-latency AHB-Lite peripheral bus.

This MCU has 1 Mbyte of flash and 564 kB of RAM, and it supports a set of DSP instructions, enabling efficient signal processing and complex algorithm execution. It also supports single and double-precision floating-point units (FPU) to speed up software development.

For onboard code implementation, STM32CubeMX and STM32CubeIDE were used for pin initialization, code generation, debugging, and download. STM32CubeMX is a graphical tool for easy configuration of STM32 microcontrollers, while STM32CubeIDE is used for onboard code debugging and download. CMSIS DSP 1.9.0 libraries were employed to perform operations in FPU, and there are functions for basic mathematical calculations and operations with complex numbers. Overall, the STM32H7 Nucleo-144 development kit is an excellent choice for various application areas, including industrial automation, medical devices, high-end audio, image, and voice processing, sensor fusion, and motor control. A complex number can be written in Euler’s formula as follows:

$$e^{\pm jx}=cosx\pm jsinx$$

(24)

Considering input numbers written as $z_1=a+jb$ and $z_2=c+jd$, a complex multiplication becomes:

$$z_1·z_2=\left(ac-bd\right)+j\left(ad+bc\right)$$

(25)

To perform the operations indicated in (25), there are two different possibilities: to perform multiplications using the library of complex functions or to perform the operations using basic mathematical functions (algebraic sums and multiplications). For the use of complex operations, it is necessary to arrange the vector to have the real value in the even places and the imaginary value in the odd places. The time to perform the CZT transformation has been measured by an oscilloscope for both cases (i.e., using complex functions and using basic mathematical operations). Less time was measured using the complex function library.

4. Experimental Results

4.1. CZT Validation Methodology

To validate the algorithm implemented on the STM32 Nucleo board, several experimental data sets were used as input for the on-board algorithm, and the CZT results were collected and compared with those of a reference virtual instrument implemented in a LabVIEW environment. As regards the experimental data, in [32], the authors described a test bench to acquire different signals, both stationary and non-stationary, to test the CZT algorithm’s performance. More precisely, the output signals from the transducers were sent to a NI (National Instrument) USB 9239 acquisition board that has the following characteristics: 24-bit simultaneous sampling, maximum sampling rate $f_{smax}=50 \mathrm{kS}/\mathrm{s}$ per channel, built-in anti-aliasing filter with automatic adjustment of the cut-off frequency as a function of sampling frequency, input voltage range ± 10 V (peak-to-peak), isolation voltage between channels 250 V. The sampling frequency was set at 10 kS/s. The acquisition board was connected to a PC, where the virtual instrument was implemented in the LabVIEW environment to perform the spectral analysis of the acquired signals. The outputs collected at the end of the algorithm were plotted, and the measurement accuracy was also calculated. CZT in the LabVIEW environment was performed using a MATLAB script. Through LabVIEW, it is possible to set the number of points of the signal portion to be analysed, the sampling frequency, and the frequency range of interest. Before elaboration, the signal has been windowed (using a Hanning window). The front panel of the tool for calculating the CZT in the LabVIEW environment is shown in Figure 9. The front panel shows the trends of the signal and its spectrum, which are used as a reference for comparisons with the algorithm implemented on the microcontroller.

The tests were carried out under different load conditions:

• A fluorescent lamp bench (80 W) with an additional resistive load (200 W);
• A vacuum cleaner;
• A dimmer;
• A hand-held electric drill.

For each test, the acquisitions were made by connecting one load. Each single load was powered by mains voltage (230 V, 50 Hz). Voltage and current were acquired, respectively, through a Hall effect differential probe (Tektronix P5200) and a shunt (resistance 0.1 Ω, accuracy 0.001%, maximum current 22 A). To obtain a non-stationary signal, an arc generator was connected in series with the above-mentioned loads.

The results produced by the virtual instrument in LabVIEW were used as reference values for comparison with those obtained with the algorithm implemented in the STM32 Nucleo. The input data for both instruments were obtained experimentally through the above-described measurement bench by acquiring the current absorbed by the different loads. For each test case, both absolute and relative errors of STM32 Nucleo CZT and harmonic measurements were evaluated with respect to the LabVIEW instrument. In detail, for the absolute errors, the differences (absolute values) between the two spectra samples were calculated and plotted; for the measurement of each harmonic component, the relative error e_h% was calculated as follows:

$$e_h\%=\frac{\left|A_{h\_STM32}-A_{h\_LabVIEW}\right|}{A_{h\_LabVIEW}}100$$

(26)

where: $A_{h\_LabVIEW}$ is the peak value of the h-harmonic lobe of the LabVIEW CZT spectrum and $A_{h\_STM32}$ is the peak value of the h-harmonic lobe of the STM32 Nucleo CZT spectrum.

The detailed results for each test case are reported in the following subsections.

4.1.1. Test Case 1: Fluorescent Lamps with Additional Resistive Load

Figure 10 shows the sampled signal of the fluorescent lamp bench with an additional resistive load, in the stationary case. This is a case with a stationary current waveform and only even harmonics in the current spectrum. Figure 11, Figure 12 and Figure 13 show the STM32 Nucleo CZT output and errors (absolute and relative, compared with the results obtained with the LabVIEW virtual instrument), respectively.

It is possible to observe that the absolute error is greater near values that are around zero. In any case, its maximum value is less than 0.01 A in absolute value (the maximum percentage error is lower than 0.35%).

Figure 14 shows the sampled signal of the fluorescent lamp bench with an additional resistive load in the nonstationary case. Figure 15, Figure 16 and Figure 17 show the STM32 Nucleo CZT output and errors (absolute and relative, compared with LabVIEW results), respectively.

The nonstationary behaviour can be visually recognized, as noise and spikes are clearly visible in the current waveform. In the nonstationary case, the current value reduced its amplitude by about 42 mA, while the odd harmonics increased their value, in particular the fifth harmonic increased by about 20 mA. As regards the STM32 Nucleo CZT results, the percentage errors on harmonic amplitudes are slightly higher for the nonstationary case than for the stationary one (maximum error of around 0.51%).

4.1.2. Test Case 2: Vacuum Cleaner

Figure 18 shows the sampled vacuum cleaner current signal in stationary conditions. This is a case where both even harmonics and inter-harmonics are present in the current spectrum (detected around fundamental and third, fifth, and seventh harmonic components). Figure 19, Figure 20 and Figure 21 show the STM32 Nucleo CZT output and errors (absolute and relative, compared with LabVIEW results).

Figure 22 shows the sampled signal of the current absorbed by the vacuum cleaner in nonstationary conditions. Figure 23, Figure 24 and Figure 25 show the STM32 Nucleo CZT output and errors (absolute and relative, compared with LabVIEW results), respectively. Although the input waveforms in both the stationary and non-stationary cases look very similar, the CZT looks very different. The third harmonic is greater than the others, and inter-harmonics are also detected in the current spectra. As for the previous case, the errors on harmonic amplitudes are slightly higher for the nonstationary case (maximum error of 0.68%) than for the stationary one (maximum error of 0.59%).

4.1.3. Test Case 3: Dimmer

Figure 26 shows the sampled signal of a dimmer current in the stationary case. This is a case where small noise is also present in stationary conditions and the current waveform is strongly partialized; in comparison with the test case 1, harmonics are higher (in percentage with respect to the fundamental) than those of the test case 1. Figure 27, Figure 28 and Figure 29 show the STM32 Nucleo CZT output and errors (absolute and relative, compared with LabVIEW results), respectively.

Figure 30 shows the sampled signal of the dimmer current signal in the nonstationary case. Figure 31, Figure 32 and Figure 33 show the STM32 Nucleo CZT output and errors (absolute and relative, compared with LabVIEW results), respectively. In comparison with the stationary case, the current waveform is even more partialized, and noise is more visible in the current waveform. The errors on harmonic amplitudes are higher for the nonstationary case (maximum error of 1.41%) than for the stationary one (maximum error of 0.25%).

4.1.4. Test Case 4: Hand-Held Electric Drill

Figure 34 shows the sampled signal of a hand-held electric drill in a stationary case. This is a case where the current waveform is strongly partialized and only the positive half-waves are present in the current waveforms. Both even and odd harmonics are present in the current spectrum, whose amplitudes are high if compared with the fundamental component (the first three harmonics amplitudes are higher than 50% of the fundamental one). Moreover, a DC component is visible in the current spectrum (higher than the fundamental component). Figure 35, Figure 36 and Figure 37 show the STM32 Nucleo CZT output and errors (absolute and relative, compared with LabVIEW results), respectively.

Figure 38 shows the sampled signal of a hand-held electric drill in a nonstationary case. Figure 39, Figure 40 and Figure 41 show the STM32 Nucleo CZT output and errors (absolute and relative, compared with LabVIEW results), respectively. Noise is visible in both the current waveform and the spectrum. As in previous cases, errors on harmonics are slightly higher in the nonstationary case than in the stationary one (with maximum values lower than 1% in both cases).

4.2. Time Analysis

Thanks to (23), it is possible to process one sample at a time. Samples are recorded in a vector, so they can be processed at maximum speed, i.e., there are no waiting times between the end of the processing of one sample and the arrival of the next. Considering the sampling frequency is equal to 10 kHz, the time interval between two samples is 100 μs. The following parameters were set for time measurements:

• Sample rate: 312.5 MS/s
• Record length: 17.1875 Mpts
• Time division: 5.5 ms/div

Figure 42 shows the measurement of the time required for performing the CZT using the MCU. To measure time, a digital I/O pin is used: when the CZT algorithm starts, the pin is set, and when the algorithm is finished, the pin is reset. It takes approximately 52.9 ms, i.e., less than the observation window of 80 ms corresponding to 4 signal periods.

4.3. Results Discussion

The results presented in the previous subsections show the feasibility of the algorithm implementation on the microcontroller board in terms of accuracy in harmonic analysis even in nonstationary conditions. In fact, in the worst case, the maximum percentage errors are lower than 1.5%. To compare such errors with some relevant standards in the field, reference can be made to IEC 61000-4-7 [34], which gives the accuracy requirements for harmonics and inter-harmonics measurements for power systems applications. This standard provides the following maximum allowable errors for current measurements: for Class I instruments, ± 5% Im for Im ≥ 3% Inom; ± 0.15% Inom for Im < 3% Inom (being Im the measured value and Inom the nominal range of the instrument); for Class II instruments: ± 5% Im for Im ≥ 10% Inom; ± 0.5% Inom for Im < 10% Inom. For the test cases herein presented, assuming as reference Inom value the amplitude of current, the following summary table can be obtained (see

Table 1. Summary of results and comparison with the maximum allowable error according to IEC 61000-4-7.

Harmonic Current Amplitudes	Test Condition (Stationary/Nonstationary)	Maximum Obtained Error in All Test Cases	IEC 61000-4-7 Maximum Allowable Error (Class I)
Im ≥ 3% Inom	Stationary	0.6% Im	±5% Im
	Nonstationary	1.1% Im
Im < 3% Inom	Stationary	0.07% Inom	±0.15% Inom
	Nonstationary	0.09% Inom

), which reports the comparison between the maximum errors (for all test cases, in stationary and nonstationary conditions, respectively) and the correspondent maximum allowable error according to IEC 61000-4-7. From Table 1, it can be seen that the standard limits (Class I) are fulfilled for all test cases.

It should be noted that the maximum allowable errors of [34] are related to stationary conditions and a DFT/FFT analysis performed on N acquired samples into an observation window of Tw = 200 ms (i.e., 10/12 cycles of the power system frequency of 50/60 Hz, respectively). As shown in Section 2, this leads to a frequency resolution of 5 Hz, independently from the choice of the sampling frequency and the consequent value of N (the higher the observation window, the higher is N, being Tw = N/fs constant). The computational burden is O(N²) or O(N logN) for the DFT or FFT algorithms, respectively. If a shorter observation window is required for the analysis of a nonstationary signal, this results in a worse frequency resolution; for example, if Tw = 80 ms (as for the test cases herein presented), the DFT/FFT frequency resolution increases to 12.5 Hz. On the other hand, with the CZT algorithm and the specified frequency band of interest (for example, 500 Hz for the test cases herein presented), the frequency resolution is improved to 0.625 Hz, according to Equation (11), with the same computational burden as the DFT algorithm.

Obviously, the computational burden will have an impact on the execution time of the spectral analysis, which will depend on both the algorithm used, and the computational resources of the board used for the implementation. For example, in [35], the authors carried out a feasibility study on the power quality metrics implementation on a different smart metering platform, i.e., the STCOMET, by STMicroelectronics [36]. The STCOMET has an integrated application core ARM^® 32-bit Cortex™-M4F CPU, with 96 MHz of maximum frequency, 640 kB or 1 MB of embedded flash, 128 kB of embedded SRAM, and 8 kB of embedded shared RAM. In comparison with the Nucleo board used in this paper, the STCOMET has more limited features in terms of both memory and computational capabilities. In [35] it was shown that this platform can be suitable for the implementation of FFT analysis and other power quality metrics. On the other hand, the implementation of the CZT algorithm, as well as the DFT one, would be not feasible on STCOMET in terms of computational time needed for the algorithm execution (estimated time for CZT of about 160 ms, with the same observation window of the test cases herein presented, which is longer than the observation window itself). As shown in Section 4.2, such an implementation is feasible on the STM32 Nucleo board, not only in terms of accuracy features but also in terms of processing time.

5. Conclusions

The increase in nonlinear loads within smart grids causes a significant increase in harmonic pollution. In particular, the increased penetration of DG, the increase in electronic loads, and the widespread presence of EV charging points contribute to amplifying any distorted waveforms. The non-stationarity of the distorted signal is one of the crucial elements in choosing the most appropriate harmonic analysis algorithm. In this paper, harmonic analysis algorithms have been studied, and the CZT algorithm has been chosen and implemented on the STM32 Nucleo microcontroller platform. The signals used for the tests were acquired in the LabVIEW environment using an NI-9239 data acquisition card with a sampling rate of 10 kS/s. The CZT of different loads were calculated: a fluorescent lamp with an additional resistive load, a vacuum cleaner, a dimmer, and an electric drill. All loads were tested in stationary and non-stationary conditions. Results were compared with the CZT algorithm performed in a LabVIEW environment. The normalized absolute percentage error was calculated. The maximum normalized error obtained is lower than 1.5%. The results allow concluding that the CZT algorithm, despite its high computational cost, achieves a very high spectral resolution, and can be implemented on a low-cost hardware platform. Besides reliability, the computational speed of the algorithm was also tested, measuring the MCU processing time. An elaboration time of 52.9 ms was obtained, which is lower than the acquisition time. This demonstrates the feasibility of the proposed solution. A possible improvement for follow-up studies can be related to the optimization of the algorithm implementation and the microcontroller resource use in order to enable the possibility of implementing the code even in lower cost devices with more limited resources and the integration of the spectral analysis algorithm with further signal processing tools for the development of a multifunction metering device to be used for both energy metering, power quality monitoring, fault detection and so on.