Algorithms for Real-Valued Noisy Damped Sinusoid Parameter Estimation

In this article, an improved version of three state-of-the-art interpolated discrete-time Fourier transform (IpDTFT) algorithms for the estimation of the frequency and the damping factor of noisy damped sinusoids is proposed. Sinusoid amplitude and phase are also estimated. Improvement is obtained by both compensating the effect on the estimated parameters of the spectral leakage from the fundamental image component and minimizing the estimator variance through a suitable selection of the interpolation points. The obtained estimates are then further improved by applying a linear signal fit (LSF) algorithm. The accuracies of the proposed algorithms are compared with each other and with the related Cramér–Rao lower bound through computer simulations. The required computational times are also briefly analyzed in order to assess the algorithms usefulness in real-time applications.


I. INTRODUCTION
D AMPED sinusoids carefully model signals encountered in many application domains, such as nuclear magnetic resonance, radar, sonar, geophysics, and biomedical signal processing [1], [2], [3]. Accurate estimates of signal parameters are often of high interest and various time-domain and frequency-domain algorithms have been proposed to that aim. Time-domain algorithms, such as the Prony [4], the MUSIC [5], the Auto Regressive [4], [5], [6], the Matrix Pencil [7], [8], [9], [10], and the Steiglitz-McBride (or iterative Prony) [7], [11] methods, or the nonlinear least-square approach [12], [13] can provide very accurate parameter estimates, but they require a high processing effort. Accurate and computationally efficient damped or pure sinusoid parameter estimators are provided by the so-called interpolated discrete-time Fourier transform (IpDTFT)-based algorithms [1], [2], [3], [7], [14], [15], [16], [17], [18], [19], [20]. In [14], the algorithms proposed by Quinn [17] and by Aboutanios and Mulgrew [18] for real-valued or complexvalued sinusoids, respectively, have been extended to the estimation of the frequency and the damping factor of complex-valued damped sinusoids. It has been shown that the extended Quinn's algorithm represents a linearization of Bertocco's algorithm [3] and a hybrid algorithm based on Quinn's and Aboutanios' estimators is proposed. Moreover, the analytical expressions for the variances of the frequency and the damping factor estimators provided by Quinn's, Bertocco's, and Aboutanios' algorithms were derived [14]. In addition, it has been shown that the estimators provided by Aboutanios' and the hybrid algorithms almost reach the related asymptotic Cramér-Rao lower bounds (CRLBs) when the number of analyzed samples is properly chosen [14]. Duda et al. [2] performed a comprehensive accuracy analysis about the frequency and the damping factor estimates returned by different IpDTFT algorithms applied to realvalued damped sinusoids. They also proposed a novel IpDTFT algorithm that leverages on the approaches proposed by Bertocco et al. [3] and Yoshida et al. [15]. An improved version of that algorithm has been later proposed in [7], in which the detrimental effect due to the spectral leakage from the fundamental image component is compensated through an iterative procedure. The obtained algorithm is more accurate than Aboutanios' one [14], but less accurate than the Steiglitz-McBride and the Matrix Pencil methods. However, the processing effort required by these latter methods is significantly higher.
In this work, at first, the accuracies of Bertocco's [3], Aboutanios' [14], and Duda's [7] algorithms are improved by compensating the spectral leakage from the fundamental image component and selecting the interpolation points that ensure the minimization of the impact of noise on the returned estimates. The obtained algorithms are called the enhanced Bertocco's, Aboutanios', and Duda's damped sinusoid parameters estimation (eBDSPE, eADSPE, and eDDSPE) algorithms, respectively. Moreover, a linearized version of the nonlinear signal fit algorithm [12] is proposed to further improve estimation accuracy while using simple closed-form formulas. The obtained procedure, called the linear signal fit (LSF) algorithm, requires a reduced additional processing effort, while the estimators root-mean-square errors (RMSEs) are very close to the related CRLB.
The remaining of this article is organized as follows. The eBDSPE, eADSPE, and eDDSPE algorithms are presented in Section II. Section III is dedicated to the presentation of the LSF algorithm. The accuracy performances of the proposed algorithms are then analyzed in Section IV. Section V concludes this article.

II. PROPOSED EBDSPE, EADSPE, AND EDDSPE ALGORITHMS
The analyzed discrete-time real-valued noisy damped sinusoid is modeled as where A, α, ν, and φ are the amplitude, the damping factor, the normalized frequency (expressed in bins), and the initial phase, respectively; e(.) is a white Gaussian noise with zero mean and variance σ 2 , and M is the acquisition length. The signal-to-noise ratio is SNR (A 2 /2σ 2 ). The normalized frequency, which represents also the number of acquired signal cycles, can be written as where l is an integer, and it corresponds to the index of the discrete Fourier transform (DFT) spectrum peak while δ is the interbin frequency location (−0.5 ≤ δ < 0.5). After simple calculations, the DTFT of (1) can be expressed as where and E(λ) is the DTFT of the wideband noise e(·). Observe that the second term in the right-hand side represents the contribution of the fundamental image component. The proposed eBDSPE, eADSPE, and eDDSPE algorithms are presented in the following.

A. EBDSPE ALGORITHM
The eBDSPE algorithm is based on the procedure proposed in [3]. It determines a first estimate of the damped sinewave parameters by means of a complex-valued interpolation based on the two largest DFT samples. The estimated values are then exploited to iteratively compensate the effect of the spectral leakage from the fundamental image component on the interpolation points. A last interpolation is finally performed using two DTFT samples located half-bin apart from the estimated frequency, so that the noise contribution on the estimated parameters is minimized [14]. The various steps of the proposed eBDSPE algorithm are detailed in the following.
Step 1: Acquire M samples of the noisy damped sinusoid x(.).
Step 2: Determine l as the normalized frequency of the acquired signal DFT spectrum peak.
Step 3: Determine a first estimate of the damped sinusoid parameters by interpolating the two largest DFT samples Step 4: Compensate the contribution of the image component to the used interpolation points , k = −1, 0, 1.
Step 5: Determine the ith estimate of the damped sinusoid parameters by using the two largest compensated DFT samples as interpolation points Step 6: Compensate the contribution of the image component on two DTFT samples located half-bin apart from the estimated frequencỹ Step 7: Determine the final estimate of the damped sinusoid parameters by interpolating the two compensated DTFT samples

B. EADSPE ALGORITHM
The eADSPE algorithm is based on the procedure proposed in [14] for complex-valued signals, but it compensates the effect of spectral leakage due to the image component on the involved interpolation points. As the eBDSPE algorithm, it determines a first estimate of the damped sinewave parameters by applying a complex-valued interpolation based on the two DFT samples X(l − 0.5) and X(l + 0.5). The estimated parameters are then exploited to iteratively compensate the contribution of the spectral leakage from the fundamental image component on the two used DFT samples. A last interpolation is finally performed using two DTFT samples located half-bin apart from the estimated frequency in order to minimize the impact of wideband noise on the returned estimates [14]. The processing steps that implement the proposed eADSPE algorithm are detailed in the following.
Step 1: Acquire M samples of the noisy damped sinusoid x(.).
Step 2: Determine l as the normalized frequency of the acquired signal DFT spectrum peak.
Step 3: Determine a first estimate of the damped sinusoid parameters by interpolating the two DFT samples X(l − 0.5) and X(l + 0.5) Step 4: Compensate the contribution of the image component to the used interpolation points Step 5: Determine the ith estimate of the damped sinusoid parameters using the compensated DFT samples as interpolation points Step 6: Compensate the contribution of the image component on two DTFT samples located half-bin apart from the estimated frequencỹ Step 7: Determine the final estimate of the damped sinusoid parameters by interpolating the two compensated DTFT samples

C. EDDSPE ALGORITHM
The eDDSPE algorithm is based on the procedure proposed in [2]. It determines a first estimate of the damped sinewave parameters by applying a complex-valued interpolation based on the three DFT samples X(l − 1), X(l), and X(l + 1). The estimated values are then exploited to iteratively compensate the contribution of the spectral leakage from the fundamental image component on the three considered DFT samples.
Simulations showed that when the magnitude of the smallest of the three considered DFT samples is maximized, the effect of noise on the estimated parameters is minimized. To that aim, a last interpolation is performed using the three DTFT samples located at the estimated frequency and one-bin apart from it. The processing steps that implement the proposed eDDSPE algorithm are detailed in the following.
Step 1: Acquire M samples of the noisy damped sinusoid x(.).
Step 2: Determine l as the normalized frequency of the acquired signal DFT spectrum peak.
Step 3: Determine a first estimate of the damped sinusoid parameters by interpolating the three DFT samples centered on the DFT spectrum peak Step 4: Compensate the contribution of the image component to the used interpolation points Step 5: Determine the ith estimate of the damped sinusoid parameters by using the compensated DFT samples as interpolation points Re ln z Step 6: Compensate the contribution of the image component on two DTFT samples located one-bin apart from the estimated frequencỹ Step 7: Determine the final estimate of the damped sinusoid parameters by interpolating the three compensated DTFT samples

III. PROPOSED LSF ALGORITHM
The accuracy of the estimators proposed above can be further improved by fitting the analyzed samples to model (1). As compared to the nonlinear signal fit algorithm [12], the LSF algorithm provides closed-form parameters estimators by exploiting the Gauss-Newton method [5], thus reducing the required processing effort [5], [21]. Simulations showed that the eADSPE and eBDSPE algorithms exhibit almost the same accuracy and outperform the eDDSPE algorithm when the integer part of the number of observed cycles l is less than or equal to 2. Conversely, when shorter observation intervals are considered, the eDDSPE algorithm returns the most accurate estimates. Thus, to determine the initial parameter estimates of the LSF algorithm, the eDDSPE algorithm is preferred when ν < 2.5 cycles, while for longer observation intervals, the eADSPE algorithm is employed since it requires a lower processing effort than the eBDSPE algorithm. The steps that implement the proposed LSP algorithm are detailed in the following.
Step 1: Acquire M samples of the noisy damped sinusoid x(.).
Step 2: Determine l as the normalized frequency of the acquired signal DFT spectrum peak.
Step 4: Determine the ith in-phase and in-quadrature signal samples and the residual signal for m = 0, 1, . . . , M − 1 Step 5: Determine the ith differences related to the damped sinusoid parameters Step 6: Determine the ith estimate of the damped sinusoid parameters Step 7: Return the estimates ν (Q) , α (Q) ,Â (Q) , and φ (Q) .
Simulations showed that the number of iterations needed to approach the best estimation accuracy decreases when the observation interval length increases (due to the smaller contribution of spectral leakage), and the SNR of the analyzed signal decreases. Fig. 1 shows the RMSEs of the parameter estimates returned by the eBDSPE, eADSPE, and eDDSPE algorithms as a function of the number of iterations when analyzing the data related to the shortest observation interval and the lowest SNR considered in this work, that is, ν = 1.5 cycles and SNR = 60 dB. The signal damping factor and phase are equal to α = 1 and φ = π /8 rad, respectively.
In Fig. 1, it can be seen that the minimum number of iterations required to approach the best accuracy are equal to 13, 10, and 6, for the eBDSPE, eADSPE, and eDDSPE algorithms, respectively. Thus, these values have been adopted when the eADSPE or eDDSPE algorithms are used in the determination of the initial parameter estimates of the LSF algorithm. Also, simulations showed that five iterations of the LSF algorithm allow us to almost achieve the best accuracy. Thus, the number of iterations specified above is adopted in the simulation results reported in the following.

IV. ALGORITHM PERFORMANCE COMPARISON
In this section, the RMSEs of both the proposed frequencydomain eBDSPE, eADSPE, and eDDSPE algorithms and the time-domain LSF algorithm are compared with each other through computer simulations. Specifically, the RMSEs are referred to the asymptotic CRLBs for unbiased estimators, so that the algorithm statistical efficiency is analyzed. It is worth remembering that, assuming s = e −(4π/M)α , the asymptotic CRLBs for the frequency and the damping factors are given by [22] while for amplitude and phase, we have [22] σ 2

A CR
Observe also that RMSEs can be consistently referred to the CRLB since simulations showed that the effect of noise on estimation biases is negligible as compared to the variance. The amplitude of the analyzed noisy damped sinusoids was set to A = 1 p.u. and for each simulation result, 10 000 runs were considered. Since the variances of the frequency and the damping factor estimators returned by the algorithms proposed in [14] are close to the related CRLBs when less than 2000 samples are analyzed, the acquisition length M was varied between 128 and 2048 samples. Fig. 2 shows the ratios RMSE/ √ CRLB for the frequency [ Fig. 2(a)], the damping factor [ Fig. 2(b)], the amplitude [ Fig. 2(c)], and the phase [ Fig. 2(d)] estimators returned by the considered algorithms as a function of the number of observed cycles ν when the damping factor is α = 1, phase is φ = π /8 rad, and SNR is 40 dB. The number of cycles ν varies in the range [1.5, 11.4] with a step of 0.1 cycles and the acquisition length is M = 512 samples. However, very similar behaviors were observed for different values of both the signal phase and the acquisition length.
As we can see, the LSF algorithm provides much more accurate estimates than the other algorithms when at least about 2 cycles are observed. Also, the returned estimates almost achieve the related √ CRLB as soon as more than about 4 cycles are analyzed. In that condition, since the accuracy loss with respect to the nonlinear signal fit algorithm (which attains the CRLBs [2]) can be hardly justified due to the lower processing effort required by the LSF algorithm.
Simulations also show that the eBDSPE and eADSPE algorithms outperform the eDDSPE algorithm. Most probably, that occurs because the latter algorithm processes three interpolation points, while the former two algorithms exploit only two suitably selected interpolation points.
As for the eBDSPE and eADSPE algorithms, they exhibit almost the same variances, whose theoretical expressions can be obtained from those derived in [14] for complex-valued damped sinusoids   [3], [14] and Duda [7] (BDSPE, ADSPE, and DDSPE) algorithms and their corresponding enhanced versions as a function of the number of observed cycles ν. The same parameters and simulation conditions adopted in Fig. 2 are considered.
As it can be seen, thanks to spectral leakage compensation, the eBDSPE, and the eADSPE algorithms well outperform the corresponding classical algorithms. Moreover, the DDSPE algorithm outperforms the other classical algorithms due to spectral leakage compensation. Anyway, the eDDSPE algorithm always returns the minimum RMSE value provided by the DDSPE algorithm due to the second interpolation performed in step 7. Fig. 4 shows the ratios RMSE/ √ CRLB for the frequency [ Fig. 4(a)], the damping factor [ Fig. 4(b)], the amplitude [ Fig. 4(c)], and the phase [ Fig. 4(d)] estimators returned by the considered algorithms as a function of the damping factor α when ν = 5.3 cycles are observed, φ = π /8 rad, and SNR = 40 dB. The damping factor α varies in the range [0.1, 2] with a step of 0.1, while the other simulation parameters are the same as in Fig. 2.
As expected, Fig. 4 shows that the LSF algorithm significantly improves the estimation accuracy of the other proposed algorithms and it is capable of almost achieving the CRLB when the damping factor α is smaller than about 1. Moreover, the accuracies of all algorithms abruptly worsen as α increases. Fig. 5 shows the RMSEs of the frequency [ Fig. 5(a)], the damping factor [ Fig. 5(b)], the amplitude [ Fig. 5(c)], and the phase [ Fig. 5(d)] estimators provided by the considered algorithms as a function of SNR when M = 512 samples, damping factor α = 1, ν = 5.3 observed cycles, and φ = π /8 rad. The SNR varies in the range [0, 60] dB with a step of 2 dB and the remaining simulation parameters are the same as in Fig. 2.
As it can be observed in Fig. 5, the LSF algorithm clearly outperforms the other analyzed algorithms in all the considered testing conditions. Moreover, the returned RMSEs steadily approach the related CRLBs when the SNR is greater than about 10 dB. Conversely, at lower SNR values, the probability that outliers occur in the DFT spectrum is no more negligible so that a threshold effect appears and the RMSE of all considered estimators abruptly increases [1], [23]. Simulations showed that the SNR threshold below which that phenomenon occurs decreases as M increases. Fig. 6 shows the average processing time required by the proposed algorithms as a function of M. The adopted platform was based on a laptop computer equipped with a Core i7-1165G7@2.80-GHz processor, 16-GB RAM, and Windows 10. The algorithms were implemented in the MATLAB environment and the timeit.m MATLAB function was exploited to determine the processing time. The integer part of the number of acquired signal cycles l was assumed known a priori, the DTFT samples were computed directly using the definition formula.
As we can see, the smallest processing time is required by the eADSPE algorithm, while the eBDSPE and eDDSPE algorithms exhibit a higher processing time due to a higher number of iterations or the use of three interpolation points.
As expected, the LSF algorithm implies the highest processing effort. Indeed, the average processing time of the LSF algorithm is about 2.00 (if M = 128) or 3.33 times (when M = 2048) higher than that of the eADSPE algorithm.

V. CONCLUSION
In this article, an improved version of the frequencydomain interpolation-based algorithms proposed by Bertocco et al. [3], Aboutanios [14], and Duda and Zielinski [7] for the estimation of the frequency and damping factor of noisy damped sinusoid have been proposed. Moreover, sinusoid amplitude and phase are also estimated. The improved algorithms-called eBDSPE, eADSPE, and eDDSPE algorithms-compensate the effect on the estimated parameters of the spectral leakage from the fundamental image component and minimize the estimator variance through a suitable selection of the interpolation points. The obtained parameter estimates are then further improved by applying an LSF algorithm. The accuracies of all the proposed algorithms have been compared with each other and with the related CRLB through computer simulations. It has been shown that the LSF algorithm is capable to almost achieve CRLB, while the eDDSPE algorithm exhibits the worst accuracy. The required computational times have also analyzed by implementing the algorithms in the MATLAB environment and running them in a laptop computer.