Design of Digital Constrained Linear Least-Squares Multiple-Resonator-Based Harmonic Filtering

: Although voiced speech signals are physical signals which are approximately harmonic and electric power signals are true harmonic, the algorithms used for harmonic analysis in electric power systems can be successfully used in speech processing, including in speech enhancement, noise reduction, speaker recognition, and hearing aids. The discrete Fourier transform (DFT), which has been widely used as a phasor estimator due to its simplicity, has led to the development of new DFT-based algorithms because of its poor performance under dynamic conditions. The multiple-resonator (MR) ﬁlter structure proposed in previous papers has proven to be a suitable approach to dynamic harmonic analysis. In this article, optimized postprocessing compensation ﬁlters are applied to obtain frequency responses of the transfer functions convenient for fast measurements in dynamic conditions. An optimization design method based on the constrained linear least-squares (CLLS) is applied. This way, both the ﬂatness in the passband and the equiripple attenuation in the stopband are satisﬁed simultaneously, and the latency is reduced.


Introduction
Signals that consist of a sum of sine waves whose frequencies are integral multiples of the lowest frequency (so-called fundamental) are said to be harmonic. Many physical signals are approximately harmonic. Examples include voiced speech and other biological signals, musical waveforms, helicopter and boat sound waves, and outputs of nonlinear systems excited by a sinusoidal input [1,2]. Many different approaches are possible; some basic underlying principles of speech enhancement techniques capitalize on the observation that waveforms of voiced sounds are periodic, with a period that corresponds to the fundamental frequency. One approach is based on comb filtering to pass the harmonics of speech but reject the frequency components [3][4][5]. In addition, there are many signal analyses and transformations that are best performed in the frequency domain [6][7][8][9][10][11][12][13][14].
Voltage and current signals in power systems are examples of true harmonic signals. The frequency spectrum of the audio signal is still in the low-frequency range and approximately 2-3 times wider than of the electric one, which states that algorithms developed for usage in electric power systems can be also applied in audio technology.
Due to the growing presence of harmonics in transmission and distribution networks, the importance of estimating harmonic phasors is an increasingly important task. Their fast measurements are very convenient for use in power system protection applications. The extension of synchrophasor algorithms to provide phasor measurement units (PMUs) with the capability of accurately estimating harmonic phasors has been recently carried out in order to facilitate the proliferation of PMU applications in active distribution networks [15][16][17][18][19][20][21]. While more and more PMUs can provide harmonic phasor values together with fundamental values, the IEC/IEEE Standard 60255-118-1:2018 [22] deals with the fundamental phasor only.
One of the primary mathematical tools used to decompose harmonic signals is the fast Fourier transform (FFT), which is an optimized discrete Fourier transform algorithm (DFT). The DFT algorithm has been widely used as phasor estimator due to its simplicity; regrettably, its poor performance under dynamic conditions [23,24] compromises compliance with the standard, leading to the development of new DFT-based algorithms such as interpolated DFT (IpDFT) [25,26]. Additionally, different techniques such as the Taylor series [27], wavelet transform [28], recursive least-squares [29], Kalman filter [15,30], and combined filter design concept [31] have also been proposed to provide solutions that improve the performance presented in traditional DTF-based techniques under dynamic conditions.
Timely measurements are critical for addressing many challenges associated with power system operation. The maximum reporting latency is perhaps one of the more restrictive requirements. It should have in mind that in many algorithms, the time reference is usually centered in the middle of the observation interval, which means an intrinsic algorithm latency is at least half of the observation interval width. In addition to that, a harmonic analysis under dynamic conditions is a growing necessity. The discrete Taylor-Fourier transform (TFT) extends and improves estimations obtained by the DFT by using a dynamic model of the signal [32]. As a result, the obtained reconstruction is more accurate than the reconstruction obtained through DFT, which is a mainstream approach in harmonic analysis. When the spectral density of narrow bandpass harmonic signals is limited to the flat-gain harmonic intervals, the dynamic model-based estimators assure good estimates of the first derivatives of complex harmonic envelopes.
Based on the multiple-resonator (MR) estimation structure, a simple recursive algorithm for dynamic harmonic analysis with some advantages over DFT/FFT, together with benefits in the reduced computational complexity, is proposed in [33]. The good sensitivity properties are assured by the infinite loop gain at the resonator frequencies. In addition, multiple zeros provide reinforcement of the required attenuations and zero-gain flatness at the harmonic components together with a high overall attenuation in the stopband(s). While the rise of the resonator multiplicity improves the amplitude response and selectivity, it adversely influences the phase response and increases the latency. Due to this feature, large values of the resonator multiplicity could be inconvenient in the control application. In order to obtain an algorithm that can be utilized in a wide range of signal dynamics in a unified way and to improve the frequency response that allows a tracking-mode harmonic estimation technique that is both accurate enough in nonstationary conditions and fast enough, a linear combination of the differentiators' outputs in the resonator cascade has been used [34][35][36]. The order of the resulting compensation filter was low and equal to the pole multiplicity. In [37], the proposed approach was generalized to any necessary order through the postprocessing compensation FIR filters applied to the output signals obtained by the basic MR structure.
A similar approach that assures equiripple frequency responses is presented in this paper. The initially formulated objective (cost) and constraining functions are nonlinear. These functions are quadratic, so convex semi-infinite programming can be applied [38]. Through the convenient linearization of the objective and constraint functions, which already was applied in various articles [35][36][37], the design problem is solved by a constrained linear least-squares (CLLS) algorithm. The difference in reference to the approach proposed in [37] is that, in this article, the equiripple-limited amplitude-frequency response in the stopband is provided.

Design Method
In this section, a design of the postprocessing compensator filters for the zeroth derivative of the harmonics for K type estimator is shown. The structure consists of two parts; see Figure 1. The first part is the basic one already described in [33,34]. The structure consists of 2M + 1 branches consisting of K + 1 resonators (having poles on the unit circle) in the cascade. Each resonator in the mth cascade is associated with complex gains g m,k , m= −M, . . . , 0, . . . , M, k = 0, 1, . . . , K. All these resonator cascades operate in parallel within a common feedback loop. At the frequencies corresponding to the resonator pole frequencies, this loop has infinite loop gain; therefore, the transfer value equals 1 independently of the other parameters of the system. The overall system order is 3N, N = 2M + 1 where Mω 1 < π; ω 1 = 2π f 1 / f S is an angular frequency of the fundamental component; f S is a sampling frequency; and M is the highest order of harmonic to be analyzed. The second part is postprocessing of the resonators' output signals. The postprocessing can be applied selectively for some particular signals, e.g., the fundamental component and dominant harmonic, whose estimation are potentially important for control applications. For dead-beat estimators, we have [33]: Calculation of gains g m,k and g m,k , k = 0 . . . K is given in [33,34,39].
In the postprocessing part, the compensation filter Q m z −1 is applied to the output signal V m z −1 = V m,0 z −1 . The resulting transfer function of the complete estimation filter (including the compensation) is where The part T m,0 z −1 = g m,0 z −1 P m z −1 originated from the basic resonator structure. It is actually the FIR filter with zeros obtained from resonator poles by the common feedback. The polynomial P m z −1 consists of all resonator poles, excluding poles z m .
Equation (4) can be written in a matrix form as follows: To assure a unity gain in the harmonic frequency, the following condition has to be satisfied: where Complex equality constraints (6) can be written as

Total Vector Gradient (TVG) Calculation
The latency of the estimator is directly related to the group delay (GD) of the filter. Due to the complexity of the estimation of the GD, in this paper a gradient of the transfer function dT Q m,0 z −1 /dz is applied, herein called a total vector gradient (TVG). It could be proved that in the flat band with small amplitude changes, TVG and GD are proportional, and optimization of one leads to the optimization of the other.
The first derivative of the transfer function T Q m,0 z −1 is where Equation (9) can be written in a matrix form as follows: where a m,2 (z) = a m,0 (z) a m,1 (z) · · · a m,N Q −1 (z) a m,N Q (z) ,

Constrainting Conditions Linearization
In the stop-and transition bands, the following condition has to be satisfied: where l m,i , i = 1, 2, · · · , N F is a settled limit in the frequency point i.
If we need to limit |TVG|, the following condition should be satisfied.
Both (11) and (12) are nonlinear and have to be linearized. For that purpose, let us define the following vectors for harmonic m and frequency point z i , depending on which function is limited: g m,i = a m,1 (z i ) or g m,i = a m,2 (z i ) Equations (11) and (12) can be linearized [35][36][37] and written in a matrix form Re g m,i cos α i2 + Im g m,i sin α i2 . . .
A number of angles L should be settled depending on the desired tolerance. By taking a sufficient number of angles, the criterion (11) can be approximated by the system (14) with the desired accuracy. Square and octagon approximations (L = 4 and L = 8, respectively) are shown in Figure 2. In this article, L = 16 is used.
Using matrix notation and collecting inequality linearization systems in all settled frequency points, (14) becomes the following linear form: where matrix A m and vector b m are given by

Sum of Squares Calculation
Let us define the following vector for harmonic m and harmonic frequency point z m , depending on which function it is: An objective is to find a minimum of the sum of squares of absolute values of h m,i q m in the assembly of the N F selected frequencies subject to the vector q m : If we apply the following equality |h m,i q m | 2 = Re 2 {h m,i q m } + Im 2 {h m,i q m } = ||C m,i , q m || 2 2 (18) where C m,i = Re{h m,i } Im{h m,i } , (17) can be written in a matrix form: where

Constrained Linear Least-Squares (CLLS) Model
The constrained linear least-squares (CLLS) is an optimization problem that deals with the maximization or minimization of a linear function called the objective function subject to linear constraints. Summarizing (19), (15) and (8), the CLLS problem is formulated as follows: min Among a variety of different optimization scenarios, here, two optimization tasks are formulated. Task 1 considers the optimization of |TVG| in the passband with a limitation of T Q m,0 z −1 in the stop-and transition bands [h m,i = a m,2 (z i ) and g m,i = a m,1 (z i )], while Task 2, similar to [37], optimizes the sum of squares of T Q m,0 z −1 in the stopband subject to the limitation of |TVG| in the passband and T Q m,0 z −1 in the transition band [h m,i = a m,1 (z i ), g m,i = a m,2 (z i ) in the passband, and g m,i = a m,1 (z i ) in the transition band].

Design Example
For example purposes, Figures 3 and 4 show frequency responses of the third harmonic's transfer function T Q 3,0 z −1 in the case of K = 2, (a) (Task 1) for different values of maximally allowed gains in the stopbands l SB m ∈ {0.1, 0.01, 0.001} corresponding to attenuations of {20, 40, 60} dB, respectively, and (b) (Task 2) for different maximum allowed values of |TVG| in the harmonic frequency z m |TVG| max ∈ {24, 32, 48}. In the inset figures at the bottom, zoomed amplitude and |TVG| characteristics around the harmonic frequency are shown. The maximum gain in the passband and transition bands is chosen to be l PB m = 1.005. Notice that a higher value of N Q provides a smaller |TVG| and a wider passband. In addition to that, it is visible from the amplitude frequency responses that sidelobes are much larger for smaller values of |TVG|, which decreases robustness against interharmonics and noise. On the other hand, the passband width is increased for smaller values of |TVG|.
In the K = 1 case, the transfer function T 3,0 z −1 has a smaller inherent |TVG| than in the K = 2 case. Hence, in this case, the total |TVG| for T Q 3,0 z −1 is smaller; see Figure 5. Due to smaller values of N Q = 16, the passband is narrower. Notice that for a requested l SB m = 0.001, the problem is infeasible.

Simulation Results
Two following tests are given to illustrate the filter's (i.e., estimator's) dynamic features, while static characteristics are obviously clear from the frequency responses. From  Figures 3 and 4, it is visible that the filters designed through Task 2 have larger sidelobes close to the passband, while filters corresponding to Task 1 have nearly equiripple amplitude responses in the stopbands. High sidelobes can cause sensibility to interharmonic disturbances. The constant fundamental frequency over time is settled, and it equals 50 Hz. The sampling rate is f S = 1600 Hz. The input signal contains the third, fifth and seventh harmonics of 30%, 10%, and 10%, respectively. In the presented examples, only the estimation of the third harmonic is shown. Estimations of the other harmonics have the same characteristics.

Amplitude Modulated Signal
A test with the amplitude modulated signal is applied to illustrate the measurement bandwidth and latency depending on the settled constraints and subject of optimization. In the test signal, the amplitude of the third harmonic is defined as follows: where X 3 is the third harmonic amplitude, the amplitude modulation depth is 0.1 (10%), and f m is the modulating signal frequency. Values of modulated signals are chosen in line with the suggestion of the PMU IEC/IEEE Standard 60255-118-1:2018-2011 [22]. The maximum modulation frequency is f m = 2 Hz. The fundamental frequency is constant over time and equals 50 Hz. X 3 = 0.3 (p.u.). The influence of modulation of the sinusoidal amplitude of the magnitude amounting to 10% of the fundamental one and the modulation frequency equal to 2 Hz is shown in Figure 6. In Figure 6, it is notable that in all cases, estimates of the amplitude follow the actual latencies, which are in line with the |TVG|s shown in Figures 3 and 4.

Amplitude Step Signal
In this test, the estimates of the third harmonic's amplitude when the signal with a step change of the amplitude of the third harmonic by 10% is applied is presented. Other harmonics are constant over time. In addition, the fundamental frequency is constant over time, and it equals 50 Hz. Figure 7 illustrates the convergence of the amplitude estimate. It is possible to observe faster responses for a higher l SB 3 for Task 1, and a smaller |TVG| max for Task 2. In Task 1 case, for N Q = 16 and l SB 3 = 0.001, the latency is higher than for the basic T 3,0 z −1 due to the high requested attenuation, while for N Q = 32, the latency for these two cases is nearly the same. For Task 2, the latency is defined by |TVG| max , and it is equal for both N Q = 16 and N Q = 32. However, the attenuation in the stopbands is higher for N Q = 32; see Figures 3b and 4b. It should be mentioned that similar estimation shapes are obtained for the phase modulation and the phase step change.

Conclusions
This article continues the research presented in the previous articles, which dealt with harmonic analysis based on the recursive MR-based parallel structure with common feedback. The CLLS optimization technique is applied to optimize frequency responses, assuring adequate flatness of the frequency responses and reduced latency in the passband and high equiripple attenuation and acceptable levels of the sidelobes in the stopband. Thus, the obtained estimates can be simultaneously accurate, resistant to the presence of harmonics/interharmonics and noise, and fast enough to allow tracking of the nonstationary signal content. For different estimators' properties, a unique compensator can be designed for each harmonic. Furthermore, different harmonic phasors can be estimated simultaneously with different estimation performances.