EVALUATION OF A STANDARDIZED SINE WAVE FIT ALGORITHM Peter H¨andel Department of Signals, Sensors and Systems Royal Institute of Technology SE-100 44 Stockholm, Sweden [email protected] ABSTRACT The IEEE Standard 1057 provides algorithms for fitting the parameters of a sine wave to noisy discrete time observations. The fit is obtained as an approximate minimizer of the sum of squared errors, i.e. the difference between observations and model output. The performance of the algorithm in §4.1.3.3 of the standard is compared both with the Cram´er-Rao bound (CRB) on accuracy, and with the performance of a non-linear least squares approach. The algorithm of IEEE-STD-1057 provides accurate estimates in most practical applications for Gaussian and quantization noise. In the Gaussian scenario it provides estimates with performance close to the CRB. In severe conditions with noisy data covering only a fraction of a period, however, it has inferior performance compared with a one-dimensional search of a concentrated cost function. 1. INTRODUCTION In testing digital waveform recorders, an important part is to fit a sinusoidal model to recorded data, and to calculate the parameters that result in the best fit (in least-squares). Several algorithms were evaluated in [1, 2], and algorithms have been standardized in IEEE Standard 1057 [3]. A survey of the standard can be found in [4]. The standard [3] was prepared by a working group that is a part of the Waveform Measurement and Analysis Technical Committee (TC-10) of the IEEE Instrumentation and Measurement Society. TC10 is currently working on a standard for analog to digital converters (IEEE-STD-1241) [5]. In IEEE-STD-1241, test methods for signal-to-noise and distortion ratio (SINAD) and effective number of bits (ENOB) rely on the sine wave fit in [3]. In this paper, the 4-parameter sine wave fit algorithm of [3, §4.1.3.3] is studied in some detail. The considered estimation problem is non-linear with respect to the parameters, and thus the 4-parameter sine fit algorithm is an iterThis work was supported in part by the Junior Individual Grant Program of the Swedish Foundation for Strategic Research.

ative method where in each iteration an updated frequency estimate is obtained based on estimated parameters. The algorithm of [3] does not exploit the fact that, although it is a non-linear optimization problem, it is linear in three out of four parameters. Thus, it can be expected that the algorithm of [3] has inferior performance compared to algorithms that utilize such a fact. In particular, one can expect that the algorithm may suffer from ill-convergence when the digital waveform recorder utilizes course quantization. In the worst case scenario with 1-bit quantization the magnitude of the sine wave amplitude is not observable from data, and thus ill-convergence may occur. Although, the recommended procedure to obtain initial estimates virtually eliminates convergence problems in most practical cases [4]. The problem under study is nonlinear in the frequency of the sine wave, and thus the main interest in this paper is to study the properties of the obtained frequency estimate. Further, we are interested in the small error properties, i.e. performance of the algorithm when it is properly initialized with initial values close to the correct ones. 2. MEASUREMENTS AND DATA MODEL Assume that the data record contains the sequence of samples y1 , . . . , yN

(1)

taken at time instants t1 , . . . , tn , . . . , tN . It is further assumed that data can be modeled by yn [A, B, C, ω] = A cos ωtn + B sin ωtn + C

(2)

where A, B are C are unknown constants. The angular frequency ω may be known, or not, leading to models with three or four parameters, respectively. Equivalently, we may write yn [x] in (2) as yn [x] = α sin(ωtn + φ) + C where A = α sin φ, B = α cos φ, and x is the vector of unknown parameters. Stressing the dependence of yn [x] on the parameter vector x turns out to be convenient for the following discussion.

Table 1: IEEE-STD-1057 for four parameter least squares fit to sine wave data using matrix operations. a) b) c) d) e) f) g)

ˆ 0 = (A0 , B0 , C0 , 0)T , ω x ˆ0, i = 0 next iteration i = i + 1 ω ˆi = ω ˆ i−1 + ∆ˆ ωi−1 ˆ i , see (12) create D ˆ Ti D ˆ i )−1 D ˆ Ti y ˆ i = (D x optional, see [3] repeat b)-f) until convergence

and sin ωtn ≈ sin ω ˆ i tn + tn cos ω ˆ i tn · ∆ωi where ∆ωi = ω − ω ˆ i . Inserting (7)-(8) into (2) gives yn [xi ] ≈ A cos ω ˆ i tn + B sin ω ˆ i tn + C −Atn ∆ωi sin ω ˆ i tn + Btn ∆ωi cos ω ˆ i tn (9) where xi is the parameter vector xi = (A, B, C, ∆ωi )T .

3. ALGORITHMS @ IEEE-STD-1057 The Standard 1057 provides algorithms both for 3-parameter (known frequency) and 4-parameter (unknown frequency) models. For easy reference, the 3-parameter algorithm is reviewed below [3, §4.1.3.1]. It is exploited by the non-linear least squares fit in Section 3.3. In Section 3.2, a derivation of the 4-parameter algorithm in [3, §4.1.3.3] is provided, a derivation that differs from the one in [3]. Thus, it is believed to be informative. 3.1. Known frequency When the frequency ω is known, estimates of the unknown parameters x = (A, B, C)T

(3)

(where T denotes transpose) are obtained by a least squares fit. Gather the data record in y = (y1 , . . . , yN )T . Then, y obeys the linear set of equations y = Dx where D is the N × 3 matrix  cos ωt1 sin ωt1  cos ωt2 sin ωt2  D= .. ..  . . cos ωtN

sin ωtN

(4)

1 1 .. . 1



  . 

(5)

Equation (4) is an overdetermined (i.e. N > 3) set of linear ˆ given by [3] equations, with the least-squares solution x ˆ = (DT D)−1 DT y. x

(6)

3.2. Unknown frequency Assuming that an estimate in iteration i, say ω ˆ i , of ω is available, a Taylor series expansion around the estimate ω ˆi gives cos ωtn ≈ cos ω ˆ i tn − tn sin ω ˆ i tn · ∆ωi

(8)

(7)

(10)

Equation (9) is still nonlinear in the parameters, but may be linearized using the observation that ∆ωi ≈ 0. Putting available estimates of A and B from previous iteration, i.e. ˆi−1 , in place of the unknown parameters in the Aˆi−1 and B two last terms in (9) results in an equation linear in the components of xi . Gathering the data record in y gives, similarly as in (4) ˆ i xi y=D

(11)

ˆ i the N × 4 augmented D matrix with an extra where D column ˆi = D  cosˆ ωi t1 sinˆ ω i t1  cosˆ ω i t2  ωi t2 sinˆ  . .. .  . . cosˆ ωi tN sinˆ ω i tN

(12)  ˆi−1 t1 cos ω 1 −Aˆi−1 t1 sin ω ˆ i t1 + B ˆ i t1 ˆi−1 t2 cos ω 1 −Aˆi−1 t2 sin ω ˆ i t2 + B ˆ i t2    .. ..  . . ˆ ˆ 1 −Ai−1 tN sinˆ ωi tN + Bi−1 tN cosˆ ω i tN

The basic idea behind the algorithm in [3, §4.1.3.3] is repeatedly to solve the linear system (11), i.e. at iteration i ˆ i in order to obtain a new set of estimates x ˆ i . The aluse D ˆ0 gorithm is summarized in Table 1. The initial estimates x and ω ˆ 0 are, for example, obtained by peak-picking the Discrete Fourier Transform (DFT) of data, followed by a prefit using the algorithm in Section 3.1. Alternative methods for finding initial estimates are discussed in [1]. The algorithm in Table 1 is an iterative process to find the parameters that minimizes the sum of squared differences, i.e. N X

(yn − yn [x])2 .

(13)

n=1

In (13), yn [x] is given by (2). In each iteration, an updated frequency estimate ω ˆ i is obtained based on ω ˆ i−1 and ∆ˆ ωi−1 , estimates that also depend on estimated values of A and B. The recommended procedure to obtain initial estimates virtually eliminates convergence problems in most practical cases [4]. However, convergence problems may occur, especially for short noisy data records or signals at low frequency (see, Section 4). An alternative to the 4parameter algorithm in [3] is derived below.

N=16, SNR=10 dB

Table 2: A non-linear least-squares fit by grid search. 50

frequency grid ωi , i = 1, . . . , M for i = 1 to M create D from ωi , see (5) gi = yT D(DT D)−1 DT y end ω ˆ = ωk where gk = max[gi |i = 1, . . . , M ] ˆ is obtained from (6) x

IEEE 1057 NLS AS CRB

40

−10LOG(MSE) (dB)

a) b) c) d) e) f) g)

30

20

10

0

3.3. Non-linear least-squares by grid search Using (4)-(5) the criterion (13) can be rewritten as (y − Dx)T (y − Dx) where x = (A, B, C)T , and D is defined in (5). The criterion can be concentrated with respect to x, and ω can be found by a one-dimensional search for the maximum of g(ω), i.e. [6, 7] ω ˆ = arg max g(ω), ω

g(ω) = yT D(DT D)−1 DT y. (14)

The parameters in x are then obtained from the least-squares fit in Section 3.1 with ω there replaced by ω ˆ given by (14). The maximization of (14) may be implemented by iterative methods, or simply by a one-dimensional grid search, as presented in Table 2. The frequency grid may be obtained by peak-picking the DFT of data and with ω1 corresponding to the frequency bin to the left of the maxima, and ωM corresponding to the bin on the opposite side. The number of grid points M is chosen depending on the desired resolution. 4. NUMERICAL EVALUATIONS 4.1. Monte Carlo simulations In order to verify the performance of the considered algorithms, the Cram´er-Rao bound (CRB) [8] is compared with the normalized sum of squared errors obtained from Monte Carlo simulations, based on 1000 independent realizations. For simplicity, we consider uniform sampling tn = n/fs with fs being the sampling frequency and n a running integer time index. The empirical mean square error (mse) is computed as mse(fˆ) =

1000 1 X ˆ (f ` − f ) 2 1000

(15)

`=1

where fˆ` denotes the frequency estimate in the `-th realization. 4.2. Data records with random initial phase In Figure 1, the 4-parameter matrix algorithm in [3] (Table 1) and the non-linear least-squares (NLS) method in Table 2

−10

−20 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

FREQUENCY

Figure 1: Mean square frequency estimation error versus frequency for noisy sinusoidal signal with random initial phase.

are evaluated for short records of N = 16 noisy samples, with sampling frequency fs = 1. A sinusoidal signal with random initial phase (uniformly distributed in [0, π]) in additive white Gaussian noise was generated, with an SNR of 10 dB. As estimator of the initial frequency a DFT based estimator was used with 4 times zero-padding, and peak-finding by triple parabolic interpolation. The grid search for the algorithm in Table 2 was performed (rather arbitrary) in M = 160 points in the symmetric interval (of length corresponding to twice the frequency resolution of the DFT) around the maxima of the DFT. The initial estimate of the nuisance parameters required for the algorithm in Table 1 was estimated using the 3-parameter fit with ω ˆ 0 in place of ω. The 4-parameter matrix algorithm was aborted after 5 iterations. Figure 1 shows the mse (15) as function of frequency f . As reference the asymptotic CRB [8] is shown. From the figure one may note an excellent performance of both algorithms for frequencies well inside (0, 0.5). For frequencies near 0 or 0.5, the algorithm in Table 2 outperforms the one in [3]. 4.3. Data records with fixed initial phase In order to compare the performance of the algorithms with the exact CRB [8] the above experiment was repeated for low frequency signals (f < 0.1). The conditions were set as in the experiment above, but now the initial phase φ was fixed and set to φ = π/7. The comparison with the exact CRB in Figure 2 reveals that the empirical mse can be seen as an empirical variance for frequencies above f = 0.03 (for NLS) and f = 0.045 (for IEEE-1057), respectively. For

N=16, SNR=10 dB

N=16, SNR=10 dB

60

60

50

50

40

−10LOG(MSE) (dB)

−10LOG(MSE) (dB)

IEEE 1057 NLS Exact CRB

30

20

40

20

10

10

0

0

−10

−10 0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

FREQUENCY

IEEE 1057 (2−BITS) NLS (2−BITS) IEEE 1057 (1−BIT) NLS (1−BIT)

30

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

0.1

FREQUENCY

Figure 2: Mean square frequency estimation error versus fre-

Figure 3: Mean square frequency estimation error versus fre-

quency for noisy sinusoidal signal with fix initial phase.

quency for quantized noisy data (1 and 2 bits, respectively).

frequencies below f = 0.03, the empirical mse increases for both methods, indicating that for (very) low frequencies both methods suffer from bias errors. For N = 16, the frequency f = 0.03 corresponds to less than half a period of a uniformly sampled sine wave.

experimental conditions, characterized by a small number of samples at a low frequency sine wave. Finally, it is worth noticing that the performance of the 4-parameter algorithm is studied in a small error context, i.e. it relies on precise initial estimates. In this paper, the initial estimates were obtained following the recommended procedure in [3]. In general, its convergence properties strongly depends on the initial estimates, whereas the NLS relies on a one-dimensional grid search, i.e. its convergence is ensured.

4.4. Quantized data records The sensitivity for course quantization is studied in this experiment. The noisy sine wave in Section 4.2, i.e. with random initial phase, is quantized with one and two bits, respectively. The results for this scenario are displayed in Figure 3. Similar results as in Figure 2 are obtained. Repeating the experiment with noise free quantized data results in similar performance as in Figure 3, except for the IEEE 1057 algorithm applied to 1-bit data where the performance coincides with the performance of NLS for f > 0.055 (an improvement compared with f > 0.065 in Figure 3). 5. CONCLUSIONS The four parameter sine wave fit algorithm in IEEE Standard 1057 [3, §4.1.3.3] has been studied in some detail. In most practical applications it seems to be an excellent method for parameter estimation, and for reconstruction of a sine wave from noisy or quantized observations. Its performance (in terms of frequency estimation error) has been compared with theoretical bounds on accuracy, and with the performance of a non-linear least squares (NLS) approach. The comparison reveals that its performance is close to optimal in a Gaussian scenario, but it also reveals that it is inferior compared with the performance of the NLS in severe

6. REFERENCES [1] J. Kuffel, T. McComb and Malewski, “Comparative evaluation of computer methods for calculating the best fit sinusoid to the high purity sine wave”, IEEE Trans. on Instrumentation and Measurement, Vol. IM-36, No. 2,1987, pp. 418-422. [2] T.R. Mccomb, J. Kuffel, and B.C. Le Roux, “A comparative evaluation of some practical algorithms used in the effective bits test of waveform recorders”, IEEE Trans. on Instrumentation and Measurement, Vol. 38, No. 1, 1989, pp. 37-42. [3] IEEE Standard for digitizing waveform recorders, IEEE Standard 1057, 1994. [4] T.E. Linnenbrink, “Waveform recorder testing: IEEE standard 1057 and You”, Instrumentation and Measurement Technology Conference, 1995. IMTC/95, pp. 241–246. [5] S.J. Tilden, T.E. Linnenbrink and P.J. Green, “Overview of IEEE-STD-1241 - Standard for terminology and test methods for analog-to-digital converters”, Proc. 16th IEEE Instrumentation and Measurement Technology Conference, IMTC/99, Vol. 3, 1999 , pp, 1498–1503. [6] P. Stoica and R. Moses, Introduction to Spectral Analysis, Prentice–Hall, Upper Saddle River, NJ, 1997. [7] S.M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory, Prentice–Hall, Upper Sadle River, NJ, 1993. [8] P. H¨andel, “Properties of the IEEE-STD-1057 four parameter sine wave fit algorithm”, Technical report IR-S3-SB-0009, Department of Signals, Sensors and Systems, KTH, Sweden, 2000.