2.1. Mathematical Description
The developed system was tested and confirmed in Greece, where the fundamental frequency of the electric grid is 50 Hz. However, the system is capable of measuring and analyzing electric current at grids where the fundamental frequency is 60 Hz, such as in the USA and Canada. This capability can easily be achieved, as the system was parametrically designed.
It is well known that an analog periodic signal can be represented by Fourier series analysis, which is the summation of a theoretically infinite number of sinusoidal signals at different frequencies (multiples of the fundamental frequency of the signal) and with different appropriate amplitudes. These sinusoidal signals are the so-called harmonics and constitute the frequency spectrum of the analog signal. On the other hand, the measurement of an unknown analog signal can be used to determine its frequency spectrum using mathematical analysis tools such as the Fourier transform (FT).
A single sinusoidal signal can be mathematically described in the time domain, through the following Equation (1) [
14]:
where the “
Re” stands for the real part of the exponential factor, “
f” stands for the frequency of the sinusoidal signal (50 Hz for the Greek electric network), and “
A” stands for the amplitude of the sinusoidal signal.
The same signal can also be presented in the frequency domain as described below [
15]:
where “
ω” is the angular frequency and “
δ” stands for the Dirac delta function.
The Fourier transform can be applied on analog and digital signals. For the case of analog signals, the continuous time FT is used, and it is calculated using the following general formula [
15,
16]:
An example of using FT on a sinusoidal signal for an electric current at 50 Hz with an amplitude of 10 A is illustrated in
Figure 3.
The previously mentioned continuous FT can be applied on continuous (analog) signals in an infinite time slot, as Equation (3) describes. However, the implementation of measuring devices with the ability of analog signal analysis in real time is a difficult process and sometimes impossible. For this reason, the majority of measuring systems initially convert the analog signal into a digital signal. For the analysis of digital signals in the frequency domain, another tool called the discrete time Fourier transform (DTFT) is used, and the conversion from the time domain to the frequency domain can be achieved using the following Equation (4) [
15]:
where
x[
n] are the samples over the analog signal
x(
t) and “
ω” stands for the angular frequency.
Even if Equation (4) is capable of accurately calculating the frequency spectrum of the signal, a specific technical problem arises. This equation needs an infinite number of samples to calculate the frequency spectrum. To overcome the previously mentioned problem, another tool is used, called the discrete Fourier transform (DFT), which is expressed by the following equation [
15,
17]:
The difference between the DTFT and DFT is the limits on the sum factor. The DFT restricts the limits in a finite number of
N samples. In Equation (5), the angular frequency is calculated using the formula [
15]:
where the domain
k Є [0~
N − 1].
The DFT mathematical tool is a strong tool in digital signal analysis with very reliable results. However, due to its complexity and the number of additions and multiplications needed to calculate the frequency spectrum, the DFT is a slow process in calculations, even on modern computers, and also difficult to implement on hardware. For this reason, some basic algorithms have been developed, resulting in a reduction in calculation complexity. These algorithms are known as fast Fourier transform (FFT) algorithms.
The FFT algorithms are based on a strategy called “divide and conquer”. This means that an
N-point DFT can be split into two smaller DFTs of
N/2 each. The calculation of the
N-point DFT would require
N2 additions and
N2 multiplications. An
N/2-point DFT requires (
N/2)
2 additions and (
N/2)
2 multiplications, which leads to reducing the calculation complexity. The two
N/2-point DFTs can also be split into
N/4-point DFTs each, etc. [
15]. In this way, the complexity of the total calculations can be reduced, and finally, the frequency spectrum of a signal can be represented with less cost in hardware and software.
The two basic algorithm categories are [
15]:
- ➢
Decimation-in-time FFT;
- ➢
Decimation-in-frequency FFT.
The usual implementation of FFT is the selection of an even number of
N and also requiring the
N to be a power of 2 (e.g.,
n = 1024 or 2
10).
Figure 4 graphically illustrates the implementation of an 8-point radix-2 decimation-in-time FFT algorithm, while
Figure 5 illustrates the 8-point radix-2 decimation-in-frequency FFT algorithm.
The signals that exist in the electricity networks present significant periodic behavior, like many other electrical signals. These signals can be analyzed as several signals consisting of basic periodic functions, commonly sines and cosines. The standard Fourier transform is the mathematical method and the tool that allows the decomposition of the measured waveforms in these components. Every measured signal that is processed should have a finite size, to reduce the complexity and the required sources for the decomposition. The size of the sampled signal area, even though it is variable based on the signal, must be limited. The harmonic analysis is of paramount importance to detect, identify and evaluate different critical incidents, through the current and voltage waveforms. The detectability of specific signals in the presence of other nearby strong ones, the resolvability of nearby signals of similar strength and shifting ones, and biases in estimating the parameters of any of the aforementioned signals are all factors to consider [
18,
19].
In
Figure 4 and
Figure 5, the factors
WNnk can be calculated using the following Equation (7) [
16]:
Τhe selection of the signal’s observation interval is of paramount importance for the correct application of the FFT transformation. When the signal observation interval is not an integral propagation of a sinusoidal signal that will be measured, this leads to an FFT result spreading out in the frequency domain, leading to a well-known phenomenon called spectral leakage. As an example, a sinusoidal signal with an amplitude of 1 volt and a frequency of 10 kHz, sampled at a frequency of 1 MHz, is considered.
Figure 6 presents this signal with five different observation time intervals, while
Figure 7 presents the corresponding FFT results. In the first case (
Figure 6a), the observation time interval is not an integer multiple of the 10 kHz period time, but this observation interval includes four full 10 kHz periods. The corresponding FFT (
Figure 7a) gives a significant spread around the 10 kHz frequency, while this spread is not symmetric around the basic signal of 10 kHz. In a second case (
Figure 6b), the observation time interval is an integer multiple of the 10 kHz signal period and includes exactly four full periods. The corresponding FFT (
Figure 7b) is a single tone at 10 kHz. In a third case (
Figure 6c), an exactly full period of 10 kHz is used as an observation interval, resulting in an also single tone as the FFT result (
Figure 7c). In a fourth case, the observation interval is increased at exactly ten full periods, resulting in a better FFT result with less symmetric spread (
Figure 6 and
Figure 7d). Finally, a case with no integer multiple of the 10 kHz signal period including ten full periods is examined (
Figure 6e). In this case, the FFT results in a non-symmetric spread of frequencies around the frequency of 10 kHz (
Figure 7e). However,
Figure 6 and
Figure 7 refer to the recognition of a specific frequency, a situation that is not the usual case when real signals are measured. In reality, such systems aim to recognize a band of frequencies, such as the band of supraharmonics described in this manuscript. That means that the spectral leakage cannot be zeroed but only limited. This is where the so-called windows are introduced in such signal-processing applications. Windows are data-weighting functions used to limit the spectral leakage caused by finite observation intervals. Many window examples are well known [
18].
This research examined the measurement and the recognition of frequencies in the supraharmonic band without using a specific window technique but only enlarging the observation time, using a new measuring and signal analysis device based on FPGAs. This approach seems to limit the spectral leakage of the FFT result but with a negative impact on the hardware throughput. This issue is expected to be solved during the further development of the measuring system, as additional functionalities and improvements will be added, including some windowing techniques and following previous work on the topic, such as the ones proposed by Sandrolini and Mariscotti [
18].
Figure 8,
Figure 9 and
Figure 10 prove the limited spectral leakage observed during the signal simulated procedures.
Figure 8 illustrates the FFT result of a 2.5 kHz sinusoidal signal observed at the same time intervals as the one used for the 10 kHz mentioned above.
Figure 9 presents the corresponding FFT results for a 150 kHz sinusoidal signal also observed at the same time intervals. From the three mentioned cases, it is obvious that, the higher the observation time interval is (for the same sampling rate of 1 Msps), the lower the spectral leakage for the supraharmonic area becomes. In these cases, the best results correspond to a total observation time interval of 1 ms (ten full periods of the 10 kHz signal). In the proposed system, a total number of 16,384 samples are captured with a sampling rate of 1 MHz, leading to a very high observation time of 16,384 ms.
Figure 10 illustrates the result of the FFT for the worst case of 2.5 kHz (the minimum frequency in the supraharmonic band) for a 16.384 ms observation time. From
Figure 8, it seems that this number of samples is sufficiently high for the aims of the designed system.
2.3. The FPGA-Based Device
The control of the ADC and the FFT implementation were developed on an FPGA-based device called MyRio. The embedded FPGA in this device is the Xilinx Zynq-7010 [
23,
24]. This device is capable of being programmed using the LabVIEW graphical language and supports the import of VHDL codes. VHDL code is a type of hardware description language. It is used for the description of hardware structure over an FPGA or a complex programmable logic device (CPLD) [
25].
Figure 13 illustrates the general structure of the measuring device. The programming structure for the developed measuring system embeds two programming parts. The first part, which is a VDHL code, controls the external ADC, and the second part, which is a LabVIEW code, implements the FFT analysis and embeds the VDHL code. The LabVIEW code runs over a processor manufactured by National Instruments (NI), which is, by default, loaded on the FPGA, while the VDHL code runs directly over the FPGA.
Figure 14 presents the LabView code that implements the proposed measurement system. The code is constituted by three frames. In the first frame, the LabView code embeds the VHDL code. In the second frame, a loop structure with a clock of 1 MHz runs continuously. During this frame, there is an exchange of commands and data to and from the ADC through a bus connection. A FIFO (first in–first out) memory manipulates the data from the ADC. The last frame implements the FFT analysis.
Figure 15 illustrates the lab setup. In this setup, a USB cable connects the PC with the MyRio device to display the FFT analysis result on its screen. The MyRio device controls the ADC chip through an SPI-based serial interface, while a pulse generator feeds the analog input of the ADC.
Figure 16 presents the manufactured custom ADC card.
The NI processor, located in the FPGA, implemented the FFT analysis. In the frame of the presented paper, the FFT calculation used the embedded FFT function that LabView offers. However, an optimized approach for the FFT implementation can be achieved, using already-proposed approaches [
25,
26].
The ADC control was described in VDHL as already mentioned. The ADC resolution, which can be up to 10 bits, was VDHL parametric, and the developed system set up used a 10-bit resolution for test purposes. That means that 5 V of measured signal represented a value of 2
10 − 1 = 1023.
Figure 17 illustrates the simulation results for the VDHL code that controlled the external ADC.