A High-Precision and Wideband Fundamental Frequency Measurement Method for Synchronous Sampling Used in the Power Analyzer

In the dedicated high-precision power quality analyzer, synchronous sampling is required to reduce the effect of spectrum leakage produced by the discrete Fourier transform process. Thus, accurate fundamental frequency measurement is urgently needed. However, due to the harmonics and noise in the power signal, it is difficult to achieve the accurate fundamental frequency measurement. Moreover, with the wide application of high-frequency programmable power supply, the fundamental frequency is gradually increasing, which requires power analyzers to have the abilities of both high precision and a wide range of the fundamental frequency measurement. To solve these issues, a new fundamental frequency measurement architecture used in synchronous sampling is proposed. This architecture consists of a small-point fast Fourier transform module, spectrum refinement algorithm, and a multimodal optimization method to calculate the accurate fundamental frequency under large harmonic conditions. In the practical hardware platform results, this architecture has a large fundamental frequency measurement range from 20 Hz to 200 kHz with a relative error which is <0.004%. The wideband fundamental frequency measurement structure proposed in this article achieves high measurement accuracy.


INTRODUCTION
With the rapid development of power electronics technology, high-power electrical equipment such as AC motors and electric vehicles, as well as new energy power generation equipment such as solar power stations and wind power stations are gradually widely used in various fields of industry and residential life. And with the access of a larger number of distributed photovoltaic power stations in the low or medium voltage distribution power grids, the harmonic issue of voltage and current generated by the inverters seriously affects the power supply safety and quality of these power networks (Haitao et al., 2018). Moreover, due to the influence of nonlinear high-power electrical loads such as AC motors, electric arc furnaces, and welding machines, harmonic components and waveform distortion are also produced. At the same time, flicker, voltage dips, supply voltage unbalance, high reactive power, and other disadvantageous effects will also be produced due to the existence of nonlinear loads, which not only leads to the decline of power supply operation efficiency but also reduces the safety and reliability of the power supply system (Jain and Singh, 2011;Repak et al., 2018). Therefore, a dedicated power quality analyzer is needed to monitor the power quality parameters in the power grid to improve the efficiency and reliability of the overall power system. The accurate acquisition of power quality parameters is based on the high-precision sampling of voltage and current signals by a power quality analyzer. Moreover, due to the development of programmable power supply and the wide application of variable frequency drives (VFDs), the fundamental frequency of power signal will be much higher or lower than 50/60 Hz in some specific application scenarios. This means that the bandwidth of the fundamental frequency of the power signal to be measured becomes wider, and the frequency of harmonics also increases. It becomes more and more difficult for the power quality analyzer to accurately capture and calculate the power quality parameters.
Because the sampling process is the core part of the power quality analyzer, the study of the sampling method is of great significance to improve the measurement accuracy of power parameters. The sampling methods commonly used in a power quality analyzer are asynchronous sampling, quasisynchronous sampling, and synchronous sampling. The asynchronous sampling and the quasi-synchronous sampling do not require the sampling rate to be an integer multiple of the signal frequency. When the discrete Fourier transform (DFT) is used to analyze the voltage and current signals sampled by these two sampling methods, it is difficult to accurately calculate the power quality parameters like inter-harmonics and harmonics because of the spectrum leakage and the fence effect (Zhou et al., 2018;Guo et al., 2019). Meanwhile, it is pointed out in the literature (Kuwalek and Otomański, 2019) that the calculation error of the total harmonic distortion (ΔTHD) can reach 4% due to the influence of spectral leakage and fence effect. And to reduce this influence on harmonic analysis and other power parameter calculations, the literature (Jin et al., 2017) processed the sampled data by windowing and interpolation operations, which significantly improved the accuracy of fundamental frequency measurement and harmonic analysis. However, this windowed interpolation method involves complex transformation formulas and heavy computation, which is difficult to apply in real-time power quality monitoring. And in the literature (Tao et al., 2010), the quasi-synchronous sampling method is used to analyze the harmonics. When the power signal frequency varies within ±10%, the relative error of harmonic amplitudes and phase angle is about 0.1%. Similarly, the quasi-synchronous sampling method increases the amount of data and computation in the interpolation section and reduces the operation efficiency.
The synchronous sampling method (Aiello et al., 2007;Jain and Singh, 2011) requires the sampling rate to be an integer multiple of the signal frequency. In theory, this integer cycle sampling of the voltage and current signal realized by the synchronous sampling method can avoid the adverse effects of spectral leakage on the calculation of power parameters (Guo et al., 2019). Song and Yang (2014) studied the harmonic analysis algorithm based on the synchronous sampling architecture. In this study, the relative measurement errors of fundamental and harmonic components are only 0.0019 and 0.003%, respectively. It shows that the synchronous sampling method can significantly improve the analysis accuracy of harmonic components and other power parameters.
To achieve synchronous sampling, high precision acquisition of fundamental frequency is a vital prerequisite. Erm et al. (2019) designed a digital phase-locked loop (PLL) to lock the frequency of the input signal, providing a frequency locking range from 177 to 222 kHz with a locking time of 12.57 us and a delay of 0.8 ns. However, when a PLL circuit is used to measure a signal with a wider bandwidth, it will result in a longer locking time. Meanwhile, the hardware circuit of the PLL is usually more complex and costly to design. When the signal is interfered with noise, harmonics, inter-harmonics, and other adverse effects, the accuracy of the PLL circuit is reduced or even the phase-locking is not possible (Tourigny-Plante et al., 2018). And another simple and direct frequency measurement method is composed of a hardware comparator and a counter, as shown in Figure 1A. In this structure, the highly stable crystal oscillator and a clock synthesizer of Channel A produce a reference signal with a period of T 0 . And the input signal from Channel B is converted to a square signal with an unknown period of T x . Then this unknown square signal is counted by the reference signal to get the specific value of the fundamental frequency. However, this hardware implementation also has difficulty in overcoming the fact that the input signal contains harmonics. As shown in Figure 1B, when the input signal contains harmonics, the actual comparison result of the zero-crossing comparator will produce additional zero-crossing points, and measured period T xe of the fundamental component is less than the real period T x of itself.
As shown in Figure 1, the harmonics caused by nonlinear electrical appliances have been a major obstacle to accurate fundamental frequency measurement. To overcome this barrier, Moreira et al. (2019) designed an online frequency estimation method based on an enhanced least mean square (LMS) method. The fundamental frequency measurement error is lower than 2 mHz under harmonic and inter-harmonic conditions. Mohan et al. (2018) proposed a data-driven method based on the dynamic mode decomposition algorithm to measure the frequency and amplitude in the power grid, and achieved a result with a relative frequency measurement error of less than −2.02e −06 in the 50 Hz power system. Li et al. (2019) adopted the frequency-shifted filtering method to estimate the fundamental frequency of the power grid, and the frequency error is <1e−4 Hz in the case of the test containing the odd harmonics ranging from 3rd to 11th. Borkowski and Kania (2019) presented new multipoint weighted interpolations of the DFT (MWIDFT) frequency estimation method based on generalized maximum sidelobe decay (GMSD) time windows in which the frequency estimation accuracy can be at the level of 10 −3 or 10 −4 Hz. The literature by Yifan et al. (2019) proposed a fundamental frequency measurement architecture implemented on an FPGA chip, and the relative error of fundamental frequency measurement is <0.2%. Also, this method extends the fundamental frequency measurement range, which is from 10Hz to 1kHz. The literature by Chen et al. (2020) designed a software-implemented fundamental frequency measurement method with absolute frequency error <0.02 Hz in the measurement frequency range from 10 to 450 Hz.
Many researchers have studied the accurate measurement strategies of the fundamental frequency under the influence of harmonics, inter-harmonics, noise, and other factors in the power grid. However, most of these studies were focused on the 50/ 60 Hz power grid system. Therefore, to satisfy the application scenarios where the fundamental frequency varies widely and to accurately measure the fundamental frequency at the condition of large harmonic contents, a new fundamental frequency measurement structure has been proposed in this study, which is based on the fast Fourier transform (FFT) method, spectral refinement method, and the improved particle swarm optimization (PSO) algorithm. Compared with the previous literature, this proposed method has a wider fundamental frequency measurement range from 20 Hz to 200 kHz and achieves far less than 0.02% relative error of fundamental frequency measurement in the case of large harmonics.
The rest of this article is organized as follows. The basic knowledge of the proposed method in this article is given in the section "Fundamental Knowledge of the Proposed Method." Then the section "Implementation of the Proposed Method" will give enough details to realize the proposed method. Afterward, the "Experiment and Verification" section shows both simulation results and the performance of this proposed method executed in a hardware platform of a power quality analyzer. Finally, the "Conclusion" section concludes this article.

FUNDAMENTAL KNOWLEDGE OF THE PROPOSED METHOD
In the actual environment of power quality testing, the harmonic components and noise of a power signal are inevitable. Therefore, in order to get over the challenge of fundamental frequency measurement caused by harmonics and noise, this article proposes a high-precision and wideband fundamental frequency measurement method. This synchronous sampling structure integrating the proposed method is shown in Figure 2.
In the synchronous sampling architecture shown in Figure 2, the three-phase four-wire mode is adopted to synchronously sample 4-channel voltage signals and 4-channel current Frontiers in Energy Research | www.frontiersin.org July 2021 | Volume 9 | Article 652386 signals, respectively. Also, the phase-A voltage signal is sampled at a fixed sampling rate by the analog-to-digital convert (ADC), ADC1. Then the fixed sampling data are transmitted to an FPGA chip and a processor chip to calculate the fundamental frequency by the proposed method. After obtaining the fundamental frequency, the synchronous sampling clock will be generated by FPGA to drive the 8-channel ADC chip ADC2 to execute the synchronous sampling process. Based on this synchronous sampling structure shown in Figure 2, the remaining part of this section will introduce the necessary knowledge to realize the proposed fundamental frequency measurement method.

Spectrum Refinement Method
Discrete Fourier transform (DFT) is the basis of spectrum analysis of signals. The DFT of a sampled digital signal sequence x(n) in which the length is N can be represented as Eq. 1: where k ∈ [0, N − 1] and W nk N e −j 2πnk N is called the rotation factor. And according to the fast Fourier transform (FFT) algorithm, the DFT of a signal can be calculated efficiently in both FPGA chips and other processors.
When the FFT algorithm is used to calculate the spectrum of the signal sequence x(n), the frequency resolution Δf is defined as Eq. 2: Δf f s N . (2) In Eq. 2, f s is the sampling rate. It can be seen that decreasing the sampling rate f s is a useful way to improve the frequency accuracy of the FFT spectrum analysis. However, according to the Nyquist sampling law, the reduction of the sampling rate f s means that the bandwidth of signal measurement is reduced as well, which is hard to meet the range of fundamental frequency measurement requirement from 20Hz to 200 kHz required in this article. On the other hand, increasing the number of sampling data used to process the FFT method is a feasible way to meet the demand of the fundamental frequency measurement resolution if the sampling rate is fixed. But the increased calculation number will cause the unacceptable time consumption of calculation. Moreover, these two possible ways still cannot radically get rid of the limitation of spectral resolution shown in Eq. 2.
Therefore, it is necessary to refine the spectrum based on the signal spectrum obtained by the FFT algorithm with a small number of calculating points N. In this way, the large frequency measurement bandwidth and the calculation efficiency can be satisfied simultaneously.
Based on the frequency resolution Δf defined by Equation 2, the frequency f k at the kth frequency point in the spectrum can be computed as Eq. 3: According to the Euler formula, Eq. 1 can be divided into the real part and the imaginary part, and the index k can be replaced by the corresponding frequency value f k calculated by Eq. 3. Then Eq. 1 can be rewritten as Eq. 4.
X R (f k ) in Eq. 4 is the real part of the spectrum of a signal, and X I (f k ) is the imaginary part. The spectrum after the transformation of Eq. 4 is still discrete. To obtain the continuous spectrum, the spectrum is needed to be infinitely refined. Since the DFT has the frequency resolution Δf , it can be considered that the real frequency value of the signal is in the neighborhood in which the central frequency is f k and the radius is Δf 2 ; that is, the real frequency f real of the sampled data is in the range of f k − Δf 2 , f k + Δf 2 . Therefore, in this neighborhood, the frequency value f k at the kth point in the spectrum obtained by the FFT method is replaced by a continuous changing frequency variable f. And Eq. 4 can be rewritten as follows: where And the amplitude of the spectrum at frequency point f can be obtained as follows: X(f ) in Eq. 6 represents the amplitude of the spectrum at a frequency point f.
Because the frequency variable f is continuously changing in the neighborhood of f k , which is the frequency at the kth spectral line obtained by the FFT method, the spectrum obtained by Eqs. 5 and 6 is also continuous on the frequency axis. In the corresponding neighborhood, the infinitely refined spectrum can be obtained by using this method, which can get rid of the frequency resolution limitation of the DFT method.

Close Neighbor Mobility Particle Swarm Optimization Method
A rough fundamental frequency value f br can be obtained by using the FFT algorithm analyzing the spectrum of the fixed-rate sampling data by ADC1 shown in Figure 2. Then the spectrum is refined according to Eqs. 5 and 6 at a neighborhood of f br . In the field of power parameter analysis, the energy of fundamental frequency is the largest compared with other harmonic components. Therefore, by adopting the frequency refinement method mentioned in the preceding subsection, the problem of fundamental frequency measurement has been transformed from the calculation of DFT to an optimization task in a certain interval. In other words, based on the rough frequency value f br calculated by the FFT method, a optimization method is used to find out a frequency point f * in the search region f k − Δf 2 , f k + Δf 2 to make Eq. 6, which is the amplitude of the DFT at this frequency point, reach its maximum value. In order to reduce the time consumption of the coarse frequency acquisition, a small-point FFT calculation will be adopted. Therefore, the frequency resolution Δf calculated by Eq. 2 will be large. When the fundamental frequency is too low, some harmonic components may fall in the search region And the maximum optimization problem becomes a multimodal optimization issue, as shown in Figure 3.
In Figures 3A and 3B, the two spectrum refinement results are based on the rough fundamental frequencies calculated by the 4,096 point FFT algorithm at the 2MSPS fixed sampling rate. The frequency resolution Δf is about 488 Hz at this configuration. Figure 3A shows the refinement result of phase-A voltage signal, where the fundamental frequency is 20 Hz, and this voltage signal is superimposed with 2nd, 3rd, 4th, and 5th harmonics. Since the frequency resolution Δf is much larger than the interval between fundamental frequency and harmonic frequencies, all the four harmonic components still appear in the refined spectrum acquired by Eq. 5. Compared with Figure 3A, the Figure 3B shows the refinement result of a signal with the fundamental frequency of 1 kHz and the same harmonic configuration. And in this condition, there will not be any harmonics in the refinement region. From this comparison, it can be seen that when the fundamental frequencies are low, obtaining the fundamental frequency by the refinement method becomes a multimodal optimization problem.
In order to solve the optimization problem, many researchers have studied this problem. The literature by Li et al. (2021) proposed a model-free optimization method based on the deep reinforcement learning algorithm which is used to optimize the control strategy for virtual synchronous generator. And the literature by Yushuai et al. (2020) proposed an improved Newton descent algorithm to the cooperative energy management task in the energy Internet. However, the optimization task in this article is based on the spectrum refinement algorithm, as shown in Eqs. 5 and 6, and is a model-based optimization task. Eqs. 5 and 6 follow the optimization model. And in some conditions mentioned above, the optimization task will convert to a multimodal optimization issue. To obtain the global maximum in the multimodal optimization problem which is caused by the frequency refinement method at low frequency and avoid falling into the local optima, this article adopts a method called the close neighbor mobility method (Zou et al., 2020) (CNMM) which is based on the PSO method. The PSO method is a swarm intelligence optimization algorithm inspired by birds' foraging behavior. Each individual in this bird population is called a particle, and the position vector of each particle represents a solution in the solution space of a certain problem. The purpose of each particle is to keep approaching the position of food which is also called the global optimal solution. In order to achieve this goal, particles are constantly dominated by the personal best position (pbest) of themselves and attracted by the global best position (gbest) of the whole population.
The dimension of the solution space is assumed as m, which means the dimension of each particle is m, too. And the size of this population is N. Then the position vector of each particle can be represented as N]. Also, the velocity vector of each particle can be represented as v i {v i1 , v i2 , . . . , v im } in the same structure. Based on the above hypothesis, the velocity and position of the i-th particle in the dimension d can be updated at the time of t + 1 by Eqs. 7 and 8.
In the velocity update Eq. 7, w is called the inertia weight, which is used to illustrate how much velocity information the current particle can retain. A large w means better global searchability; however, the local search ability will be declined. The c 1 and c 2 are called the self-learning coefficient and the group-learning coefficient, respectively. The particle is greatly attracted by its personal best position when the self-learning factor c 1 is large. Conversely, the particle is more easily affected by the global best position when c 2 is large. And the parameters r 1 and r 2 are two random numbers in the range [0, 1]. The personal best position pbest i (t) of a particle records the best position of itself from initialization to current time t. And the global best position gbest(t) records the best position of the whole population until now. And the fitness which is calculated by the aim function is used to judge how best the particle is. In this study, Eq. 6 is regarded as the aim function or also called the fitness function. And the larger fitness of a particle means this particle is good and has a larger probability of becoming the global optimal.
And to improve the multimodal optimization ability of the basic PSO method described by Eqs. 7 and 8, the following strategies are needed to be described.

Elite Selection Strategy
The main purpose of the elite selection strategy is to find out the particles which have excellent convergence ability, adaptability, and distribution, and to retain these particles which have good performance. And then, these elite particles will guide the search behavior and evolution of the remaining particles.
The procedure of the elite selection strategy is shown in Figure 4A. To begin this process, the particles in the whole population P are arranged in a descending order according to these fitness values, which are calculated by Eq. 6. And the collection P sort is created to represent the sorted population. After the arrangement, the particle which has the largest fitness value is added to the elite population S. And its fitness value is recorded as fb. Then the remaining particles in the sorted set P sort where the differences between their fitness value and fb that are less than a threshold e are considered as the candidate elite particles. Then the Euclid distances defined by Eq. 9 are calculated between each candidate elite particle and each particle in the elite set S. The candidate particle in which the Euclid distance is larger than the threshold r will be added to the elite population S.
In Eq. 9, d represents the dimension of a particle.
To illustrate the elite selection strategy more vividly, we assume there are 9 particles in the population P, and the aim function is a function with three identical peak points. Through the elite selection strategy, Figure 4B can be obtained. In Figure 4B, the circles represent all the particles. It can be found that particle F has the largest fitness, so this particle F is added to the elite population S. Its fitness value is recorded as fb. And particles B, H, and I are the candidate elite particles because their differences of the fitness value to fb is less than threshold e. Then the Euclid distances defined by Eq. 9 between particles in set S which is particle F and candidate particles which are particles B, H, and I are calculated. From Figure 4B, particles H and I are too close, which means particle H is dominated by particle I. And distances among particles B, C, F, and I are larger than the distance threshold r. Therefore, the final elite population S contains the four particles B, C, F, and I. And these four particles have a greater fitness value which means they are good particles and more nearer to optima than others. Also, they are not too close to each other and that indicates the four particles are independent.

Close Neighborhood Strategy
The close neighborhood strategy is based on the concept of a niche in biology. And in biological evolution, organisms usually live and breed with their own species, and form multiple subpopulations. And this condition increases the diversity of the population and improves the ability of local searching.
To improve the local searching ability of the PSO method, a close neighborhood strategy shown in Figure 5A is used instead of the global best position gbest(t) used in the velocity update Eq. 7. For any particle I, this strategy first calculates the Euclid distance between this particle I and the remaining particles according to Eq. 9. Then the remaining particles will be sorted at the ascending order of Euclid distance, and L particles closest to particle I are selected. Finally, the particle with the largest fitness value in these L particles will be found, and this particle is used as the global best particle of particle I and represented as gbest i (t), (i I). Then the velocity update function 7 can be rewritten as Eq. 10.
In this equation, the global best particle gbest(t) is replaced by the neighbor global best particle gbest i (t) of particle I.  As shown in Figure 5B, we assume there are 9 particles in the population. And in the current iteration, particle F is the particle with the largest fitness value. And according to the traditional PSO method, all particles are attracted by particle F in the next iteration. In the close neighborhood strategy, taking particle D as an example, we assume that L equals 3. Then particles B, C, and E are found, which are three particles closest to particle D. And from Figure 5B, it can be seen that particle B has the largest fitness value in this neighborhood. So B is regarded as the neighbor global best particle gbest i (t) to particle D, and particle D needs to converge to particle B in the following iteration.
Unlike the classical PSO method, where all particles converge to one global best position, the close neighborhood strategy allows the particles to continuously move to the best position in their neighborhood. And this strategy helps the PSO method avoid converging to one local optimum in a multimodal situation.

Differential Evolution Strategy
The differential evolution (DE) strategy uses mutation, crossover, and selection operations to generate a new population from the current population. By these operations, the PSO method can expand the search region and improve the global optimization ability.
First, the mutation operation is realized using Eq. 11. v where i ∈ [1, N] and N is the size of the particle population.
Parameter v i is the mutation result of the i-th particle. The three particles x r1 , x r2 , x r3 are randomly selected from the population (r1, r2, r3 ≠ i), which do not belong to the elite set S. And in Eq. 11, parameter F is a positive real number which controls the strength of the mutation operation. Based on mutation result v i of each particle i obtained from Eq. 11, the crossover operation is performed by Eq. 12.
where j ∈ [1, m] is the dimension of particle i. And the rand in Eq. 12 is a random number with the uniform distribution from 0 to 1. And CR is a crossover control parameter. And when the random number rand ∼ U(0, 1) is less than the parameter CR, particle u i at dimension j will be replaced by the mutation result at the same dimension. The idea of this operation is to probabilistically exchange some dimensions between mutation result v i and current particle x i to form a new particle u i . By creating new particles, the diversity of the particle population and the search area is increased. Thus, the probability of converging to the global optimum is improved. After finishing the mutation and crossover operation, a selection process is used to form the final new population according to Eq. 13.
In Eq. 13, x * i is the i-th particle in the new population. The selection operation is processed between the crossover new particle u i and current particle x i , and the particle which has a larger fitness value will be selected to the new population of the next generation.
This differential evolution strategy generates the new population by mutation, crossover, and selection operation on the current population. This strategy improves the abundance of individuals in the population and expands the search scope.
Based on the above strategies, the procedure of the complete CNMM PSO algorithm is shown in Figure 6. At the beginning of this procedure, the initial population P is generated randomly, and each particle in P has its own velocity vector v i and position vector x i . Then the initial personal best position is assigned by their initial position x i . After the initialization stage, the outer loop is entered as the maximum iteration number is not reached, and the outer loop represents the evolution times. At the beginning of the outer cycle, the elite set S is found by the elite selection strategy. Then the close neighborhood strategy is used to find out the global best positions of each particle to form set G. Furthermore, the inner particle update loop is entered to update the velocity and position vectors of each particle in population P. If the current particle is an elite particle in S, the particle will keep their characters and not be updated. Otherwise, particles will be used to calculate the fitness value, and the personal best position pbest(i) will be updated. According to Equations 10 and 8, the new velocity and position vectors of current particle can be obtained. Before recording these new vectors, the boundary check process must be done. If the new vectors satisfy the boundary requirement, the new velocity vectors are used to update the position vector of particles. Otherwise, the boundary will be used to replace the undesired velocity vector. And after updating all the particles in population P, if the number of generations reaches a specific number, the DE strategy will be performed. When the number of iteration reaches the setting value, the CNMM algorithm is completed. Finally, the position vectors x i of the whole population P are output as the final solution.

IMPLEMENTATION OF THE PROPOSED METHOD
Based on the spectrum refinement method and the CNMM PSO algorithm described above, the high precision and wide range fundamental frequency measurement architecture proposed in this article are shown in Figure 7.
In Figure 7A, the structure consists of an FPGA chip and a processor. The phase A voltage signal is sampled by ADC1 at the sampling rate of 2MSPS. Then in the FPGA chip, the sampled digital signal is divided into two channels, one is the FFT coarse frequency measurement branch and the other is the down-sampling branch. In the FFT channel, a Blackman window is used to preprocess the sampled data. The peak amplitude of the sidelobe of this Blackman window is only −57 dB, which can effectively reduce the spectrum leakage at the fixed sampling rate condition. After this preprocessing operation, an FFT module with N 4,096 points is used to obtain the discrete spectrum sequence of the sampled digital signal. Since the power signal has the prior information that the energy of the fundamental frequency component is greater than that of other harmonics, the spectrum index k corresponding to the frequency component with the largest amplitude can be found through a simple maximum value searching module called the "Max Search" in Figure 7A. Therefore, the coarse fundamental frequency f k can be obtained by Eq. 3.
And another channel is the down-sampling channel which reduces the sampling rate. Then the decimated digital signal is transferred to the processor to complete the spectrum refinement method and the CNMM PSO method. The target of this fundamental frequency measurement method proposed in this article is from 20 Hz to 200kHz, and the analog-to-digital converter ADC1 works at the fixed sampling rate mode. Therefore, according to the Nyquist sampling law, the sampling rate of ADC1 needs to be twice the maximum signal frequency at least. However, the high sampling rate is unnecessary and redundant for the low fundamental frequency measurement conditions. Therefore, it is necessary to reduce the sampling rate of the sampled digital signal by the decimation method to improve the calculation efficiency when the ADC1 is working at the fixed sampling rate. And according to the rough frequency measurement result f k , the different decimation coefficient M is selected. Also, to avoid the spectral aliasing caused by decimation, the low-pass finite impulse response (FIR) filter is used as the anti-aliasing filter which is implemented in the FPGA chip. Since different decimation coefficients require different coefficients of the FIR filter to form low-pass filters with different cutoff frequencies, the down-sampling factor M is classified into several groups shown in Table 1 according to the spectrum index k calculated by the FFT module.  In terms of the cutoff frequencies required for the FIR low-pass filter listed in Table 1, the corresponding filter coefficients can be designed by some software design tools, like Matlab. And through the spectrum index k found by the "Max search" module, the corresponding FIR filter coefficients and decimation factor M are found in Table 1 and brought into the "FIR filter" and "down sampling" modules shown in Figure 7 to complete the decimation stage.
After the decimation process, the decimated digital signal and the rough fundamental frequency f k are sent to the processor chip from FPGA by the PCIe interface. In this processor, the search region [f k − Δf , f k + Δf ] is calculated first, and it is found that the search region is a little larger than the theoretical one mentioned in the last section. And the frequency resolution Δf is calculated by Eq. 2, and f s is the sampling rate of ADC1 which is not decimated and equals to 2MSPS in this article. After determining the search region, the CNMM PSO algorithm with specific iteration times will be executed. And the fitness value of each particle is calculated according to Eqs. 6 and 5. And in the fitness computing process, the sampling rate f s in Eq. 5 is the sampling rate after the decimation, rather than the original sampling frequency of ADC1. When the requisite number of the iteration of the CNMM PSO method has reached, the fitness values of the final population P will be calculated again. The position vector x i of the particle in the final population which has the largest fitness value is regarded as the result of the fundamental frequency measurement. The complete flow of the proposed method is shown in Figure 7B.

Simulation
In this section, the proposed fundamental frequency measurement method was performed in the Matlab programming environment. A simulation voltage signal was used to test the functionality and correctness of this proposed method.
The simulation signal is constructed according to Eq. 14.
The simulation signal formed by Eq. 14 consists of n frequency components, of which A 1 and f 1 is the amplitude and frequency of the fundamental component, and the remaining part, i ∈ [2, n], are the amplitude and frequency of the harmonic components, respectively.
First, the fundamental frequency f 1 is set to 50 Hz to verify the frequency measurement accuracy of the most commonly used power frequency. And assuming that the amplitude A 1 of the fundamental frequency component is 1V, n equals 5, which means the simulation signal contains harmonics from the second order to the fifth order. And the amplitude of each harmonic component is 0.5V, which is a test environment with more severe harmonic content. Then according to the simulation structure implemented on a computer shown in Figure 8, the test of 50 Hz fundamental frequency measurement is executed. In Figure 8, the maximum velocity of each particle is defined as 10, which balances the convergence speed and the convergence accuracy. And based on the parameter settings in Figure 8, the proposed fundamental frequency measurement method is verified for the signal with a fundamental frequency of 50 Hz and containing harmonics from the 2nd order to the 5th order. And the simulation result is shown in Figure 9. Figure 9A shows the coarse frequency measurement result after the 4,096 point FFT calculation of the simulation signal. Since the frequency resolution of this FFT module is about 488.3Hz, the spectrum calculated by the FFT module has the maximum spectrum amplitude at k 0, which means the fundamental component and the harmonics are mixed into one spectral line. After obtaining k 0, the rough fundamental frequency f k can be calculated by Eq. 3, which is 0Hz. The frequency range needed to be refined can be calculated, which is in the range from −244.15 to  Figure 9B, it can be seen that as the frequency resolution is larger than the frequency interval between harmonics, the refined spectrum still has 3 harmonic components. And in this situation, it represents the problem of the multimodal optimization. Based on the refined spectrum, the convergence result of the CNMM PSO method is shown in Figure 9C. The red circles in Figure 9C are particles in the population. It can be seen that after 100 iterations, all the 150 particles converge to the region very close to 50 Hz. And the particles do not fall into the local optimal solution which represents the harmonic components. Then the fitness values of these final 150 particles are calculated again. The position x i of the particle which has the largest fitness value is considered as the fundamental frequency. Thus, the final calculated fundamental frequency is 49.99976 Hz, and the absolute error is only −0.00024 Hz.
After verifying the capability of the proposed algorithm for 50 Hz fundamental frequency measurement, we examine the influence of the signal-to-noise rate (SNR) of the input power signal on the fundamental frequency measurement performance. In this test scenario, the SNR of the simulation signal constructed by Eq. 14 changes from 3 to 60 dB. A group of simulation signals with different fundamental frequencies are tested. The number of harmonic components of each simulation signal is 4, which contains harmonics from the second order to the fifth order. And the amplitude of the fundamental component and harmonics of each simulation signal are 1 and 0.5 V, respectively. Figure 10A demonstrates the SNR performance of this proposed fundamental frequency measurement method. The performance is evaluated by the absolute value of the relative error, which is calculated by Eq. 15.
In this Eq. 15, the error means the absolute value of the relative error, and f measure is the fundamental frequency measured by the proposed method. The parameter f real is the setting value of the simulation signal.
In the SNR range from 3 to 60 dB, the absolute values of the relative error of each fundamental frequency are <0.02%, as shown in Figure10A. Even if the SNR of the simulation signals is 3 dB, the maximum relative error is only 0.0114%. These results shown in Figure 10A indicate that this proposed method has the ability to overcome the changing SNR of the input power signal.
In this next test case, the fundamental frequency measurement accuracy of the proposed method under different harmonic amplitudes is evaluated. In this set of simulation signals, each signal contains 4  Figure 10B, it can be seen that the relative error of each fundamental frequency is lower under the small harmonic amplitude. The relative errors increase with the increase of the harmonic amplitudes. This is because the attraction of local optima increases with the increase of harmonic amplitude. Taking the fundamental frequency of 50 Hz as an example, as shown in Figure 9C, the amplitude of the global optimum, the 50 Hz, is fixed. And the amplitudes of other local optima, which are 100, 150, and 200 Hz, vary from 0.1 to 0.7 V. When the amplitudes of harmonics are small, the amplitudes of local optima are small, which means the fitness values of particles which are in the position close to the local optima are small, too. In this situation, the global optimum is very attractive for particles. However, when the amplitudes of harmonics are gradually increasing, the fitness values of particles near the local optima are gradually as well. The attractiveness of local optima increases. And the speed of particles moving to global optimum will decrease. After a fixed number of evolution, the particles are not as close to the global optimum as the case of small harmonic amplitude. The distances between global optimum and particles are the frequency measurement error. Although the relative errors increase when the harmonic amplitudes are large, the relative error is <0.02%.
In the third test scenario, the relationship between the number of harmonic components and measurement accuracy is examined. The amplitude of the fundamental component is set to 1 V. And the amplitude of each harmonic component is 0.5 V. Then a group of simulation signals for which the number of harmonic varies from 2 to 30 is used to evaluate the performance of the proposed method. From Figure 10C, it can be seen that the relative error is <0.01% in the whole range of the number of harmonic components varying from 2 to 30. Based on the above examination, the simulation signals with a fundamental frequency in the range of 20 Hz to 200 kHz and with the harmonics of the 2nd order to the 5th order are tested. The amplitudes of both fundamental component and harmonic components are configured in the same manner as the simulation signal with a fundamental frequency of 50 Hz. And the test results are shown in Figure 11. In this Figure 11, the two red lines are absolute error boundaries which are calculated by the relative error requirement of 0.02%*fb. And fb is the theoretical true value of the fundamental frequency measured. From Figure 11, in the measurement range from 20 Hz to 200 kHz, all the absolute errors are within these boundaries, which means the measurement accuracy of this proposed method is not affected by the measurement bandwidth and satisfies the design goal mentioned in the Introduction section.

Hardware Verification
The hardware verification platform based on the proposed high precision and wide range fundamental frequency measurement method is shown in Figure 12. And Figure 12A shows the power analyzer used for the verification of the proposed method in this article. And in this power analyzer, four identical acquisition boards as shown in Figure 12B are used to realize the current and voltage signal acquisition test of the three-phase four-wire power system. And in Figure 12B, the acquisition board of phase-A consists of a pair of voltage signal input ports, the analog signal conditioning circuit of the phase-A voltage signal, the 16-bits ADC1 working at a fixed sampling rate of 2MSPS, and the FPGA chip. And Figure 12B shows the connection relationship of the hardware verification platform between the signal source and the power quality analyzer. In this connection relationship diagram, the red lines represent the signal path of this proposed fundamental frequency measurement method. After completing the calculation of the coarse spectrum for the FPGA chip using the FFT module, and the decimation by the FIR filter module and down sampling module, the processed data will be transferred to the industrial personal computer (IPC) through the PCIe interface. Then another processor module in the IPC will execute the CNMM PSO method to accurately calculate the fundamental frequency through these decimated data.
And based on the hardware verification platform, the power signals with fundamental frequency varying from 20 Hz to 200 kHz are used to test the performance of this proposed fundamental frequency measurement method. The SNR of the test signal is 30 dB. The configuration of all other parameters is  the same as that of the simulation subsection. The typical actual measurement results are shown in Table 2.
From Table 2, it can be seen that the largest relative error is 0.003950% at the fundamental frequency setting value of 20 Hz. And all the hardware verification results show that the relative error of this proposed method is far less than the target value mentioned above. This proposed fundamental frequency measurement method in this study effectively improves the frequency measurement accuracy and has great significance in synchronous sampling applications.
After verifying the measurement accuracy of this proposed method on a practical hardware platform, Table 3 compares several fundamental frequency measurement methods used in the power quality analysis. Method 1 by Li et al. (2019), method 2 by Chen et al. (2020), and method 3 by Yifan et al. (2019) are used as benchmarks to illustrate the performance of this proposed method. From Table 3, it can be seen that method 1 has the least relative error of fundamental frequency measurement. However, this method only focuses on the application scenarios of 60 Hz power grid standard. And the amplitudes of the harmonic components are just 10% of that of the fundamental component in test process, which is far <70% of that of the proposed method. Methods 2 and 3 have extended the measurement range of the fundamental frequency in varying levels. These two methods are validated only in test cases with fewer harmonic components. And the maximum amplitude of harmonics is only 10% of that of the fundamental frequency. The proposed method has been examined in the environment in which the amplitudes of harmonics are 70% of the amplitude of fundamental frequency or the number of harmonics is 30. And the worst relative error of the proposed method is no greater than 0.0114%. Nevertheless, the worst relative error is obtained in the situation that the SNR of the input phase-A voltage signal is only 3 dB, which can hardly be seen in the actual test environment. In the same SNR condition, for instance, 50 dB, the relative error of this proposed method is not >0.002%, as shown in Figure 10A, which is better than that of other methods.

CONCLUSION
This article sets out to improve the accuracy of fundamental frequency measurement under large harmonic conditions to meet the requirement of the synchronous sampling architecture used in the power quality analyzer. And to achieve this goal, in this article, the spectrum refinement strategy is studied first to overcome the disadvantage of the limited frequency resolution of the DFT algorithm. Then a modified PSO method is analyzed to resolve the multimodal optimization problem in fundamental frequency acquisition under large harmonics situation. After researching these two main components, the entire fundamental frequency measurement structure has been designed. And the measurement convergence ability and accuracy of this proposed method have been assessed.
This method has a high fundamental frequency measurement range, which is from 20 Hz to 200 kHz. And this wide-bandwidth measurement ability helps to test the performance of some programmable power source and variable frequency drivers. And in the premise of a wide measurement range, this proposed method achieves high precision that the largest relative error of fundamental frequency measurement in the practical hardware platform test is only 0.00395%. This minimal relative error of fundamental frequency measurement provides a basis for the high-performance synchronous sampling method. Also, in the situation of low SNR, a large number of harmonic components, or large amplitudes of harmonic components, this method can maintain high accuracy of fundamental frequency measurement.