Perspective Chapter: Cascaded-Resonator-Based Recursive Harmonic Analysis

Abstract

It is well known that recursive algorithms for harmonic analysis have better characteristics in terms of monitoring the change of the spectrum in comparison to methods based on the processing of blocks of consecutive samples, such as, for example, discrete Fourier transform (DFT). This property is particularly important when applying spectral estimation in real-time systems. One of the recursive algorithms is the resonator-based one. The approach of the parallel cascades of multiple resonators (MR) with the common feedback has been generalized as the cascaded-resonator (CR)-based structure for recursive harmonic analysis. The resulting filters of the CR structure can be finite impulse response (FIR) type or the infinite impulse response (IIR) ones as a computationally more efficient solution, optimizing the frequency responses of all harmonics simultaneously. In the case of the IIR filter, the unit characteristic polynomial present in the FIR filter is replaced with an optimized characteristic polynomial of the transfer function. Such a change does not lead to an increase in computing requirements and changes only the resonator gain values. By using a conveniently linearized iterative algorithm for stability control purpose, based on the Rouche’s theorem, the iterative linear-programming-based or the constrained linear least-squares (CLLS) optimization techniques can be used.

1. Introduction

In recent years, a lot of various algorithms for harmonic analysis have been proposed in the literature. Good surveys of some techniques are presented in Refs. [1, 2]. The discrete Fourier transform (DFT)-based method, as a mainstream approach, is widely used for harmonic analysis, thanks to its low computational burden, especially with the fast Fourier transform (FFT). However, errors arise when the power system is operating at off-nominal frequency, especially under dynamic conditions. Harmonic estimates under oscillating conditions were recently proposed in several studies. A huge volume of papers has been written on harmonics tracking in power systems. The focus of recent literature has been on preprocessing and postprocessing methods for fixed-sample-rate algorithms surrounding a core DFT (or similar) analysis with a fixed number of samples [3, 4].

Idea of considering a dynamic model to better estimate the fundamental and harmonic phasors has been emerging in Refs. [4, 5, 6, 7, 8, 9], and its importance has been pointed out in Ref. [10]. In Ref. [9], the discrete Taylor-Fourier transform (TFT) was proposed as an extension of the full DFT. The TFT by using a dynamic model of the signal extends and improves estimations obtained by DFT [9, 11]. This transformation corresponds to an FIR filter bank with a maximally flat frequency response. Each filter in the bank has maximum flat gain around the harmonic frequency and near-ideal attenuation around the other harmonics. This results in less distortion of the signal and less influence of disturbances present in the signal. In this way, the periodicity restriction assumed by the Fourier analysis is mitigated. As result, so obtained reconstruction is more accurate than the reconstruction obtained through DFT. When harmonics are narrow-band pass signals with spectral density confined into the flat-gain harmonic intervals, the coefficients of the TFT provide good estimates of the first derivatives of their complex envelopes. The digital TFT formulation in a matrix form that facilitates its implementation with the FFT to reduce the computational complexity of its straightforward implementation has been given in [11].

The multiple-resonator (MR)-based recursive estimators have been introduced in Ref. [12]. In Ref. [13], the MR-based observer structure is proposed for the implementation of TFT. Their good properties are provided by their parallel form, a recursive implementation, and good sensitivity properties assured by the infinite loop gain at the resonator frequencies [14]. Multiple zeros also provide reinforcing of the required attenuations and zero-gain flatness at the harmonic components with a high overall attenuation in the stopbands. For the known frequency of the periodic signal, the estimator based on resonators with common feedback enables the estimation of Fourier components even in cases when the sampling rate is not synchronized with the signal frequency. Also, this harmonic analyzer shows robustness in real conditions where there is noise and nonlinearity of the analog part of the equipment. MR-based harmonic analysis provides better performances of the spectral estimation than the single-resonator-based observer that corresponds to the classical DFT estimator.

This approach has been generalized as the cascaded-resonator (CR)-based structure for harmonic analysis. In Ref. [15], the cascaded-dispersed-resonator-based (CDR-based) structure for harmonic analysis is proposed. Although the design objectives in Refs. [15, 16] are different, the design technique is the same in both cases. In Ref. [16], the task is to replace multiple resonators with a cascade of single resonators. In this way, for the design purpose, it is possible to use the classic Lagrange interpolation technique instead of the more complex Hermitian interpolation. The condition that the poles are distributed in a narrow band around the resonant frequencies, as close as possible to each other, which is however limited by numerical accuracy. In Ref. [15], the task is to arrange the poles in the cascade in such a way as to enable optimal attenuation in the entire range around the harmonic frequencies. Practically, the only difference is in the arrangement of the resonator poles around the harmonic frequencies. The frequency deviation issue can be resolved by adaptive estimators based on the actual frequency feedback. This approach has drawbacks as a stability issue, due to an internal delay. Instead of that, usage of the external module for the fundamental frequency estimation is proposed in Ref. [17].

2. Cascaded-resonator-based structure harmonic analysis

Figure 1 shows the block diagram of the $K$-type CR-based harmonic analyzer. The structure includes $\left(K+1\right)\left(2M+2\right)$ resonators with poles $\left\{z_{m,k},m=-M,\dots ,0,\dots ,M+1,k=0,2,\dots K\right\}$, placed in the $2M+2$ cascades each of them having $K+1$ cascaded complex poles on the unit circle around related harmonic frequency [13, 15, 16]. Each resonator has its belonging complex gains $g_{m,k}$. A complete set of resonator cascades is connected in parallel in a common feedback loop. The gains of the transfer functions at the resonator frequencies are equal to unity due to the infinite loop gain at these frequencies. The number of cascades (for the coherent sampling) is $2M+2$ ($M$ is a number of harmonics) and depends on ${\omega }_1$, because the condition $M{\omega }_1<\pi$ has to be satisfied. The ${\omega }_1=2\pi f_1/f_S$, $f_1$ and $f_S$ are the nominal fundamental frequency and sampling rate. The overall system order is $\left(K+1\right)\left(2M+2\right)$.

Figure 1.

Block diagram of the $K$-type CR-based harmonic analyzer.

We have in every mth channel of the structure, as an internal transfer function

$$H_m\left(z\right)=\frac{V_m^F\left(z\right)}{E\left(z\right)}=z^{-1}\sum _{k=0}^K\frac{g_{m,k}^{\prime }}{{\prod }_{i=0}^{k-1}\left(1-z_{m,i}{z}^{-1}\right)},g_{m,k}^{\prime }=\prod _{i=k}^K{g}_{m,i}$$E1

where $m=-M,\dots ,0,\dots ,M+1$, $k=0,1,\dots ,K$, and ${\prod }_{i=0}^{k-1}\left(1-z_{m,i}{z}^{-1}\right)=1$ for $k=0$. $V_m^F\left(z\right)$ is the total feedback signal corresponding to $m$th channel, composed as the linear combination of the output of each resonator, i.e., each channel contributes to the filter output with $K+1$ complex weights $\left\{g_{m,k},k=0,1,\dots ,K\right\}$.

The closed transfer function for every channel $m$ has the form of [15, 16].

$$T_{m,0}\left(z\right)=\frac{V_{m,0}\left(z\right)}{V_1\left(z\right)}=g_{m,0}^{\prime }\frac{z^{-1}{P}_m\left(z\right)}{A\left(z\right)},P_m\left(z\right)=\prod _{\begin{array}{c}n=-M\\ n\ne m\end{array}}^{M+1}\prod _{i=0}^K\left(1-z_{n,i}{z}^{-1}\right),$$E2

$$T_{m,k}\left(z\right)=\frac{V_{m,k}\left(z\right)}{V_1\left(z\right)}=g_{m,k}^{\prime }\frac{z^{-1}{P}_m\left(z\right){\prod }_{i=0}^{k-1}\left(1-z_{m,i}{z}^{-1}\right)}{A\left(z\right)}$$E3

$$A\left(z\right)=\prod _{n=-M}^{M+1}\prod _{i=0}^K\left(1-z_{n,i}{z}^{-1}\right)+z^{-1}\sum _{n=-M}^{M+1}\left[P_n\left(z\right)\sum _{k=0}^K\left(g_{n,k}^{\prime }\prod _{i=0}^{k-1}\left(1-z_{n,i}{z}^{-1}\right)\right)\right]$$E4

It can be seen that all poles of the resonators are mapped to the zeros of the transfer function $T_{m,0}\left(z\right)$ due to the common feedback, with the exception of the poles belonging to the cascade of the harmonic $m$, which are automatically canceled by the poles that generated them. In differentiators transfer functions $\left\{T_{m,k}\left(z\right),k=1,\dots ,K\right\}$, $\left(k\ne 0\right)$, zeros $\left\{z_{m,i},i=0,\dots ,k-1\right\}$, originated from poles in mth channel, exist providing zero gain.

From Eq. (2), it is obvious that the filter corresponding to $m$th-channel provides the maximally flatness property in the stop band around the remaining harmonic frequencies. For small pole displacements, $∆f$ frequency response reshaping is negligible in comparison to the multiple-resonator case [16]. It is important to mention that the lower border of $∆f$ is limited by the computational accuracy. On the other hand, avoiding of the multiple poles allows design by the direct usage of the classical Lagrange interpolation formula rather than the Hermite one.

Although the characteristic polynomial of the transfer functions can be chosen in different ways, under some conditions it is possible to choose one so that the error is driven to zero in exactly $\left(K+1\right)\left(2M+2\right)$ samples. This is provided by what is called a dead-beat observer, for which the coefficients are calculated from the condition that the observer has deadbeat settling, i.e., it finds the unknown state within at most $\left(K+1\right)\left(2M+2\right)$ steps. That leads to FIR filters (with $A\left(z\right)=1$) in each channel. This way, although the structure is realized by resonators, which are IIR filters, the resulting filters in each channel are FIR type. Figure 2 shows frequency responses of $T_{1,0}\left(z\right)=g_{1,0}^{\prime }{z}^{-1}{P}_1\left(z\right)$ (corresponding to the fundamental component) of the dead-beat observer for the first up to the sixth order of resonator multiplicity ($K=0,1,\dots ,5$). It is observed that (quasi) MR structures with a higher order multiplicity of poles provide smaller sidelobes and thus ensure a lower sensitivity to noise and to harmonic and interharmonic disturbances. The negative effect is the increase in the order of the filter, which increases the group delay and response time of the filter, as well as the numerical complexity. Due to this feature, large values of the resonator multiplicity could be inconvenient in the control application. The case $K=0$ corresponds to a classic DFT estimator, while the cases $K>0$ correspond to the TFT.

Figure 2.

Frequency responses for $T_{1,0}\left(z\right)$ (the zeroth differentiator of the first harmonic) for $K=0,1,\dots ,5$ (the first- to sixth-order resonator structure).

Frequency responses of zeroth-, first-, and second-order differentiator discrete FIR filters corresponding to the transfer functions $T_{m,k}\left(z\right)$ related to the fundamental component ($m=1$, $k=0,1,2$), for $K=2$ and $K=3$ (the third- and forth-order resonator structure) are given in Figure 3.

Figure 3.

Frequency response of the zeroth-, first-, and second-order discrete FIR filters corresponding to $T_{m,k}\left(z\right)$, ($m=1$ and $k=0,1,2$), related to the fundamental component, for $f_s=800\mathit{Hz}$ and $f_1=50\mathit{Hz}$, for (a) $K=2$ estimator (the third-order resonator structure) and (b) $K=3$ estimator (the fourth-order resonator structure).

In order to obtain wider flatness intervals in the pass band, the feedback signals $V_m^F\left(z\right)$ could be used for harmonic estimation instead of $V_{m,0}\left(z\right)$ [13]. The global transfer function of the feedback loop is a sum of the transfer functions of all $K+1$ differentiators $T_m^F\left(z\right)=V_m^F\left(z\right)/V\left(z\right)={\sum }_{k=0}^K{T}_{m,k}\left(z\right)$. The frequency responses of the estimation of the fundamental component obtained by $V_{m,0}\left(z\right)$ and estimation obtained by $V_m^F\left(z\right)$ are given together in Figure 4a. The good properties of the filters corresponding to the transfer functions $T_m^F\left(z\right)$ are related to the phase responses. Frequency responses have a zero phase response in the frequency bands around the harmonic frequencies, which means that in those frequencies the group delay is equal to zero. Bad properties are high resonant gains at the edges of bandwidths and high sidelobes. The zero flat gains in the stop band are preserved, although their intervals are narrowed. It should be mentioned that the peaks of the interharmonic gains and the side lobes increase by the multiplicity of the resonators (Figure 4b).

Figure 4.

Frequency responses of (a) $T_{1,0}\left(z\right)$ and $T_1^F\left(z\right)$ for $K=2$ estimator (the third-order cascade) and (b) $T_1^F\left(z\right)$ for $K=1$, $K=2$ and $K=3$, related to the fundamental component, for $f_s=800\mathit{Hz}$ and $f_1=50\mathit{Hz}$.

2.1 Optimization problem statement

In order to adapt the achieved digital differentiators to their ideal frequency responses around the harmonic frequencies, it is possible to modify the filters transfer functions. An optimization technique is utilized to reshape frequency responses of the filters transfer functions, avoiding resonant frequency peaks and reducing a group delay simultaneously, that can be rather important in control applications. The optimization task can also be different, e.g., maximization of the selectivity.

The transfer function of the extended structure, including the compensation FIR filter $B\left(z\right),$ for the mth component channel is as follows:

$$T_{m,0}^{\mathit{AB}}\left(z\right)=\frac{V_{m,0}\left(z\right)}{V_1\left(z\right)}=B\left(z\right)T_{m,0}\left(z\right)=\left[g_{m,0}^{\prime }{z}^{-1}{P}_m\left(z\right)\right]\frac{B\left(z\right)}{A\left(z\right)}=\left[g_{m,0}^{\prime }{z}^{-1}{P}_m\left(z\right)\right]\frac{{\mathit{q}}_{\mathit{B}}{\mathit{x}}_{\mathit{B}}}{1+{\mathit{q}}_{\mathit{A}}{\mathit{x}}_{\mathit{A}}}$$E5

where

$$\begin{array}{c}B\left(z\right)=b_0+b_1{z}^{-1}+\cdots +b_{N_B-1}{z}^{-\left(N_B-1\right)}+b_{N_B}{z}^{-N_B},\\ A\left(z\right)=1+a_1{z}^{-1}+\cdots +a_{N_A-1}{z}^{-\left(N_A-1\right)}+a_{N_A}{z}^{-N_A},\end{array}$$

$${\mathit{x}}_{\mathit{B}}={\left[b_0{b}_1\cdots \begin{array}{ccc}& b_{N_B-1}& b_{N_B}\end{array}\right]}^{\mathrm{T}}\mathrm{,}{\mathit{x}}_{\mathit{A}}={\left[\begin{array}{cc}{a}_1& a_2\end{array}\cdots \begin{array}{cc}{a}_{N_A-1}& a_{N_A}\end{array}\right]}^{\mathrm{T}}\mathrm{,}\mathit{x}\mathit{=}{\left[\begin{array}{cc}{{\mathit{x}}_{\mathit{B}}}^{\mathrm{T}}& {{\mathit{x}}_{\mathit{A}}}^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}}\mathrm{,}$$

$${\mathit{q}}_{\mathit{B}}\mathit{=}\left[\begin{array}{ccc}1& {\mathit{z}}^{-1}& {\mathit{z}}^{-2}\end{array}\begin{array}{ccc}\mathit{\cdots }& {\mathit{z}}^{-\left({\mathit{N}}_{\mathit{B}}-1\right)}& {\mathit{z}}^{-{\mathit{N}}_{\mathit{B}}}\end{array}\right]\mathrm{,}{\mathit{q}}_{\mathit{A}}=\left[\begin{array}{cc}{\mathit{z}}^{-1}& {\mathit{z}}^{-2}\end{array}\begin{array}{ccc}\mathit{\cdots }& {\mathit{z}}^{-\left({\mathit{N}}_A-1\right)}& {\mathit{z}}^{-{\mathit{N}}_{\mathit{A}}}\end{array}\right]\mathit{.}$$

The polynomial $A\left(z\right)$ does not cause any additional computation and only the polynomial $B\left(z\right)$ represents an additional numerical burden. Even more, in some cases, it is possible to choose $B\left(z\right)=b_0$ which causes only one additional multiplication. Nevertheless, for the purposes of design, we will consider the IIR filter $B\left(z\right)/A\left(z\right)$ as a common compensation for the total set of FIR filters $P_m\left(z\right)$, $m=-M,\dots ,0,\dots ,M+1$.

With a given weighting function $W\left(\omega \right)$, the weighted Chebyshev error between the desired and actual frequency responses is defined as follows:

$$J^1\left(\mathit{x}\right)=\underset{\mathrm{\omega }\in \mathrm{\Omega }}{max}W\left(\omega \right)\left|T^{\mathit{AB}}\left(e^{\mathit{j\omega }}\right)-H^d\left(e^{\mathit{j\omega }}\right)\right|$$E6

where $H^d\left(e^{\mathit{j\omega }}\right)$ is the desired frequency response in angular frequency $\omega$ specified in the frequency region $\Omega$ (or a union of several compact frequency bands) of the interests and 0 ≤ ω ≤ π.

The sum of squares of absolute values of errors in N_F angular frequencies as follows:

$$J^2\left(\mathit{x}\right)=\sum _{i=1}^{N_F}W\left(z_i\right)\left[T^{\mathit{AB}}\left(z_i\right)-H^d\left(z_i\right)\right]$$E7

where $z_i=e^{j{\omega }_i}$.

In order to minimize the error $J^1\left(\mathit{x}\right)$ defined in Eq. (6), a new variable $\delta$ can be introduced and the problem reformulated as follows:

$$\begin{array}{c}minimize\delta \\ \text{subject to}\left|E\left({\omega }_i\right)\right|\le \delta ,{\omega }_i\in \Omega ,i=1,2,\cdots ,N_F\end{array}$$E8

where $E\left(z_i\right)=W\left(z_i\right)\left[T^{\mathit{AB}}\left(z_i\right)-H^d\left(z_i\right)\right]$, $z_i=exp\left(j{\omega }_i\right)$, for the total number $N_F$ of points defined in $\Omega$.

Further, it is:

$$\frac{W\left(z_i\right)}{\left|A\left(z_i\right)\right|}\left|T\left(z_i\right)B\left(z_i\right)-A\left(z_i\right)H^d\left(z_i\right)\right|\le \delta$$E9

In Ref. [18], a suitable method has been described to linearize the error function $J^1\left(\mathit{x}\right)$ such that the design problem can be solved by the linear programming (LP) method. However, this method neglects the denominator part $\left|A\left(z_i\right)\right|$. In Ref. [19], the performance of the LP method was improved by eliminating the above drawback by using the following iterative constraints scheme:

$$\frac{W\left(z_i\right)}{\left|A^{\left(k-1\right)}\left(z_i\right)\right|}\left|T\left(z_i\right)B\left(z_i\right)-A^{\left(k\right)}\left(z_i\right)H^d\left(z_i\right)\right|\le \delta$$E10

For the sake of notational simplicity, we denote

$$\frac{W\left(z_i\right)}{\left|A^{\left(k-1\right)}\left(z_i\right)\right|}\left|\left[{\left.T\left(z_i\right){\mathit{q}}_{\mathit{B}}\right|}_{z=z_i},{\left.-H^d\left(z_i\right){\mathit{q}}_{\mathit{A}}\right|}_{z=z_i}\right]\mathit{x}\mathbf{-}{H}^d\left(z_i\right)\right|\le \delta$$E11

The vector of unknown coefficients $\mathit{x}$ is expanded with an additional variable δ, so that the expanded vector of unknowns is obtained:

$${\mathit{x}}^{\mathit{\delta }}={\left[\begin{array}{cc}\mathit{x}& \delta \end{array}\right]}^T.$$E12

The constraints defined by inequality (11) refer to the frequency ranges in which the error optimization is performed. In addition, sometimes it is necessary to keep the error within predefined limits, such as for example the gains in the stopbands and/or the transition bands:

$$\frac{1}{\left|A^{\left(k-1\right)}\left(z_i\right)\right|}\left|\left[T\left(z_i\right){\left.{\mathit{q}}_{\mathit{B}}\right|}_{z=z_i},{\left.-H^d\left(z_i\right){\mathit{q}}_{\mathit{A}}\right|}_{z=z_i}\right]\mathit{x}-H^d\left(z_i\right)\right|\le l_i$$E13

$z_i=exp\left(j{\omega }_i\right)$, $i=1,2,\cdots ,N_G$

$l_i$, $i=1,2,\cdots ,N_G$, are fixed borders of the absolute values of the set of complex error $E\left(z\right)$ along assemblies $z_i=exp\left(j{\omega }_i\right)$, $i=1,2,\cdots ,N_G$.

If one wants to ensure unity gain in harmonic frequencies, the following condition must be met:

$${\left.\frac{B\left(z\right)}{A\left(z\right)}\right|}_{z=z_m}=1,\mathrm{i}.\mathrm{e}.B\left(z_m\right)=A\left(z_m\right).$$E14

where $m=-M,\dots ,0,\dots ,M+1.$

In a matrix form, it can be written as follows:

$${\mathit{q}}_{\mathit{m}}\mathit{x}=1.$$E15

where $\mathit{q}=\left[\begin{array}{cc}{\mathit{q}}_{\mathit{B}}& \mathbf{-}{\mathit{q}}_{\mathit{A}}\end{array}\right]$, ${\mathit{q}}_{\mathit{m}}={\left.\mathit{q}\right|}_{z=z_m}$.

Complex equality constraints (15) can be written as follows

$$\left[\begin{array}{c}Re\left\{{\mathit{q}}_{\mathit{m}}\right\}\\ Im\left\{{\mathit{q}}_{\mathit{m}}\right\}\end{array}\right]\mathit{x}=\left[\begin{array}{c}1\\ 0\end{array}\right],m=-M,\dots ,0,\dots ,M+1.$$E16

2.2 Linearization of constraints

The inequalities (11) and (13) are nonlinear. The convex semi-infinite programming can be applied [20], thanks to the quadratic property of the functions. Furthermore, a convenient approximation of these inequalities by the system of the linear ones [21, 22, 23, 24, 25] allows us to solve this constrained optimization problem through the LP or the constrained linear least-squares (CLLS) optimization technique.

It is valid:

$$\left|E\left(z_i\right)\right|=\left|E\left(z_i\right)\right|\left({cos}^2{\alpha }_i+{sin}^2{\alpha }_i\right)=Re\left\{E\left(z_i\right)\right\}cos{\alpha }_i+Im\left\{E\left(z_i\right)\right\}sin{\alpha }_i$$E17

where ${\alpha }_i=arg\left\{E\left(z_i\right)\right\}$.

Since ${\alpha }_i$ is not known a priory, the nonlinear constraints in Eq. (8) can be approximated by the system of linear constraints:

$$Re\left\{E\left(z_i\right)\right\}cos{\mathit{\alpha }}_{\mathit{i},\mathit{j}}+Im\left\{\mathit{E}\left({\mathit{z}}_{\mathit{i}}\right)\right\}sin{\mathit{\alpha }}_{\mathit{i},\mathit{j}}\le \mathit{\delta }$$E18

where $i=1,2,\cdots ,N_F$. If we choose $L$ equidistantly distributed angles then it is ${\alpha }_{i,j}={\alpha }_{i,0}+\left(j-1\right)2\pi /L$, where $j=1,2,\cdots ,L$. Figure 5 shows that approximations by square and octagon ($L=4$ and $L=8$, respectively) allow only rough approximations. A higher accuracy is obtained by increasing $L$. Herein, $L=32$ is used.

Figure 5.

Approximation of a cycle with a square and an octagon.

Let us define as following for frequency point $i$:

$$W^A\left(z_i\right)=\frac{W\left(z_i\right)}{\left|A^{\left(k-1\right)}\left(z_i\right)\right|},{\mathit{g}}_i=\left[{\left.T\left(z_i\right){\mathit{q}}_{\mathit{B}}\right|}_{z=z_i},{\left.-H^d\left(z_i\right){\mathit{q}}_{\mathit{A}}\right|}_{z=z_i}\right].$$E19

Hence, Eq. (13) can be linearized and written in a matrix form

$$W^A\left(z_i\right){\mathit{A}}_i^{\prime }\mathit{x}\le W^A\left(z_i\right){\mathit{b}}_i^{\prime }+l_i{\mathbf{1}}_{L\times 1},i=1,2,\cdots ,N_G$$E20

where matrix ${\mathit{A}}_i^{\prime }$ and vector ${\mathit{b}}_i^{\prime }$ are given by

$${\mathit{A}}_i^{\prime }=\left[\begin{array}{c}Re\left\{{\mathit{g}}_i\right\}cos{\alpha }_{i1}+Im\left\{{\mathit{g}}_i\right\}sin{\alpha }_{i1}\\ Re\left\{{\mathit{g}}_i\right\}cos{\alpha }_{i2}+Im\left\{{\mathit{g}}_i\right\}sin{\alpha }_{i2}\\ ⋮\\ Re\left\{{\mathit{g}}_i\right\}cos{\alpha }_{\mathit{iL}}+Im\left\{{\mathit{g}}_i\right\}sin{\alpha }_{\mathit{iL}}\end{array}\right],{\mathit{b}}_i^{\prime }=\left[\begin{array}{c}Re\left\{H^d\left(z_i\right)\right\}cos{\alpha }_{i1}+Im\left\{H^d\left(z_i\right)\right\}sin{\alpha }_{i1}\\ Re\left\{H^d\left(z_i\right)\right\}cos{\alpha }_{i2}+Im\left\{H^d\left(z_i\right)\right\}sin{\alpha }_{i2}\\ ⋮\\ Re\left\{H^d\left(z_i\right)\right\}cos{\alpha }_{\mathit{iL}}+Im\left\{H^d\left(z_i\right)\right\}sin{\alpha }_{\mathit{iL}}\end{array}\right].$$

and $l_i$ is a constraint limit of the error in the point $z_i$.

Using matrix notation, and collecting inequality linearization systems in all settled frequency points, (20) becomes the following linear form:

$${\mathit{A}}^{\prime }\mathit{x}\le {\mathit{b}}^{\prime }\text{or}\left[\begin{array}{cc}{\mathit{A}}^{\prime }& {\mathbf{0}}_{\left(N_{G}L\right)\times 1}\end{array}\right]{\mathit{x}}^{\mathit{\delta }}\le {\mathit{b}}^{\prime }$$E21

where matrix ${\mathit{A}}^{\prime }$ and vector ${\mathit{b}}^{\prime }$ are given by

$${\mathit{A}}^{\prime }=\left[\begin{array}{c}{W}^A\left(z_1\right){\mathit{A}}_1^{\prime }\\ W^A\left(z_2\right){\mathit{A}}_2^{\prime }\\ ⋮\\ W^A\left(z_{N_G}\right){\mathit{A}}_{N_G}^{\prime }\end{array}\right],{\mathit{b}}^{\prime }=\left[\begin{array}{c}{W}^A\left(z_1\right){\mathit{b}}_1^{\prime }+l_1{\mathbf{1}}_{L\times 1}\\ W^A\left(z_2\right){\mathit{b}}_2^{\prime }+l_2{\mathbf{1}}_{L\times 1}\\ \mathbf{⋮}\\ W^A\left(z_{N_G}\right){\mathit{b}}_{N_G}^{\prime }+l_{N_G}{\mathbf{1}}_{L\times 1}\end{array}\right].$$

In case of the (11), we have:

$$\left[\begin{array}{cc}{W}^A\left(z_i\right){\mathit{A}}_i^{\prime }& -{\mathbf{1}}_{L\times 1}\end{array}\right]{\mathit{x}}^{\mathit{\delta }}\le W^A\left(z_i\right){\mathit{b}}_i^{\prime },i=1,2,\cdots ,N_F$$E22

or in a matrix notation

$$\left[\begin{array}{cc}{\mathit{A}}^{\prime }& \mathbf{-}{\mathbf{1}}_{\left(N_{F}L\right)\times 1}\end{array}\right]{\mathit{x}}^{\mathit{\delta }}\le {\mathit{b}}^{\prime }$$E23

where

$${\mathit{A}}^{\prime }=\left[\begin{array}{c}{W}^A\left(z_1\right){\mathit{A}}_1^{\prime }\\ W^A\left(z_2\right){\mathit{A}}_2^{\prime }\\ ⋮\\ W^A\left(z_{N_F}\right){\mathit{A}}_{N_F}^{\prime }\end{array}\right],{\mathit{b}}^{\prime }=\left[\begin{array}{c}{W}^A\left(z_1\right){\mathit{b}}_1^{\prime }\\ W^A\left(z_2\right){\mathit{b}}_2^{\prime }\\ \mathbf{⋮}\\ W^A\left(z_{N_F}\right){\mathit{b}}_{N_F}^{\prime }\end{array}\right].$$

2.3 Design (optimization) approach 1: CLLS minimization

An objective is to find a minimum of the sum of squares of absolute values of $\mathit{hx}\mathit{-}d$ in the assembly of the $N_F$ selected frequencies subject to the vector $\mathit{x}$

$$\underset{x}{min}\sum _{i=1}^{N_F}{\left|{\mathit{h}}_{\mathit{i}}\mathit{x}\mathbf{-}{d}_i\right|}^2$$E24

where ${\mathit{h}}_{\mathit{i}}\mathbf{=}W\left(z_i\right){\mathit{g}}_{\mathit{i}}/\left|A^{\left(k-1\right)}\left(z_i\right)\right|$ and $d_i=W\left(z_i\right)H^d\left(z_i\right)/\left|A^{\left(k-1\right)}\left(z_i\right)\right|$.

If we apply the following equality

$${\left|{\mathit{h}}_{\mathit{i}}\mathit{x}\mathbf{-}{d}_i\right|}^2={Re}^2\left\{{\mathit{h}}_{\mathit{i}}\mathit{x}\mathbf{-}{d}_i\right\}+{Im}^2\left\{{\mathit{h}}_{\mathit{i}}\mathit{x}\mathbf{-}{d}_i\right\}={‖{\mathit{C}}_{\mathit{i}}\mathit{x}\mathbf{-}{\mathit{d}}_i‖}_{\mathbf{2}}^{\mathbf{2}}$$E25

where ${\mathit{C}}_{\mathit{i}}=\left[\begin{array}{c}Re\left\{{\mathit{h}}_{\mathit{i}}\right\}\\ Im\left\{{\mathit{h}}_{\mathit{i}}\right\}\end{array}\right]$, ${\mathit{d}}_{\mathit{i}}=\left[\begin{array}{c}Re\left\{d_i\right\}\\ Im\left\{d_i\right\}\end{array}\right]$, (24) can be written in a matrix form:

$$\underset{x}{min}{‖\mathit{Cx}-\mathit{d}‖}_{\mathbf{2}}^{\mathbf{2}}.$$E26

where $\mathit{C}$ and $\mathit{d}$ include ${\mathit{C}}_{\mathit{i}}$ and ${\mathit{d}}_i$, respectively, $i=1,2,\cdots ,N_F$.

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 (16), (21), and (26), the CLLS problem is formalized as follows:

$$\begin{array}{l}\underset{\mathit{x}}{min}\frac{1}{2}{‖\mathit{Cx}-\mathit{d}‖}_{\mathbf{2}}^{\mathbf{2}}\text{subject to}\mathit{Ax}\le \mathit{b}\left(\text{and}\left[\begin{array}{c}Re\left\{{\mathit{q}}_{\mathit{m}}\right\}\\ Im\left\{{\mathit{q}}_{\mathit{m}}\right\}\end{array}\right]\mathit{x}=\left[\begin{array}{c}1\\ 0\end{array}\right]\right),\\ \left(m=-M,\dots ,0,\dots ,M+1\right)\end{array}$$E27

where $\mathit{A}$ and $\mathit{b}$ include ${\mathit{A}}^{\prime }$ and ${\mathit{b}}^{\prime }$, respectively, defined in (21), for all frequency points in which the constraints are defined.

2.4 Design (optimization) approach 2: minimax optimization

The LP optimization problem can be formalized in the following way:

$$\begin{array}{l}minimize\mathbf{c}{\mathit{x}}^{\mathit{\delta }}\mathrm{s}\text{ubject to}\mathit{A}{\mathit{x}}^{\mathit{\delta }}\le \mathit{b}\left(\text{and}\left[\begin{array}{cc}Re\left\{{\mathit{q}}_{\mathit{m}}\right\}& 0\\ Im\left\{{\mathit{q}}_{\mathit{m}}\right\}& 0\end{array}\right]{\mathit{x}}^{\mathit{\delta }}\mathbf{=}\left[\begin{array}{c}1\\ 0\end{array}\right]\right),\\ \left(m=-M,\dots ,0,\dots ,M+1\right)\end{array}$$E28

where $\mathbf{c}=\left[\begin{array}{cc}{\mathbf{0}}_{\mathbf{1}\times \left({\mathit{N}}_{\mathit{A}}+{\mathit{N}}_{\mathit{B}}+\mathbf{1}\right)}& 1\end{array}\right]$, and $\mathit{A}$ and $\mathit{b}$ include ${\mathit{A}}^{\prime }$ and ${\mathit{b}}^{\prime }$, respectively, for all frequency points in which the constraints or objective functions are defined in (21) and (23), respectively.

3. IIR cascaded-resonator-based harmonic analysis

In accordance with the prevailing trends in works dealing with this issue, in the initial works [13, 16, 22, 23, 25, 26] the resulting filters of CR structures were of the FIR type. Later, in [27], IIR filters were used, which represent a computationally more efficient solution [28, 29]. The unit characteristic polynomial of the transfer function is replaced by the optimized one. Such a change does not lead to an increase in the volume of numerical calculations and only requires a change in the gain values associated with the resonators. Since the optimization of frequency characteristics for all harmonics is carried out at the same time, it is possible to obtain frequency responses of the same shape. By using a linearized iterative scheme [30] based on Rouche’s theorem with the aim of stability control, it is possible to use iterative optimization techniques based on LP or CLLS.

3.1 Problem statement

The task of optimization is to design a filter $B\left(z\right)/A\left(z\right)$ where the order of the characteristic polynomial $A\left(z\right)$ is $N_A=\left(K+1\right)\left(2M+2\right)$ and the polynomial $B\left(z\right)$ is of order $N_B$. We seek to find a causal stable rational function $T_{m,0}^{\mathit{AB}}\left(z\right)=\left[g_{m,0}^{\prime }{z}^{-1}{P}_m\left(z\right)\right]B\left(z\right)/A\left(z\right)$ for $m=-M,\dots ,0,\dots ,M+1$ that best approximates $H_m^d\left(e^{\mathit{j\omega }}\right)$.

In order to make the notation as simple and short as possible, let us form a virtual transfer function so that in each bandwidth centered in $m{f}_1$with width of $f_1$, i.e., for $f\in \left(m{f}_1-f_1/2\right.,\left.m{f}_1+f_1/2\right]$, it corresponds to the transfer function belonging to the harmonic $m$. It follows:

$$T^A\left(z\right)=\frac{T\left(z\right)}{A\left(z\right)},T\left(z\right)=g_{m,0}^{\prime }{z}^{-1}{P}_m\left(z\right)$$E29

for $f\in \left(m{f}_1-f_1/2\right.,\left.m{f}_1+f_1/2\right],m=-M,\dots ,0,\dots ,M+1$.

In addition, we define a unique transfer function

$$T^{\mathit{AB}}\left(z\right)=B\left(z\right)T^A\left(z\right)=B\left(z\right)\frac{T\left(z\right)}{A\left(z\right)}$$E30

Similarly, a virtual unique desired transfer function in an angular frequency $\omega$ has the following form [27]:

$$H^d\left(e^{\mathit{j\omega }}\right)=\left\{\begin{array}{c}{e}^{-j2\mathit{\pi \tau }\left(f-m{f}_1\right)/f_S},\text{for}f\in \left[m{f}_1-f_{\mathit{PB}}\right.,\left.m{f}_1+f_{\mathit{PB}}\right]\\ 0,\text{for}f\in \left[m{f}_1-f_1/2\right.,\left.m{f}_1-f_{\mathit{SB}}\right]\bigcup \left[m{f}_1+f_{\mathit{SB}}\right.,\left.m{f}_1+f_1/2\right]\end{array}\right.$$E31

where pass and stop bands are defined by $f_{\mathit{PB}}$ and $f_{\mathit{SB}}$, respectively. A desired group delay in the passband is denoted as $\tau$.

3.2 Stability constraint

To obtain a stable IIR filter $T\left(z\right)$, stability constraint must be imposed on the coefficient vector ${\mathit{x}}_{\mathit{A}}$. In [30], the more convenient stability condition which is based on Rouche’s theorem was proposed.

Rouche’s Theorem. If$f\left(z\right)$and$g\left(z\right)$are analytic inside and on a closed contour$C$, and$\left|g\left(z\right)\right|<\left|f\left(z\right)\right|$on$C$, then$f\left(z\right)$and$g\left(z\right)+f\left(z\right)$have the same number of zeros inside$C$.

Let

$$f\left(z\right)=z^{N_A}A\left(z\right)=z^{N_A}+a_1{z}^{N_A-1}+\cdots +a_{N_A-1}z+a_{N_A},$$E32

$$g\left(z\right)=z^{N_A}∆\left(z\right)={\delta }_0{z}^{N_A}+{\delta }_1{z}^{N_A-1}+\cdots +{\delta }_{N_A-1}z+{\delta }_{N_A},$$E33

where $∆\left(z\right)$ is the update of the characteristic polynomial $A\left(z\right)$ of the transfer function at each iteration step. Since the functions $f\left(z\right)$ and $g\left(z\right)$ are analytic, except at z = ∞, and have the same zeros as $A\left(z\right)$ and $∆\left(z\right)$, according to Rouche’s theorem, if the polynomial $A^{\left(k-1\right)}\left(z\right)$ in the iteration step, $k-1$ has all its zeros inside a circle of radius $\rho \left(0<\rho <1\right)$ with the center at the origin of the complex plane, then also the polynomial in the iteration step $k$ given by [30]

$$A^{\left(k\right)}\left(z\right)=A^{\left(k-1\right)}\left(z\right)+\alpha {∆}^{\left(k\right)}\left(z\right),0<\alpha <1$$E34

will retain the zeros within this circle provided that in step $k$ the following condition satisfied

$$\left|{∆}^{\left(k\right)}\left(z\right)\right|\le \left|A^{\left(k-1\right)}\left(z\right)\right|,\left|z\right|=\rho .$$E35

If (34) is included in (35), we get

$$\left|{\overline{A}}^{\left(k\right)}\left(z\right)-{\overline{A}}^{\left(k-1\right)}\left(z\right)\right|\le \alpha \left|A^{\left(k-1\right)}\left(z\right)\right|$$E36

where ${\overline{A}}^{\left(k\right)}\left(z\right)=A^{\left(k\right)}\left(z\right)-1$, ${\overline{A}}^{\left(k-1\right)}\left(z\right)=A^{\left(k-1\right)}\left(z\right)-1$, or in a matrix notation:

$$\left|\left[{\mathbf{0}}_{\mathbf{1}\times \left({\mathit{N}}_{\mathit{B}}+\mathbf{1}\right)},{\left.{\mathit{q}}_{\mathit{A}}\right|}_{z=z_i}\right]\mathit{x}-{\overline{A}}^{\left(k-1\right)}\left(z_i\right)\right|\le \alpha \left|A^{\left(k-1\right)}\left(z_i\right)\right|$$E37

As for the initial value of the vector $\mathit{x}$, it is simplest to take ${\mathit{x}}^{\left(0\right)}=\mathbf{0}$, when all the roots of the polynomial $A^{\left(0\right)}\left(z\right)$ lie within the circle of radius $\rho$ (0< $\rho$ <1).

If constraint (37) is applied to a sufficiently dense set of points lying on a circle of radius $\rho \left(\left|z\right|=\rho \right)$, of total length $N_S$, we get

$${\mathit{A}}^{\mathrm{S}}\mathit{x}\le {\mathit{b}}^{\mathrm{S}}$$E38

where matrix ${\mathit{A}}_i^{\mathrm{S}}$ and vector ${\mathit{b}}_i^{\mathrm{S}}$ are given by

$${\mathit{A}}_i^S=\left[\begin{array}{c}{\mathbf{0}}_{\mathbf{1}\mathbf{\times }\left({\mathit{N}}_{\mathit{B}}\mathbf{+}\mathbf{1}\right)},Re\left\{{\left.{\mathit{q}}_{\mathit{A}}\right|}_{z=z_i}\right\}cos{\alpha }_{i1}+Im\left\{{\left.{\mathit{q}}_{\mathit{A}}\right|}_{z=z_i}\right\}sin{\alpha }_{i1}\\ {\mathbf{0}}_{\mathbf{1}\mathbf{\times }\left({\mathit{N}}_{\mathit{B}}\mathbf{+}\mathbf{1}\right)},Re\left\{{\left.{\mathit{q}}_{\mathit{A}}\right|}_{z=z_i}\right\}cos{\alpha }_{i2}+Im\left\{{\left.{\mathit{q}}_{\mathit{A}}\right|}_{z=z_i}\right\}sin{\alpha }_{i2}\\ ⋮\\ {\mathbf{0}}_{\mathbf{1}\mathbf{\times }\left({\mathit{N}}_{\mathit{B}}\mathbf{+}\mathbf{1}\right)},Re\left\{{\left.{\mathit{q}}_{\mathit{A}}\right|}_{z=z_i}\right\}cos{\alpha }_{\mathit{iL}}+Im\left\{{\left.{\mathit{q}}_{\mathit{A}}\right|}_{z=z_i}\right\}sin{\alpha }_{\mathit{iL}}\end{array}\right],$$

$${\mathit{b}}_i^S=\left[\begin{array}{c}Re\left\{{\overline{A}}^{\left(k-1\right)}\left(z_i\right)\right\}cos{\alpha }_{i1}+Im\left\{{\overline{A}}^{\left(k-1\right)}\left(z_i\right)\right\}sin{\alpha }_{i1}+\alpha \left|A^{\left(k-1\right)}\left(z_i\right)\right|\\ Re\left\{{\overline{A}}^{\left(k-1\right)}\left(z_i\right)\right\}cos{\alpha }_{i2}+Im\left\{{\overline{A}}^{\left(k-1\right)}\left(z_i\right)\right\}sin{\alpha }_{i2}+\alpha \left|A^{\left(k-1\right)}\left(z_i\right)\right|\\ ⋮\\ Re\left\{{\overline{A}}^{\left(k-1\right)}\left(z_i\right)\right\}cos{\alpha }_{\mathit{iL}}+Im\left\{{\overline{A}}^{\left(k-1\right)}\left(z_i\right)\right\}sin{\alpha }_{\mathit{iL}}+\alpha \left|A^{\left(k-1\right)}\left(z_i\right)\right|\end{array}\right].$$

Thus, the set of constraints (38) is added to the set of the above constraint conditions. In this way, the iterative methods mentioned above solve the LP or CLLS problem by taking into account constraints (38) in each iteration step. A fixed step size $\alpha$ can be used, while a gradual decrease (e.g., exponential) can help the convergence of the solution.

3.3 Resonators’ gains calculation

After the polynomial $A\left(z\right)$ having been determined, the direct usage of the Lagrange interpolation formula provides the closed-form formulas [26]. It should be taken into account that these formulas are valid only in the case of single resonators. If a quasi-MR-based analyzer is designed, the resonator poles connected to the same harmonic should be arranged close enough to each other with a minimum distance that is limited by numerical precision. Its lower border depends on the resonator multiplicity and the sampling rate. A chosen displacement of $0.1\mathrm{Hz}$ allows a fair approximation for sampling frequencies up to 6.4 kHz ($M=63$ for $f_1=50\mathrm{Hz}$) and $K=5$ [16].

A generalized closed-form formula for gains calculation for any $K$ for previously chosen polynomial $A\left(z\right)$ is given as follows:

$$g_{m,k}^{\prime }={\left.\frac{A\left(z\right)-z^{-1}{P}_m\left(z\right){\sum }_{j=0}^{k-1}\left[g_{m,j}^{\prime }{\prod }_{i=0}^{j-1}\left(1-z_{m,i}{z}^{-1}\right)\right]}{z^{-1}{P}_m\left(z\right)\prod _{i=0}^{k-1}\left(1-z_{m,i}{z}^{-1}\right)}\right|}_{z=z_{m,k}}$$E39

As a final result, the designed resonator gains are as follows:

$$g_{m,k}=g_{m,k}^{\prime }/g_{m,k+1}^{\prime },\left(g_{m,K+1}^{\prime }=1\right),m=-M,\dots ,0,\dots ,M+1,k=0,1,\dots ,K.$$

It should be mentioned that polynomial $B\left(z\right)$ can be conveniently implemented by adding its roots as poles in additional parallel channels to the basic structure (see Figure 1). In this case, the existing formulas for gains calculation are valid only for the identical structure of the extension cascades (they have to consist $K+1$ resonators). Otherwise, the formulas are not valid and need a completely new derivation.

3.4 Design example

In the next section, three demonstration examples, with frequency responses and pole-zeros maps, of the designed $K=2$ type CR-based harmonic analyzer are shown, for $f_S=800\mathrm{Hz}$ and $f_1=50\mathrm{Hz}$. For a clear readability, lower values of $f_S=800\mathrm{Hz}$ and $M=7$ are selected. The following parameters are prescribed in all three examples: $f_{\mathit{SB}}=17\mathrm{Hz}$, $l_{\mathit{TB}}=1.005,W\left(z\right)=1,{\alpha }_0=0.5,\rho =0.95,N_A=\left(K+1\right)\left(2M+2\right)$. ${\mathit{x}}^{\left(0\right)}=\mathbf{0}.$ Other parameters were variated, depending on the chosen optimization criteria. It should mention that a large variety of optimization scenarios is possible, allowing design spectrum analyzers for a wide scope of different applications.

3.4.1 Example 1: flat-top passbands

The filters with a wider flatness in the pass band allow better signal tracking in the dynamic conditions. Since it is difficult task to provide the tracking of the parameters changes together with good attenuation in the stopband, a relatively high order of $N_B=\left(K+1\right)\left(2M+2\right)$ is settled. The desired group delay (in samples) in the passband is $\tau$=0.9$\left(K+1\right)\left(2M+2\right).$$f_{\mathit{PB}}=1.7\mathrm{Hz}.l^{\mathit{PB}}=0.01$. Obtained frequency responses (Figure 6) show that the passband flatness is not derogated, while the selectivity and attenuation in the stopbands are increased thanks to the zeros of the polynomial $B\left(z\right)$ which are located between the existing multiple zeros of the resonator structure that had been obtained through the common feedback. A cost is an increased total group delay which causes higher latency.

Figure 6.

(a) Frequency responses of $g_{3,0}^{\prime }{z}^{-1}{P}_3\left(z\right)$, $B\left(z\right)/A\left(z\right)$ and $T_{3,0}^{\mathit{AB}}\left(z\right)$and (b) Pole-zero map of $T_{3,0}\left(z\right)$ and $B\left(z\right)$ (for the third harmonic), for $N_A=\left(K+1\right)\left(2M+2\right)$, $N_B=\left(K+1\right)\left(2M+2\right)$, and $\rho =0.95,f_{\mathit{PB}}=1.7\mathrm{Hz}$ and $\tau$=0.9$\left(K+1\right)\left(2M+2\right).$

3.4.2 Example 2: narrow selective passbands

In this example, the requests for passband and transition bands are omitted, which decrease a numerical burden. To obtain high selectivity, $N_B=\left(K+1\right)\left(2M+2\right)$ is kept. The obtained frequency responses (see Figure 7) show that selectivity and attenuation in the stopbands are increased. This is achieved thanks to the poles of the transfer function located on a circle of radius $0.95$ ($\left|z\right|=0.95$) very close to the zeros located in the harmonic frequencies, as well as the additional zeros of the polynomial $B\left(z\right)$. Such high selectivity caused a large increase in a group delay.

Figure 7.

(a) Frequency responses of $g_{3,0}^{\prime }{z}^{-1}{P}_3\left(z\right)$, $B\left(z\right)/A\left(z\right)$ and $T_{3,0}^{\mathit{AB}}\left(z\right)$and (b) Pole-zero map of $T_{3,0}\left(z\right)$ and $B\left(z\right)$ (for the third harmonic), for $N_A=\left(K+1\right)\left(2M+2\right)$, $N_B=\left(K+1\right)\left(2M+2\right)$, and $\rho =0.95,f_{\mathit{PB}}=0\mathrm{Hz}$.

3.4.3 Example 3: numerically cost-effective solution

This example is very similar to the previous one with different that now is $N_B=0,$ which means that there is no extension to the existing resonator structure with common feedback. Obtained frequency responses (see Figure 8) show that selectivity and attenuation in the stopbands are smaller than in the previous case, however, with a smaller group delay too.

Figure 8.

(a) Frequency responses of $g_{3,0}^{\prime }{z}^{-1}{P}_3\left(z\right)$, $B\left(z\right)/A\left(z\right)$ and $T_{3,0}^{\mathit{AB}}\left(z\right)$and (b) Pole-zero map of $T_{3,0}\left(z\right)$ and $B\left(z\right)$ (for the third harmonic), for $N_A=\left(K+1\right)\left(2M+2\right)$, $N_B=0$, and $\rho =0.95,f_{\mathit{PB}}=0\mathrm{Hz}$.

4. FIR cascaded-resonator-based harmonic phasor estimation

Instead of the common simultaneous compensation of the frequency responses for all harmonics through the compensating filter $B\left(z\right)$ placed in the front of the parallel resonator structure, a more flexible solution is shown in Figure 9 with postprocessing by the set of compensators $Q_m\left(z\right)$ designed particularly for each harmonic $m$. In this case, $A\left(z\right)=1$ is chosen to allow the use of linear optimization techniques such as LP and CLLS.

Figure 9.

Block diagram of the $K$-type CR-based harmonic analyzer with postprocessing.

In order to obtain an algorithm that can be utilized in a wide range of signal dynamics in a unified way and improve the frequency response, a linear combination of the differentiators’ outputs in the cascade can be used [22, 25, 26]. The goal of this compromised solution was to propose a tracking-mode harmonic estimation technique. In Ref. [31], it is shown that this estimation technique exhibiting maximally flat frequency responses can be efficiently used for implementation of P-Class Compliant PMU in accordance with IEC/IEEE Standard 60255-118-1:2018 for harmonic phasors estimation. In this approach, the order of the resulted compensation filter was low and equals to the pole multiplicity. In Refs. [21, 23, 24], the proposed approach was generalized to any necessary order through the postprocessing compensation FIR filters applied to the output signals obtained by the CR structure. The drawback of this approach is that we have to use as many postprocessing FIR filters as there are harmonic phasors that we need to estimate (one estimator per one harmonic phasor). On the other hand, the advantage is that it is possible to obtain a filter bank, surrounding the core CR structure, with a set of different compensation filters corresponding to different signal dynamics.

The transfer function for every mth channel has the form of

$$T_{m,0}^Q\left(z\right)=\frac{V_{m,0}^Q\left(z\right)}{V\left(z\right)}=T_{m,0}\left(z\right)Q_m\left(z\right)=\left[g_{m,0}^{\prime }{z}^{-1}{P}_m\left(z\right)\right]{\mathit{q}}_{\mathit{Q}}{\mathit{x}}_{\mathit{Q}}\mathbf{.}$$E40

where $Q_m\left(z\right)=q_{m,0}+q_{m,1}{z}^{-1}+\cdots +q_{m,N_Q-1}{z}^{-\left(N_Q-1\right)}+b_{m,N_Q}{z}^{-N_Q}$. ${\mathit{q}}_{\mathit{Q}}=\left[\begin{array}{ccc}1& z^{-1}& z^{-2}\end{array}\begin{array}{ccc}\cdots & z^{-\left(N_Q-1\right)}& z^{-N_Q}\end{array}\right]$, ${\mathit{x}}_{\mathit{Q}}={\left[\begin{array}{ccc}{q}_0& q_1& \cdots \end{array}\begin{array}{ccc}& q_{N_Q-1}& q_{N_Q}\end{array}\right]}^{\mathrm{T}}$.

Eq (40) has the same form as Eq. (5) with the following constraints: $A\left(z\right)=1$, $N_A=0,$${\mathit{x}}_{\mathit{A}}\mathit{=}\left[\right]$, ${\mathit{q}}_{\mathit{A}}\mathit{=}\left[\right]$, $\mathit{x}\mathit{=}{\mathit{x}}_{\mathit{Q}}$, $\mathit{q}\mathit{=}{\mathit{q}}_{\mathit{Q}}$. ${\mathit{q}}_{\mathit{B}}$ is replaced by ${\mathit{q}}_{\mathit{Q}}$, and ${\mathit{x}}_{\mathit{B}}$ by ${\mathit{x}}_{\mathit{Q}}\mathit{.}$ Since $N_A=0$, the stability constraints are not present. The desired frequency response is related only to the actual harmonic m and does not consider the frequency responses of the other ones.

4.1 Total vector gradient (TVG) calculation

The response time and delay of the estimator are directly correlated with the group delay (GD) of the filter. Due to the more complex calculation of GD, it is possible to use the gradient of the transfer function $d{T}_{m,0}^Q\left(z\right)/\mathit{dz}$, which is called the total vector gradient (TVG) here. In a flat range with small amplitude changes, TVG and GD are proportional, and optimization of one leads to optimization of the other. The first derivative of the transfer function $T_{m,0}^Q\left(z\right)$ is as follows:

$$\begin{array}{l}d{T}_{m,0}^Q\left(z\right)/\mathit{dz}=g_{m,0}^{\prime }\left[-z^{-2}{Q}_m\left(z\right)P_m\left(z\right)+z^{-1}{P}_m\left(z\right)d{Q}_m\left(z\right)/\mathit{dz}+z^{-1}{Q}_m\left(z\right)d{P}_m\left(z\right)/\mathit{dz}\right]\\ =g_{m,0}^{\prime }{z}^{-1}\left[P_m\left(z\right)d{Q}_m\left(z\right)/\mathit{dz}+{\Psi }_m\left(z\right)Q_m\left(z\right)\right]\end{array}$$E41

where

$$\frac{d{Q}_m\left(z\right)}{\mathit{dz}}=-q_{m,1}{z}^{-2}-2q_{m,2}{z}^{-3}-\cdots -N_Q{q}_{m,N_Q}{z}^{-\left(N_Q+1\right)}$$

$${\Psi }_m\left(z\right)=-z^{-1}{P}_m\left(z\right)+d{P}_m\left(z\right)/\mathit{dz}$$

$$\frac{d{P}_m\left(z\right)}{\mathit{dz}}=P_m\left(z\right)\sum _{i=-Mi\ne m}^{M+1}\sum _{k=0}^K\frac{z_{i,k}{z}^{-2}}{1-z_{i,k}{z}^{-1}}.$$

Eq. (41) can be written in a matrix form as follows:

$$d{T}_{m,0}^Q\left(z\right)/\mathit{dz}={\mathit{\psi }}_{\mathit{m}}{\mathit{x}}_{\mathit{Q}\mathbf{,}\mathit{m}}.$$E42

where

$${\mathit{\psi }}_{\mathit{m}}=\left[\begin{array}{ccc}{\psi }_{m,0}\left(z\right)& {\psi }_{m,1}\left(z\right)& \cdots \end{array}\begin{array}{ccc}& {\psi }_{m,N_Q-1}\left(z\right)& {\psi }_{m,N_Q}\left(z\right)\end{array}\right],$$

$${\psi }_{m,n}\left(z\right)=g_{m,0}^{\prime }{z}^{-\left(n+1\right)}\left[{\Psi }_m\left(z^{-1}\right)-n{z}^{-1}{P}_m\left(z^{-1}\right)\right],n=0,1,\dots ,N_Q.$$

4.2 Optimization criteria

Selections of the object and functions and constraints can be very different depending on the optimization criteria scenario. Herein will be considered three criteria summarized in Table 1 [21, 24]. In the first Criterion 1, the cost function in which absolute values are minimized is the transfer function $T_{m,0}^Q\left(z\right)$ in the stop bands. Weighting function $W_m\left(\omega \right)$ is settled to 1. The function which absolute values are kept under settled limits is an error in the passband $E_m\left(z_i\right)={\left.T_{m,0}^Q\left(z\right)\right|}_{z=z_i}-T_m^d\left(z_i\right)$ where $T_m^d\left(z_i\right)=exp\left(-\mathit{j\tau }\left({\omega }_i-{\omega }_m\right)\right).$ In addition to this, a control of overshoots in the transition bands is performed.

Criteria	Object function	Desired values	Frequency range	Constrained functions	Reference values	Frequency range
Criterion 1	$T_{m,0}^Q\left(z\right)$	0	Stopbands	$T_{m,0}^Q\left(z\right)$	$T_m^d\left(z\right)=e^{-\mathit{j\tau }\left(\omega -{\omega }_m\right)}$	Passband
				$T_{m,0}^Q\left(z\right)$	0	Transition band
Criterion 2	$T_{m,0}^Q\left(z\right)$	0	Stopbands	$d{T}_{m,0}^Q\left(z\right)/\mathit{dz}$	0	Harmonic frequency
				$T_{m,0}^Q\left(z\right)$	0	Transition band
Criterion 3	$d{T}_{m,0}^Q\left(z\right)/\mathit{dz}$	0	Harmonic frequency	$T_{m,0}^Q\left(z\right)$	0	Stopband
				$T_{m,0}^Q\left(z\right)$	0	Transition band

Table 1.

Considered design criteria.

Criterion 2 is similar with Criterion 1 with the difference that the absolute value of the TVG in the harmonic frequency is limited, that is ${\left|d{T}_{m,0}^Q\left(z\right)/\mathit{dz}\right|}_{z=z_m}\le {\left|\mathit{TVG}\right|}_{\mathit{max}}$, where ${\left|\mathit{TVG}\right|}_{\mathit{max}}$ is a maximally allowed absolute total vector gradient. Like in Criterion 1, the limitation of overshoots in the transition bands is necessary.

Criterion 3 minimizes the absolute value of the TVG in the harmonic frequency subject to the limitation of the gain in the stopband. Similarly with the previous cases, the limitation of overshoots in the transition bands is necessary.

4.3 Design example

In order to illustrate the described algorithms, examples overtaken from [21] are shown for Criteria 2 and 3 defined in Table 1. Figures 10 and 11 show the frequency responses of the transfer function of the third harmonic $T_{3,0}^Q\left(z\right)$ in the case of $K=2$. For Criterion 3, the maximum allowed gains in the stopband was selected as $l_m^{\mathit{SB}}\in \left\{0.1,0.01,0.001\right\}$, which corresponds to attenuations of $\left\{20,40,60\right\}\mathrm{dB}$. For Criterion 2, the maximum value of $\left|\mathit{TVG}\right|$ in the harmonic frequencies $z_m$, ${\left|\mathit{TVG}\right|}_{\mathit{max}}\in \left\{24,32,48\right\}$ was selected. Zoomed amplitude and $\left|\mathit{TVG}\right|$ characteristics around the harmonic frequency are shown in the inset figures at the bottom of the figures. It can be seen that the higher value of $N_Q$ gives a smaller value of $\left|\mathit{TVG}\right|$ and wider bandwidth. It is also visible that for smaller values $\left|\mathit{TVG}\right|$ sidelobes are larger, which reduces robustness to interharmonics and noise. In addition, the bandwidth increases for smaller values of $\left|\mathit{TVG}\right|$.

Figure 10.

Frequency responses for the basic ($T_{3,0}\left(z\right)$) and reshaped ($T_{3,0}^Q\left(z\right)$) transfer function for $K=2$, for $f_s=1.6\mathrm{kHz}$, $N_Q=16$ and (a) $l_3^{\mathit{SB}}\in \left\{0.1,0.01,0.001\right\}$ and (b) ${\left|\mathit{TVG}\right|}_{\mathit{max}}\in \left\{24,32,48\right\}$.

Figure 11.

Frequency responses for the basic ($T_{3,0}\left(z\right)$) and reshaped ($T_{3,0}^Q\left(z\right)$) transfer function for $K=2$, for $f_s=1,6\mathrm{kHz}$, $N_Q=32$ and (a) $l_3^{\mathit{SB}}\in \left\{0.1,0.01,0.001\right\}$ and (b) ${\left|\mathit{TVG}\right|}_{\mathit{max}}\in \left\{24,32,48\right\}$.

Figure 12 shows the frequency responses of the transmission functions $T_{3,0}\left(z\right)$ and $T_{3,0}^Q\left(z\right)$ for different values of given parameters (a) $l_m^{\mathit{SB}}$ and (b) ${\left|\mathit{TVG}\right|}_{\mathit{max}}$, in the case of $K=1$. In this case, the total $\left|\mathit{TVG}\right|$ is smaller than in the case of $K=2$. The order of the compensation filter $N_Q=16$ is smaller, so the bandwidth is narrower. In the case of $l_3^{\mathit{SB}}=0.001$, the optimization problem has no solution.

Figure 12.

Frequency responses for the basic ($T_{3,0}\left(z\right)$) and reshaped ($T_{3,0}^Q\left(z\right)$) transfer function for $K=1$, for $f_s=1,6\mathrm{kHz}$, $N_Q=16$ and (a) $l_3^{\mathit{SB}}\in \left\{0.1,0.01,0.001\right\}$ and (b) ${\left|\mathit{TVG}\right|}_{\mathit{max}}\in \left\{16,24,32\right\}$.

5. Conclusions

CR-based algorithms for harmonic analysis and estimation of harmonic phasors are described in this chapter. The resulting filters for extracting harmonic signals can be of the FIR or IIR type. Algorithms for the optimization of frequency responses are presented and corresponding examples of synthesis are given. Linearized mathematical models were used, which enabled the use of linear optimization methods such as LP and CLLS. When designing the IIR analyzers, a linearized iteration scheme based on the Rouche’s theorem was used to control the stability of the system. As for the optimization algorithms, they can potentially be improved by various modifications such as, for example, by nesting optimization loops related to different constraint conditions and/or objectives, and adaptation of iteration steps. It is notable that approximating result could be obtained heuristically thanks to the characteristic position of the pole and zeros. In addition, it seems that closed-form calculation expressions derivation could be possible. On the other hand, the FIR-type algorithm particularly optimizes frequency responses through the postprocessing compensation FIR filters applied to the output signals obtained by the CR structure. This approach allows the usage of a set of compensation filters corresponding to different signal dynamics.

Conflict of interest

The author declares no conflict of interest.