A High-Resolution Algorithm for Supraharmonic Analysis Based on Multiple Measurement Vectors and Bayesian Compressive Sensing

: Supraharmonics emitted by electrical equipment have caused a series of electromagnetic interference in power systems. Conventional supraharmonic analysis algorithms, e.g., discrete Fourier transform (DFT), have a relatively low frequency resolution with a given observation time. Our previous work supplied a signiﬁcant improvement on the frequency resolution based on multiple measurement vectors and orthogonal matching pursuit (MMV-OMP). In this paper, an improved algorithm for supraharmonic analysis, which employs Bayesian compressive sensing (BCS) for further improving the frequency resolution, is proposed. The performance of the proposed algorithm on the simulation signal and experimental data show that the frequency resolution can be improved by about a magnitude compared to that of the MMV-OMP algorithm, and the signal frequency estimation error is about 20 times better. In order to identify the signals in two adjacent frequency grids with one resolution, a normalized inner product criterion is proposed and veriﬁed by simulations. The proposed algorithm shows a potential for high-accuracy supraharmonic analysis.


Introduction
With the popularity of power electronic equipment, such as photovoltaic (PV) inverters, electric vehicle chargers, switching power supplies, new lighting devices and household devices, more and more supraharmonics [1] in a frequency range of 2 kHz to 150 kHz are emitted into the power systems. Supraharmonics can cause a series of electromagnetic interference issues [2], e.g., it can affect the power line communications, such as automatic meter readers, lower the measurement accuracy of a power meter [3], generate undesired noise in electrical equipment, etc. In recent years, more and more research has been focusing on supraharmonics. For example, the IEC SC 77A, a technical committee for standardization in the field of electromagnetic compatibility with regards to low frequency phenomena, formed four working groups engaged in compatibility, emission, immunity and measuring method for supraharmonics. The European Committee for Electrotechnical Standardization (CENELEC) has gathered many electromagnetic interference cases in this frequency range. Related international conferences have set up innovation forums or special sessions for supraharmonics, such as the International Conference on Electricity Distribution (CIRED) and the IEEE Power and Energy Society Conference. Undoubtedly, supraharmonics have become a new research hot-spot.
In order to research the propagation property, interaction mechanism, emission limits, and mitigation for supraharmonics, it is necessary to measure them accurately. In the newly revised standard IEC 61000-4-30 [4], annex C, three supraharmonic measurement methods are recommended. One method is the International Special Committee on Radio Interference (CISPR) 16-1-2 [5]. This method uses a frequency-scan strategy. It assumes the supraharmonic distribution does not change during the measurement, which, however, in reality is not always true and could mislead the measurement result. The other two methods are DFT-based methods. This kind of method has the fastest calculation speed, but is limited by the spectrum leakage, the DFT-based method typically has a low frequency resolution and measurement accuracy. The pros and cons of the three methods are given in Table 1. Table 1. The pros and cons of the three methods.

Method Pros Cons
CISPR 16-1-2 high accuracy long measurement interval IEC 61000-4-7 gapless low frequency resolution and measurement accuracy IEC 61000-4-30 fast measurement speed low frequency resolution and measurement accuracy In order to improve the frequency resolution and mitigate the spectrum leakage, [6] proposed a supraharmonic algorithm based on a single measurement vector (SMV) model [7][8][9] and OMP algorithm [10][11][12], which can improve the frequency resolution by an order of magnitude for a single measurement vector at a time. To fulfill the multiple measurement occasions, ref. [13] proposed an enhancement supraharmonic algorithm combined MMV and OMP. The MMV-OMP [14,15] algorithm can achieve high-resolution analysis for multiple spectra vectors simultaneously and the frequency resolution can be improved by one order of magnitude. However, the MMV-OMP algorithm needs to estimate the sparsity of the measured signal in advance [13], and it requires the orthogonality between the basis functions in any subset of the sensing matrix. As a result, when high frequency resolution is desired, the orthogonality becomes difficult to be satisfied, and the interpolation factor is not much larger than an order of 10 [10].
A recent study on the compressive sensing technique shows a potential for accurate signal frequency analysis. Shihao Ji et al. applied the relevant vector machine (RVM) [16] and Sparse Bayesian Learning (SBL) [16] to the compressive sensing reconstruction, a Bayesian compressive sensing (BCS) [17,18] reconstruction method is proposed from the perspective of Bayesian statistical learning. Furthermore, by exploiting the joint sparsity between multiple tasks, a multi-task BCS algorithm for the need to reconstruct multiple measurement vectors is proposed also [19]. However, the BCS algorithms are only suitable for reconstructing the sparse real-valued signals, and are not suitable for analyzing the sparse complex signals, such as the spectrum.
In recent years, some new BCS sparse reconstruction algorithms for MMV problems have been proposed, i.e., in [20], to recover the joint sparse MMV signals with impulsive noise, Zhuhong Zhang et al. proposed a variational Bayesian method with two three-level hierarchical Bayesian models to characterize the joint sparsity and the impulsive noise respectively. In [21,22], Mohammad Shekaramiz et al. proposed a new Sparse Bayesian Learning (SBL) method with a binary vector, which was modeled as a hierarchical Bernoulli-Beta distribution, to characterize the supports for solving the MMV joint sparse reconstruction with unknown clustering patterns. In [23], Jian A. Zhang et al. proposed a MMV based implementation for SMV-SBL problems by converting a SMV model to a MMV one to reduce the computational complexity. In [24], Xiaolin Huang et al. proposed a re-weighted Bayesian method for multiple vectors, and a block coordinate decent technique was used to reconstruct the joint sparse vectors sequentially. In [25], Hang Xiao et al. introduced a MMV pattern-coupled sparse Bayesian learning algorithm to reconstruct block-sparse signals simultaneously, in which a pattern-coupled hierarchical Gaussian prior model was proposed to characterize both the block-sparsity of the coefficient and the statistical dependency between neighboring coefficients.
RVM is the simplified model for Bayesian model selection and has a relatively high computation efficiency [17]. In this paper, a Bayesian compressive sensing algorithm for multiple DFT vectors is proposed, denoted as MMV-BCS, which uses the extended RVM model [19] and can achieve high-resolution analysis for multiple spectrum vectors of supraharmonics simultaneously from the Bayesian viewpoint. Specifically speaking, by introducing an interpolation factor F, a MMV compressive model for multiple DFT vectors is constructed. The likelihood function is considered as a Gaussian distribution and a Gaussian-Gamma hierarchical prior model, with hyperparameters introduced to promote sparsity for the multiple high-resolution spectrum vectors. Then, a robust Student-t posterior distribution, according to the Bayesian rule, is obtained. Accordingly, an iteration algorithm is proposed to reconstruct multiple high-resolution spectrum vectors simultaneously, which can improve the DFT frequency resolution by F times while significantly suppressing the spectrum leakage. The MMV-BCS algorithm can provide some distinct advantages over the MMV-OMP algorithm, such as providing the posterior distribution of high-resolution spectrum vectors, automatic determination of joint sparsity and automatic estimation of model parameters, etc. The analysis shows that the proposed MMV-BCS algorithm has a good robustness with high interpolation factors. The frequency resolution and the accuracy of estimation can be improved at least one order of magnitude compared to the MMV-OMP algorithm, and at least two magnitudes compared to DFT-based methods.
The rest of the paper is organized as follows: In Section 2, the principle of the MMV-BCS algorithm is presented. Section 3 introduces numerical simulations to check the performance of the proposed algorithm in analyzing supraharmonics. The evaluation is compared to our previous work, i.e., MMV-OMP method. In Section 4, we discuss the presented algorithm based on experimental data. In Section 5, the conclusion is drawn.

MMV Compressive Sensing Model
According to IEC 61000-4-30, for a 200 ms power signal containing supraharmonics, a digital high-pass filter is used to remove the fundamental component (50 Hz or 60 Hz) and conventional harmonics below 2 kHz, and then the filtered data is uniformly divided into 400 sets of small data blocks with ∆T = 0.5 ms by applying a rectangular window function. Supraharmonic components in a small data block can be expressed as a multi-tone signal, i.e., where A sh , f sh and θ sh are respectively the magnitude, frequency and initial phase of each supraharmonic component, and in this case stationary parameters are assumed. T s = 1/ f s is the sample interval and the data length is N = f s ∆T. The DFT of x(n) can be expressed as where frequency resolution is defined as ∆ f = f s /N, and D N (·) is the Dirichlet kernel function, i.e., In order to improve the frequency resolution, an interpolation factor F is introduced. Therefore, the frequency resolution is refined to ∆ f = ∆ f /F and the total spectrum grid is changed from N to M = NF . The final spectrum expression in the new frequency resolution can be written as where f sh T s ≈ r/M. Equation (4) is a SMV compressive sensing model, and it can be written in the form of vectors, i.e., y ≈ Da + w, (5) where y is the original DFT vector with a dimension N × 1, D is the sensing matrix of dimension N × M, a is a K sparse high-resolution spectrum vector to be reconstructed of dimension M × 1, w is an independent-identically distributed (i.i.d.) and zero-mean Gaussian noise vector N(0, δ 2 ). Because Equation (5) is under-determined, by exploiting the sparsity of a, it can be solved by an l 1 -regularized formulation [8], i.e., where the parameter κ(κ > 0) controls the importance of the Euclidian error and the sparseness term in (6). Equation (6) is the basic framework of several CS inversion algorithms.

Sparse Bayesian Model
In order to analyze multiple spectrum vectors simultaneously, L sets of DFT vectors are assembled into a spectrum matrix {y i } i=1,L , and in such a way the MMV compressive sensing model is constructed, i.e., where y i denotes the i-th DFT vector, a i the i-th high-resolution spectrum vector to be reconstructed, w i the corresponding noise vector with an unknown noise inverse variance β 0 = δ −2 . Note that {a i } i=1,L is joint sparse, and {w i } i=1,L is a zero-mean Gaussian noise. Through the central limit theorem, the likelihood function can be obtained as a multi-parameter Gaussian likelihood model, i.e., For the prior distribution of {a i } i=1,L , [17] chose a Laplace prior function, and it is demonstrated to comply with Equation (6) by using the maximum a posteriori (MAP) estimation for a SMV in {a i } i=1,L . The major drawback of using Laplace prior is that it is not conjugate to the Gaussian likelihood and hence the analytical solutions cannot be obtained. In [17], the RVM was adopted, a hierarchical Gaussian-Gamma prior is applied to a SMV in {a i } i=1,L , which works for Gaussian likelihood estimation and also promotes the sparsity of {a i } i=1,L . The hierarchical Gaussian-Gamma prior, however, requires a large scale of calculations during the iteration. In order to accelerate the calculation, in [19], a fast algorithm is introduced, but the model is sensitive to initialization of the hyperparameter β 0 , and an inappropriate initial value could significantly lower the algorithm accuracy. Aiming to solve the above issues, a fast marginal likelihood maximization for the sparse Bayesian model, proposed in [19], is applied to extend the RVM in this paper. The prior distribution of {a i } i=1,L is defined as a Gaussian distribution function with double hyperparameters β 0 and β while the hyperparameter β 0 of the noise is defined as a Gamma distribution and is called a hyper-prior. The main advantage of the model is that the distribution function is conjugate to the posterior function and likelihood function, β 0 can be integrated out conveniently, and has better robustness.
The zero-mean Gaussian prior distribution with double hyperparameters is defined for each element of {a i } i=1,L , and the Gamma hyper-prior is defined for the noise inverse variance β 0 [19,21], as follows.
where β = {β j } j=1,M is composed of M independent hyperparameters, and each hyperparameter independently controls the corresponding elements of the vectors {a i } i=1,L . a and b are the shape parameter and scale parameter of the Gamma distribution, respectively. Since the hyper-parameter β = {β j } j=1,M is shared by all L sets of vectors, {y i } i=1,L will contribute to the estimation of hyperparameters β, meaning the inter-relationships among the multiple spectrum vectors can be fully utilized. According to Equations (8)-(10), the MMV hierarchical Bayesian model is shown in Figure 1, and the dependence among different parameters in the probability model is also shown.
with the mean and covariance of {a i } i=1,L given by with B = diag{β 1 , β 2 , · · · , β M }. (·) H denotes the conjugate transpose of a matrix or vector. As expressed in Equation (7), the noise is modeled as a zero-mean Gaussian noise with inverse variance β 0 . By employing a Gamma prior on the inverse variance β 0 and then marginal integrating out the hyper-parameter β 0 , the noise structure in Equation (7) is not changed. After integrating out the hyperparameter β 0 , the posterior distribution function of {a i } i=1,L is transformed from a multi-parameter Gaussian distribution into a multi-parameter Student-t distribution. Compared with the Gaussian distribution, the Student-t distribution is a heavy-tailed distribution that can capture the data trend of variables accurately and eliminate the outlier interference.

Marginal Likelihood Maximization
The unknown hyperparameters β in Equations (12) and (13) need to be estimated. Here, a strategy for seeking hyperparameters β to maximize the logarithm joint marginal likelihood is employed, i.e., When the iterative method is used to estimate the hyperparameters β directly, it requires an O(M 3 ) operation per iteration, because a matrix with a dimension of M × M needs to be inverted in solving Equation (13). In order to promote sparsity and to decrease the computation amount, a fast marginal likelihood maximization algorithm is derived. The C i in Equation (15) is decomposed as where C i,−j is the result of C i , excluding the contribution of β j . In order to simplify Equation (14), the determinant and inverse matrix of C i are rewritten as Substituting Equations (17) and (18) into Equation (14), L(β) is rewritten as , q i,j and g i,j are respectively defined as At the maximum value of L(β), the differential form of the joint marginal likelihood with respect to β j should be zero, i.e., To simplify the zero-finding of (21), β j << s i,j is assumed (typically s i,j > 20β j [19]). Then the first and second derivatives of L(β) with respect to β j are can be solved as L(β) has a maximum value, where β j is solved as Otherwise, β j → ∞. In fact, it can be easier if we maintain and update values of the following variables during iteration, where d j is the j-th basis function, and D is the basis function set contained in the present calculation. In such a way, the number of basis functions included in the iteration is relatively small, so the computation amount is greatly reduced and a faster algorithm can be achieved. The computation complexity of the MMV-BCS algorithm is about O(LMK 2 ) [19].

Flow for MMV-BCS Algorithm
The procedure for the MMV-BCS algorithm is summarized in Table 2. Table 2. Procedure of the MMV-BCS algorithm. Update s i,j , q i,j and g i,j according to (26) 13. end while Step 4, the hyperparameter β j calculated in each iteration should result in the greatest increase in the joint marginal likelihood L(β) in (19), which can yield a faster convergence. The updates of {Σ i } i=1,L , {μ i } i=1,L , S i,j , Q i,j and G i in the add, re-estimate and delete operations are given in the Appendix A. In order to realize the calculation of a complex spectrum, and the matrix power in all the equations are changed into the Hadamard product of the original matrix and its conjugate matrix, i.e., . Note that β 0 can be estimated using the maximum margin likelihood method, however, due to the compressive sensing reconstruction being a under-determined problem, the estimation of β 0 is unreliable at early iterations and an inappropriate estimate value of the hyperparameter β 0 could affect the algorithm accuracy significantly. Therefore, the value of the hyperparameter β 0 here is fixed during iterations, i.e., β 0 = 10 2 /std(abs(y)) 2 [17,19]. In order to make the mean of the Gamma prior Γ(β 0 |a, b) consistent with the fixed value of β 0 and to promote sparsity of priors, based on experiments, the parameters a and b in the Gamma prior Γ(β 0 |a, b) are set a = 10 2 /std(abs(y)) 2 and b = 1 respectively [19,21].
The supraharmonic frequency is obtained as The mean values {μ i } i=1,L are assigned to the estimated values {â i } i=1,L . Since the estimated values of {â i } i=1,L are complex, the magnitude matrix of supraharmonics are obtained by Note that as the phase of supraharmonic is independent of the fundamental [26], the phase accuracy of the supraharmonic analysis is not considered in this paper.
The algorithm does not need to know the sparsity of {a i } i=1,L in advance, which, in fact, is done automatically during the iteration. This can be considered as an advantage over the greedy and convex compressive sensing algorithm. In addition, because the algorithm does not require the orthogonality between the basis functions in any subset of the sensing matrix D, there is no limit on the interpolation factor F, and frequency resolution can hence be improved sufficiently.
As given in [10,13], the minimum frequency interval between two supraharmonic components that the MMV-OMP algorithm can reconstruct accurately is about 1.5 times the DFT frequency resolution ∆ f . When the frequency of a supraharmonic component is not integer multiples of the improved frequency resolution ∆ f , it will be reconstructed on the closest frequency grid by the MMV-OMP algorithm, and hence a frequency estimation error, within ±(∆ f /2) [10], is caused in the MMV-OMP algorithm.
In order to overcome the above issue of the MMV-OMP algorithm, the MMV-BCS algorithm employs a reconstruction technique, which allows representation of the component not aligned to the integer multiples of ∆ f . The magnitude decomposition diagram for the supraharmonic components that are not on the spectrum grid is shown in Figure 2.
In the analysis, it is found that the sum of the magnitude of the two spectrum grids is approximately equal to the true magnitude of the supraharmonic component. Furthermore, the magnitude of the two spectrum grids are inversely proportional to their distance to the true frequency. Therefore, the magnitude and frequency of the supraharmonic component can be estimated by the following two equations, i.e., where ∆ f is the improved frequency resolution, Mag m and Mag m+1 are the corresponding magnitudes of the m-th grid and (m + 1)-th grid in the high-resolution average spectrum, Mag and f are the estimated magnitude and frequency of the true supraharmonic component, that is not on the spectrum grid.

Numerical Verification
It is found in [13] that the PV inverter emits supraharmonics at the switching frequency of around 10 kHz and its integer multiples. Besides, the magnitude of each supraharmonic component has a periodic fluctuation. In order to check the performance of the MMV-BCS algorithm proposed in this paper, two simulation models, a supraharmonic stationary model and a supraharmonic magnitude-modulated model, are constructed similar to the supraharmonic spectrum characteristic of the PV inverter.

Supraharmonic Stationary Model
The supraharmonic stationary model is set as the form in Equation (1). As listed in Table 3, the signal contains 6 supraharmonic components, ranging from 10 kHz to 150 kHz. The initial phase of each supraharmonic component is set randomly in [0, 2π). The sampling frequency is set to 500 kHz and then a 0.5 ms rectangular window is applied to the 200 ms sampled data to generate 400 data blocks. All the data blocks are transformed to the spectrum vectors by DFT with a 2 kHz frequency resolution. In order to improve the frequency resolution of the 400 sets of DFT spectrum by 10 times and compare the results with the MMV-OMP algorithm [13], the interpolation factor F is set to 10. According to Equation (7), a MMV model is constructed. Before high-resolution analysis for this MMV model, some parameters should be initialized, i.e., the noise standard deviation σ is set to 0.02 a.u. and initial values of a and b in the Gamma prior Γ(β 0 |a, b) are set to a = 10 2 /std(abs(y)) 2 and b = 1 respectively, such that the Gamma prior will promote sparsity [19], and the convergence threshold ε is set to 10 −8 .
By applying the proposed MMV-BCS algorithm, 400 sets of supraharmonic spectrum vectors are reconstructed with an enhancement frequency resolution of ∆ f =200 Hz, 10 times better than the DFT analysis. The frequency and magnitude of each supraharmonic are solved following (28) and (29). The average spectrum of the 400 sets of supraharmonic spectra is shown in Figure 3. From Figure 3a, it can be seen that the spectrum represented by DFT has a considerable spectrum leakage. The spectrum leakage can be mitigated significantly using the MMV-OMP algorithm and MMV-BCS algorithm and the frequency resolution is enhanced by an order of magnitude. As discussed in Section 2, the MMV-OMP algorithm cannot recognize the two supraharmonic components of 10 kHz and 10.2 kHz. Instead, the two supraharmonic components are combined into a new magnitude-modulated signal with a 10 kHz frequency, as shown in Figure 3c. The magnitude fluctuates between 0.5 a.u. and 1.5 a.u., and the mean magnitude is 1.060 a.u. The MMV-BCS algorithm proposed in this paper can well reconstruct the two supraharmonic components of 10 kHz and 10.2 kHz.
The 70.2 kHz and 90.2 kHz supraharmonic components are integer multiples of the reconstruction frequency resolution ∆ f = 200 Hz, and both the MMV-OMP algorithm and the MMV-BCS algorithm can reconstruct them on grids with high accuracy. However, for the 29.9 kHz and 50.1 kHz supraharmonic components, which are not integer multiples of the reconstruction frequency resolution ∆ f = 200 Hz, the MMV-OMP algorithm reconstructs them on a 30.0 kHz grid and 50.2 kHz grid, respectively, with a 100 Hz frequency estimation error, and the mean magnitudes are 0.597 a.u. and 0.796 a.u., respectively. The MMV-BCS algorithm can represent the two components on the nearest two grids with one ∆ f : the 29.9 kHz supraharmonic component is reconstructed on 29.8 kHz and 30.0 kHz grids, and the 50.1 kHz supraharmonic component is reconstructed on 50.0 kHz and 50.2 kHz grids.
As we know, in the spectrum obtained by the MMV-BCS algorithm, the 29.8 kHz and 30.0 kHz grids belong to the 29.9 kHz supraharmonic component, and the 50.0 kHz and 50.2 kHz grids belong to the 50.1 kHz supraharmonic component. The four grids are processed by Equations (30) and (31), and the estimated frequencies and magnitudes for the 29.9 kHz and 50.1 kHz supraharmonic components are shown in Figure 3b. From Figure 3b, it can be seen that the estimated frequencies for the 29.9 kHz and 50.1 kHz supraharmonic components are the same as the set values, and the mean magnitudes for them are 0.603 a.u. and 0.803 a.u., respectively. Compared with the MMV-OMP algorithm, the frequency estimation errors are improved by a factor of 100.
The estimation error is evaluated by ∆y = |y −ŷ|, y is the set value,ŷ is the estimated average value. The error bar obtained by the MMV-OMP and MMV-BCS algorithm is shown in Figure 3d. It can be seen from Figure 3d that the frequency estimation error obtained by the MMV-OMP algorithm is approximately within ±(∆ f /2) [10], and for the supraharmonic components that the MMV-OMP algorithm cannot recognize, the frequency estimation error reaches to one ∆ f at least, because the 10.2 kHz supraharmonic component is combined into a new magnitude-modulated signal with the 10 kHz supraharmonic component. For the MMV-BCS algorithm proposed in this paper, the frequency estimation error for all six supraharmonic components is zero. As for the magnitude estimation error, the MMV-BCS algorithm is also much better than the MMV-OMP algorithm. When the frequency interval between the 10 kHz supraharmonic component and the 10.2 kHz supraharmonic component is only one ∆ f , the MMV-BCS algorithm can also reconstruct them, which improves the accuracy by a factor of 15 when the interpolation factor is set to 10.
As discussed above, when one supraharmonic component is not at the high-resolution spectrum grid, the supraharmonic component is reconstructed on the nearest two spectrum grids by the MMV-BCS algorithm with one frequency interval ∆ f . However, when the frequencies of two supraharmonic components are integer multiples of the improved frequency resolution ∆ f and with just one frequency interval ∆ f under a certain interpolation factor, the MMV-BCS algorithm can also reconstruct them with high accuracy. How to distinguish the two types of spectrum grids is important. Here, the normalized inner product (NIP)-based decision criterion is introduced. The NIP coefficient is defined as where a = [a n,1 , a n,2 , · · · , a n,L ] and b = [a n+1,1 , a n+1,2 , · · · , a n+1,L ] are the spectrum vectors on the left and right grids respectively with one frequency interval ∆ f . When the NIP coefficient is bigger than 0.5, the two spectrum vectors belong to one supraharmonic component, whose frequency is not integer multiples of the improved frequency resolution ∆ f . Otherwise, the two spectrum vectors belong to two independent supraharmonic components that are on the spectrum grids with just one frequency interval ∆ f . Using (32), it is found in this example that, the NIP coefficient between the 10 kHz spectrum vector and the 10.2 kHz spectrum vector is 0.0152. The NIP coefficient between the 29.8 kHz spectrum vector and the 30 kHz spectrum vector is 0.9963, and the NIP coefficient between the 50 kHz spectrum vector and the 50.2 kHz spectrum vector is 0.9978. The results support the above NIP criterion well.
From the above discussion, under the interpolation factor F set to 10, the MMV-BCS algorithm can improve the frequency resolution by an order-of-magnitude compared with DFT, and can mitigate the spectrum leakage significantly. The minimum frequency interval between two supraharmonic components that the MMV-BCS algorithm can reconstruct accurately is reduced from 1.5∆ f to just one frequency resolution interval ∆ f , i.e., 200 Hz. The frequency estimation error obtained by the MMV-BCS algorithm is zero and is independent of the supraharmonic component frequency location. The magnitude error obtained by the MMV-BCS algorithm is also much better than that of the MMV-OMP algorithm.
In order to further check the performance of the MMV-BCS algorithm, the reconstruction errors versus the number of spectrum vectors L and noise standard deviation δ at different b values are shown in Figure 4. Note the calculation is made with a fixed frequency of 10 kHz. The reconstruction error is evaluated by y −ŷ 2 / y , y is the set vector,ŷ is the reconstructed vector obtained by the MMV-BCS algorithm.  From Figure 4a, when the number of spectrum vectors increases, the reconstruction error is stable, and the reconstruction error is proportional to b. It can be seen from Figure 4b that the reconstruction error increases with the noise standard deviation δ. When the value of b is set below 10, the reconstruction error will not be affected by the parameter b.
As discussed, the MMV-OMP algorithm requires the orthogonality between the basis functions in any subset of the sensing matrix D. As a consequence, increasing the interpolation factor to improve the frequency resolution will lead to non-orthogonality and hence an incorrect signal reconstruction. Therefore the interpolation factor F cannot be much larger than an order of 10. To further verify the frequency resolution improvement capability of the MMV-BCS algorithm, the stationary model is modified correspondingly. As the interpolation factor F increases, the frequency of the 10 kHz suparharmonic component is changed to 10 kHz plus ∆ f under different interpolation factor F. Therefore, the suparharmonic components are always on the grid in the high-resolution spectrum. The frequency of the 50.1 kHz suparharmonic component is randomly assigned in (50 kHz, 50 kHz + ∆ f ) under different interpolation factor F, therefore the suparharmonic component is off the grid in the high-resolution spectrum. The reconstruction errors, estimation errors are shown in Figure 5.
From Figure 5a, with the interpolation factor F increasing to 100, the on-grid reconstruction error for the 10 kHz plus ∆ f suparharmonic component is basically unchanged. For the off-grid suparharmonic component randomly in (50 kHz, 50 kHz + ∆ f ) in the high-resolution spectrum, the reconstruction error decreases as the interpolation factor F increases. The off grid reconstruction error is slightly larger than the on grid reconstruction error. This proves that by increasing the interpolation factor F to 100, the MMV-BCS algorithm can improve the frequency resolution by a factor of 100 compared to the DFT, and by an order of magnitude compared to the MMV-OMP algorithm. The magnitude estimation errors and the frequency estimation errors for the off-grid suparharmonic component randomly in (50 kHz, 50 kHz + ∆ f ) are shown in Figure 5b. The magnitude estimation errors and the frequency estimation errors are relatively small, i.e., the maximum magnitude estimation error is 0.0112 a.u., the minimum magnitude estimation error is zero. For the frequency estimation error, the maximum error is 6 Hz, the minimum error is 0 Hz, which is about 20 times better than the MMV-OMP algorithm.

Supraharmonic Magnitude-Modulated Model
In this subsection, a magnitude-modulated model is constructed with six supraharmonic components, written as where A sh , f sh and θ sh are the magnitude, frequency and initial phase of each supraharmonic component, respectively. k m is the magnitude modulation factor and f m is the modulation frequency. In this model, the frequencies of the six supraharmonic components are 10 kHz, 20.2 kHz, 29.8 kHz, 50.2 kHz, 70.2 kHz and 90.2 kHz, and the magnitude of each supraharmonic component corresponds to the magnitude listed in Table 3. The modulation frequency is set to f m = 300 Hz and the magnitude modulation factor is set to k m = 0.1. To simulate the intermittent characteristics of a PV inverter, the 50.2 kHz supraharmonic component outputs only one quarter of the time in 200 ms.
As normal, the 200 ms signal is divided into 400 small data blocks, and then all the data blocks are transformed into spectrum vectors using DFT with the frequency resolution of 2 kHz.
The interpolation factor F is set to 10. The noise standard deviation δ is set to 0.02 a.u., and the parameters a and b are set as a = 10 2 /std(abs(y)) 2 , b = 1. The convergence threshold ε is set to 10 −8 . The MMV-BCS algorithm is applied to enhance the frequency resolution for the 400 sets of spectrum vectors and to mitigate spectrum leakage. The reconstructed frequency and magnitude of each supraharmonic component are calculated by (28) and (29). The reconstructed 3-D spectrogram and magnitude fluctuating characteristics of three supraharmonic components with time are shown in Figure 6.   Figure 6c. Although the frequency resolution is 5 Hz, it is impossible to reflect the fluctuation characteristics of the supraharmonic components, and there is also a serious spectrum leakage. Since the spectrum calculates the average value of 200 ms, for the 50.2 kHz supraharmonic component, the magnitude is a quarter of the set magnitude in Figure 6b, only about 0.2 a.u. Obviously, if DFT analysis is used, it could yield a misleading supraharmonic amplitude estimation when the signal is not continuous. This could cause unexpected interference and malfunction for the electrical device if the supraharmonic source is not located and filtered. The MMV-BCS algorithm can address such issues successfully.

Applications
Taking a typical supraharmonic emission source-PV inverter as an example, the MMV-BCS algorithm proposed in this paper is used to perform supraharmonic analysis on the output signal. The switching frequency of the PV inverter is about 10 kHz, and the PV inverter emits supraharmonic components at the switching frequency and its integer multiples. The current of the PV inverter is recorded by a YOKOGAVA DL850E recorder and a CAE3N current clamp with a sampling frequency of 500 kHz.
To avoid contributions from the fundamental component of the current signal in the PV inverter, the intercepted 200 ms sampling signal is filtered to remove the fundamental component with a Chebyshev II high-pass digital filter recommended by IEC 61000-4-30, as shown in Figure 7, which is the magnitude frequency response characteristic curve of the Chebyshev II high-pass digital filter. The cut-off frequency is set to 2 kHz, and the stop-band attenuation is −60 dB. Using this filter, 99.9 % of the fundamental component is filtered out. The filtered signal is divided into 400 small data blocks with 0.5 ms, and then an array of 400 sets of spectrum is obtained using the DFT. In order to balance frequency resolution and calculation accuracy, the interpolation factor is set to 20. Then the DFT frequency vectors are analyzed using the MMV-BCS algorithm proposed in this paper. The spectrogram of the PV inverter is shown in Figure 8.  Figure 8a is the 3-D spectrogram of the supraharmonics obtained by the MMV-BCS algorithm with a 100 Hz frequency resolution. The spectrum leakage can be mitigated significantly. It can be seen that the PV inverter emits supraharmonics at the switching frequency of 10 kHz and its integer multiples, and the magnitude of each supraharmonic component still has a common periodic fluctuating characteristic, which is a resonance phenomenon caused by an LC filter at the output unit of the PV inverter. From the 3-D spectrogram of the displayed 200 ms, it can be seen that the modulation frequency is about 300 Hz, and the magnitude modulation factor is about 0.1. Figure 8b is the average spectrum of 400 sets of high-resolution spectra. The frequency and magnitude of the main supraharmonic components obtained by the MMV-BCS algorithm are listed in Table 4. It can be seen that the PV inverter emits supraharmonics at the switching frequency of 10 kHz and its integer multiples. Furthermore, the magnitude of the supraharmonic components gradually decreases as the frequency increases. Besides that, there are also some side-band components around each supraharmonic component.
Here, we use another supraharmonic emission source, a wireless electric vehicle (WEV) charger. The input current of the WEV charger is recorded by a HIOKI MR8875 recorder and CT9692 current clamp with a sampling frequency of 100 kHz. Similarly, a high-pass filter is used to remove the fundamental and harmonic components below 2 kHz, and the signal after filtering is divided into 100 small data blocks with 2 ms data length. All the data blocks are transformed to a spectrum domain by DFT with a frequency resolution of 500 Hz. The interpolation factor F is set to 20, and hence the frequency resolution is improved to 25 Hz. The parameters of the MMV-BCS algorithm are set the same as above PV example. By applying the MMV-BCS algorithm to the multiple DFT spectrum vectors, 100 sets of high-resolution supraharmonic spectrum vectors are reconstructed with an enhancement frequency resolution of 25 Hz. The spectrogram and the average spectrum are shown in Figure 9.  Figure 9a is the 3-D spectrogram of the supraharmonics obtained by the MMV-BCS algorithm with a 25 Hz frequency resolution, and it can be seen that the spectrum leakage can be suppressed significantly compared to the DFT algorithm. There are two main supraharmonic components, and the magnitude of each supraharmonic component has a fluctuating characteristic. Figure 9b is the average spectrum of the reconstructed 100 sets of high-resolution spectra. Two major peaks are obtained: at 15.97 kHz, the magnitude of the supraharmonic component is 0.7134 A, and at 31.93 kHz, the magnitude of the supraharmonic component is 0.7925 A. Some adjacent small supraharmonic components are also successfully detected.

Conclusions
In this paper, using the RVM and a hierarchical Bayesian model, a new supraharmonic high-resolution reconstruction algorithm, MMV-BCS, is proposed. The new algorithm can simultaneously achieve high-resolution reconstruction for multiple spectrum vectors.
Experiments and applications show that the MMV-BCS algorithm can improve the frequency resolution by two order-of-magnitude compared with DFT and by one order-of-magnitude compared with the MMV-OMP algorithm without increasing the observation time. The MMV-BCS algorithm has zero frequency estimation error when the supraharmonics are on the high-resolution spectrum grid. While the supraharmonics are off-grid in the high-resolution spectrum grid, the algorithm can also reconstruct them with high accuracy. The frequency estimation error is about 20 times less than the frequency estimation error of the MMV-OMP algorithm. The minimum frequency interval between two supraharmonic components that the MMV-BCS algorithm can reconstruct with high accuracy is only one frequency resolution interval, which is at least 15 times better than the MMV-OMP algorithm. Besides, tests show the MMV-BCS algorithm can also handle the analysis of a magnitude-modulated supraharmonic signal.
The MMV-BCS algorithm can measure supraharmonic accurately and shows good application prospects in troubleshooting, emission limit setting and mitigation of supraharmonics. Since the iterative process of the algorithm is still complicated, it is the next research direction to study a more efficient and fast Bayesian compressed sensing reconstruct algorithm. Q i,l = Q i,l − µ i,j (d H l e i,j ) (A5) where Σ i,(jj) = (β j + S i,j ) −1 is the j th diagonal element of Σ i , µ i,j = Σ i,(jj) Q i,j and e i,j = d j − DΣ i D H d j