Synchronized Measurement of the Fundamental Voltage and Harmonic, Interharmonic, and Subharmonic Components of the Electrical Grid Using an Adaptive Kalman Filter

Abstract

The effects of harmonics, interharmonics, and subharmonics on low-voltage distribution networks, leading to a deterioration in electrical power quality, have become more evident in recent years. The main harmonic sources are power electronic devices due to their implicit nonlinearity. Interharmonic and subharmonic components are mainly caused by a lack of synchronization between the grid frequency and the switching frequency of the power converters. This can be caused by asynchronous modulated devices, or more commonly by fluctuations in the fundamental grid frequency. Interharmonic currents cause interharmonic voltage distortions that affect grid-synchronized or frequency-dependent systems. The IEC-61000-4-7 proposes a general guide on harmonics, interharmonic measurements, and instrumentation in current supply systems. However, the techniques proposed in the standard are intended for measurement and do not enable a precise identification of the interharmonic components in a signal. This work proposes new definitions for the spectral energy aggrupation to improve signal component detection for the IEC standard. Furthermore, an adaptive Kalman filter algorithm is developed that enables the exact identification in real time of the frequency, amplitude, and phase of these components. The proposed system will become the basis for the implementation of a new range of measurement systems that provide improved accuracy and real-time operation. The work is supported by simulated results analysing various scenarios (including transients after changes in harmonic content in the grid voltage) that demonstrate the effectiveness of the proposed method.

1. Introduction

In recent decades, attention has been directed towards the harmful effects of harmonics, interharmonics, and subharmonics in the distribution network, and so the number of standards and regulations concerning the quality of electrical energy has grown [1]. The IEC 61000-2-1:1990 standard [2] provides definitions for the frequency components of a signal, including interharmonics. In [3,4], a precise definition of what constitutes a subharmonic component is formulated. Frequency component definitions are shown in

Table 1. Frequency components of a signal.

Component	Frequency
DC	0 Hz
Fundamental component	$f_1$
Subharmonic	$0 < f < f_1$
Harmonic	$\left.\begin{array}{c}f=h\ast f_1\\ f\ne h\ast f_1\end{array} \right\}$	h = integer > 1
Interharmonic

. These disturbances cause operational problems in grid-connected systems, especially when dealing with microgrids [5], weak low-power grids, and in residential distribution networks [6]. Problems such as errors in metering equipment, malfunctioning of relays or differential protections, and vibrations in induction motors also occur. Interharmonics can also cause synchronization problems in equipment requiring grid synchronization. Grid-connected power electronic systems are the main cause of these disturbances [7,8,9,10,11], and this results in an overall decrease in the power quality of the electrical system.

Power electronic converters are present in most domestic and industrial electrical and electronic equipment, and grid-connected photovoltaic generation systems can produce harmonics and interharmonics [8,12] under different operating conditions [10]. In the same way, it has been demonstrated that wind-generation systems are also responsible for the generation of harmonics, interharmonics, and subharmonics [11]. In the case of industrial applications, AC/DC conversion [7], motor control [13], and arc furnaces [9] stand out among the major generators of interharmonics and subharmonics.

Due to the occurrence of these phenomena, new power-quality measurement standards have been developed [1] to achieve better electrical power measurements. The IEEE Std. 1459-2010 [14] provides definitions of electrical power to quantify the flow of electrical energy in single-phase and three-phase circuits under sinusoidal, non-sinusoidal, balanced, and unbalanced conditions. IEC 61000-4-7 [15] and IEC 61000-4-30 [16] define specifications to be met by electrical measurement devices and instrumentation intended for measuring harmonics and interharmonics. In turn, new measurement systems employing different estimation methods and techniques have been developed [17]. Some of these systems are intended for renewable energy-generation systems, as in [12,18,19], while others are oriented to more general applications, as in [17,20,21]. Measurement systems based on the discrete Fourier transform (DFT) [8,18], which is the method recommended in [15], have the disadvantage of being very sensitive to picket fence and spectral leakage effects [22,23], and these are described in [24]. In addition, the error between the actual and the computed frequency value of a harmonic or interharmonic spectral component depends on the acquisition window. The IEC standard [15] does not identify the actual interharmonics existing between two consecutive harmonic frequencies. However, by grouping the spectral components between them, it provides a single rms value quantifying the interharmonic spectral information and indicating the presence of interharmonics in this range of frequencies.

To improve the identification of interharmonics, new grouping techniques are proposed in [3] to achieve accurate measurements of interharmonic component values. These new techniques use longer times for the acquisition window [3], which makes them more sensitive to picket fence and leakage effects while also increasing the computational cost. However, these methods do not solve the main problem of providing accurate values for the frequency and rms values of the interharmonic components. Techniques based on the wavelet transform [25], hybrid models, including interharmonic estimation techniques and filter banks [26], or, more recently, those that incorporate artificial intelligence (AI) require considerable processing time [17,27]. AI-based approaches will become an alternative to traditional harmonic analysis in the near future, although some challenges must be faced first. To date, their effectiveness is low when dealing with unknown signals containing unknown harmonics and interharmonics [27]. Another technique that has become widespread is the linear Kalman filter (LKF). This is a well-known technique based on a recursive algorithm that can be used to accurately obtain the amplitude, frequency, and phase angle of a signal [28,29]. The main concern with the application of LKF in measuring non-sinusoidal signals (voltages or currents) is the need for prior knowledge of the real number of harmonic and interharmonic components present in the voltage. Without the correct dimension of its matrices, the algorithm may fail to converge or may produce inaccurate results [30].

Intended for the future development of new power-quality measurement equipment, this work proposes a new real-time adaptive algorithm that improves the frequency component identification and quantification of an unknown electrical signal. New definitions are proposed for spectral energy grouping, allowing for the identification of the real number of harmonic and interharmonic components. Therefore, an accurate model of the signal is obtained, enabling the implementation of an LKF with adequate dimensions to perform a precise tracking of the signal components. Applied to network voltage, it results in a highly accurate identification method for harmonic and interharmonic voltage components. The proposed algorithm also enables for fast convergence when the harmonic spectrum of the network voltage changes.

This article is structured as follows: Section 2 presents a review of the IEC 61000-4-7 standard. Section 3 details the operation of the LKF in obtaining real-time values of the individual components present. Section 4 delineates the developed method. In Section 5, the presented method is validated through simulations. Finally, Section 6 presents the conclusions.

2. IEC 61000-4-7 Recommendations for the Measurement of Harmonics and Interharmonics

The IEC 61000-4-7 standard [15] specifies a series of groupings of spectral components obtained after applying a DFT to the grid voltage signal in order to identify and quantify the harmonics and interharmonics present in the electrical power network. The standard specifies that the DFT must be applied over the signal samples acquired during a synchronized time window ($T_w$) of 10 or 12 cycles of the fundamental frequency for 50 Hz or 60 Hz networks, respectively. The acquired values are weighted using a rectangular weighting. However, the Hanning weighting is allowed if there is synchronization loss. After this process, the DFT is applied, and the spectral information obtained is grouped according to the definitions presented below (for 50 Hz networks). The distribution of the spectral components is shown in Figure 1, where k is the spectral component order, being k = 1 the dc component of the signal. The components are spaced every 5 Hz. In this way, the fundamental component n = 1 and its harmonics are associated with the spectral components of order $k=(10·n)+1$.

2.1. Harmonic Group

The harmonic group “n” represents the sum of the spectral components around the harmonic “n”, as presented in Equation (1) for 50 Hz networks (only the expressions for 50 Hz networks are presented). It attempts to represent the value of the harmonic “n” following the concept of energy conservation according to Parseval’s theorem. The rms value of the harmonic group of order “n” ($G_{g,n}$) is calculated in (1) and shown in Figure 2 for the harmonic group n + 1.

$$G_{g,n}^2={\displaystyle \frac{C_{k-5}^2}{2}}+{\sum }_{i=-4}^4{C}_{k+i}^2+{\displaystyle \frac{C_{k+5}^2}{2}}$$

(1)

where g denotes a harmonic group, n represents the harmonic order, and $C_{k+i}$, is the rms value of the spectral component k + i, with $k=(10·n)+1$, as indicated previously.

2.2. Interharmonic Group

An interharmonic group represents a regrouping of the spectral components in the interval between two consecutive harmonic components, as presented in Equation (2). The rms value of the interharmonic group of order n ($C_{ig,n}$) is calculated in (2) and is shown in Figure 2 for the interharmonic group n = 3 (where k + i will be in the range of 32–40).

$$C_{ig,n}^2=\sum _{i=1}^9{C}_{k+i}^2$$

(2)

where ig denotes an interharmonic group.

The result obtained in (2) quantifies the combined rms value of the possible interharmonics present between two consecutive harmonics. However, with this aggrupation, it is not possible to establish the number of interharmonics present in this range of spectral components. Consequently, neither their frequencies nor their rms values can be known.

2.3. Harmonic Subgroup

The harmonic subgroup is the grouping of a harmonic with its adjacent spectral components. It is intended to regroup the energy of a harmonic that may have spread among the nearest components due to the effects of spectral leakage [18,19] and the picket-fence effect [20]. The rms value of the harmonic subgroup n ($G_{sg,n}$) is calculated in Equation (3) and shown in Figure 3 for the harmonic subgroup n = 2.

$$G_{sg,n}^2=\sum _{i=-1}^1{C}_{k+i}^2$$

(3)

where sg denotes a harmonic subgroup.

2.4. Interharmonic Subgroup

Excluding the spectral components adjacent to the harmonics n, n + 1, etc., the sum of the remaining spectral components found between two consecutive harmonics form an interharmonic subgroup. The rms value of the interharmonic subgroup n ($G_{sg,n}$) is calculated in Equation (4) and shown in Figure 3 for the interharmonic subgroup n = 3.

$$C_{isg,n}^2=\sum _{i=2}^8{C}_{k+i}^2$$

(4)

where isg denotes an interharmonic subgroup.

The IEC standard provides values for harmonics and interharmonics, whether as harmonic groups (1) or interharmonic groups (2), as an association of several spectral components in order to maintain the principle of energy conservation in Parseval’s theorem. The same applies to harmonic subgroups (3) and interharmonic subgroups (4). These values are produced in a manner that does not enable the distinguishment of whether one or more spectral components are present within the same interharmonic space or the exact frequency value, as the provided value is a single value at the central frequency position. Therefore, in applications requiring a more precise measurement of the components forming the grid signal, it is necessary to employ techniques that offer more accurate solutions. In this way, new definitions to improve the spectral energy grouping will be presented in Section 4. The use of the proposed expressions will result in a precise identification of the real components of the voltage signal, which is crucial to improving the power-quality measurement equipment accuracy.

3. Linear Kalman Filter and Derivates

The linear Kalman filter (LKF) [31] is a recursive algorithm that performs real-time estimation of the state variables of a system. It is based on minimizing the estimation error between the state real value $x_k$ and the estimated value ${\widehat{x}}_k$ using the least-squares estimation method. The algorithm operates recursively in two stages: estimation and correction, as shown in Figure 3. The mathematical model description of a sinusoidal and non-sinusoidal signal is presented in Appendix A.

Among many other applications, the LKF is used to estimate the amplitude and phase of the fundamental frequency component of the voltage in the electrical power network [28]. The LKF assumes a stochastic model for the estimation of the state variables. This model should be defined before starting the recursive estimation process. The equations that define the system model to be estimated are defined as follows:

$$x_{k+1}={\Phi }_k{x}_k+w_k$$

(5)

$$z_k=H_k{x}_k+{\upsilon }_k$$

(6)

where

$x_k$ corresponds to the value of the state vector at the current instant (k).

${\Phi }_k$ corresponds to the value of the state transition matrix at the current instant (k).

$x_{k+1}$ corresponds to the value of the state vector at the next sampling instant (k + 1).

$w_k$ corresponds to the value of the system’s intrinsic error at the current instant (k).

$z_k$ corresponds to the value of the system’s output at the current instant (k).

$H_k$ corresponds to the value of the system’s output matrix at the current instant (k).

$v_k$ corresponds to the value of the error introduced and measured by the measurement system at the current instant (k).

$w_k \mathrm{a}\mathrm{n}\mathrm{d} {\upsilon }_k$ are uncorrelated white noises, where $w_k$ is the intrinsic noise of the system, and ${\upsilon }_k$ is the noise introduced by the measurement system. Their known covariances are:

$$E\left\{{w_k{w}_i}^T\right\}=\left\{\begin{array}{c}{Q}_k;i=k\\ 0,i\ne k\end{array}\right.$$

(7)

$$E\left\{{v_k{v}_i}^T\right\}=\left\{\begin{array}{c}{R}_k;i=k\\ 0,i\ne k\end{array}\right.$$

(8)

$$E\left\{{w_i{v}_j}^T\right\}=0,{\forall }_{k,i}$$

(9)

where,

$E\left\{.\right\}$ represents the mathematical expectation.

$Q_k$ is the covariance of the intrinsic system noise.

$R_k$ is the covariance of the noise measured in the output.

It is assumed that matrices ${\Phi }_k$ and $H_k$ are known. The state error in the estimation stage is defined as:

$$e_k^{-}=x_k-{\widehat{x}}_k^{-}$$

(10)

where $x_k$ is the actual state value at time k, and ${\widehat{x}}_k^{-}$ is the estimation of the stage at time instant k, which is obtained by model prediction from previous information at instant k − 1. The superscript (−) indicates that the calculation for this time instant (k) was performed using the values obtained or estimated at the previous instant (k − 1), and the caret symbol (^) indicates that the displayed value is an estimation rather than a real value.

Equations (11) to (15) constitute the recursive LKF model. Equation (11) shows the correction phase, where the value of the state is estimated based on the prediction state estimation ${\widehat{x}}_k^{-}$ and the measured ($z_k$), which leads to the error between the output measured value $\left(z_k\right)$ and its estimated value ($H_k{\widehat{x}}_k^{-}$). The Kalman filter gain ($K_k$) is the weighting factor of the error:

$${\widehat{x}}_k={\widehat{x}}_k^{-}+K_k\left(z_k-H_k{\widehat{x}}_k^{-}\right)$$

(11)

The optimal error weighting factor ($K_k$), also known as the Kalman gain matrix, is calculated using Equation (12).

$$K_k=P_k^{-}{H}_k^T{\left(H_k{P}_k^{-}{H}_k^T+R_k\right)}^{-1}$$

(12)

where $P_k^{-}$ is the covariance of the prediction stage estimation error defined as:

$$P_k^{-}=E\left\{e_k^{-}·{\left(e_k^{-}\right)}^T\right\}$$

(13)

The update of the error covariance matrix ($P_k$), that is the variance of the correction stage, is calculated as in (14).

$$P_k=\left(I-K_k{H}_k\right)P_k^{-}$$

(14)

where I is the identity matrix.

The value of the state at the future instant, using the current values (${\widehat{x}}_{k+1}^{-}$), is obtained in Equation (15) as follows:

$${\widehat{x}}_{k+1}^{-}={\Phi }_k{\widehat{x}}_k$$

(15)

Finally, the estimation of the error value that will be obtained ($e_{k+1}^{-}$) is calculated in (16). This value is used to estimate the covariance of the future error ($P_{k+1}^{-}$) in (17).

$$e_{k+1}^{-}=x_{k+1}-{\widehat{x}}_{k+1}^{-}=\left({\Phi }_k{x}_k+w_k\right)-\left({\Phi }_k{\widehat{x}}_k^{-}\right)={\Phi }_k{e}_k+w_k$$

(16)

$$P_{k+1}^{-}=E⌈e_{k+1}^{-}\ast {e_{k+1}^{-}}^T⌉=E\left[\left({\Phi }_k{e}_k+w_k\right)\ast {\left({\Phi }_k{e}_k+w_k\right)}^T\right]={\Phi }_k{P}_k{\Phi }_k^T+Q_k$$

(17)

Once the equations and matrices are defined, the recursive algorithm of the LKF can be initiated. Figure 4 shows the flow chart of the recursive algorithm.

${\Phi }_k$ is the system’s state matrix, and it must match the system model, meaning that to adequately size this matrix, the number of harmonic components must be known, corresponding to Equation (5). The accuracy of the estimation then depends on the LKF model [30]. The LKF will converge if the model contains a number of states equal to the harmonic components of the signal to be estimated.

The LKF will attempt to adjust the developed model to the real value of the signal, aiming to adjust the output of each component so that the result is as close as possible to the actual value. This approach may result in the overall output value of the LKF being exactly equal to the desired signal, but the obtained components may differ from the actual components [30].

4. Description of the Developed Method

The developed method obtains the number of harmonic components present in the signal under study through an algorithm that modifies Equations (3) and (4) from [15]. The resulting outcome is used to set the structure of the ${\Phi }_k$ matrix of the Kalman filter (5), as prior knowledge of the number of harmonic components in the signal to be analysed is required for the implementation of the LKF. Without this prior knowledge, the identification will contain errors [30] and the results obtained will be invalid for an accurate measurement of the network signal. The proposed method is detailed below.

4.1. Data Acquisition and DFT

The data acquisition follows the IEC 61000-4-7 standard. The sampling period is 20,480 Hz and the samples are stored in a buffer of 4096 data points, corresponding to an acquisition window of 200 ms (T1), to achieve a DFT base frequency of 5 Hz when the fundamental harmonic is 50 Hz. As shown in Figure 4, data acquisition occurs at the start of the algorithm, and whenever the system fails to meet the error condition, the signal obtained with the LKF does not correspond to the input signal.

The acquired signals are stored in an input data buffer. A Hanning weighting window [23] is applied to the buffer values, and the result is used to compute the DFT. The DFT provides a combination of spectral components consisting of real values, at frequencies close to the spectral components actually forming the signal to be measured, and a set of component values that do not actually exist, originating from leakage and picket-fence effects.

4.2. Harmonics Selection

The method shown below aims to identify, among all the values obtained after applying a DFT to the input signal, the values that are real and those that result from leakage and picket-fence effects. The algorithm for performing the identification is based on [15], with some equations from [15] being modified. The modifications made are detailed below. The rms value of the harmonic subgroup for the fundamental harmonic, at 50 Hz, is calculated in (18).

$$G_{sg,1}^2=\sum _{i=-2}^2{C}_{11+i}^2$$

(18)

where $C_{11}$ is the effective value of the spectral component corresponding to 50 Hz. It is the same value as $K_{11}$ in Figure 1. They are equivalent, maintaining the nomenclature $“C_x”$ to preserve similarity with the standard. Compared to [15], two additional spectral components are added to the calculation of the effective value of the harmonic subgroup, as shown in Figure 5.

By adding these two bands, part of the energy dispersed by the spectral leakage and the picket-fence effect are included in the 50 Hz harmonic subgroup, preventing them from affecting the calculation of adjacent harmonic subgroups. Hence, they are included in Equation (18). The rms value of the harmonic subgroups other than the fundamental is computed following (19):

$$G_{sg,n}^2=\sum _{i=-1}^1{C}_{k+i}^2$$

(19)

where n indicates the order number of the harmonic for which the harmonic subgroup is being calculated, n > 1.

The rms value of the first interharmonic subgroup, corresponding to the group between 5 and 50 Hz (subharmonics), using the developed Equation (20):

$$C_{isg,1}^2=\sum _{i=4}^8{C}_{k+i}^2$$

(20)

Rms value of the second interharmonic subgroup (21):

$$C_{isg,2}^2=\sum _{i=3}^8{C}_{k+i}^2$$

(21)

Rms value of the remaining interharmonic subgroups is calculated according to Equation (22).

$$C_{ig,n}^2=\sum _{i=2}^8{C}_{k+i}^2$$

(22)

The minimum value of the first harmonic subgroup (MVFHS) is defined as:

$$MVFHS=K_{LS}\ast G_{sg,1}=0.003\sqrt{\left(\sum _{i=-2}^2{C}_{11+i}^2\right)}$$

(23)

The constant $K_{LS}$ relates the sensitivity and error of the data-acquisition system to the results obtained in the DFT. After observing the datasheets of numerous sensors, the authors have determined that the value of $K_{LS}$ = 0.003 effectively enables distinguishing the results that a sensing system could not acquire. In the next stage, all values in the DFT output vector that are less than the MVFHS value are removed, as shown in Figure 6a.

For the next stage of the method, the minimum value of the interharmonic subgroups (MVIS) is defined as:

$${MVIS}_n=K_{RS}\ast C_{isg,n}=0.1\ast C_{isg,n}$$

(24)

The value of $K_{RS}$ corresponds to a constant value, similar to $K_{LS}$, as a design parameter that determines the minimum value a spectral component must have to be employed in the calculation of the interharmonic subgroup. All spectral component values lower than the MVIS of that interharmonic subgroup will be removed from the DFT output as they are considered products of peak fence and/or leakage effects. Empirically, a value of $K_{RS}=0.1$ has been found to provide optimal results in most cases.

Figure 6 represents how the method works with the DFT output to define the spectral components considered correct. Figure 6a shows that the application of the algorithm makes the first selection among the spectral components based on MVFHS. Figure 6b shows the second selection of the spectral components based on MVIS. Note that the equations related to every subprocess have been included in its corresponding block.

4.3. Obtaining Harmonic, Interharmonic, and Subharmonic Components in Real Time

After discarding spectral components that do not exceed the MVIS values, a grouping of spectral bars is performed to estimate which constitutes a true spectral component.

As previously stated, the aim is not to find the exact values of the spectral components present in the signal under study but rather to identify the number of spectral components (including harmonics, subharmonics, and interharmonics), along with their approximate effective values and frequencies. These approximate values serve as “seed” or starting values for developing the KF-PLL structure.

The grouping criterion follows the approach outlined in [15,32,33], which groups adjacent spectral components to obtain a single result. The algorithm employs the following criteria for grouping:

(a) Results calculated with non-zero value harmonic subgroups are considered as harmonics.

(b) When a spectral component with a value exceeding MVIS is detected in an interharmonic subgroup, it checks for adjacent components with values exceeding MVIS. If there are no adjacent components, the first component is considered a single interharmonic, with its value and frequency identical to that of the component.

(c) If there are two or three adjacent components, they are considered a single interharmonic. In the case of two components, the output value in magnitude and frequency is that of the first component. For three consecutive components, the values of the central component are considered valid.

(d) In the case of four-to-six consecutive components, it is considered that there are two consecutive interharmonic components. These two components will have the following values:

(d1) In the case of four, the central value of the first three is considered valid, and the second value is that of the fourth component.

(d2) For five components, the first is calculated in the same way as four components, and the second value is adopted from the fourth component.

(d3) With six components, the values are obtained as if they were two sets of three components, taking the central values of the two groups of three. Figure 6 depicts the result of applying the proposed method to the output of a hypothetical discrete Fourier transform of an unknown signal.

Figure 7 illustrates the spectral component grouping process to determine which components are real. In the interharmonic subgroup (m), the minimum value of the interharmonic subgroup m (${MVIS}_m)$ is shown as a dashed red line. The algorithm groups the spectral components with values exceeding ${MVIS}_m$ that are adjacent, as seen in the case of those situated at $f_{m_{k+2}}$ and $f_{m_{k+3}}$. From this grouping, the algorithm determines that the spectral component at $f_{m_{k+3}}$ is a real component. The spectral component located at $f_{m_{k+7}}$, being solitary, is considered a real spectral component since the algorithm identifies any solitary spectral components as real, with the MVIS of each interharmonic group defining the limit. This process continues until all output spectral bars of the DFT are analysed.

4.4. Linear Kalman Filter

With the results obtained in the previous section, the number of harmonic components composing the input signal is determined. The number of components is used to determine the number of variables necessary to implement the equations of the LKF (5–17). The obtained voltage and frequency values of each component are used to provide initial values for the LKF variables. Once the LKF is implemented, during the first 200 ms of its execution, the values of each of the spectral components are considered invalid, as it is assumed that the LKF is still in its transient stage. After the initial instants, the error in the estimation of the input signal is compared with the maximum allowable error (Figure 8).

4.5. Evaluation of the Estimation and Restart Condition of the Calculation

For evaluating the estimation performed with the LKF and determining the need to reevaluate the number of variables used in the LKF, two values are employed.

The maximum allowable error (MAE) represents the maximum value that can occur in the LKF estimation and is 10% of the value of the main spectral component obtained in the DFT.

$$\mathrm{M}\mathrm{A}\mathrm{E}=0.1\ast C_{11}$$

(25)

where $C_{11}$ is the value of the 50 Hz spectral component obtained in the DFT.

The root mean square error (RMSE) represents the average value of the error committed by the estimation during the previous 10 ms. It is calculated using the equation:

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{{\displaystyle \frac{\sum _{t=1}^T{\left({\widehat{y}}_t-y_t\right)}^2}{T}}}$$

(26)

where,

T = number of samples used in the calculation.

${\widehat{y}}_t$ is the value of the estimation at time t made with the LKF.

$y_t$ is the value of the input signal at time t.

The proposed method compares the output of the LKF with the maximum allowable error. If the result is greater, the RMSE is calculated, and the maximum allowable error is compared again. If the result is negative, the method determines that there has been a disturbance and the system has returned to its steady state, so it outputs the values of the estimations made by the LKF. If the result is positive, it determines that there has been a sufficiently high change in the input values to justify performing a new DFT and running the entire algorithm again to implement the LKF with a new number of variables. Figure 8 depicts the stages and phases of the method. The equations related to every process have been included in its corresponding block.

5. Results

Two simulations have been conducted in Matlab 2021b–Simulink to demonstrate the validity of the method. A third simulation, including noise, has been conducted in order to analyse its effect on the proposed system. The simulations begin by forming an input signal with a defined number of spectral components. This signal is then subjected to the developed method, which involves the stages outlined in Figure 8. The first simulation involves applying the method to a signal with 11 spectral components, which, at a certain time, changes to 5 spectral components. The second simulation starts with a signal of 5 spectral components that is modified to 11 components. Despite these similarities, the frequencies of the spectral components used in the two simulations differ so that various aspects of the method can be tested. In the third case, the second simulation is repeated, adding a white noise to the input signal in order to evaluate the behavior of the developed system. The flowchart of the step-by-step research methodology is presented in Figure 9.

5.1. Simulation 1: Signal with 11 Harmonic Components Reduced to 5 Components

The simulation begins with the generation of the signal, which is analysed by the proposed method. The input signal is composed of 11 spectral components, as shown in

Table 2. Values of the spectral components contained in the input signal.

No. Harmonic	Voltage rms (V)	Frequency (Hz)	Frequency (Rad/s)
1	23	22	138.226
2	230	49.9	307.870
3	46	71.9	451.761
4	34.5	122	766.548
5	69	149.7	940.592
6	23	171.7	1078.822
7	23	249.5	1567.654
8	23	271.8	1707.769
9	23	449.1	2821.778
10	23	471.1	2960.008
11	23	548.9	3448.840

.

In Figure 10, the time evolution of the input signal, consisting of the 11 spectral components, is shown over 10 fundamental cycles (200 ms) required for the calculation of the DFT [15].

In Figure 11a, the result of the DFT performed on the input signal is shown. In Figure 11b, the DFT result is displayed in greater detail to better illustrate the effects of leakage and the picket-fence effect in the DFT output. In Figure 11b, it is clearly observed how a multitude of spectral components appear in the DFT output that do not actually exist. The main objective of the method shown is to differentiate between the real components of the DFT output and those that are not, as shown in Figure 11b.

The signal shown in Figure 10 and Table 2 is applied to the measurement system. The system computes the discrete Fourier transform (DFT) of the input signal and applies the algorithm described in the preceding sections. The obtained results are presented in

Table 3. Values obtained from DFT using the proposed method.

No. Harmonic	Voltage (Peak) Measured	Frequency (Hz) Measured	Frequency (Hz)
1	14.63	20	22
2	199.2	50	49.9
3	29.65	70	71.9
4	21.97	120	122
5	59.75	150	149.7
6	15.13	170	171.7
7	19.9	250	249.5
8	14.97	270	271.8
9	19.9	450	449.1
10	15.79	470	471.1
11	13.93	550	548.9

. It is observed how the first part of the method selects 11 spectral components from the output values shown in Figure 11a. Additionally, Table 3 provides a comparison between the frequency values obtained in the selection and the actual values of the components. In the comparison, it is observed that the obtained values approximate the real values.

The values obtained in Table 3 are used to construct the state transition matrix (5) of the LKF. Once the state matrix is constructed, the system initiates the LKF and the results are depicted in Figure 12 and Figure 13.

Figure 12 shows the signal estimated by the LKF compared to the input signal. LKF starts at t = 0.3 s after the DFT and the proposed harmonic selection algorithm are completed. The evolution of the error signal is shown in Figure 13. As the LKF structure is established with the appropriate dimensions, the algorithm converges and the estimated signal matches the real signal after 300 ms, showing minimum error.

This initial phase of Simulation 1 shows how the developed method calculates the number of spectral components forming the input signal. Additionally, with this information, it appropriately implements the equations constituting the Kalman filter (5–17), causing the Kalman filter output to converge to the true value (Figure 12), thereby minimizing the error between the actual and estimated signals (Figure 13).

Reduction of the Number of Spectral Components in the Input Signal

To validate the operation of the proposed method, the input signal is modified by varying its composition at simulation time t = 0.8 s. At that moment, it ceases to be composed of 11 spectral components and instead becomes composed of 5 components. The waveform of the input signal is shown in Figure 14, and

Table 4. Spectral components contained in the input signal from t = 0.8 s.

No. Spectral Component	Voltage rms (V)	Frequency (Hz)	Frequency (Rad/s)
1	23	22	138.23
2	230	50.1	307.876
3	46	71.9	451.761
4	34.5	122	766.548
5	69	149.7	940.592

displays the value of the components comprising it.

In Figure 15, it is shown how the LKF estimation was accurate in the moments prior to the change in the value of the input signal (t = 0.8 s). With the change, the LKF attempts to estimate the new signal with the number of variables from the previous signal, resulting in an error between the real signal and the estimation (Figure 16).

Figure 16 shows the error between the real signal and the estimation made by the LKF. It is observed that in moments prior to the change in the value of the input signal, this variation is practically non-existent, while once the change occurs, the error increases. This variation in the error signal value causes the method to restart, reacquiring the data for the calculation of the DFT (Figure 8). Figure 17, Figure 18 and Figure 19 depict the evolution of some of the estimations made for the spectral components in the period between t = 0.75 and t = 1 s. During this period, the estimations are made with the number of states prior to the change in the input signal.

Figure 17, Figure 18 and Figure 19 show the evolution of the estimation (red) of three out of the five spectral components composing the input signal. It is shown how the estimation is erroneous, and the obtained result is invalid for the application. These errors occur because the LKF attempts to minimize the error (Figure 15 and Figure 16) between the input signal and the total estimation of the model, thereby modifying the values of all the internal variables it possesses. In this case, it is observed how spectral components that do not exist in the real signal (Figure 19, blue) are assigned an estimated value to minimize the error between the estimation value and the real signal (Figure 19, red). These results confirm that when there is a change in the number of spectral components composing a signal, a new set of equations must be implemented to account for the actual number of spectral components [34]. Otherwise, as demonstrated, it may result in the estimation of the overall signal, which, when compared with the input signal, has minimal error, but the estimations of the variables are far from their real values.

The proposed algorithm is applied again at the end of the DFT at t = 1 s, and the output values are reassessed by the method to determine which values are real and which are the result of the undesired effects of the DFT (Figure 5, Section 4). In the next stage, the equations of the LKF (5–17) are rewritten, adapting the Kalman filter to the actual components of the input signal so that they only contain equations corresponding to the true spectral components. The simulation continues, and the new LKF performs the estimation of the different spectral components.

In Figure 20, the evolution of the estimation made with the LKF (red) and the actual input signal (blue) is shown. It can be observed how the difference between them decreases as the simulation progresses until it is minimized. In Figure 21, the evolution of the error is shown.

In Figure 21, it can be observed how the error occurs between the estimation made with the LKF using the new equations implemented after estimating the number of spectral components and the actual input signal. It shows how the error minimizes over the course of the simulation and reaches a minimum (less than 0.1% of the maximum value of the main harmonic) after 260 ms from the beginning of the new estimation. A detail of the estimation error is shown in Figure 22.

Figure 23 illustrates the evolution of the estimation of the 22 Hz spectral component. It shows from t = 0.95 s, moments prior to the start of the estimation with the new variables of the LKF, but still with the variables for the previous 11 spectral component signals. It demonstrates how the estimation, after implementing new equations of the LKF according to the results provided by the method (Table 4), tends to converge to the value of the real signal. Figure 24 shows the convergence of the estimation of the frequency value for all the harmonic components composing the input signal. The convergence of the frequency measures the time it takes to occur. For this harmonic component (22 Hz), the LKF requires about 400 ms.

In this second phase of Simulation 1, it is demonstrated how the Kalman filter can produce erroneous results when the variables to be estimated do not correspond to the actual variables of the input signal (Figure 16, Figure 17 and Figure 18). It has been shown how the developed method addresses this issue by recalculating the DFT of the input signal and, through this method, determining which DFT values are real and which are not. With this information, the method appropriately implements the equations constituting the Kalman filter (5–17), causing the filter’s output to converge to the true value (Figure 19), thereby minimizing the error between the actual and estimated signals (Figure 21 and Figure 22). Figure 23 illustrates how the method enables the estimation of the value of each harmonic component present in the input signal. In this case, the 22 Hz component is shown. The same applies to the other components, which are not displayed to avoid redundancy. Figure 24 shows that another piece of information provided by the method is the temporal evolution of the frequency value of each spectral component forming the input signal.

In summary, these results have demonstrated how the method resolves the problem that occurs when using the LKF in measuring spectral components that vary over time (Figure 17, Figure 18 and Figure 19). The simulation has also shown the necessity of determining, among the output values of the DFT (Figure 10), which values correspond to components that really exist and shows how the presented method is suitable for identification.

5.2. Simulation 2: Input Signal with 5 Spectral Components That Increases to 11 Components at t = 0.8 s

The following simulation can be divided into two parts: the initial phase, during which a signal with five spectral components (values in

Table 5. Spectral components comprising the input signal.

No. Spectral Component	Voltage rms (V)	Frequency (Hz)	Frequency (Rad/s)
1	23	33.9	213.032
2	230	51.1	321.070
3	46	75.15	472.181
4	34.5	121.75	764.977
5	69	150.45	945.305

) is introduced into the system, and various figures obtained from the simulation are shown to demonstrate the efficacy of the method in detecting these components. The second phase starts at 800 ms of the simulation when the input signal varies and is formed by 11 components. Among these components are two interharmonics that share the same interharmonic subgroup (in the frequency range between 205 and 245 Hz). These components are introduced to assess the validity of the method in identifying spectral components.

The value of the spectral components contained in the input signal is shown in Table 5. Its waveform is presented in Figure 25.

The method begins at t = 0.1 s when the DFT is initiated. Upon completion of the DFT at t = 0.3 s, the harmonic filtering process presented in Figure 6a,b and the selection of spectral components described in Section 4.3 is applied. The obtained result is shown in

Table 6. Spectral components obtained after applying the proposed method.

No. Spectral Component	Voltage rms (V)	Frequency (Hz)
1	11.24	30
2	199	50
3	15.79	70
4	22.9	120
5	59.75	150

.

With the results obtained from Table 6, the matrices and equations of the LKF (Equations (5)–(17)) are implemented, initiating the prediction of the different spectral components contained in the input signal. Figure 26 illustrates the evolution of the input signal (blue) and the estimation performed with the LKF (red).

Figure 27 displays the temporal evolution of the error between the input signal and the estimation made. The error evolution reveals the time taken by the LKF to reach a steady state, thus providing estimated spectral component values close to the actual ones. In this case, the algorithm requires 300 ms to reach a steady state at t = 0.6 s. Figure 28 is an enlargement of Figure 27, providing a more detailed view of the time period during which steady state is achieved (t = 0.55 s to t = 0.7 s) and the error becomes minimal.

Figure 29 shows the evolution of the subharmonic (red waveform) present in the input signal. The estimation begins at t = 200 ms, reaching a steady state at t = 500 ms. This result is similar across all estimations of the spectral components that form the input signal. Thus, they have been omitted to avoid redundancy. The evolution of the subharmonic is analogous to the overall estimation of the input signal in Figure 26.

As in the first simulation, the Kalman filter provides the temporal evolution of the frequency estimation of the input signal, as shown in Figure 30. These values can be useful in multiple applications for measuring power quality.

Increase to 11 Components at t = 0.8 s

In the second part of the simulation, the evolution of the system is shown when the input signal varies and transitions from being composed of 5 components to 11 spectral components.

Table 7. Harmonic components present in the input signal at the instant t = 800 ms.

No. Spectral Component	Voltage rms (V)	Frequency (Hz)	Frequency (Rad/s)
1	23	33.9	213
2	230	51.1	321.0707
3	46	75.15	472.1813
4	34.5	121.75	764.9778
5	69	150.45	945.3052
6	23	214.3	1346.4866
7	23	229	1438.8494
8	23	250.5	1573.9379
9	23	450.7	2831.8316
10	23	497	3122.7431
11	16.1	549	3449.4687

displays the new components forming the input signal, and Figure 31 illustrates the temporal evolution of the new signal.

The input signal changes at t = 0.8 s. As shown in Figure 32, there is a rapid increase in the error produced between the real signal and the signal estimated by the LKF. Figure 33 illustrates the evolution of the input and estimated signals.

In response to the increase in estimation error, the method starts again, performing the DFT and applying the algorithms described in Section 4. This leads to a new identification of the number of components in the real signal, and the obtained result is presented in

Table 8. Harmonic components identified in the input signal after the increase of components.

No. Spectral Component	Voltage rms (V)	Frequency (Hz)
1	11.26	30
2	199.4	50
3	15.14	70
4	22.55	120
5	59.74	150
6	9.851	210
7	15.87	230
8	19.91	250
9	19.91	450
10	19.57	500
11	13.93	550

.

Figure 34 depicts the evolution of the input signal to the algorithm and the estimation made by the LKF (in red). The estimation of the input signal begins at t = 1 s, which corresponds to the completion of the DFT calculation, implementing the LKF with the number obtained from the spectral components forming the signal.

In Figure 35, the evolution of the error signal between the input signal to the algorithm and the estimation made by the LKF (Figure 34, red) is shown. The estimation of the input signal begins at t = 1 s, with the error being higher at the start of the estimation. It is observed that the error signal decreases, reaching a steady state at t = 1.3 s of simulation (300 ms after the estimation begins). Figure 36 provides a closer look at this approximation. It is noteworthy that at t = 1.25 s, the error is less than 0.5% of the value of the main harmonic and less than 1% of that provided by most commercial sensors.

Figure 37, Figure 38, Figure 39, Figure 40 and Figure 41 present the results obtained for the spectral components of 214.3 Hz, 229 Hz, and 250.3 Hz. These results show how the proposed algorithm allows for the distinguishment between two interharmonics within the same interharmonic group (214.3 Hz and 229 Hz), demonstrating one of the advantages of the proposed method compared to the formulation suggested by the standard [15]. The results for the remaining spectral components are not shown to avoid redundancy.

Figure 41 illustrates the evolution of the spectral components between Harmonics 4 and 5, specifically within the range of 200 to 250 Hz. The method effectively distinguishes the presence of the two interharmonic components (Table 8) and constructs an LKF encompassing all of the harmonics present. Consequently, the LKF provides accurate estimations of all the spectral components, both in terms of rms value and frequency, enabling a precise determination of the instantaneous value of each component, as demonstrated in Figure 37, Figure 38, Figure 39, Figure 40 and Figure 41. This identification is particularly significant when compared to the IEC standard [15] as it accurately identifies the two interharmonic spectral components within the 200 to 250 Hz range, rather than grouping them into a single interharmonic value with a frequency of 225 Hz. Therefore, this method enables better identification of the interharmonic components present in the signal under study compared to [15].

5.3. Simulation 3: Noise-Effect Analysis

To study the effect of noise on the results obtained by the proposed algorithm, a white noise signal (30% of the fundamental component, Figure 42) has been generated in order to be added to the input signal. For this analysis, the signal of Figure 42 is added to the input signal used in Point 5.2 for the first part of Simulation 2 (Figure 25), composed of the spectral components summarized in Table 5. The resulting input signal is presented in Figure 43.

Table 9. Harmonic components identified in the input signal.

No. Spectral Component	Voltage rms (V) (w/o Noise)	Frequency (Hz) (w/o Noise)	Voltage rms (V) (w/Noise)	Frequency (Hz) (w/Noise)
1	11.24	30	11.15	30
2	199	50	199.7	50
3	15.79	70	15.99	70
4	22.9	120	22.81	120
5	59.75	150	59.98	150

summarizes the spectral components, as well as their corresponding frequency and magnitude values, obtained once the first part of the algorithm is applied to the output of the DFT. In order to facilitate the comparison, the results obtained with the input signal without noise are presented in Table 6 and are also included in Table 9. The results demonstrate that the effect of the noise does not affect the proposed identification algorithm. The same number of spectral components, as well as the same frequencies, are obtained. Only neglectable differences can be observed in the identified rms voltage values.

Once the first identification is done, LKF starts the estimation of the five components at t = 0.3 s. Following the same procedure as in the previous cases, Figure 44 shows the input signal (in blue) compared with the signal estimated by the Kalman filter (in red), and Figure 45 shows the error as the difference between the two signals. It can be observed that the error signal basically corresponds to the large noise of the input signal. Therefore, to better appreciate the effect of noise on Kalman performance, it is necessary to show the estimated components individually. In this way, Figure 46 shows the estimation (in red) obtained for the first spectral component of the input signal (in blue, 23 Vrms and 33.9 Hz). Figure 47 and Figure 48 present the frequency estimation for the 33.9 Hz component and for all components, respectively.

Comparing these results with those presented in Figure 29 and Figure 30, it can be clearly seen how noise affects the signals estimated by the LKF. In terms of frequency, the estimated signals show an error of less than 1% for all the components studied.

Figure 49 shows the estimation of the rms voltage value of the subharmonic component. It can be seen how, after the LKF stationary is reached at t = 0.5 s, the estimated rms value varies around the real value of the input signal (23 V) with a maximum deviation of 2.9 V (20.1 V) at t = 0.64 s. This result shows an instantaneous maximum error of 12.6%. However, as explained in the previous section, the error condition is computed by using the RMSE (26). In this case, RMSE is 0.65 V for the last 100 ms (t = 0.7 s to t = 0.8 s), which implies an error of 2.84%.

Figure 50 shows the time evolution of the estimated rms voltage value using the LKF of the spectral components that form the input signal.

As a result,

Table 10. Measurement of error from estimation with noise.

No. Spectral Component	Abs Max Error (%)	RMSE (V)	RMSE (%)
1	6.48	0.66	2.84
2	0.22	0.20	0.09
3	4.27	1.01	2.19
4	7.63	0.88	2.56
5	4.02	1.20	1.74

shows the absolute maximum error and the RMSE measured during the last 100 rms of the simulation.

Since noise is not a component that can be modeled in the Kalman filter structure, noise is translated into Kalman states because the recursive Kalman structure attempts to reproduce the original signal. Frequency estimations show fluctuations within a band very close to the real value. The rms estimated values show a higher deviation, with RMSE errors up to 2.84%, which is good enough considering the high noise added to the input signal. Regarding the estimations made for the instantaneous voltage values of the spectral components, they lack good precision. However, Figure 44 demonstrates good results in global estimation, with instantaneous error never reaching the maximum allowable error (MAE) that could have forced the algorithm to restart the identification process unnecessarily.

The problems raised above can be corrected by studying the noise signal and calculating an optimal R matrix for high noise values. On the other hand, a Kalman filter with greater immunity to noise, such as the extended Kalman filter or some other type of hybrid Kalman filter, can be used to improve the noise immunity and the estimation precision.

6. Conclusions

The article presents a method that extracts the frequency and instantaneous voltage values of the spectral components in an unknown alternating voltage signal in real time. The most immediate application is for systems measuring and managing the quality of electrical energy.

This method addresses inaccuracies encountered by the linear Kalman filter when the observed signal model does not match reality, using an algorithm developed by the authors and inspired by the clustering of results of the Discrete Fourier Transform, as shown in the IEC 61000-4-7 standard. The proposed algorithm enables the determination of which values provided by the discrete Fourier transform should be considered real and which are the result of picket fence and spectral leakage effects. This information is used for constructing the appropriate state matrix of the linear Kalman filter. Therefore, it is synchronized with all the harmonic components present in the signal under analysis. The linear Kalman filter provides synchronized output for each harmonic component in real time, supplying the instantaneous voltage and frequency values of each spectral component forming the input signal.

The results obtained by applying the proposed method to various electrical network signals show how it corrects errors in the estimations made with the linear Kalman filter when the number of states used in the construction of the state matrix does not coincide with the number of components that make up the real signal to be estimated. Additionally, the method provides information on the number of spectral components in the interharmonic groups, enabling a better identification of the signals from the discrete Fourier transform performed on the input signal and providing more information than the IEC 61000-4-7 standard. Using the linear Kalman filter enables the individual extraction of each spectral component present in the signal in real time, so it continuously provides the voltage and frequency values of the instantaneous value of that component. The results confirm the validity of the method in systems with harmonic, interharmonic, and subharmonic signals in power distribution networks. The limitations in applicable systems arise from the need for Fourier transform decomposition and the time required for the linear Kalman filter to converge to the real value of the observed states. It has been verified that noisy signals affect the LKF estimations but not the identification algorithm. The results obtained validate the developed method.