Portal de Periódicos da CAPES

Spectrum estimation using frequency shifting and decimation

2019; Institution of Engineering and Technology; Volume: 14; Issue: 3 Linguagem: Inglês

10.1049/iet-spr.2019.0306

1751-9683

Raymundo Albert, Cecilia G. Galarza,

Blind Source Separation Techniques

IET Signal ProcessingVolume 14, Issue 3 p. 134-141 Research ArticleFree Access Spectrum estimation using frequency shifting and decimation Raymundo Albert, Corresponding Author Raymundo Albert araymundo@fi.uba.ar Facultad de Ingeniería, Universidad de Buenos Aires, Paseo Colón 850, Buenos Aires, Argentina Centro de Simulación Computacional, CONICET, Godoy Cruz 2390, Buenos Aires, ArgentinaSearch for more papers by this authorCecilia Galarza, Cecilia Galarza Facultad de Ingeniería, Universidad de Buenos Aires, Paseo Colón 850, Buenos Aires, Argentina Centro de Simulación Computacional, CONICET, Godoy Cruz 2390, Buenos Aires, ArgentinaSearch for more papers by this author Raymundo Albert, Corresponding Author Raymundo Albert araymundo@fi.uba.ar Facultad de Ingeniería, Universidad de Buenos Aires, Paseo Colón 850, Buenos Aires, Argentina Centro de Simulación Computacional, CONICET, Godoy Cruz 2390, Buenos Aires, ArgentinaSearch for more papers by this authorCecilia Galarza, Cecilia Galarza Facultad de Ingeniería, Universidad de Buenos Aires, Paseo Colón 850, Buenos Aires, Argentina Centro de Simulación Computacional, CONICET, Godoy Cruz 2390, Buenos Aires, ArgentinaSearch for more papers by this author First published: 01 May 2020 https://doi.org/10.1049/iet-spr.2019.0306AboutSectionsPDF ToolsRequest permissionExport citationAdd to favoritesTrack citation ShareShare Give accessShare full text accessShare full-text accessPlease review our Terms and Conditions of Use and check box below to share full-text version of article.I have read and accept the Wiley Online Library Terms and Conditions of UseShareable LinkUse the link below to share a full-text version of this article with your friends and colleagues. Learn more.Copy URL Share a linkShare onFacebookTwitterLinkedInRedditWechat Abstract Parametric spectral estimation techniques are widely used to estimate the parameters of sums of complex sinusoids corrupted by noise. In this work, the authors show that the numerical stability of the estimated frequencies not only depends on the size of the amplitudes associated with the real frequencies, but also to the distance among frequencies. Therefore, for closely spaced frequencies, the estimates are vulnerable to large deviate from their true values. To overcome this problem, they propose a strategy to artificially increase the frequency separation by downsampling the baseband equivalent of the noisy signal before applying a spectral estimation technique. This methodology significantly improves the estimation performance especially in the low signal to noise ratio regime. The performance of the technique is assessed in terms of the root mean square error and it is compared to results obtained in previous publications. 1 Introduction Problems dealing with a mixture of closely spaced damped sinusoidal signals in a noisy environment are regularly observed in engineering. This is a problem of major significance in applications such as seismic exploration, communication, radar, and sonar, to name but a few. In those cases, the signal of relevance contains unknown parameters such as amplitude, phase, and frequency of each mode. This path led in the past to developing parametric spectral estimation techniques, in particular, the family of subspace methods. These algorithms generate frequency component estimates of a given signal based on the decomposition of an observation vector space into two subspaces: one associated with the signal and the other with the noise [1, 2]. A major concern when developing these techniques was to cope with noise and model uncertainties, in particular, the unknown order of the model that corresponds to the dimension of the signal subspace. Different strategies for obtaining this magnitude have been worked in the past [3, 4]. Lately, novel optimisation techniques have been proposed to estimate the model order. This method, which relies on Kronecker's theorem for Hankel operators, was used to formulate an appropriate non-linear least squares problem in terms of a rank constraint on the Hankel matrix associated with the observed signal [5]. Given the model order, the frequencies making up the observed signal are estimated on a second step. A largely used technique, known as ESPRIT (Estimation of Signal Parameters by Rotational Invariance Technique) [6], exploits the structure of the signal space and it formulates a generalised eigenvalue problem to obtain the frequency estimates. An alternative approach, known as the Matrix Pencil Method (MPM), has its roots in the Vandermonde decomposition of the Hankel matrix built from the observed signal [7]. A more recent approach based on the Hankel matrix associated with the observation was proposed in [8], where the parameter estimation is performed by minimising the nuclear norm of the Hankel matrix. A different path was followed in [9], where the authors considered sums of undamped complex exponentials and they exploited a sparse description for them. In this case, using a previously defined set of discrete frequencies, compressed sensing techniques are used to estimate the unknown parameters. An improvement to this approach was presented in [10, 11]. There, the solution is obtained by minimising the atomic norm of the estimated signal. It was shown that, in this case, the exact recovery of the observed signal is guaranteed when the frequencies are adequately separated [12]. Unfortunately, these methods are not suitable for a sum of damped complex exponentials. The numerical conditioning of the estimates when recovering a sum of undamped complex exponentials was analysed in [13]. In this work, an upper bound of the condition number of the related Vandermonde matrix was provided, and a decimative approach for a single cluster of frequencies was proposed. As it was pointed out in [14] the performance of spectrum estimation methods such as ESPRIT, MPM, etc., degrades for closely spaced complex exponentials. To overcome this problem, the frequency separation was artificially increased by decimating the signal and then applying ESPRIT. In [15], a decimation approach is applied, and a total least squares problem is solved. However, in [14, 15], the downsampling factor is chosen to make sure that no aliasing is introduced. To overcome this restriction, the original signal was oversampled. On a different angle, the authors of [16] analyse the spectral content of geophysical time series by frequency shifting and later filtering the observed signal. By doing so, they were able to identify the resonant frequencies contained in the observed signal. Nevertheless, accurate estimation of spectral parameters for sums of damped complex exponentials is performed via the solution of generalised eigenvalue problems. These eigenvalue problems are usually ill-conditioned and non-square. Solving generalised eigenvalue problems of these characteristics is a challenge from an algorithmic perspective. This difficulty turns to be harder when the number of frequencies in the original mixture is very large and the frequencies are clustered in small regions of the complex plane [17, 18]. Our focus is to analyse the behaviour of the frequency estimation step when the observed data is subject to small perturbations. In particular, we study the numerical stability of the generalised eigenvalue problem that is built from the Hankel matrix associated to the observed signal. A related problem has been addressed before [19, 20]. In particular, it has been observed that the sensitivity of each eigenvalue is inversely proportional to the inner product of the left and right generalised eigenvectors weighted by the Hankel matrix. This inner product is proportional to the amplitude of the damped oscillation corresponding to the analysed eigenvalue. In this paper, we extend these results and we show that the first-order approximation of the perturbed eigenvalues not only depends on the amplitudes associated with the eigenvalues, but also on the distance among eigenvalues. Therefore, the eigenvalues that are close to each other are prone to exhibit large deviations from their actual values, even when the observed signal is only lightly perturbed. Similar results were obtained for a different context by the authors of [21]. Based on the analytical results, we propose a strategy to pre-processed the observed signal using multi-rate processing. This has the ability to accommodate the complex frequencies for more stable computation. The idea is validated with numerical experiments and compared to the results obtained in previous publications. The paper is organised as follows. In Section 2, we make a brief introduction to subspace methods, in particular, ESPRIT and MPM. In Section 3, we study the numerical stability of the non-square generalised eigenvalue problem and we obtain a new bound for the first-order derivative of the generalised eigenvalue. A method to obtain better numerical stability and improved frequency estimation is presented in Section 4. Finally, numerical experiments are presented in Section 5 and conclusions are given Section 6. Throughout the paper, we use the standard notation: the lower case for scalars, boldface lower case for vectors, uppercase bold face for matrices. The k th entry of a vector is denoted . Given a matrix , we denote its transpose, Hermitian, and Moore–Penrose pseudo-inverse as , , , respectively. refers to the element in row i and column j of a matrix . refers to the Frobenius norm of the matrix , is reserved for the induced 2-norm. Finally, the notation is used for the identity matrix. 2 Spectrum estimation methods 2.1 Background Spectrum estimation techniques are usually treated as a two-step procedure. In the first step, the model order is estimated and the observed signal is cleaned up from the perturbing noise. In the second step, the values for the complex frequencies and their amplitudes are obtained [2]. This paper is focused on the analysis of this second step and how sensitive the estimated complex frequencies are when the observed signal is subject to small perturbations. To begin with, the description of the traditional techniques, consider a noiseless discrete-time signal composed by a sum of damped complex exponential as follows: (1)where n is the number of exponential terms, is a complex resonant frequency, and the amplitude associated to it. We assume that are all mutually distinct and for all i. In the first approach, we assume that n is known. This assumption is revised in Section 2.2. Given , define the following -complex Hankel matrix (2)and the complex vector . The solution to the linear system of equations (3)is a complex vector whose entries are the coefficients of the polynomial . This simple setup is the starting point for classical Prony's method [22]. Unfortunately, is an ill-conditioned matrix, and solving (3) for could be very sensitive to noise in the data. Moreover, obtaining as the roots of is also a challenging task to perform when (1) is perturbed by noise. More robust approaches have been proposed by exploiting the Vandermonde decomposition of . Let us define the -Vandermonde matrix as (4)Then, the Vandermonde decomposition of is (5)where . Consider the singular value decomposition of , where the columns of span the range of and they are identified with the signal subspace in . Since are distinct and , , and from (5), we see that the columns of also span the signal subspace. Then, there exists an invertible matrix such that (6)Define now the following complex matrices: (7)Matrices and are obtained from by deleting its first and last rows, respectively. Analogously, we define and from . Using these definitions, it is easy to verify the rotational invariance property that says that (8)where . By replacing (6) in (7), we have that and . Therefore, we obtain from (8) that (9)This equation leads to a generalised eigenvalue problem that obtains the unknown frequencies and characterises the algorithm known as ESPRIT [6]. A related technique exploits the relationship between and , which are defined in accordance with (7) using the Hankel matrix . Notice that . Then, using (5) and (8), we obtain (10)Now for any , we have that (11)Since , , and are all full rank matrices, the unknown diagonal entries of are the solutions to the generalised eigenvalue problem (12)Equation (12) is known as the Matrix Pencil Method [7] because the set constitutes a matrix pencil. Notice that in the noiseless case, when the model order n is known, (9) and (12) are equivalent, and MPM and ESPRIT obtain the same result. General matrix pencils do not always have eigenstructures that satisfy . However, the pencil defined in (9) or (12) will always have a solution as long as , are all distinct and , . This becomes clear by observing (11). Notice that for all , , , where is the i th unit vector. Let be the solution to an under-determined system of equations . Clearly, satisfies the matrix pencil in (12). Similarly, we find the set of left eigenvectors by solving . 2.2 Signals with noise When working with noisy signals, we deal with the following model: (13)where is given in (1) and is an additive noise term. We assume that is a circularly symmetric complex Gaussian process with identically distributed zero-mean uncorrelated real and imaginary parts whose variance is . Generally, the order model n is unknown and it needs to be estimated from the collected data before solving the estimation problem. Even when n is known, the observed data needs to be processed to mitigate the effect of on . For that, several procedures have been proposed in the past to denoise before proceeding to the spectral estimation step. When n is unknown, a first Hankel matrix is built using . A first simple approach searches for the dominant singular values of and defines n as their number [23]. Then, to perform the spectral estimation procedure, is replaced by its truncation to the first n singular directions. Unfortunately, the resulting approximation does not preserve the Hankel structure. However, Kronecker's theorem [24] states that there is a one-to-one correspondence between a linear combination of n complex exponentials and a Hankel matrix with rank n. Then, accuracy is lost when a simple truncation of the Hankel matrix is performed. This problem was observed in [18, 25] where an iterative procedure that is performed in two steps was proposed. For each iteration, the first step consists of performing the matrix truncation to the first n singular directions of . In the second step, the truncated matrix is modified by forcing a Hankel structure on it. Both steps are repeated until a stopping criterion is satisfied. An alternative solution performs a search on , the space of Hankel matrices [26]. For a given n, it is possible to formulate the following optimisation problem: (14)In [27, 28], this problem was tackled by solving a total least square problem. Recently, a new denoising technique was proposed in [5] by revising Kronecker's theorem, which states that the Hankel matrix generated by a signal has a rank n, if and only if coincides with a function that is a linear combination of n exponential functions. Let be the Hankel matrix obtained from the sequence . Now, the following optimisation problem returns an approximate signal that satisfies: (15)Due to the constraint on the Hankel matrix, this is a non-convex problem that may be reformulated as (16)where is an indicator function for matrices that are defined such that if and 0 otherwise. In [5], the authors propose to solve the problem using the Alternative Direction Method of Multiplier (ADMM) [29]. This is an efficient technique for solving non-convex optimisation problems. The formulation of ADDM associates an augmented Lagrangian associated with the original problem (16) (17)where is the Lagrange multiplier matrix and is a constant penalty parameter. The solution is obtained performing the following steps until convergence: (18)The final estimate is obtained by averaging the anti-diagonal terms of the matrix . We notice that the goal of these procedures is always to obtain an approximated sequence and its associated Hankel matrix . Then, for estimating , , a second step is performed by solving (9) or (12). Even when the first step that mitigates the effect of noise is successful, the computation of the generalised eigenvalues in (9) or (12) may be unstable. This problem is particularly acute when the number of oscillation modes is very large or the complex frequencies are close to each other in the plane. 3 Numerical stability 3.1 Derivatives of generalised eigenvalues Suppose that is the i th generalised eigenvalue of the pair with associated right and left eigenvectors and , respectively. We assume further that and are full rank and are all distinct. By definition (19)Now, consider smooth variations of the matrices and with a real parameter in a neighbourhood of the origin, i.e. . Assume that there exist differentiable functions , , , so that (20)where , , , . Theorem 1.Given that satisfy (20) for all , we have that (21)Here, and are the first derivatives with respect to of and . Proof.See the Appendix. □ The eigenvalue will be ill-posed, i.e. small perturbations of and produce large deviations of from its unperturbed value, if is small. This is the case when is close to perpendicular to . The product has been associated with the numerical condition of [30, ch. 7] 3.2 Perturbations on When analysing (12), the perturbed matrices and are obtained from , where is a scalar parameter. In that case, we can formulate the following corollary: Corollary 1.Let and be differentiable functions that represent perturbed versions of the pencil . Assume that for each , the eigenvalues , and the left and right eigenvectors and exist and they satisfy (20). Then (22)where , is a polynomial in z of degree , and . Proof.Applying Theorem 1 to the pencil , we obtain (23)where and are left and right eigenvectors of the pencil .Recalling the decomposition in (10), we have that and satisfy (24)where is the i th unitary vector in . Moreover (25)Since , with , one may choose to satisfy (24). On the other hand, and since is an invertible matrix, (24) implies that , is the i th column of . Then, using the result given in [31], we have that for (26)where is a combination of n - k elements taken from . Given a choice for the left eigenvector, . Then (27)Notice that is the k th coefficient of the polynomial . Re-arranging terms, we have that Then, the Parseval's theorem between the sequence and its discrete-time Fourier transform states that (28)where to compute the integral, we have evaluated in the unit circle, . Replacing this equation in (27), and its result together with (25) into (23), we obtain (29) Corollary 1 shows that when performing spectral estimation, the estimation of is sensitive to but also to . The constant represents the residue associated with and it is related to the observed strength of the complex frequency. As one may suspect, small lead to an associated eigenvalue that can be easily perturbed. On the other hand, when an eigenvalue is located inside a cluster of eigenvalues, may take a large value. Then, the sensitivity to numerical errors of the estimated frequency is affected not only by its strength in the mixture, but also by its proximity to the other complex frequencies. In the next section, we present a technique to cope with this problem. 4 Downsampling the base-band equivalent Consider that in (13) is obtained after uniformly sampling a continuous-time version where . In this case, , where is the sampling period, and . Clearly, the location of in the complex plane may be controlled by judiciously adjusting the sampling time. However, the selection of responds to several constraints when designing the signal acquisition system, and adding additional requirements may not be feasible. Nevertheless, the location of may also be changed by the decimation of the observed signal [32]. Suppose that the observed signal has its spectral content concentrated on , i.e. for (30)According to Section 3, one would like to 'zoomed in' before performing the estimation of to improve the performance of the estimation algorithm. For that, we proceed as follows. Let and define The signal is known as the baseband equivalent of . Consider a scalar Q such that . Then, using Fig. 1, we construct after filtering and decimating . To avoid aliasing, is filtered before decimation with a linear phase finite impulse response filter [33]. Now, we apply the frequency estimation procedure to and shift back the estimated modes. The final structure is shown in Fig. 1. We refer to this estimation strategy as a shift-and-zoom technique. Fig. 1Open in figure viewerPowerPoint Shift-and-zoom scheme for estimating complex modes Using the decimated baseband signal , we construct a Hankel matrix as in (2) and estimate the complex frequencies as in Section 2. Let be the frequency estimates obtained from . Clearly When decimating or downsampling by a factor Q, the complex frequencies move in the complex plane. For , the frequencies that were clustered in may spread apart, diminishing the value of , and making the spectral estimation more accurate. It is important to remark that the low-pass filter included in Fig. 1 modifies the amplitudes related to the complex modes , affecting the value of the derivative of as seen in (22). Now, suppose that , for all , where are disjoints intervals. In this case, are separated in L different clusters in the complex plane. To improve performance, the procedure outlined above could be repeated for each cluster. In particular, since are intervals of unequal lengths, the decimation factor could be different for each l. The complete procedure is summarised in Fig. 2. Fig. 2Open in figure viewerPowerPoint Block diagram for processing the signal compound of L clusters 5 Numerical experiments In this section, we evaluate the performance of the shift-and-zoom procedure outlined in the previous section. For that, we analyse two examples where sums of complex exponentials are simulated for different signal-to-noise ratios (SNRs). In both cases, we estimated the complex modes and compared the results using the procedure presented in [5]. To appraise the performance of each algorithm, we repeat the experiments K times and we compute the root mean square errors (RMSEs) for the estimated frequencies and damping factors (31)where and are the estimations obtained after the k th Monte Carlo run. 5.1 Weak modes As a first example, we consider the sum of four modes described in Table 1. This example is taken from [5]. We simulate using a sampling interval . Table 1. Parameters for the example taken from [5] i 1 2 3 4 −7.68 39.68 40.96 99.84 −0.274 −0.150 0.133 −0.221 0.4 1.2 1.0 0.9 −0.93 −1.55 −0.83 0.07 When applying the shift-and-zoom technique, we consider three disjoint frequency intervals, , , and . Clearly, , , and . Here, we will assume that the cluster locations are known a priori. When this information is not available, a preliminary run of the estimation algorithm can be performed in order to obtain rough estimates of the frequency bands. Notice that all three segments have the same bandwidth. Although this is not a necessary condition for the algorithm, it simplifies its implementation. Using this fact, we use the same low-pass filter and downsampling coefficient for all three frequency segments. In particular, we design an FIR linear phase filter using a rectangular window of 16 coefficients, with bandwidth 15.39 Hz. The downsampling coefficient is . For we compute and for each SNR. The results are shown in Fig. 3 as functions of the SNR. In Fig. 4, we also show the RMSEs of the amplitude and phase. Fig. 3Open in figure viewerPowerPoint RMSE as a function of SNR: dashed lines correspond to [5]; solid lines correspond to shift-and-zoom (a) , (b) Fig. 4Open in figure viewerPowerPoint RMSE as a function of SNR: dashed lines correspond to [5]; solid lines correspond to shift-and-zoom (a) Amplitude RMSE, (b) Phase RMSE For high SNR, shift-and-zoom has a similar performance to [5]. However, for small SNR, the better conditioning of becomes critical, resulting on remarkably better performance of shift-and-zoom when compared to a more traditional approach. However, it is important to notice that when we downsample the signal, there is a tradeoff between the number of samples remaining after decimation and how far apart the complex modes result. To evaluate our procedure, we compute the Cramér Rao bound (CRB) using the expressions in the Appendix [34]. The CRB is a lower bound on the variance of any unbiased estimator. Although we have not proved that the proposed estimators are unbiased, we will use this bound as an indicator of the margin for improvement in the estimation procedure. Notice that the CRB depends on and not on . Therefore, when downsampling , the CRB changes. In Fig. 5, we compared the ratios and for the estimations obtained through shift-and-zoom and through the procedure in [5]. Fig. 5Open in figure viewerPowerPoint Ratio between RMSE and the CRB. Solid lines correspond to shift-and-zoom and dashed lines to [5] (a) Ratio between RMSE and the CRB for , (b) Ratio between RMSE and the CRB for Notice that the downsampling factor has an upper bound determined by the number of samples in the original signal. As was pointed out in [34], the CRB depends not only on the distance among frequencies but also on the number of samples. In fact, it is pointed out that the CRB depends on the product of the number of samples and the separation among frequencies. The length of the observed signal is also significant when estimating the damping factors. This issue becomes critical when the estimated model order is larger than the actual one. In this case, erroneous phantom damping factors could affect the estimation of the actual parameters. A solution in [35] tackles this problem by using . To analyse the effect of the downsampling factor Q, we have computed the residual RMSE between the estimated signal and the true signal for different values of Q. The results are shown in Table 2. The minimum RMSE is obtained for . As we increase Q, the RMSE becomes higher because fewer samples are used. On the other hand, as Q decreases, even though a longer signal record is analysed, we lose the effect of frequency separation. Table 2. Residual RMSE as a function of the oversampling factor Q 1 2 4 7 9 14 15 RMSE 15.27 15.26 7.82 11.25 12.09 14.54 15.82 5.2 Clustered modes In the second example, we consider the sum of 25 damped oscillations with different complex amplitudes each. This example that corresponds to data obtained from magnetic resonance spectroscopy (MRS) was presented in [17]. The sampling interval is . The parameters are described in Table 4 where frequencies, damping factors, phases and amplitudes are provided for each complex mode. The authors grouped the 25 modes in 6 disjoint clusters. In Fig. 6, we have plotted the locations of , . Fig. 6Open in figure viewerPowerPoint Location of the modes described in Table 4 when For the purpose of using the shift-and-zoom procedure, we consider different frequency intervals. Interval limits and decimation factors are found in Table 3. Table 3. Frequency segments and decimation factors , Hz , Hz , Hz −349.45 −390.41 −308.49 4 −130.29 −191.73 −68.85 3 14.95 4.71 25.19 4 140.95 100.61 166.45 5 436.75 416.27 457.23 6 We processed with the shift-and-zoom algorithm and compared the estimations with those obtained with no preprocessing (shifting and decimation). Figs. 7 and 8 show the RMSE of the estimated frequencies and damping factors for the different clusters using both techniques. For a better presentation of the results, the Cluster IV is separated into two. Cluster IV (a) with four different frequencies and Cluster IV (b) with three frequencies. Fig. 7Open in figure viewerPowerPoint RMSE for as a function of SNR for all clusters. Dashed line corresponds to [5]; solid lines correspond to shift-and-zoom (a) Cluster I, (b) Cluster II, (c) Cluster III, (d) Cluster IV (a), (e) Cluster IV (b), (f) Cluster VI Fig. 8Open in figure viewerPowerPoint RMSE for as a function of SNR for all clusters. Dashed line corresponds to [5]; solid lines correspond to shift-and-zoom (a) Cluster I, (b) Cluster II, (c) Cluster III, (d) Cluster IV (a), (e) Cluster IV (b), (f) Cluster VI At low SNR, shift-and-zoom performs better than a more traditional approach. Since shift-and-zoom increases the frequency separation, improving makes the algorithm more resilient to noise disturbances. This observation is even more relevant when the amplitude is small. This is the case of modes to in Cluster IV. In this example, it was critical to decimate instead of directly. Should we have done the latter, the maximum decimation factor would have been in order to avoid aliasing. Such value for Q results in a poor estimation of some modes, in p