Design of Compressive Sensing Adaptive Taylor-Fourier Comb Filters for Harmonic Synchrophasor Estimation

In modern power systems, phasor measurements are expected to deal with challenging conditions, e.g., fast dynamics and high distortion levels. Taylor-Fourier Multifrequency models represent a promising solution, but their performance is strongly related to the accurate extraction of the signal spectral support. In this context, this paper proposes an enhanced method for support recovery that exploits the inherent block-sparsity properties of electrical signals. The proposed method is fully characterized in diverse and distorted test conditions, inspired by reference standards and real-world scenarios. The comparison against another Compressive Sensing based approach confirms the significant improvement in terms of both recovered support exactness and synchrophasor measurement accuracy.


I. INTRODUCTION
M ODERN power systems are rapidly evolving towards active power infrastructures, characterized by an ever-increasing penetration of renewable energy sources (RESs) [1], [2]. In this context, system operators are expected to deal with time-varying conditions to be properly tracked and monitored, and Phasor Measurement Units (PMUs) might represent a promising solution thanks to their remarkable accuracy and responsiveness [3], [4].
PMUs were originally conceived for transmission networks, where synchronous generators guarantee high inertia and dampened dynamics. In compliance with such a scenario, the reference IEEE/IEC Std. 60255-118-1 (briefly, IEC Std) defines signal models, test conditions and performance requirements for synchrophasor measurements [5]. In particular, the signal under investigation is approximated by a slowly-varying fundamental tone and stationary narrow-band spectral components [6], [7].
The recent literature has proposed to refine the signal model through a Taylor series expansion of the frequency components: the constant term corresponds exactly to the only one considered by the usual DFT approach, whereas the higher-order terms account for time-varying contributions [8]. According to this approach, estimates are computed by convolving the sampled signal with Taylor-Fourier Filters (TFFs), in turn obtained by inverting the aforementioned signal model. Among several different implementations [9], [10], [11], adopting the TFFs for PMU measurements allows for suitably addressing both static and dynamic conditions [12].
In order to minimize the interference coming from spurious contributions, the Compressive Sensing Taylor-Fourier Multifrequency (CS-TFM) algorithm was proposed [13]. First, a CS routine identifies the spectral support of the signal, so that it defines a parametric signal model that includes these components. Then, the corresponding TFF obtained from a least squares (LS) fitting allows for estimating the dynamic phasor associated with the fundamental. In this context, the algorithm estimation accuracy primarily depends on a successful support recovery, which dictates how well the resulting TFM model is able to match the actual signal. In particular, if it contains all the significant signal components, the corresponding TFF introduces properly located zeros mitigating the interference with harmonics and interharmonics, thus it guarantees a noticeable level of accuracy even when off-nominal test conditions are accompanied by high distortion levels [14]. Otherwise, if the support is incomplete or the included frequency components are not near those actually present in the signal, the estimation accuracy is heavily deteriorated and even the numerical stability of the estimation process is no longer guaranteed.
A performance enhancement of the CS-TFM algorithm is possible only if the support recovery procedure is made more stable and capable of locating components with a finer frequency resolution. To this end, in the present paper, we investigate the application of the recent concept of block-sparsity to regularize the CS-based super-resolution technique [15]. In particular, we formulate the TFM model by embedding additional a priori information about electrical signals in ac power systems: the constant ratio between fundamental and harmonic frequencies. This property promises to be particularly relevant in modern power systems, as RESs are typically interconnected via dedicated inverters whose switching circuitry produces significant harmonic pollution [16], [17]. Traditionally, harmonic estimation is performed by power quality meters that operate with long observation intervals (e.g., 200 ms) and timestamp accuracy in the order of milliseconds [18], [19]. In the Smart Grid scenario, it is desirable that PMU-like devices (i.e., synchronized, low-latency, high-reporting rate instruments) will not be limited to the analysis of the fundamental component, but they will be required to estimate the dynamic synchronized phasors associated with a restricted set of harmonic components [20], [21], [22]. Therefore, in this paper we develop a novel support recovery method, named CS-COMB, that yields a significant performance enhancement with respect to the previous formulation in terms of both estimation accuracy and robustness. CS-COMB relies on a dictionary of harmonic blocks, specifically designed to maximize the spectral resolution while reducing the coherence between adjacent blocks. In turn, the resulting TFFs (for fundamental and harmonics) contain a comb of evenly spaced zeros to suppress spectral interference from fundamental and harmonics.
In order to validate the proposed method, we carry out a thorough performance characterization under diverse and distorted test conditions. In more detail, we reproduce a more realistic operating scenario by combining the IEC Std test conditions with harmonic distortion levels compatible with Standard EN 50160 [23]. For this analysis, we compare the CS-COMB based algorithm against CS-TFM algorithm in terms of support recovery and harmonic phasor estimation.
Furthermore, accuracy has been also compared to that obtained with another high-performance method enabling harmonic measurements that was recently proposed in the literature, based on a completely different approach [24]. The obtained results prove a significant enhancement in all the considered test conditions and confirm the potential of block-sparsity constraints in power system monitoring and tracking applications.
The paper is organized as follows. Section II recalls the principles of spectral support recovery and the inherent resolution limits of the CS-TFM algorithm. In Section III we introduce the concept of block-sparsity, and discuss its possible application to a novel support recovery method, thus to a new harmonic synchrophasor estimation algorithm. Section IV presents a numerical validation of the proposed algorithm, and Section V provides some closing remarks.

A. THE TAYLOR-FOURIER MODEL
Let us introduce a generic electrical signal x(t) in an ac power system having rated frequency f 0 (corresponding to the angular frequency ω 0 ), acquired with sampling time T s . Let us suppose that a symmetric interval centered around time instant t r is observed, thus resulting in an odd number N w of samples. Excluding highly dynamic events, an accurate representation of the collected data is obtained by assuming that the considered signal segment can be decomposed into a set of N frequency components having slowly-modulated magnitudes and phase angles. In terms of equations, we have is the instantaneous phasor of the lth component using ω l as reference angular frequency. In order to set a one-to-one relationship between signal and model, we assume that ω l corresponds to the actual frequency at the reporting instant. This implies that the derivative of ϕ ω l (t) in t r is zero. Adopting the Taylor-Fourier approach,X ω l (t) is approximated with its Taylor expansion around t r truncated to the K ω l th order, thus leading to the signal model is the kth order derivative of the lth phasor in t = t r . Usually, the real part operator in (2) is removed by introducing complex-conjugated couples of modulated rotating exponentials. Using trigonometric expressions, (2) can be also rewritten as a linear combination of modulated sine-cosine pairs as the real and imaginary part ofX (k) ω l (t r ), respectively. When passing to the discrete-time domain, the N w samples of x(t) can be arranged in a column vector x(t r ). Adopting the signal model (3), we have . . .
K = H ω 1 · · · H ω N defines the model structure, while p = g ω 1 · · · g ω N contains the model parameters. In general, they are functions of the reporting instant, but it is omitted in the following for a lighter notation. The lth submatrix or subvector refers to the lth spectral component and they can be partitioned as ω l refer to the kth order contribution in the Taylor expansion of the lth spectral component. Reminding (3), H (k) ω l is a N w × 2 matrix whose first and second column are respectively

B. DICTIONARY, SUPPORT RECOVERY AND ESTIMATION
When looking at (4), the vector of the samples x(t r ) is approximated by a linear combination of vectors, namely the columns of K. Each of them belongs to a broad set of vectors (atoms) defined by (6) that represents a dictionary of candidate atoms. Typically, since x(t r ) is obtained from a segment of a power system signal, it can be accurately modeled considering a concise set of components, namely its support that, for this reason, can be considered as sparse, thus framing the component identification and measurement problem from a CS point of view. Moreover, vectors appear in blocks of amplitude-modulated cosines and sines sharing the same frequency; in this respect, matrix H ω l represents the generic block and thus x(t r ) is said to be block sparse. This feature is exploited by the CS-TFM algorithm (from here on CS-TFM) for synchrophasor, frequency and ROCOF measurements. A grid of possible angular frequency values { d } d∈{1,...,N D } is generated a priori by selecting a resolution δω, typically much lower than the DFT resolution corresponding to the observed time interval and thus called super-resolution. Furthermore, a Taylor expansion order is adopted, which could be different for each frequency. These choices define a finite-dimensional dictionary, thus the set of N D matrices {H d } d∈{1,...,N D } corresponding to the blocks that potentially could be combined to compose the matrix K associated with the recovered support. Once it has been obtained, it leads to a TFM signal model. An estimatep(t r ) of its parameters can be computed by minimizing the norm of the residual, that is the mismatch between measured samples and their reconstruction from the model, thereforê with · denoting the Euclidean norm. Assuming thatK has full column rank, we havê whereK † is the Moore-Penrose inverse ofK. Let us suppose that subvectorĝ d is that whose reference frequency d is the closest to the fundamental. Assuming that an expansion order ≥ 2 has been selected for this block, fundamental synchrophasor, frequency and ROCOF estimates can be computed.
The key challenge to be faced is obtainingK, namely retrieving the set of dictionary blocks that, for a given cardinality, defines the TFM model resulting in the best LS approximation of x(t r ). This corresponds to finding the spectral support, i.e., the main components of the signal in the considered window, once having adopted a suitable expansion order for each component.
Considering practical applications, finding the best spectral support means solving an NP-hard optimization problem and mostly results in huge computational burden that does not comply with reporting rates and maximum latency requirements for PMU applications. For this reason, CS-TFM adopts Orthogonal Matching Pursuit (OMP) for support recovery [13]. On the one hand, it may lead to a suboptimal solution, on the other hand it enables a drastic reduction of computational complexity. OMP is an iterative procedure, so that before iteration number n we have a model with n − 1 spectral components, a recovered matrixK [n−1] , an estimate p [n−1] of the model parameters and the corresponding residual r [n−1] . At the first step, n = 1,K [0] is empty while During the generic nth step, residual is projected on all the zeroth order frequency components H (0) d present in the dictionary. We obtain: From a geometric point of view, u d, [n] and v d, [n] are obtained by projecting the signal samples on a rectangular coordinate system rotating with angular speed d . According to the OMP paradigm, the block whose zeroth order terms are able to capture the highest energy is added to the support: the reason is that it is likely to produce the highest residual energy reduction. Since it can be easily shown that (c Once having updated the support, signal parameters and the corresponding residual are updated using (8). The process stops when a convergence criterion has been reached, e.g., when a maximum number of blocks has been identified (each one corresponds to an iteration), or when the residual energy falls below a predetermined threshold, thus meaning that the samples have been modeled with satisfactory accuracy, also taking into account the signal to noise ratio. Even if CS-TFM is a computationally intensive algorithm, a real-time implementation is feasible through a high performance FPGA [25]. Finally, it is worth noting that the method is generally presented by considering a dictionary made of modulated rotating and counter-rotating exponentials. This alternative formulation, instead, reaches identical results, with the advantage of dealing only with real-valued quantities.

III. ENHANCED SUPPORT RECOVERY THROUGH COMB DICTIONARY BLOCKS
The performance obtained with CS-TFM is strictly connected with the capability of properly retrieving the spectral support of the signal. It is rather intuitive to understand that it depends on the degree of similarity between the zeroth order blocks to be projected on the residual. In fact, similar blocks result in similar projections, thus there is higher chance to select a non-optimal block to be included into the support. The recent literature about CS has discussed several methods to quantify the performance of different block configurations to retrieve the correct signal support. In particular, inter-block coherence [15] has been introduced, namely a number between 0 and 1 that measures their similarity level. It is evident that the resolution of the frequency grid (summarized by δω if it is regularly spaced) employed to generate the dictionary has a strong impact on inter-block coherence, and thus on performance. It comes out from a critical tradeoff, since a denser dictionary potentially results in a better approximation of the actual spectral support of the signal. However, there is a higher risk to make mistakes during OMP-based recovery.

A. NEW SIGNAL MODEL AND COMB BLOCKS
A possible idea to improve support recovery consists in introducing further a priori knowledge, specifically about the peculiar spectral content of electrical signals in ac power systems. By definition, excluding abrupt transients, these signals are quasi periodic, hence most of the energy is confined around the fundamental frequency and its harmonics. Subharmonics and interharmonics could be present, but their energy content is expected to be much lower. According to these considerations, the signal model can be rewritten as: where: while: Moving to the discrete time domain, according to (11) we can write x(t r ) ≈ x c (t r ) + x nh (t r ), thus assuming that the spectral support of the samples consists of a "comb" of frequencies multiples of ω 1 , plus (few) other components. Equation (4) still applies, but we have: where: The new vector of the signal parameters becomes being This suggests introducing a new class of blocks in the dictionary, whose structure resembles that of C ω 1 . They can be generated by selecting a set of possible fundamental frequencies { F,c } c∈{1,...,N DF } , a maximum harmonic order N H , and a Taylor expansion order for each one. In addition, the dictionary also comprises another set of blocks devoted to represent other components that may be present in the signal, defined by their angular frequencies { I,i } i∈{1,...,N DI } and the corresponding expansion orders.

B. SUPPORT RECOVERY AND ESTIMATION PROCEDURE
In order to exploit the new dictionary defined in the previous paragraph, support recovery should be divided into two stages. The target of the first one is finding the comb block that allows the best modeling of the contribution x c in x. The second stage looks for non-harmonic components, basically following the procedure already discussed in Section II-B.
Going back to the first step, the OMP approach could be adopted also to identify the best comb block. From each comb block in the dictionary, we extract just the zeroth order terms, thus obtaining {C The target is selecting the comb whose projection has the highest total energy. However, since in general the components of C are not orthogonal, we cannot compute it simply as e c 2 . A possible solution is precomputing orthonormal zeroth order comb blocks through economy size QR decomposition: Signal is then projected on the set of orthonormal blocks {Q is added to the support. The approach significantly reduces the number of iterations with respect to the classic CS-TFM method, since all the N H components included in the comb block are added to the support in a single step. Preliminary signal parameters are estimated, together with the corresponding residual. If its energy is below a predetermined threshold, the process stops, otherwise the second stage of the procedure is executed. The target is including in the TFM model significant interharmonics or harmonics above the N H th order. Algorithm 1 reports the pseudocode of the CS-COMB support recovery.
It is worth highlighting that, during the first recovery stage, the adoption of high-cardinality comb blocks instead of H (0) d has the clear advantage that the larger number of atoms in the first decreases their mutual coherence, thus resulting in a more effective discrimination between them. On the other hand, we are introducing a stiffer underlying signal model that may capture more noise when some frequency components are not present in the signal. These considerations should be taken into account when selecting N H , which determines the size of the comb blocks. A block composed only of most probable harmonic orders might also be considered.
As a final remark, a possible application enabled by an accurate support recovery is harmonic synchrophasor measurement. In this case, they are easily obtained by combining residual update: r = x −K † x end return the elements of the estimated signal parameters that correspond to the zeroth order harmonic terms captured by the comb.

IV. ALGORITHM VALIDATION
In this section we validate the algorithm proposed in the previous section (also called CS-COMB from here on for the sake of brevity) under realistic operating conditions simulated in MATLAB programming environment. We design the test conditions in order to reproduce a synchrophasor measurement scenario, inspired by not only the compliance requirements of the IEEE Std [5], but also considering plausible distortion levels in distribution networks, as inferred by EN50160 specification [23] and real-world scenarios published in the recent literature [26], [27], [28]. The aim is to assess the support recovery capability of the CS-COMB and, as a consequence, its fundamental and harmonic synchrophasor estimation performance. CS-TFM is used as benchmark since it is based on the same CS principle and it may exploit the same basic dictionary (just not grouped into comb blocks), thus enabling to carefully analyze the brought improvements.
In all the tests discussed in the following, we consider the following base configuration: 1) Sampling rate f s = 5 kHz and window length N w = 401. The window does not include an integer number of cycles at nominal frequency (f 0 = 50 Hz, thus leading to about 4 cycles per window) since we want to explore a condition where the super-resolved grid does not exactly fit the signal characteristics. Indeed, it would be helpful to choose a matched condition in terms of grid frequencies and N w , but here we aim at stressing the method. 2) Fundamental is assumed to be in the range (f 0 − 5 Hz, f 0 + 5 Hz) (a slightly broader range is used in practice for completeness). Frequency combs made of N H = 5 components have been generated by considering 1 Hz resolution at the highest harmonic; the corresponding frequency components have been added to the dictionary. It is worth noting that the resulting frequency separation is inversely proportional to the harmonic order and at fundamental frequency becomes 1 N H Hz. 3) 1 Hz frequency resolution has been used for generating the sub-dictionary for the components not included into a comb. 4) Order 2 (i.e., K ω 1 = 2) is used for the fundamental (which requires ROCOF computation), while order 1 is adopted for all the other components. 5) The reporting rate is set to 100 frames/s, i.e., an estimate every 10 ms, and test duration is 5 s.
The configuration of the method is extremely flexible and different choices are also allowed, e.g., to include additional prior information about the signal (limited frequency ranges, forbidden bands, etc.). The sampling rate has to be selected according to the maximum harmonic order to be measured and considering computational constraints. In practical applications, a proper low-pass filter is typically employed to avoid aliasing artifacts and suppress high-frequency disturbances. The length of the time window has major impact on support recovery and thus on achieved performance. For a given frequency separation between dictionary blocks, their coherence increases as the observed time interval is reduced. The four cycle length has been selected as a favorable tradeoff between accuracy and latency, as typically happens for PMU algorithms.
Before discussing the results, it is extremely interesting to do some preliminary considerations about the blocks adopted by the two methods, obtained from the same set of dictionary atoms. During the first recovery stage, CS-COMB uses zeroth-order, orthonormal comb blocks, and their maximum coherence is 0.099. This is rather impressive if we consider that the maximum coherence between blocks representing single spectral components is 0.994. Therefore, CS-COMB is expected to provide a much more robust support recovery for fundamental and harmonics.
Accuracy of the estimates provided by the two methods is quantified in terms of percent total vector error (TVE), together with absolute frequency and ROCOF error (|FE|, |RFE|).
The test signals contain 4 harmonics (orders 2 to 5) with random phases and amplitudes within the limits of EN50160. In particular, the results with harmonic levels 1.1 %, 6.1 %, 0.5 % and 4.9 % for h = 2, . . . , 5, respectively, are reported. Fundamental frequency is varied between 45.05 and 54.95 Hz with 0.1 Hz step. First of all, it is interesting to analyze the support recovery performance. Figure 1 shows the detection rate, namely, the percentage of correct support recovery (considering all the reporting instants in the test), for the fundamental and the harmonic components as function of the fundamental frequency. A component is considered as correctly identified if the closest frequency available in the dictionary has been added to the support. It is worth noting that, because of the selected fundamental frequency values, synchronous sampling never occurs. The results show that CS-COMB always picks the correct frequencies from the grid at every harmonic order, whereas CS-TFM shows correct detection rates down to 14 % for the fundamental and down to 0 % and 38 % for second and third harmonic, respectively. Higher harmonics are less influenced by the fundamental tone leakage and thus also CS-TFM properly detects them in most of the cases. The above result is quite self-explanatory, but we do not want it to appear misleading or overreaching. Identifying a suboptimal set of frequencies does not mean that they could not be estimated, but it clearly reduces accuracy due to a model that is less capable of fitting the signal, thus introducing definitional uncertainty. From a different point of view, a mismatch between the actual frequency component in the signal and the counterpart embedded into the model produces scalloping loss when that component has to be extracted and spectral interference when evaluating the others. To better study this aspect, Fig. 2 reports, for both methods, the root mean square (RMS) frequency distance of each retrieved component with respect to the optimal value belonging to the dictionary as function of the fundamental frequency. Its RMS is computed among all the measurements in a test and gives further insight into the goodness of support recovery.  The results follow the conclusions of Fig. 1, but reveal that incorrect identifications can be significant (the RMS is up to 0.33, 0.43 and 0.47 Hz for h = 1, 2 and 3, respectively), thus making support recovery less effective. The larger the distance (it is worth reminding that it is an averaged value) the more likely estimation performance is jeopardized. Figure 3 shows an example of the estimation results for the fundamental component. With CS-COMB, |FE| is reduced by one order of magnitude and the same happens to TVE, but the results are not reported here for the sake of space and because, for both methods, TVE is always below 0.01 % and thus the synchrophasor is accurately measured.
|RFE| goes from a maximum value of 0.058 Hz/s (corresponding to f = 51.45 Hz) with CS-TFM to 0.0046 Hz/s with CS-COMB, thus showing that the proposed support recovery allows reducing error of about ten times also for the most sensitive PMU measurement. CS-COMB |RFE| is thus always below 0.006 Hz/s, that is compliant even with the strictest requirement reported in [29].
Even more interesting are the results of harmonic synchrophasor estimations. Figure 4 shows an example of the results. CS-COMB reduces the TVE at the second harmonic by more than 64 % (up to 95 %), thus achieving a TVE below 0.05 % for all the considered fundamental frequency values. This is somehow expected since Figs. 1 and 2 already highlight that the second harmonic is matched during the first step of the CS-COMB algorithm. However, the results are much more interesting, as it can be deduced from Fig. 5, where TVE values are reported for h = 4. The fourth harmonic is indeed correctly identified also by CS-TFM during support recovery (see Figs. 1 and 2) but CS-COMB provides much more accurate estimates. This reveals the side effect of an optimal identification of the support: all the components are located in the best available position on the grid and, thus, the bank of filters that correspond to the rows of the pseudoinverse in (8) is optimally tuned both in terms of passband (i.e., properly centered on the harmonic frequency of interest) and in terms of zeros (located near the other harmonic frequencies to be cancelled). For this reason, with CS-TFM the fourth harmonic estimation is affected by the spillover of the fundamental and other harmonics that are not fully rejected.
Similar results can be found for all the estimated synchrophasors, summarized in Table 1. Both maximum and average TVE, |FE|, and |RFE| values across all the fundamental frequencies are reported. TVEs are available for all the components, because they correspond to the estimated harmonic synchrophasors, while there are unique frequency and ROCOF values.
To assess the robustness of the above results in the presence of wideband noise, tests with different signal to noise  ratios (SNRs) have been performed adding white and uniform noise (AWUN) to the signal of the previous tests. In particular, the same harmonic orders and levels are considered. Figure 6 shows the FE results of the two methods as function of the SNR (from 60 dB to 90 dB with a 5 dB step) when the fundamental frequency is 50.65 Hz. RMS values are used to prevent outliers due to random extractions, thus allowing a more significant comparison.
The difference between the two methods is small for the lowest SNR values, because noise tends to mask the other effects. However, the advantage brought by CS-COMB is always significant and at higher SNRs it reflects the enhanced support recovery capability as discussed above. Table 2 reports the results at three different SNRs in terms of TVE, FE and RFE for the fundamental, and TVE for the two lowest harmonics. The results underline again the advantages of the proposed method. The improvement can be partially concealed for some indexes (TVE and ROCOF) and specific components by higher noise (SNR = 60 dB), i.e., when the key factor is not the algorithm itself but the equivalent noise bandwidth of the filter used to extract the  component. Clearly the impact can be reduced using a higher sampling rate or an elementary prefiltering stage.
To further stress the disturbance immunity, the interharmonic rejection performances achieved by CS-TFM and CS-COMB have been compared. In this respect, a 1 % interharmonic component has been superimposed to the base signal employed in the previous tests, with fundamental frequency equal to 50.65 Hz and containing harmonics up to the 5th. Interharmonic frequency has been swept within the band prescribed by [5] for 50 fps reporting rate. When considering the fundamental synchrophasor, both the methods result in very low TVEs. However. the values obtained by CS-COMB are considerably smaller, in particular when the interharmonic frequency is within 75 − 90 Hz, hence when there is a strong spectral interference between interharmonic and second harmonic. Differences are more significant when looking at the RMS value of FE and RFE: they may exceed 2 mHz and 0.3 Hz/s with CS-TFM, while they remain always below 0.5 mHz and 0.06 Hz/s when CS-COMB is adopted. This performance improvement occurs thanks to the enhanced support recovery that characterizes CS-COMB: also in this case 100 % success rate is reached for fundamental and harmonics. This is reflected also on harmonic phasor estimates: for example, the RMS TVE value at the 4th harmonic exceeds 1 % with CS-TFM, while it is lower than 0.23 % with CS-COMB.
Then, we tested also the tracking capability of the method when harmonics move in the spectrum. In particular, frequency ramp in presence of the same harmonic pattern as before has been adopted. The ROCOF is 1 Hz/s and the fundamental frequency evolves from 48 Hz to 52 Hz in 4 seconds. Harmonic frequencies are tied to the fundamental and are thus scaled-up replicas according to their harmonic orders. Table 3 reports the maximum TVEs for fundamental and harmonics along with maximum |FE| and |RFE|. CS-COMB still performs better for all the considered estimates, even if for some quantities the dynamic effects start to be the main error source. The table also shows the detection rate and the RMS frequency distance after support recovery, thus highlighting that the proposed approach always finds the best set of frequencies, while CS-TFM may have an average frequency distance up to about 5/4 of the grid resolution (e.g., for the fundamental frequency).
Finally, we evaluated the performance of CS-COMB under real-world operating conditions. Reference waveforms have been synthesized starting from a current sampled in the medium voltage distribution network of the EPFL campus in Lausanne, Switzerland (further details about the measurement setup in [21]). Due to the relevant penetration of inverter-connected resources, the signal is characterized by a significant harmonic distortion. Consistently with the previous tests, this analysis considers only the first five harmonic terms whose normalized amplitude is equal to 0.008, 0.004, 0.002, and 0.054, respectively; frequency is around 50.15 Hz. As benchmark, we consider the CS-TFM and the Iterative Interpolated DFT (i-IpDFT) method [24]. In the first case, we determine the improvement introduced by the enhanced support recovery. In the second case, instead, we compare the CS-COMB against an inherently static estimator that minimizes spectral leakage contributions. In this regard, Table 4 reports the maximum TVE, |FE|, and |RFE| values. In all the considered error metrics, the CS-COMB outperforms the two benchmarks. In particular, it is interesting to observe that the enhanced support recovery guarantees a TVE not larger than 0.1 % for all the considered harmonic orders.

V. CONCLUSION
Due to the ever-increasing penetration of RESs, power systems are characterized by faster dynamics and higher distortion levels. In such a challenging scenario, CS-TFM proved to be a viable and robust solution for accurate PMUbased measurements, but its performance is dependent on the definition of the signal spectral support.
The CS theory allows for super-resolved spectral estimation, without affecting the measurement reporting interval and latency. By exploiting the inherent block-sparsity properties of electrical signals, this paper proposed an enhanced support recovery method. A thorough characterization in diverse and distorted conditions confirmed the significant performance improvement with respect to traditional OMP-based approaches and proved a remarkable estimation accuracy not only in the recovered spectral support, but also in the dynamic phasors associated with fundamental and harmonic components.