Frequency Estimation Techniques Peter J. Kootsookos p.kootsookos@ieee.org

1. Frequency Estimation Techniques Peter J. Kootsookos p.kootsookos@ieee.org comp.dsp conference

2. Frequency Estimation Techniques Talk Summary Some acknowledgements What is frequency estimation? What other problems are there? Some algorithms Maximum likelihood Subspace techniques Quinn-Fernandes Associated problems Analytic signal generation Kay / Lank-Reed-Pollon estimators Performance bounds: Cramér-Rao Lower Bound comp.dsp conference

3. Frequency Estimation Techniques Some Acknowledgements Eric Jacobson – for his presence on comp.dsp and for his work on the topic. Andrew Reilly – for his presence on comp.dsp and for analytic signal advice. Steven M. Kay – for his books on estimation and detection generally, and published research work on the topic. Barry G. Quinn – as a colleague and for his work the topic. I. Vaughan L. Clarkson – as a colleague and for his work on the topic. CRASys – Now defunct Cooperative Research Centre for Robust & Adaptive Systems. comp.dsp conference

4. Frequency Estimation Techniques Talk Summary Some acknowledgements What is frequency estimation? What other problems are there? Some algorithms Maximum likelihood Subspace techniques Quinn-Fernandes Associated problems Analytic signal generation Kay / Lank-Reed-Pollon estimators Performance bounds: Cramér-Rao Lower Bound comp.dsp conference

5. Frequency Estimation Techniques What is frequency estimation? Find the parameters A, w, f, and s2 in y(t) = A cos [w(t-n) + f)] + e(t) where t = 0..T-1, n = T-1/2 and e(t) is a noise with zero mean and variance s2. q is used to denote the vector [Awfs2 ]T. comp.dsp conference

6. Frequency Estimation Techniques What other problems are there? y(t) = A cos [w(t-n) + f)] + e(t) What about A(t) ? Estimating A(t) is envelope estimation (AM demodulation). If the variation of A(t) is slow enough, the problem of estimating w and estimating A(t) decouples. What about w(t)? This is the frequency tracking problem. What’s e(t) ? Usually assumed additive, white, & Gaussian. Maximum likelihood technique depends on Gaussian assumption. comp.dsp conference

7. Frequency Estimation Techniques What other problems are there? [continued] Amplitude-varying example: condition monitoring in rotating machinery. comp.dsp conference

8. Frequency Estimation Techniques What other problems are there? [continued] Frequency tracking example: SONAR Thanks to Barry Quinn & Ted Hannan for the plot from their book “The Estimation & Tracking of Frequency”. comp.dsp conference

9. Frequency Estimation Techniques What other problems are there? [continued] p m=1 • Multi-harmonic frequency estimation • y(t) = S Am cos [mw(t-n) + fm)] + e(t) • For periodic, but not sinusoidal, signals. • Each component is harmonically related to the fundamental frequency. comp.dsp conference

10. Frequency Estimation Techniques What other problems are there? [continued] p m=1 • Multi-tone frequency estimation • y(t) = S Am cos [wm (t-n) + fm)] + e(t) • Here, there are multiple frequency components with no relationship between the frequencies. comp.dsp conference

11. Frequency Estimation Techniques Talk Summary Some acknowledgements What is frequency estimation? What other problems are there? Some algorithms Maximum likelihood Subspace techniques Quinn-Fernandes Associated problems Analytic signal generation Kay / Lank-Reed-Pollon estimators Performance bounds: Cramér-Rao Lower Bound comp.dsp conference

12. Frequency Estimation Techniques The Maximum Likelihood Approach The likelihood function for this problem, assuming that e(t) is Gaussian is L(q) = 1/((2p)T/2Ö|Ree|)exp(–(Y –Ŷ(q))TR-1ee(Y –Ŷ(q))/ 2) where Ree = The covariance matrix of the noise e Y = [y(0) y(1) … y(T-1)]T Ŷ = [A cos(f) A cos(w + f) … A cos(w(T-1) + f)]T Y is a vector of the date samples, and Ŷ is a vector of the modeled samples. comp.dsp conference

13. Frequency Estimation Techniques The Maximum Likelihood Approach [continued] Two points to note: The functional form of the equation L(q) = 1/((2p)T/2Ö|Ree|)exp(–(Y –Ŷ(q))TR-1ee(Y –Ŷ(q))/ 2) is determined by the Gaussian distribution of the noise. If the noise is white, then the covariance matrix R is just s2I – a scaled identity matrix. comp.dsp conference

14. Frequency Estimation Techniques The Maximum Likelihood Approach [continued] Often, it is easier to deal with the log-likelihood function: ℓ (q) = –(Y –Ŷ(q))TR-1ee(Y –Ŷ(q)) where the additive constant, and multiplying constant have been ignored as they do not affect the position of the peak (unless s is zero or infinite). If the noise is also assumed to be white, the maximum likelihood problem looks like a least squares problem as maximizing the expression above is the same as minimizing (Y –Ŷ(q))T(Y –Ŷ(q)) comp.dsp conference

15. Frequency Estimation Techniques The Maximum Likelihood Approach [continued] If the complex-valued signal model is used, then estimating w is equivalent to maximizing the periodogram: P(w) =|Sy(t) exp(-i w t) |2 For the real-valued signal used here, this equivalence is only true as T tends to infinity. T-1 t=0 comp.dsp conference

16. Frequency Estimation Techniques Talk Summary Some acknowledgements What is frequency estimation? What other problems are there? Some algorithms Maximum likelihood Subspace techniques Quinn-Fernandes Associated problems Analytic signal generation Kay / Lank-Reed-Pollon estimators Performance bounds: Cramér-Rao Lower Bound comp.dsp conference

17. Frequency Estimation Techniques Subspace Techniques The peak of the spectrum produced by spectral estimators other than the periodogram can be used for frequency estimation. Signal subspace estimators use either PBar(w) = v*(w) RBar v(w) or PMV(w) = 1/( v*(w) RMV-1 v(w) ) where v(w) = [ 1 exp(iw) exp(i2w) .. exp(I(T-1)w)] and an estimate of the covariance matrix is used. Note: If Ryy is full rank, the PBar is the same as the periodogram. ^ ^ comp.dsp conference

18. Frequency Estimation Techniques Subspace Techniques - Signal Bartlett: RBar = Sl k e k e*k Minimum Variance: RMV -1 = S 1/l k e k e*k Assuming there are p frequency components. p ^ k=1 p ^ k=1 comp.dsp conference

19. Frequency Estimation Techniques Subspace Techniques - Noise Pisarenko: RPis -1 = e p+1 e*p+1 Multiple Signal Classification (MUSIC): RMUSIC -1 = S e k e*k Assuming there are p frequency components. Key Idea: The noise subspace is orthogonal to the signal subspace, so zeros of the noise subspace will indicate signal frequencies. ^ While Pisarenko is not statistically efficient, it is very fast to calculate. M ^ k=p+1 comp.dsp conference

20. Frequency Estimation Techniques Quinn-Fernandes The technique of Quinn & Fernandes assumes that the data fits the ARMA(2,2) model: y(t) – by(t-1) + y(t-2) = e(t) – ae(t-1) + e(t-2) Set a1 = 2cos(w). Filter the data to form zj (t) = y(t) + ajzj (t-1) – zj(t-2) Form bj by regressing ( zj (t) + zj (t-2) ) on zj (t-1) bj = St( zj (t) + zj (t-2) ) zj (t-1) / Stzj2(t-1) If |aj - bj | is small enough, set w = cos-1(bj / 2), otherwise set aj+1 = bj and iterate from 2. comp.dsp conference

21. Frequency Estimation Techniques Quinn-Fernandes [continued] The algorithm can be interpreted as finding the maximum of a smoothed periodogram. comp.dsp conference

22. Frequency Estimation Techniques Talk Summary Some acknowledgements What is frequency estimation? What other problems are there? Some algorithms Maximum likelihood Subspace techniques Quinn-Fernandes Associated problems Analytic signal generation Kay / Lank-Reed-Pollon estimators Performance bounds: Cramér-Rao Lower Bound comp.dsp conference

23. Frequency Estimation Techniques Associated Problems Other questions that need answering are: What happens when the signal is real-valued, and my frequency estimation technique requires a complex-valued signal? Analytic Signal generation How well can I estimate frequency? Cramer-Rao Lower Bound Threshold performance comp.dsp conference

24. Frequency Estimation Techniques Associated Problems: Analytic Signal Generation Many signal processing problems already use “analytic” signals: communications systems with “in-phase” and “quadrature” components, for example. An analytic signal, exp(i-blah), can be generated from a real-valued signal, cos(blah) , by use of the Hilbert transform: z(t) = y(t) + i H[ y(t) ] where H[.] is the Hilbert transform operation. Problems occur if the implementation of the Hilbert transform is poor. This can occur if, for example, too short an FIR filter is used. comp.dsp conference

25. Frequency Estimation Techniques Associated Problems: Analytic Signal Generation [continued] Another approach is to FFT y(t) to obtain Y(k). From Y(k), form Z(k) = 2Y(k) for k = 1 to T/2 - 1 Y(k) for k = 0 0 for k = T/2 to T and then inverse FFT Z(k) to find z(t). Unless Y(k) is interpolated, this can cause problems. Makes sure the DC term is correct. comp.dsp conference

26. Frequency Estimation Techniques Associated Problems: Analytic Signal Generation [continued] If you know something about the signal (e.g. frequency range of interest), then use of a band-pass Hilbert transforming filter is a good option. See the paper by Andrew Reilly, Gordon Fraser & Boualem Boashash, “Analytic Signal Generation : Tips & Traps” IEEE Trans. on ASSP, vol 42(11), pp3241-3245 They suggest designing a real-coefficient low-pass filter with appropriate bandwidth using a good FIR filter algorithm (e.g. Remez). The designed filter is then modulated with a complex exponential of frequency fs/4. comp.dsp conference

27. Frequency Estimation Techniques Kay’s Estimator and Related Estimators If an analytic signal, z(t), is obtained, then the simple relation: arg( z(t+1)z*(t) ) can be used to find an estimate of the frequency at time t. See this by writing: z(t+1)z*(t) = exp(i (w(t+1) + f) ) exp(-i (wt + f) ) = exp(i w) comp.dsp conference

28. Frequency Estimation Techniques Kay’s Estimator and Related Estimators [continued] What Kay did was to form an estimator w = arg( w(t) z(t+1)z*(t) ) where the weights, w(t), are chosen to minimize the mean square error. Kay found that, for very small noise w(t) = 6t(T-t) / (T(T2-1)) which is a parabolic window. T-2 S t=0 ^ comp.dsp conference

29. Frequency Estimation Techniques Kay’s Estimator and Related Estimators [continued] If the SNR is known, then it’s possible to choose an optimal set of weights. For “infinite” noise, the rectangular window is best – this is the Lank-Reed-Pollon estimator. The figure shows how the weights vary with SNR. comp.dsp conference

30. Frequency Estimation Techniques Associated Problems: Cramer-Rao Lower Bound The lower bound on the variance of unbiased estimators of the frequency a single tone in noise is var(w) >= 12s2 / (T(T2-1)A2) ^ comp.dsp conference

31. Frequency Estimation Techniques Associated Problems: Cramer-Rao Lower Bound [continued] The CRLB for the multi-harmonic case is: var(w) >= 12s2 / (T(T2-1) m2Am2) So the effective signal energy in this case is influenced by the square of the harmonic order. p S m=1 ^ comp.dsp conference

32. Frequency Estimation Techniques Associated Problems: Threshold Performance Key idea: The performance degrades when peaks in the noise spectrum exceed the peak of the frequency component. Dotted lines in the figure show the probability of this occurring. comp.dsp conference

33. Frequency Estimation Techniques Associated Problems: Threshold Performance [continued] For the multi-harmonic case, two threshold mechanisms occur: the noise outlier case and rational harmonic locking. This means that, sometimes, ½, 1/3, 2/3, 2 or 3 times the true frequency is estimated. comp.dsp conference

34. Frequency Estimation Techniques Talk Summary Some acknowledgements What is frequency estimation? What other problems are there? Some algorithms Maximum likelihood Subspace techniques Quinn-Fernandes Associated problems Analytic signal generation Kay / Lank-Reed-Pollon estimators Performance bounds: Cramér-Rao Lower Bound comp.dsp conference

35. Frequency Estimation Techniques Thanks! Thanks to Lori Ann, Al and Rick for hosting and/or organizing this get-together. comp.dsp conference

36. Frequency Estimation Techniques Good-bye! comp.dsp conference