Abstract:
This is a method and associated apparatus for accurately and quickly estimating the amplitude, frequency and phase of a signal of interest. The method comprises the steps of, inputting the signal of interest; generating a reference signal with adjustable amplitude, frequency and phase at an output thereof; mixing the signal of interest with the reference signal and a signal 90° out of phase with the reference signal to provide a pair of quadrature sample signals comprising respectively a difference between the signal of interest and the reference signal and a difference between the signal of interest and the signal 90° out of phase with the reference signal; using the pair of quadrature sample signals to compute estimates of the amplitude, frequency, and phase of an error signal comprising the difference between the signal of interest and the reference signal employing a least squares estimation; adjusting the amplitude, frequency, and phase of the reference signal from the numerically controlled oscillator in a manner which drives the error signal towards zero; and, outputting the estimates of the amplitude, frequency, and phase of the error signal in combination with the reference signal to produce a best estimate of the amplitude, frequency, and phase of the signal of interest. The preferred method includes the step of providing the error signal as a real time confidence measure as to accuracy of the estimates wherein the closer the error signal is to zero, the higher the probability that the estimates are accurate. A matrix in the estimation algorithm provides an estimate of the variance of the estimation error.

Description:
ORIGIN ON THE INVENTION 
     The invention described herein was made in the performance of work under a NASA contract, and is subject to the provisions of Public Law 96-517 (35 USC 202) in which the Contractor has elected not to retain title. 
    
    
     This application is a continuation of application Ser. No. 523,692, filed May 15, 1990, now abandoned. 
    
    
     TECHNICAL FIELD 
     The invention relates to signal analysis apparatus and methods and, more particularly, a method of accurately and quickly estimating the amplitude, frequency and phase of a signal of interest. 
     BACKGROUND ART 
     The problem of estimating the parameters of a sinusoidal signal has received considerable attention in the literature and in various prior art patents. Such a problem arises in diverse engineering situations such as carrier tracking for communications systems and the measurement of Doppler in position location, navigation, and radar systems. 
     A variety of techniques have been proposed in the literature to solve such problems including, to mention a few, the application of the Fast Fourier Transform (FFT), one and two-dimensional Kalman filters based on a linearized model, a modified Kalman filter that results in a phase locked loop, and a digital phase locked loop derived on the basis of linear stochastic optimization. 
     U.S. Pat. No. 4,686,532 of McAulay describes a system for determining the strength of distributed scatterers (say P in number) that are possibly close to the sonar or radar relative to the array dimension of size I. Assuming that P scatterers rediate signals at some of the L Doppler frequencies known in advance, the method uses a least squares approach to estimate the complex amplitude of these P scatterers at each of these L possible frequencies. It is also assumed that various attenuations, time delays and angles between the P scatterers and I receivers are available via some other means. Since, of course, the Dopplers are not known a-priori, the method would essentially segment the entire Doppler range into intervals of delta-f and measure the complex signal amplitude of each of these frequencies and at each of the I receivers by some prior known method. Essentially, the method obtains the knowledge of the spectral contents of the scatterers from the known spectral contents of the received signals. It does not address the problem of how to obtain the spectral contents of the received signal. This latter problem is the subject of the present invention. Moreover, the framework of the scheme is also very computationally intensive as it involves an exhaustive search in the Doppler domain. It also does not take into consideration the time-varying Doppler--a clear possibility in a near field situation. Thus, in this framework, a discrete scatterer with time-varying Doppler would be identified by L complex amplitudes that are functions of time rather than identifying it as a single sinusoidal signal with time-varying amplitude, phase and frequency as would be more desirable. 
     U.S. Pat. No. 4,533,918 of Virnot describes a device working with radar equipment to obtain the location of a target. This is achieved by a Kalman filter type of circuit that operates on the samples of range measurements from the ship to a plurality of known landmarks. The range measurement is done by conventional pulse radar. It does not describe any scheme for Doppler measurement, which is part of the area of interest addressed by this invention. 
     U.S. Pat. No. 4,179,696 of Quesinberry contains a Kalman estimator which essentially transforms the line-of-sight (LOS) measurement data (i.e. range, angle, etc.) to obtain the target position as referenced to a stable coordinate system. This again is a quite different problem. It involves no scheme for a Doppler measurement and is essentially concerned with the transformation of one set of data (LOS coordinates) into a different set (fixed or stable coordinates) as is the case with the Virnot patent mentioned above. 
     The patent of Wu (U.S. Pat. No. 4,471,357) describes a method of two-dimensional processing of SAR images which is an entirely different problem than that of primary concern herein. 
     The following general commentary may be helpful in placing a proper perspective on the present invention as will be described hereinafter in relation to least squares algorithms and the Kalman filter. The basic idea of minimizing the sum of squares of the residuals goes back to Gauss in the 19th century who applied it to the study of problems in celestial mechanics. Since then, as is evidenced by a very large body of published scientific research, this generic idea has permeated through many scientific, engineering, and other fields. In this regard, the celebrated Kalman filter is also a specific implementation of the least squares algorithm; but, limited to a very structured problem of linear systems with Gaussian noise. For the more general nonlinear case (as is true with the environment of principal concern herein) there is no single scheme of a general nature. One possible solution, although not always satisfactory, is to use a linear approximation and apply the Kalman algorithm to this linearized model. A sinusoidal function can hardly be considered to be linear except over a very small segment. Thus, the above-referenced algorithms, although using least squares and Kalman filters, do not solve this problem. 
     Indeed, the fact that there are so many different techniques to solve the problem indicates the importance of the problem. This, however, also implies that there is no single technique superior to all others in all possible situations and/or with respect to different criteria such as computational complexity, statistical efficiency, etc. 
     STATEMENT OF THE INVENTION 
     Accordingly, it is the object of this invention to provide a novel and elegant least squares scheme for the nonlinear estimation problem of estimating the phase, frequency, and amplitude of a sinusoid, very rapidly and precisely without making any linearizing approximation. 
     Other objects and benefits of the invention will become apparent from the detailed description which follows hereinafter when taken in conjunction with the drawing figures which accompany it. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a simplified functional block diagram of a system employing the method of the present invention to estimate the frequency, phase, and amplitude of an incoming signal of interest. 
     FIG. 2 is a conceptual aid to the understanding of the overall algorithm of the present invention. 
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     The methods and apparatus to be described hereinafter are based on a new least squares algorithm which can be employed to advantage for fast frequency and phase acquisition of sinusoids in the presence of noise. The algorithm is a special case of more general adaptive parameter estimation techniques. The advantage of this particular algorithm is its conceptual simplicity, flexibility, and application to general situations. For example, the frequency to be acquired can be time varying, and the noise can be non-Gaussian, non-stationary, and colored. As the algorithm can be made recursive in the number of observations, it is not necessary to have a-priori knowledge of the received signal-to-noise ratio or to specify the measurement time. This, of course, would be required for batch processing techniques such as the Fast Fourier Transform (FFT). The algorithm also improves the frequency estimate on a resursive basis as more and more observations are obtained. When the algorithm is applied in real time, it has the extra advantage that the observations need not be stored. The algorithm also yields a real time confidence measure as to the accuracy of the estimator. 
     The inventor herein has authored a paper entitled &#34;Fast frequency acquisition via adaptive least-squares algorithm&#34; (hereinafter referred to as &#34;the Paper&#34;) which describes his early work with respect to his algorithm. The Paper was published in the Proceedings of the International Telemetering Conference, Las Vegas, pp. 91-101, Oct., 1986 and later re-published in extended and modified form in the IEE PROCEEDINGS, Vol. 136, Pt.F, No.4, Aug. 1989. Therefore, the contents of the Paper in its 1986 version only are prior art to the present invention which comprises a novel improvement to that prior work thereby solving certain problems associated therewith. For simplicity, a copy of the extended and modified form of the Paper is attached hereto as Appendix A and the contents thereof are incorporated herein by reference. Only those aspects which are an improvement over the disclosure of the Paper will be addressed in detail hereinafter. 
     The algorithm as disclosed herein differs from the present inventor&#39;s earlier published work in the Paper in several non-obvious and important ways which will now be delineated. 
     To keep the truncation errors small, the dimension n of the state vector x k  in equations (4) and (5) of the Paper must satisfy the following bound: 
     
         n≧ΩT 
    
     where Ω is the initial frequency uncertainty in rad/s and T is the observation interval. Thus, for example, if Ω=2π×100 rad/s and T=0.1s, then the dimension n must be at least 60. If the observation interval is increased to 1s, n must be greater than 600. Moreover, as the computational requirements (in one implementation) for processing each signal sample increase as n 2 , one would require in the order of 3,600 and 360,000 operations respectively for these two examples. Obviously, this is a very high computational requirement and in actual practice one has to simplify the algorithm perhaps by some simple ad hoc procedure to keep the computational requirements manageable. However, such a possible ad hoc simplification would result in the loss of optimality in terms of the estimation errors. 
     Related to the foregoing is the fact that the variance of the frequency estimation error given by equation (16) of the Paper contains a constant K. For low values of the state dimension n, the constant K is equal to 4 and then the results correspond to the bound for optimum variance derived in the literature for the case of high signal-to-noise ratio (SNR). For high value of n, however, K may be higher. The reason for this is that equation (16) is derived when ω d  is estimated via equation (12) and is based on just two major elements of the parameter matrix θ. For a low SNR case and/or high value of n, to keep the optimality, the final estimate of ω d  should be based on the complete matrix. Thus, a procedure is required for optimally combining the elements of θ so as to yield an optimum estimate for ω d  for all conditions. The present invention accomplishes this. 
     Since the algorithm is based on minimizing the error between the estimated frequency and phase at any particular instant in time and the present value of a standard against which the incoming signal of interest is compared, the computational requirements of the algorithm as presented in the Paper tended to increase continually as the number of terms to be computed also continued to increase. This, of course, is not a desirable attribute under actual working conditions. That is, the present inventor&#39;s prior work provided to be suitable for its intended purpose under short term laboratory conditions; but, would impose an inordinate computational load on actual field equipment employing the algorithms thereof. The present invention eliminates that problem. 
     Finally, it was not obvious or apparent from the disclosure of the Paper how the basic algorithm could be extended to the important case of M sinusoids. The present invention discloses this aspect for the first time. 
     Thus, although the present inventor&#39;s prior publications introduced a novel concept, they did not contain a simple realizable algorithm except in some limited cases as for the case of small values of ΩT product. The extension to the more general problem of estimating several sinusoid was not attempted as the way for such a generalization was not clear at the time. The present invention provides elegant solutions to all of the above-described problems as will now be discussed in detail. FIG. 1 shows the manner in which the present invention can be implemented. The system 10 inputs the signal of interest at 12. The signal 12 is mixed with the output of a numerically controlled oscillator (NCO) 14 and a signal 90° out of phase with the output of the NCO 14 at the mixers 16, 16&#39;, respectively. As is well known, when two frequencies are mixed together, there are two outputs from the mixing process, their sum and their difference. The lower possible frequency (i.e. the difference), of course, also represents the error between the two frequencies. Thus, in the system 10 of FIG. 1, the sum frequency signals are removed at 18 and the difference frequency signals are sampled at 20 to provide quadrature sample signals at 22 and 22&#39;. The two quadrature sample signals 22,22&#39; are then employed at 24,26, and 28 to compute estimates of the amplitude, frequency, and phase of the error signal generally according to the techniques of the Paper. At 30 the NCO 14 is adjusted in a manner which will tend to drive the error towards zero. The error amount is what acts as a real time confidence measure as to the accuracy of the estimator as mentioned earlier herein; that is, as the computation progresses over multiple samples and the estimation closes in on the actual amplitude, frequency and phase, the error approaches zero. Thus, the closer the error is to zero, the higher the probability that the estimate is accurate. A matrix implemented as part of the estimation algorithm actually provides an estimate of the variance of the estimation error. 
     In the algorithm according to the present invention, the dimension of the state vector n is determined from equation (18) of the Paper and is given by: ##EQU1## and thus n is made independent of the observation interval T. Thus, for Ω=200 πrad/s and a typical value of P/N O  24 dB-Hz(251), the required value of n is equal to 8 as against n&gt;60 required before for the case of T=0.1s. This reduces the computational requirements by a factor of more than 50. For the case of T=1s, the improvement factor is 2,500. The important thing to realize is that n is no longer a function of T. This is achieved by segmenting the total observation period into T O  =nΩ -1  s intervals, applying the algorithm to successive intervals as described in the Paper, and at the end of each interval applying a frequency and phase correction to the NCO 14 of FIG. 1. The overall effect of this method is that the algorithm effectively seems to be operating on a smaller than ω d0  initial frequency ω d1  between (0, 2T O ), an even smaller frequency ω d2  between (0, 3T O ) and so on without losing any optimality in terms of the estimation accuracy. This is because, according to the present invention, the frequency and phase of the NCO 14 are constantly being adjusted towards the actual incoming signal of interest such that the error is driven towards zero and, correspondingly, the number of terms to be computed decreases as the error also decreases rather than increasing as in the prior art approach described in the Paper. This is undoubtedly the major contribution of the present invention as it turns what would otherwise have been a very limited use algorithm into a most viable approach to solving the problem at hand on a commercial and practical basis. 
     The present invention can employ either of two methods which combine the elements of the matrix θ so as to produce the optimal estimate of ω d  rather than deriving it on the basis of just two elements of θ via equation (12) as was previously done at the cost of losing optimality, especially for the cases of large n and low SNR. According to one method, by simply adding the appropriate terms of θ one can estimate Sin ω d  and Cos ω d . From a well known theorem in estimation theory, optimality of θ also implies the optimality of the estimates of Sin ω d  and Cos ω d . From these estimates ω d  is then uniquely estimated. The alternate method is to apply Prony&#39;s method [see N. Mohanty, Random Signal Estimation and Identification, Van Nostrand Reinhold Company, New York, 1986] to optimally estimate ω d  from the matrix θ as shown in the Paper. FIG. 2 is included herewith as a conceptual aid in the understanding of the overall algorithm of the present invention. 
     A detailed description by the inventor herein relative to extending the adaptive least-squares algorithm of the present invention to multiple sinusoids is attached hereto as Appendix A and the contents thereof are incorporated herein by reference. 
     APPENDIX A 
     Fast Frequency Acquisition Via Adaptive Least-Squares Algorithm 
     Consider the problem of estimating an unknown frequency w d  from the measurement y k , z k  below: ##EQU2## Here, {y k , z k  } represent the samples of the inphase and quadrature components of a received signal s(t) obtained by demodulating s(t) by a carrier reference signal r(t) and its 90° phase-shifted version respectively, i.e. 
     
         s(t)=A sin (w.sub.o t+φ.sub.o)+n(t) 
    
     
         r(t)=2 sin (w.sub.c t+φ.sub.c);φ=φ.sub.o -φc, w.sub.d =w.sub.o -w.sub.c                                         (2) 
    
     with n ik  and n qk  denoting the inphase and quadrature components of white noise n(t) with variance σ 2 . The algorithm can be easily extended to the case where n(t) is a colored noise. 
     The measurement equations can be written in alternative forms as follows: 
     
         y.sub.k =A sin (w.sub.d t.sub.k) cos φ+A cos (w.sub.d t.sub.k) sin φ+n.sub.ik 
    
     
         z.sub.k =A cos (w.sub.d t.sub.k) cos φ-A sin (w.sub.d t.sub.k) sin φ+n.sub.qk                                            (3) 
    
     or, with a power series expansion for the sine and cosine functions, as ##EQU3## 
     In the above approximation the terms of the order (w d  t k ) n  /n! and smaller order have been ignored (assuming here that w d  t k  &lt;n). With obvious definitions, the measurement equation can be written in a form `linear in parameters`: 
     
         Z.sub.k =θ&#39;x.sub.k +n.sub.k                          (5) 
    
     In the above, &#39; denotes transpose, Z&#39; k  =[y k  z k  ], n&#39; k  =[n ik  n qk  ], x&#39; k  denotes the observable state vector [1 t k  t k   2  . . . t k   n-1  ] and θ&#39; is the unknown parameter matrix. A standard least-squares algorithm can be applied to estimate the unknown parameter matrix θ&#39; from the sequence of noisy observations Z k , k=1, 2, . . . , N. 
     The parameter matrix θ&#39; can be estimated by either a recursive or nonrecursive form. We consider the nonrecursive form. The estimate of θ on the basis of measurement Z k , k=1,2, . . . , N, denoted θ N , is given by ##EQU4## where 0&lt;λ&lt;1 is the exponential data weighting factor. From θ N , the estimates of A, w d  and φ can be obtained. 
     The algorithm of Equation 6 requires an inverse of a symmetric (n×n) matrix once, requiring order n 2  computations. It may appear that the computation of each x j  x&#39; j  term requires n 2  computations. However, detailed examination shows that only 2n computations are required. Thus, the total number of computations is equal to 6nN+O(n 2 ). In practice, the matrix inverse can be precomputed, thus reducing the data-dependent computations to only 2nN+(n 2  /2). 
     The matrix  -1 Δ  (Σ j=1   N  x j  x&#39; j  λ N-j ) in Equation 6 has a very special structure as can easily be seen by explicit computation of the term x j  x&#39; j  of the summand. Thus, ##EQU5## 
     Each of the matrices x j  x&#39; j  and P -1  is a Hankel matrix. That is, all the elements of each crossdiagonal are the same. The structure of Hankel matrix is very similar to that of a Toeplitz matrix wherein the elements along the various subdiagonals are equal. The well-known fast algorithm disclosed in Kumar, &#34;A Fast Algorithm for Solving a Toeplitz System of Equations,&#34; IEEE Transactions, 1985, ASSP-33 (1), pp. 254-267 for the solution of a Toeplitz system of equations can be slightly modified so as to become applicable to the present problem. Thus, Equation 6 can be solved in order n (log 2  n) 2  computations, resulting in a considerable reduction in the requirement for large values of n. 
     If the matrix inverse is precomputed then, with the well-known fast algorithm referenced above, the solution for θ N  can be obtained in approximately 6n log 2  n operations. In the implementations above, it is sufficient to store only the first row and column of P or P -1 . 
     In this case only the measurements {y k  } are available and the parameter matrix θ&#39; is of dimension n×1. In such an implementation, however, there may result an ambiguity of π radians in the phase estimate if the sign of w d  is also unknown. 
     Assuming that the model Equation 5 is exact (the dimension n of the parameter matrix is Equation 4 is sufficiently high), then the substitution of Equation 5 in Equation 6 yields ##EQU6## 
     A simple manipulation of Equation 8 yields the estimation error ##EQU7## As the state vector x j  is deterministic and n j  is a zero mean process, θ N  has its mean equal to zero. The error covariance matrix of θ N  can also be evaluated in a straightforward manner. Postmultiplying Equation 9 by the transpose of θ N  and taking expected values of both sides gives ##EQU8## considering the case of λ=1 and recalling that {n j  } is a white noise sequence ##EQU9## 
     A simple approximate expression can also be obtained for the frequency estimation error when the amplitude A is known and uniform sampling is used. The frequency estimate w d  can be obtained as 
     
         w.sub.d,N ={(θ.sub.N.sup.21).sup.2 +(θ.sub.N.sup.22).sup.2 }.sup.1/2 A.sup.-1                                        (12) 
    
     When the amplitude A is also unknown, it can be replaced by its estimate given by 
     
         A.sub.N ={(θ.sub.N.sup.11).sup.2 +(θ.sub.N.sup.12).sup.2 }.sup.1/2                                                 (13) 
    
     In the above expressions, θ N   i ,j denotes the (i,j)th element of the parameter matrix θ N . The error variance of these elements of interest is given by ##EQU10## where K approaches a constant with the increase in the numbers of observations N. 
     For relatively small errors, the frequency estimation error w d ,N =w d  -w d ,N has variance approximately ##EQU11## For the case of uniform sampling t j  =jT s , where T s  is the sampling period. Substituting for t j  and letting T=NT s  denote the observation period, ##EQU12## 
     In terms of the unsampled system, if the additive noise process has one-sided noise spectral density N 0 , then σ 2  =2N 0  /T s . Thus, ##EQU13## where P=A 2  /2 is the received signal power and K has value approximately equal to 4 for low values of n. 
     We note here that in the derivation of Equation 16, the approximation error in Equation 5 has been ignored. It is difficult to estimate the error due to such finite approximation. However, from a few computer simulations, it appears that for n&gt;w d  T=(w d  T s )N, such error is small. 
     The estimator Equation 12 for w d  from the parameter matrix θ, though simple from the computational viewpoint, may not provide the minimum possible estimation error variance* for w as it ignores the other components of θ. In theory, if the number of terms in the parameter matrix is sufficiently large, one may first estimate, by direct summation, the terms A cos φ sin w d , A cos φ cos w d , A sin φsin w d , and A sin φ cos w d  without losing any optimality in the estimation, as these terms are related linearly to θ. From these estimates one obtains the estimates for φ and w d . Provided that |w d  |&lt;π/2, such an estimate for w d  is unique. The bound on w d  is ensured via a simple normalization described subsequently in Section 6 of the paper. It may be noted that the errors associated with various elements of θ are necessarily correlated. 
     Alternatively, one may resort to one of the several possible ad hoc procedures to exploit the additional information contained in the estimates w d   n , n&gt;1. In one such procedure, in estimating the magnitude of w d , one may simply compute the average of the nth roots of |w d   n  | for n=1, 2 . . . Perhaps a very systematic procedure to obtain w d  from θ is the application of the least square method to this problem. According to this method, with Y i  denoting the estimate of A|w d  | i  for i=1, 2, . . . , n-1, the estimate of |w d  | is ##EQU14## 
     To obtain w d  (including sign) we repeat the above procedure with Y i  replaced by the estimate of one of A (sin φ) w d   i  A (cos φ) w d  or an average of these two terms. As the estimate of |w d  | is independent of the phase φ, it may be more robust and thus it may be worthwhile to compute both |w d  | and w d  at the cost of some redundancy. 
     It may however be noted that the improvements provided by such modifications may not be significant, if the expectations of w d   n , are relatively small in relation to their variances, as would be the case for small values of w d . 
     To keep the dimension n of the parameter matrix small, the following estimation method is proposed. Dividing both sides of Equation 16 by Ω 2 , where Ω denotes the maximum possible value of w d , and substituting T=n/Ω, one obtains ##EQU15## 
     Selecting a value of 1/36 for the left hand side of the above equation allows one to express the maximum frequency uncertainty that can be resolved by the algorithm has a function of n. Thus, ##EQU16## The rationale for selecting the value of 1/36 for E[w d   2  ]/Ω 2 , is as follows. Since the additive noise has a Gaussian distribution, one may assume that the frequency estimation error has a Gaussian distribution with its standard deviation denoted by σ(w d ). The above selection thus ensures that 3σ(w d )&lt;Ω/2. 
     Now, for any given value of n (related to the algorithmic complexity permitted) and the P/N 0  ratio, the value of maximum possible initial frequency uncertainty Ω is calculated from the expression 18 above. The total estimation period is then divided into intervals of T 0  =nΩ -1  seconds and the numerically controlled oscillator (NCO) frequency and phase are adjusted at the end points of these intervals in the following manner. 
     Denoting by w d0  the difference frequency (frequency of the signal at the phase detector output) at t=0+, its estimate w d0  is obtained on the basis of observations y k , z k  ; 1&lt;k&lt;K 0  where K 0  =T 0  /T s  is the number of samples (assumed to be an integer) in any T 0  seconds interval and t k  =(k-1)T s . At the end of the first T 0  -seconds interval, the NCO frequency is adjusted by w d0  such that the resulting difference frequency w d1  in the time interval (T, 2T) is simply equal to w d0  =w d0 . We also note that from the above selection of T 0 , max (w d1 )&lt;0.5 max (w d0 ). The cumulative phase φ(t) of the signal at the phase detector output is then ##EQU17## Thus, if in addition to the frequency correction a phase correction Δφ=w d0  T 0  is introduced to the NCO (over an interval of T 0  -T s  to T 0  seconds), then the resulting phase is simply equal to w d1  t+φ 0  for T 0  &lt;t&lt;2T 0 . Thus, the estimation of w d1  can be obtained by processing observations y k , z k , K 0  &lt;k&lt;2K 0  using the least-squares algorithm with its initial estimate set equal to zero. It should be noted that the knowledge about the error variance of the previous estimates must also be included as a priori information in the least-squares algorithm. 
     With θ 1  denoting the parameter matrix obtained by replacing w d  by w d1  in the matrix θ defined earlier in Equations 4 and 5, its estimate θ 1 ,2K.sbsb.0 obtained on the basis of observations y k , Z k  ; 1&lt;k&lt;2K 0  is given by ##EQU18## In Equation 20 θ 1 ,2K.sbsb.0 is obtained by replacing w d  by 0 (the a priori estimate of w d1 ) in θ 0 ,K.sbsb.0 which denotes the estimate of θ 1  obtained on the basis of observations y k , z k  ; 1&lt;k&lt;K 0 . We note that the maximum product of frequency uncertainty and observation period given by max w d1 . 2T 0  is smaller than ΩT 0  &lt;n and thus the required dimension n of the parameter matrix need not be increased for the increased time interval 2T 0 . 
     The above procedure can be generalized to subsequent T 0  -seconds intervals in a straightforward manner, while the bound on frequency uncertainty and the length of observation period product remains satisfied for all t. This is so because, in view of Equation 16 above, the product E[w 2  ]T 2  in fact decreases with the increase in total observation period T. 
     From the development above it is apparent that the estimation error variance is still given by Equation 16 providing that the peak frequency uncertainty is within three times the RMS frequency estimation error. If this condition is violated, increased truncation errors may result. Also, by selecting the data weighting coefficient λ to be strictly less than 1, the algorithm can track a slowly varying input signal frequency w 0 . In this case, for example, w d1  =w d0  +Δw, where Δw represents the change in the input signal frequency w 0  over the first T 0  -seconds intervals. Such a variation is taken into account in the error covariance matrix in Equations 11 and 20 by adding an appropriate `dynamic noise covariance matrix` Q to it at each T 0  -seconds interval. 
     In an equivalent implementation, θ and hence w d  may be made recursive in each measurement y k , z k , using a standard recursive least-squares algorithm and thereby also eliminating the need for a matrix inversion. Detailed calculations of the number of computations per sample (a function of n, T, K 0  etc.) would dictate the more specific implementation so as to minimize the computational requirements. 
     In an alternative approach to keep the value of n fixed and small with an increase in the total observation period, instead of resetting the frequency reference (making a correction in the NCO frequency) the time reference is reset to zero. Subsequent observations in Equation 1 are now with respect to a different phase-reference, say φ denoted φ. At this stage the observations in the second T 0  interval are processed to have the phase reference φ and are then combined with the first set of observations. Equivalently, it is required to postmultiply the second sum on the right hand side of Equation 6, obtained for the second subinterval of T 0 , by the following matrix: ##EQU19## where: Δφ=φ-φand add the result to the corresponding sum for the first T 0  interval. The first sum on the right hand side of Equation 6 is simply multiplied by a factor of 2. This procedure is extended in an appropriate manner to subsequent intervals, so as to obtain a final estimate for w d  and φ based on the complete set of observations.