Phase estimation using sine fitting with low number of samples is biased in the presence of additive noise

It is demonstrated here that the initial phase estimation of a sine wave corrupted by additive noise using the least squares sine fitting algorithm is not biased, even for low number of samples, contrary to what happens in the case of the amplitude estimation. Monte Carlo simulations are presented which corroborate the conclusion reached. This result is of particular importance when studying the uncertainty of measurement methods based on phase estimation like ultrasonic non-destructive testing


Introduction
Many signals encountered in innumerous applications is sinusoidal.This can be due to the engineer choice of a sinusoidal stimulus signal for system characterization, due to normal operation, like the power grid, of due to the bandpass nature of the phenomena which can be considered quasi-sinusoidal.
An example of the first case is the measurement of distance using ultrasound waves, like that used as parking assistant in cars.A sinewave is applied to the ultrasonic emitter and the signal obtained in the receiver located next to it will be sinusoidal, due to the reflected acoustic wave sensed [1].Information about the distance can be obtained from the time of flight which can be estimated, for example, by fitting sinewave to the generated and received signals.The time of flight will lead to a phase difference between those two sinewaves that need to be estimated [2].Many other applications use the same idea as sonar [3], underwater surveillance, liquid velocity [4]- [6], geophysical exploration [7], etc.
Other examples where the stimulus signal is chosen to be sinusoidal is in defect detection in conductive materials [8] and in impedance measurement [9].A sinusoidal voltage is applied to the unknown impedance and the current, which turns out to be also sinusoidal due to the linear behavior of the impedance, is measured.The difference in the initial phase of the two sinewaves fitted to the voltage and current waveforms gives information about the argument of the complex number which is the impedance, and which is related with the reactive part of the electronic circuit under test, that is, with the amount and flow of energy stored.
A third example of choosing sinusoidal stimulus signals is analog-to-digital converter (ADC) characterization where different test methods use sinewaves to estimate the linearity, gain or hysteresis of the ADC, for example [10], by fitting data points to sinewaves.
An example where the normal operation of a system involves the use of sinewave is the power grid.Information regarding the phase of the current flowing relates with the capacitive or inductive behavior of the loads [11].
There are also situations where the signals involved are naturally bandpass and can thus be considered approximately sinusoidal.A common reason for heart failure is insufficient pump function which is usually diagnosed with magnetic resonance imaging (MRI).The characterization of this pathology can be done from a sequence of images of the crosssection of the heart obtained using MRI [12].Displacement over time of the heart walls is mapped over time using optical flow techniques.The sequence of images obtained is bandpass filtered around the first harmonic.This allows the modeling of the image by a sine wave with spatial variation of phase [13].Sine fitting techniques are used to estimate that variation.
In the following we derive the bias of the initial phase estimated using the three-parameter sine fitting algorithm as described in [14]- [15] for any number of samples.At the end we conclude that this bias is null.Note that this estimator is a function of the estimated in-phase and in-quadrature amplitudes, as it will be shown in section II and these two estimators are the maximum likelihood estimators which leads to the amplitude and initial phase being asymptotically unbiased, that is, unbiased when an infinite number of samples is used ( [16], chapter 7).Here we show that the phase estimation is unbiased even for a low number of samples.
Considering the mathematical approximations made, it is advisable to substantiate this claim using numerical simulations using a Monte Carlo procedure.This will confirm the initial findings made.This work addresses the presence of additive noise in the measurement system [17]- [19].There are, however many other factors that also have an influence, like sampling jitter [20]- [22], phase noise [23], noise in the power supply [24], quantization error [25]- [27], frequency error [28]- [29] and harmonic distortion [30]- [31].
Additive noise affects also the precision of other estimators like, for example, the histogram test of ADCs [32]- [35] or procedures that determine the amount of noise present in a setup [36]- [37].Here we considered exclusively the case where sinusoidal shaped signals are involved, there are, however, other types of stimulus signals that could be used [38]- [39].

Sinewave Fitting
Consider M data points acquired at instants ti with voltages zi, given by   =  +  ⋅ cos(    + )……….. (1) where  is the amplitude (greater than 0),  is the offset,  is the initial phase (between − and ) and   is the angular frequency (2  ).This data is affected by additive voltage white Gaussian noise,   , with null mean and standard deviation     =   +   ,…………….. (2) where is the sample voltage in the absence of additive noise.
We wish to estimate the sine wave that best fits, in a least square error sense, to these M points.The estimates of the sine wave are obtained, in a matrix form, with [14] [ where ],……….. (5) and where a is the angular frequency of the sinusoid we are trying to adjust to the data.The estimated sine wave amplitude is given by and the initial phase by Here we will assume that the number of samples acquired (M) covers exactly an integer number of periods (J) of the sine wave we are trying to fit to the data.This means that the sine wave frequency (fa), sampling frequency (fs) and number of samples satisfy Note that J and M should be mutually prime so that the M different samples acquired at M different time instants, correspond to M different sine wave phases.If not, you will have less than M different phases which will increase the uncertainty in the estimation of the sine wave parameters.In the case that J is a multiple of M/2, the sampling instants will correspond to only 2 different sine wave phases and matrix D T D will be singular and hence not invertible (you cannot estimate the 3 sine wave parameters with only two data points).
Note that the sampling instants are given by is t i f  . The assumption is reasonable because we can choose whatever values we want for those frequencies and the number of samples.In practice, however, due to instrument inaccuracies, the actual value of those frequencies may not be exactly the values chosen and which satisfy (8) but are close enough considering typical frequency errors smaller than 100 ppm.If a non-integer number of periods is acquired a bias will affect the estimator.In this work, however, we will not consider this scenario.
If the samples cover an integer number of sine wave periods, we have .………. (13) The sine wave parameters can thus be estimated with ,……….( 14)

In-Phase and In-Quadrature Amplitude Mean
The in-phase amplitude is thus given by The expected value of the in-phase amplitude estimation is ,………… Using ( 2) and ( 3) we can write In this work we consider that the frequency of the sine wave is known, and we use its value in the sine fitting algorithm that is, we make a = x.We thus have .

…………(19)
Since we it is also assumed that the acquisition is done during an integer number of periods we have For the expected value of the in-quadrature amplitude, given by the reasoning is similar, and the result is

In-Phase and In-Quadrature Amplitude Second Moment
From (15)

………….
( The expected value of the product of yi and yj is, using (2),

In-Phase and In-Quadrature Amplitude Variance
The variance of both in-phase and in-quadrature amplitudes can be obtained from the mean and the second moment using where x is any random.The variance of the in-phase amplitude is thus, using ( 20) and (33), For the case of the in-quadrature amplitude we have, using ( 22) and ( 34),

In-Phase and In-Quadrature Amplitude Covariance
The covariance between the in-phase and in-quadrature components is obtained from The expected value of the product of the two amplitudes is, using ( 15) and ( 21 Following the same procedure as before it is possible to reach the conclusion that this term is also zero.We thus have which proves that the in-phase and in-quadrature amplitudes have a null covariance.This is going to be used next when we determine the expected value of the estimated initial phase.

Estimated Phase
To obtain the expected value of the estimated sine wave initial phase from the moments of the in-phase an in-quadrature amplitudes, we are going to use a Taylor series approximation given in [14], namely As seen in the previous section, the covariance of the two amplitudes in the presence of additive noise is zero.Also, it was seen that their variances are the same.We thus have The term in the square brackets turns out to be 0 and thus we have which shows that the phase estimator is unbiased for any number of samples.

Monte Carlo Analysis
To validate the approximations made in the mathematical derivations, a numerical simulation, using a Monte Carlo procedure was carried out where a sinusoidal stimulus signal with amplitude , null offset, ( = 0) and frequency   was created, a set of data points were sampled at a rate   covering an integer number of periods of the stimulus signal.The parameters of the best fit sinusoid were obtained using the least squares procedure studied here, and the value of the resulting estimated initial phase was plotted as a function of additive noise standard deviation (Figure 1), initial phase value ( Figure 2) and number of samples ( Figure 3).In Figure 1 we plot the average value of the estimated sinewave as a function of the additive noise standard deviation (circles).The vertical bars represent the confidence interval obtained from the standard deviation of the values obtained for a normal distribution considering a probability of the actual value being on the interval of 99.9% (for all simulations).

Conclusion
The least squares procedure of fitting a sinewave to a set of data points was studied from a mathematical point of view.
The focus here was on the estimation of initial phase.The derivation of its bias lead to the result that this particular estimator is not biased in the case where just additive random noise with a normal distribution is present.
This proves useful in many applications that use sine fitting to estimate different physical quantities, as exemplified in the introduction, and the knowledge of which needs to be unbiased.This complements existing work in the literature that addressed the bias and precision of the estimates of sinusoidal amplitude which is also important for many applications.
In the future we will address the issue of precision of the initial phase estimation of sinewave fitted to a set of data points.Note that many quantities that are derived from the sinewave amplitude and initial phase, like total harmonic distortion (THD), signal to noise ratio (SNR), signal to noise and distortion ratio (SINAD) and noise floor, for example, could also be the subject of studies relating to their bias and precision in different situations like the presence of additive noise, phase noise, jitter, frequency error in the signal or the sampling frequency.

Figure 1
Figure 1 Expected value of the estimated sine wave initial phase as a function of the additive noise standard deviation.The circles represent the values obtained with the Monte Carlo analysis.The confidence intervals for a confidence level of 99.9% are represented by the vertical bars The sine fitting procedure was repeated many times () in the same conditions and the average value of estimated initial phase was computed and the difference to the actual value was determined -the phase estimation error, represented by   ̂ .That is what is shown in the vertical axis of the charts.

Figure 2 Figure 2
Figure 2 Expected value of the estimated sine wave initial phase as a function of the initial phase.The circles represent the values obtained with the Monte Carlo analysis.The confidence intervals for a confidence level of 99.9% are represented by the vertical bars