Multiharmonic Tracking Using Sigma-Point Kalman Filter
MULTIHARMONIC TRACKING USING SIGMA-POINT KALMAN FILTER Sunghan Kim1 , Anindya S. Paul2, Eric A. Wan2 , and James McNames1 1
2
Biomedical Signal Processing Laboratory, Portland State University, Portland, Oregon, U.S.A. OGI School of Science and Engineering, Oregon Health & Science University, Beaverton, Oregon, U.S.A.
ABSTRACT Several groups have proposed the state-space approach to track time-varying frequencies of multi-harmonic quasi-periodic signals contaminated with white Gaussian noise. We compared the extended Kalman filter (EKF) and sigma-point Kalman filter (SPKF) algorithms on this problem. On average, the SPKF outperformed the EKF and more accurately tracked the instantaneous frequency over a wide range of signal-to-noise (SNR) ratios.
(EKF) because they use statistical sampling to estimate the propagation of the noise covariance through system nonlinearities. The EKF relies on local linearizations. In this paper we compare the performance of these two approaches to the multi-harmonic frequency tracking problem [7]. 2. METHODOLOGY
The Kalman filter recursively estimates the state of a linear Index Terms— Extended Kalman filter, Sigma-Point Kalman stochastic process minimizing the mean squared error (MSE) filter, instantaneous frequency (IF), normalized frequency mean- [12]. The Kalman filter could not be applied directly to our application because our statistical model is nonlinear. square-error (NFMSE), and local minima issue The extended Kalman filter (EKF) is a generalized version of the Kalman filter for the case of a nonlinear state-space 1. INTRODUCTION model that use a local linear approximation of the model. The Sigma-Point Kalman filter (SPKF) is another generalization Many natural signals contain nearly periodic rhythms with to nonlinear state-space models that includes the Unscented slowly varying morphologies. Biomedical signals with this Kalman Filter (UKF) [13], Central Difference Kalman Filter property include tremor, speech, and the electrocardiogram. (CDKF) [14], and their square-root variants [15]. In many applications the instantaneous frequency of these sig-
nals contains useful diagnostic information. Many signal processing methods have been applied to the problem of frequency tracking in quasi-periodic signals. Examples include phase locked loops [1], peak tracking in spectrograms [2], linear filters [3], comb filters [4], Hilbert transforms [5], and hidden Markov models [6]. For this paper we have adopted a Fourier series representation of these signals in which the amplitudes, phases, and frequencies are allowed to change slowly over time. The application of state space methods to track the amplitudes, phases, and frequencies was pioneered in [7] with many subsequent investigations [8, 9]. Recently there have been several proposed methods based on particle filters [10, 11]. Sigma-point Kalman filters (SPKF) are generally accepted as producing better performance than extended Kalman filters This work was supported in part by the Thrasher Research Fund. Corresponding Author: Sunghan Kim Anindya S. Paul and Eric A. Wan (Professor) are with OGI School of Science and Engineering at Oregon Health and Science University Sunghan Kim and James McNames (director) are with Biomedical Signal Processing Laboratory at Portland State University Email Address: [email protected], [email protected], [email protected], [email protected] Tel.: 503.725.5390, (503) 748-1164 Fax: 503.725.3807
2.1. Filter Design 2.1.1. State Space Model We used boldface notation to denote random processes, normal face for deterministic parameters, upper case letters for matrices, lower case letters for vectors and scalars, and subscripts for time indices. The observed signal is denoted as yn where n = 0, . . . , N is the independent variable representing discrete time. We used the following observation model which is called a rectangular model [7], yn = y¯n +
m X
ak,n cos (kθn ) + bk,n sin (kθn ) + vn (1)
k=1
where m is the number of the harmonics that is assumed to be known, θn is the instantaneous phase, ak,n and bk,n are the amplitudes of the k th harmonic sinusoidal components, y¯n is the trend of yn , and vn is a white noise process with zero-mean and variance r. We modeled fluctuations in the instantaneous phase θn as a first-order approximation of an integral of the instantaneous
• Measurement-update equations:
frequency (IF) fn , θn+1 =
mod 2π {θn + 2πTs fn }
(2)
ak,n+1 = ak,n + ua,n bk,n+1 = bk,n + ub,n y¯n+1 = y¯n + uy¯,n
ˆ n|n−1 ) yˆn|n−1 = h(x ` ´ ˆ n|n = x ˆ n|n−1 + Kn yn − yˆn|n−1 x Pn|n = Pn|n−1 − Kn re,n KnT ˛ ∂fn (x) ˛˛ Fn = ∂x ˛ ˆ n|n x=x
• Time-update equations: Pn+1|n = Fn Pn|n FnT + Q ˆ n+1|n = f (x ˆ n|n ) x
ˆ n|n and These recursions produce both the filtered estimates x ˆ n|n−1 of all the state variables. the predicted estimates x 2.1.3. SPKF frequency tracker Recursions
(4)
where u·,n are mutually uncorrelated white noise processes and the variance of these processes determine how quickly the parameters are expected to change over time. The state vector xn is defined as, T xn = θn fn a1,n . . . am,n b1,n . . . bm,n y¯n .
−1 Kn = Pn|n−1 HnT re,n
(3)
where α is an autoregressive (AR) process coefficient, and un is a white noise process with zero mean and variance qf . A value of α = 1 results in a random walk model and α = 0 results in a white noise model. The coefficients ak,n , bk,n , and y¯n were modeled as random walk processes as follows,
n|n−1
re,n = r + Hn Pn|n−1 HnT
where Ts = 1/fs is the sampling interval. The modulus operator, mod 2π , has no effect on the model mathematically, but keeps θn+1 bounded to 0 ≤ θn ≤ 2π and reduces roundoff error. The instantaneous frequency (IF) fn is modeled as a firstorder autoregressive (AR) process whose mean is f¯, fn+1 = α fn − f¯ + f¯ + un
˛ ∂hn (x) ˛˛ ∂x ˛x=xˆ
Hn =
(5)
The SPKF represents the probability density of the state error covariance by a set of carefully chosen deterministic sample points. These sigma points are then propagated through the true nonlinear system, with the posterior mean and covariance calculated using simple weighted averaging. This simple approach captures the posterior mean and covariance accurately to the 2nd order1 for all nonlinearities. The basic sigma-point propagation scheme is used to form a recursive estimation approach leading to the SPKF. The SPKF recursions are, • Measurement-update equations:
Then, the state-space model can be written as follows, xn+1 =f (xn ) + un mod 2π {x1,n + 2πTs x2,n } u1,n u2,n α(x2,n − f¯n ) + f¯n u3,n x 3,n = + .. .. . . x2m+2,n u2m+2,n x2m+3,n u2m+3,n
´ ` h χxn|n−1 , χv
γn|n−1
=
Py˜n|n−1
=
2L X 2L X
` ´` ´T c wij γi,n|n−1 γj,n|n−1
=
2L X 2L X
´` ´T c ` x wij χi,n|n−1 γi,n|n−1
i=0 j=0
Px n|n−1 y n|n−1
i=0 j=0
Kn
=
Px n|n−1 y n|n−1 Py˜−1 n|n−1
(6)
yˆn|n−1
=
2L X
where f and h are the state transition and observation functions.
ˆ n|n x
=
´ ` ˆ n|n−1 + Kn yn − yˆn|n−1 x
Pn|n
=
yn =h(xn ) + vn
wim γi,n|n−1
i=0
Pn|n−1 − Kn Py˜n|n−1 KnT
• Calculate sigma-points: 2.1.2. EKF frequency tracker Recursions χan|n =
The filtered and predicted state estimates can be computed directly from the well-known EKF recursions,
1 Third
h
ˆ an|n x
ˆ an|n + ζ x
q
a Pn|n
ˆ an|n − ζ x
order is achieved for symmetric distributions.
q
a Pn|n
i
i=0 j=0
χxn+1|n ˆ n+1|n x
= =
i=0
• Parameters: xa
=
xT
χa
ˆ
=
ˆ
(χx )T
a Pn|n
=
Pn|n 4 0 0
uT
vT
˜T
(χu )T 0 Q 0
0.2
0.4 0.6 Time (s)
0.8
(χv )T
3 0 0 5 R
0.4 0.2 0.2
0.4 0.6 Time (s)
0.8
1
1
0.6 0.4 0.2 0.2
0.4 0.6 Time (s)
0.8
1
0.8 0.6 0.4 0.2 0 0
0.2
0.4 0.6 Time (s)
0.8
1
(c) Spectrogram of an estimated sig- (d) Spectrogram of a residual with nal with SPKF frequency tracker SPKF frequency tracker
˜T
Fig. 1. (a) EKF Estimated Signal (b) EKF Estimation Residual (c) SPKF Estimation Signal (d) SPKF Estimation Residual
2.2. Synthetic Multi-Harmonic Signals We generated synthetic signals with a sample rate fs = 2 kHz, mean frequency f¯ = 100 Hz, and duration 1 s using (1)–(4). At a given SNR ( dB) the variance of added noise vn was determined as, r = var [h(xn )] × 10−SNR/10
0.6
0 0
1
0.8
0 0
where ζ is scaling parameter that determines the spread of the sigma-points around the prior mean, L is the dimension of the augmented state, Q is the process-noise covariance, R is the c observation-noise covariance, and wim and wij are the scalar weights.
(7)
Tables 1 and 2 lists the user-specified parameters that we used for the results and examples in this paper. These values were selected based on empirical results obtained during the development of the trackers. 2.3. Performance Measures
0.2
1
wim χxi,n+1|n
2
0.4
0.8
(a) Spectrogram of an estimated sig- (b) Spectrogram of a residual with nal with EKF frequency tracker EKF frequency tracker
` ´ f χxn|n , χu
2L X
0.6
Frequency (kHz)
=
` x ´` ´T c wij χi,n+1|n χxj,n+1|n
1
0.8
0 0
Frequency (kHz)
Pn+1|n
2L 2L X X
1 Frequency (kHz)
• Time-update equations:
Table 2. Synthetic signal generation parameters. Name Symbol Value AR Coefficient α 0.99999 Frequency process noise variance qf 0.2 Amplitude process noise variance qs 10−5 Mean frequency f¯ 100 Hz
Frequency (kHz)
Table 1. Summary of user-specified design parameters. Name Symbol Value AR Coefficient α 0.9987 Phase process noise variance qθ 10−5 Ts Frequency process noise variance qf 10 Ts Amplitude process noise variance qs 0.0002 Ts Average process noise variance qy¯ 0.001 Ts Measurement noise variance r 0.1 Mean frequency f¯ 100 Hz where Ts is a sample interval
(ANMSE) of all harmonics’ amplitudes. The NFMSE measures the accuracy of tracking the instantaneous frequency (IF). The second metric, NMSE, measures the accuracy of estimating the “denoised” signal h(xn ). They can be expressed as follows, PN
(fi,n − fˆi,n )2 ¯2 n=1 (fi,n − f ) 2 PN ˆ n|n−1 ) n=1 h(xn ) − h(x NMSE = PN ¯n )]2 n=1 [h(xn ) − y
NFMSE = Pn=1 N
(8)
where fi,n is the true IF, fˆi,n is the estimated IF, f¯ is the mean IF, and N is a number of samples. NFMSE has a natural scale ranging from 0 to 1. A value NFMSE = 1 means that the average accuracy of the estimated IF is no better than simply using the mean IF as an estimate. Values of NFMSE > 1 indicate poor frequency tracking and those of NFMSE ≪ 1 indicate accurate frequency tracking.
We used two metrics to measure the EKF and SPKF frequency trackers’ performance. They are normalized frequency meanWe calculated the NFMSE and NMSE over 300 simulasquare-error (NFMSE), conventional normalized mean-squaretions. error (NMSE), and averaged normalized mean-square-error
1
1
1 EKF SPKF
0.8
0.6
0.6
NMSE
NFMSE
0.8
0.4
0.2
EKF SPKF Averaged NMSE of Amplitudes
EKF SPKF
0.4
0.2
0 0
2
4
6
8
10 12 SNR (dB)
14
16
18
20
0 0
2
(a) NMSE versus SNR
4
6
8
10 12 SNR (dB)
14
(b) NFMSE versus SNR
16
18
20
0.8
0.6
0.4
0.2
0 0
2
4
6
8
10 12 SNR (dB)
14
16
18
20
(c) ANMSE versus SNR
Fig. 2. (a) NMSE versus SNR (b) NFMSE versus SNR. (c) NFMSE versus SNR. In the plots the shaded regions represent the 5th and 95th percentile ranges of the NMSE,NFMSE, and averaged NMSE, respectively. 3. RESULTS AND DISCUSSION Four plots in Fig. 1 show the spectrograms of estimated signals (two on the left) and one-step-ahead prediction residuals (two on the right) using the EKF and SPKF frequency trackers. There isn’t much difference between two spectrograms of the SPKF and EKF estimation residuals.
1
Frequency (kHz)
0.8
0.6
0.4
0.2
0 0
0.2
0.4
0.6
0.8
1
Time (s)
(a) Spectrogram and Estimated Frequency with EKF 1
Frequency (kHz)
0.8
0.6
In many cases the tracker will lock on to a “subharmonic” or harmonic of the fundamental frequency [7]. Fig. 3 (SNR= 0 dB) shows an example of this where the EKF frequency tracker locks on to the wrong frequency while the SPKF frequency tracker does not.
0.4
0.2
0 0
Fig. 2(a) shows NMSE versus SNR of the EKF and SPKF frequency trackers. It demonstrates that the power of the residuals of two frequency trackers are comparable in terms of NMSE. Fig. 2(b) depicts NFMSE versus SNR of two frequency trackers. The SPKF frequency tracker performed better than the EKF frequency tracker on average over the entire range of SNR. This is probably due to a better approximation of the state error covariance matrix with the sampling approach of the SPKF, as compared to the local linearization approach of the EKF. Fig. 2(c) illustrates the averaged NMSE of all harmonics’ amplitudes. Although there wasn’t much performance difference between two filters in terms of the NMSE, the averaged NMSE of all harmonics’ amplitudes demonstrates that the SPKF can estimate the amplitude of each harmonic component more accurately than the EKF.
0.2
0.4
0.6
0.8
1
Time (s)
(b) Spectrogram and Estimated Frequency with SPKF
Fig. 3. (a) Erratic EKF frequency tracker (b) Successful SPKF frequency tracker. In both plots the white lines are the estimated harmonics.
4. CONCLUSION We designed the EKF and SPKF frequency trackers which can track harmonically related frequencies and their amplitudes using a state-space model proposed in [7]. The SPKF performed better than the EKF on average over a wide range of SNR.
5. REFERENCES [1] D. R. Polk and S. C. Gupta, “An approach to the analysis of performance of quasi-optimum digital phase-locked loops,” IEEE Transactions on Communications, vol. 21, no. 6, pp. 733–738, Jun 1973. [2] Robert J. McAulary and Thomas F. Quatieri, “Speech analysis/synthesis based on a sinusoidal representation,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. ASSP-34, no. 4, pp. 744–754, Aug 1986. [3] Petr Tichavsk´y and Arye Nehorai, “Comparative study of four adaptive frequency trackers,” IEEE Transactions on Signal Processing, vol. 45, no. 6, pp. 1473–1484, Jun 1997. [4] Arye Nehorai and Boaz Porat, “Adaptive comb filtering for harmonic signal enhancement,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 34, no. 5, pp. 1124–1138, Oct. 1986. [5] Boualem Boashash, “Estimating and interpreting the instantaneous frequency of a signal—Part 1: Fundamentals,” Proceedings of the IEEE, vol. 80, no. 4, pp. 520– 538, Apr 1992. [6] J Tabrikian, S Dubnov, and Y Dickalov, “Maximum a-posteriori probability pitch tracking in noisy environments using harmonic model,” IEEE TRANSACTIONS ON SPEECH AND AUDIO PROCESSING, vol. 12, no. 1, pp. 76–87, Jan 2004. [7] Philip J. Parker and Brian D.O. Anderson, “Frequency tracking of nonsinusoidal periodic signals in noise,” Signal Processing, vol. 20, no. 2, pp. 127–152, Jun 1990. [8] Barbara F. La Scala, Robert R. Bitmead, and Barry G. Quinn, “An exteded Kalman filter frequency tracker for high-noise environments,” IEEE Transactions on Signal Processing, vol. 44, no. 2, pp. 431–434, Feb. 1996. [9] Sergio Bittanti and Sergio M. Savaresi, “On the parameterization and design of an extended Kalman filter frequency tracker,” IEEE Transactions on Automatic Control, vol. 45, no. 9, pp. 1718–1724, Sep 2000. [10] E. Fischler and B. Z. Bobrovsky, “Mean time to loose lock of phase tracking by particle filtering,” Signal Processing, vol. 86, pp. 3481–3485, 2006. [11] Corentin Dubois and Manuel Davy, “Joint detection and tracking of time-varying harmonic components: A flexible bayesian approach,” IEEE Transactions on Audio, Speech, and Language Processing, vol. 15, no. 4, pp. 1283–1295, May 2007.
[12] R. E. Kalman, “A new approach to linear filtering and prediction problems,” Transactions of the ASME - Journal of Basic Engineering, vol. 82, pp. 35–45, 1960. [13] S. J. Julier and J. K. Uhlmann, “Unscented filtering and nonlinear estimation,” in Processings of The IEEE, March 2004, vol. 92, pp. 401–422. [14] M. Norgaard, N. Poulsen, and O. Ravn, “New developments in state estimation for nonlinear systems,” Automatica, vol. 36, pp. 1627–1638, Nov 2000. [15] R. van der Merwe and E. Wan, “The square-root unscented kalman filter for state and parameter estimation,” in Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), May 2001, vol. 6, pp. 3461–3464.
2
Biomedical Signal Processing Laboratory, Portland State University, Portland, Oregon, U.S.A. OGI School of Science and Engineering, Oregon Health & Science University, Beaverton, Oregon, U.S.A.
ABSTRACT Several groups have proposed the state-space approach to track time-varying frequencies of multi-harmonic quasi-periodic signals contaminated with white Gaussian noise. We compared the extended Kalman filter (EKF) and sigma-point Kalman filter (SPKF) algorithms on this problem. On average, the SPKF outperformed the EKF and more accurately tracked the instantaneous frequency over a wide range of signal-to-noise (SNR) ratios.
(EKF) because they use statistical sampling to estimate the propagation of the noise covariance through system nonlinearities. The EKF relies on local linearizations. In this paper we compare the performance of these two approaches to the multi-harmonic frequency tracking problem [7]. 2. METHODOLOGY
The Kalman filter recursively estimates the state of a linear Index Terms— Extended Kalman filter, Sigma-Point Kalman stochastic process minimizing the mean squared error (MSE) filter, instantaneous frequency (IF), normalized frequency mean- [12]. The Kalman filter could not be applied directly to our application because our statistical model is nonlinear. square-error (NFMSE), and local minima issue The extended Kalman filter (EKF) is a generalized version of the Kalman filter for the case of a nonlinear state-space 1. INTRODUCTION model that use a local linear approximation of the model. The Sigma-Point Kalman filter (SPKF) is another generalization Many natural signals contain nearly periodic rhythms with to nonlinear state-space models that includes the Unscented slowly varying morphologies. Biomedical signals with this Kalman Filter (UKF) [13], Central Difference Kalman Filter property include tremor, speech, and the electrocardiogram. (CDKF) [14], and their square-root variants [15]. In many applications the instantaneous frequency of these sig-
nals contains useful diagnostic information. Many signal processing methods have been applied to the problem of frequency tracking in quasi-periodic signals. Examples include phase locked loops [1], peak tracking in spectrograms [2], linear filters [3], comb filters [4], Hilbert transforms [5], and hidden Markov models [6]. For this paper we have adopted a Fourier series representation of these signals in which the amplitudes, phases, and frequencies are allowed to change slowly over time. The application of state space methods to track the amplitudes, phases, and frequencies was pioneered in [7] with many subsequent investigations [8, 9]. Recently there have been several proposed methods based on particle filters [10, 11]. Sigma-point Kalman filters (SPKF) are generally accepted as producing better performance than extended Kalman filters This work was supported in part by the Thrasher Research Fund. Corresponding Author: Sunghan Kim Anindya S. Paul and Eric A. Wan (Professor) are with OGI School of Science and Engineering at Oregon Health and Science University Sunghan Kim and James McNames (director) are with Biomedical Signal Processing Laboratory at Portland State University Email Address: [email protected], [email protected], [email protected], [email protected] Tel.: 503.725.5390, (503) 748-1164 Fax: 503.725.3807
2.1. Filter Design 2.1.1. State Space Model We used boldface notation to denote random processes, normal face for deterministic parameters, upper case letters for matrices, lower case letters for vectors and scalars, and subscripts for time indices. The observed signal is denoted as yn where n = 0, . . . , N is the independent variable representing discrete time. We used the following observation model which is called a rectangular model [7], yn = y¯n +
m X
ak,n cos (kθn ) + bk,n sin (kθn ) + vn (1)
k=1
where m is the number of the harmonics that is assumed to be known, θn is the instantaneous phase, ak,n and bk,n are the amplitudes of the k th harmonic sinusoidal components, y¯n is the trend of yn , and vn is a white noise process with zero-mean and variance r. We modeled fluctuations in the instantaneous phase θn as a first-order approximation of an integral of the instantaneous
• Measurement-update equations:
frequency (IF) fn , θn+1 =
mod 2π {θn + 2πTs fn }
(2)
ak,n+1 = ak,n + ua,n bk,n+1 = bk,n + ub,n y¯n+1 = y¯n + uy¯,n
ˆ n|n−1 ) yˆn|n−1 = h(x ` ´ ˆ n|n = x ˆ n|n−1 + Kn yn − yˆn|n−1 x Pn|n = Pn|n−1 − Kn re,n KnT ˛ ∂fn (x) ˛˛ Fn = ∂x ˛ ˆ n|n x=x
• Time-update equations: Pn+1|n = Fn Pn|n FnT + Q ˆ n+1|n = f (x ˆ n|n ) x
ˆ n|n and These recursions produce both the filtered estimates x ˆ n|n−1 of all the state variables. the predicted estimates x 2.1.3. SPKF frequency tracker Recursions
(4)
where u·,n are mutually uncorrelated white noise processes and the variance of these processes determine how quickly the parameters are expected to change over time. The state vector xn is defined as, T xn = θn fn a1,n . . . am,n b1,n . . . bm,n y¯n .
−1 Kn = Pn|n−1 HnT re,n
(3)
where α is an autoregressive (AR) process coefficient, and un is a white noise process with zero mean and variance qf . A value of α = 1 results in a random walk model and α = 0 results in a white noise model. The coefficients ak,n , bk,n , and y¯n were modeled as random walk processes as follows,
n|n−1
re,n = r + Hn Pn|n−1 HnT
where Ts = 1/fs is the sampling interval. The modulus operator, mod 2π , has no effect on the model mathematically, but keeps θn+1 bounded to 0 ≤ θn ≤ 2π and reduces roundoff error. The instantaneous frequency (IF) fn is modeled as a firstorder autoregressive (AR) process whose mean is f¯, fn+1 = α fn − f¯ + f¯ + un
˛ ∂hn (x) ˛˛ ∂x ˛x=xˆ
Hn =
(5)
The SPKF represents the probability density of the state error covariance by a set of carefully chosen deterministic sample points. These sigma points are then propagated through the true nonlinear system, with the posterior mean and covariance calculated using simple weighted averaging. This simple approach captures the posterior mean and covariance accurately to the 2nd order1 for all nonlinearities. The basic sigma-point propagation scheme is used to form a recursive estimation approach leading to the SPKF. The SPKF recursions are, • Measurement-update equations:
Then, the state-space model can be written as follows, xn+1 =f (xn ) + un mod 2π {x1,n + 2πTs x2,n } u1,n u2,n α(x2,n − f¯n ) + f¯n u3,n x 3,n = + .. .. . . x2m+2,n u2m+2,n x2m+3,n u2m+3,n
´ ` h χxn|n−1 , χv
γn|n−1
=
Py˜n|n−1
=
2L X 2L X
` ´` ´T c wij γi,n|n−1 γj,n|n−1
=
2L X 2L X
´` ´T c ` x wij χi,n|n−1 γi,n|n−1
i=0 j=0
Px n|n−1 y n|n−1
i=0 j=0
Kn
=
Px n|n−1 y n|n−1 Py˜−1 n|n−1
(6)
yˆn|n−1
=
2L X
where f and h are the state transition and observation functions.
ˆ n|n x
=
´ ` ˆ n|n−1 + Kn yn − yˆn|n−1 x
Pn|n
=
yn =h(xn ) + vn
wim γi,n|n−1
i=0
Pn|n−1 − Kn Py˜n|n−1 KnT
• Calculate sigma-points: 2.1.2. EKF frequency tracker Recursions χan|n =
The filtered and predicted state estimates can be computed directly from the well-known EKF recursions,
1 Third
h
ˆ an|n x
ˆ an|n + ζ x
q
a Pn|n
ˆ an|n − ζ x
order is achieved for symmetric distributions.
q
a Pn|n
i
i=0 j=0
χxn+1|n ˆ n+1|n x
= =
i=0
• Parameters: xa
=
xT
χa
ˆ
=
ˆ
(χx )T
a Pn|n
=
Pn|n 4 0 0
uT
vT
˜T
(χu )T 0 Q 0
0.2
0.4 0.6 Time (s)
0.8
(χv )T
3 0 0 5 R
0.4 0.2 0.2
0.4 0.6 Time (s)
0.8
1
1
0.6 0.4 0.2 0.2
0.4 0.6 Time (s)
0.8
1
0.8 0.6 0.4 0.2 0 0
0.2
0.4 0.6 Time (s)
0.8
1
(c) Spectrogram of an estimated sig- (d) Spectrogram of a residual with nal with SPKF frequency tracker SPKF frequency tracker
˜T
Fig. 1. (a) EKF Estimated Signal (b) EKF Estimation Residual (c) SPKF Estimation Signal (d) SPKF Estimation Residual
2.2. Synthetic Multi-Harmonic Signals We generated synthetic signals with a sample rate fs = 2 kHz, mean frequency f¯ = 100 Hz, and duration 1 s using (1)–(4). At a given SNR ( dB) the variance of added noise vn was determined as, r = var [h(xn )] × 10−SNR/10
0.6
0 0
1
0.8
0 0
where ζ is scaling parameter that determines the spread of the sigma-points around the prior mean, L is the dimension of the augmented state, Q is the process-noise covariance, R is the c observation-noise covariance, and wim and wij are the scalar weights.
(7)
Tables 1 and 2 lists the user-specified parameters that we used for the results and examples in this paper. These values were selected based on empirical results obtained during the development of the trackers. 2.3. Performance Measures
0.2
1
wim χxi,n+1|n
2
0.4
0.8
(a) Spectrogram of an estimated sig- (b) Spectrogram of a residual with nal with EKF frequency tracker EKF frequency tracker
` ´ f χxn|n , χu
2L X
0.6
Frequency (kHz)
=
` x ´` ´T c wij χi,n+1|n χxj,n+1|n
1
0.8
0 0
Frequency (kHz)
Pn+1|n
2L 2L X X
1 Frequency (kHz)
• Time-update equations:
Table 2. Synthetic signal generation parameters. Name Symbol Value AR Coefficient α 0.99999 Frequency process noise variance qf 0.2 Amplitude process noise variance qs 10−5 Mean frequency f¯ 100 Hz
Frequency (kHz)
Table 1. Summary of user-specified design parameters. Name Symbol Value AR Coefficient α 0.9987 Phase process noise variance qθ 10−5 Ts Frequency process noise variance qf 10 Ts Amplitude process noise variance qs 0.0002 Ts Average process noise variance qy¯ 0.001 Ts Measurement noise variance r 0.1 Mean frequency f¯ 100 Hz where Ts is a sample interval
(ANMSE) of all harmonics’ amplitudes. The NFMSE measures the accuracy of tracking the instantaneous frequency (IF). The second metric, NMSE, measures the accuracy of estimating the “denoised” signal h(xn ). They can be expressed as follows, PN
(fi,n − fˆi,n )2 ¯2 n=1 (fi,n − f ) 2 PN ˆ n|n−1 ) n=1 h(xn ) − h(x NMSE = PN ¯n )]2 n=1 [h(xn ) − y
NFMSE = Pn=1 N
(8)
where fi,n is the true IF, fˆi,n is the estimated IF, f¯ is the mean IF, and N is a number of samples. NFMSE has a natural scale ranging from 0 to 1. A value NFMSE = 1 means that the average accuracy of the estimated IF is no better than simply using the mean IF as an estimate. Values of NFMSE > 1 indicate poor frequency tracking and those of NFMSE ≪ 1 indicate accurate frequency tracking.
We used two metrics to measure the EKF and SPKF frequency trackers’ performance. They are normalized frequency meanWe calculated the NFMSE and NMSE over 300 simulasquare-error (NFMSE), conventional normalized mean-squaretions. error (NMSE), and averaged normalized mean-square-error
1
1
1 EKF SPKF
0.8
0.6
0.6
NMSE
NFMSE
0.8
0.4
0.2
EKF SPKF Averaged NMSE of Amplitudes
EKF SPKF
0.4
0.2
0 0
2
4
6
8
10 12 SNR (dB)
14
16
18
20
0 0
2
(a) NMSE versus SNR
4
6
8
10 12 SNR (dB)
14
(b) NFMSE versus SNR
16
18
20
0.8
0.6
0.4
0.2
0 0
2
4
6
8
10 12 SNR (dB)
14
16
18
20
(c) ANMSE versus SNR
Fig. 2. (a) NMSE versus SNR (b) NFMSE versus SNR. (c) NFMSE versus SNR. In the plots the shaded regions represent the 5th and 95th percentile ranges of the NMSE,NFMSE, and averaged NMSE, respectively. 3. RESULTS AND DISCUSSION Four plots in Fig. 1 show the spectrograms of estimated signals (two on the left) and one-step-ahead prediction residuals (two on the right) using the EKF and SPKF frequency trackers. There isn’t much difference between two spectrograms of the SPKF and EKF estimation residuals.
1
Frequency (kHz)
0.8
0.6
0.4
0.2
0 0
0.2
0.4
0.6
0.8
1
Time (s)
(a) Spectrogram and Estimated Frequency with EKF 1
Frequency (kHz)
0.8
0.6
In many cases the tracker will lock on to a “subharmonic” or harmonic of the fundamental frequency [7]. Fig. 3 (SNR= 0 dB) shows an example of this where the EKF frequency tracker locks on to the wrong frequency while the SPKF frequency tracker does not.
0.4
0.2
0 0
Fig. 2(a) shows NMSE versus SNR of the EKF and SPKF frequency trackers. It demonstrates that the power of the residuals of two frequency trackers are comparable in terms of NMSE. Fig. 2(b) depicts NFMSE versus SNR of two frequency trackers. The SPKF frequency tracker performed better than the EKF frequency tracker on average over the entire range of SNR. This is probably due to a better approximation of the state error covariance matrix with the sampling approach of the SPKF, as compared to the local linearization approach of the EKF. Fig. 2(c) illustrates the averaged NMSE of all harmonics’ amplitudes. Although there wasn’t much performance difference between two filters in terms of the NMSE, the averaged NMSE of all harmonics’ amplitudes demonstrates that the SPKF can estimate the amplitude of each harmonic component more accurately than the EKF.
0.2
0.4
0.6
0.8
1
Time (s)
(b) Spectrogram and Estimated Frequency with SPKF
Fig. 3. (a) Erratic EKF frequency tracker (b) Successful SPKF frequency tracker. In both plots the white lines are the estimated harmonics.
4. CONCLUSION We designed the EKF and SPKF frequency trackers which can track harmonically related frequencies and their amplitudes using a state-space model proposed in [7]. The SPKF performed better than the EKF on average over a wide range of SNR.
5. REFERENCES [1] D. R. Polk and S. C. Gupta, “An approach to the analysis of performance of quasi-optimum digital phase-locked loops,” IEEE Transactions on Communications, vol. 21, no. 6, pp. 733–738, Jun 1973. [2] Robert J. McAulary and Thomas F. Quatieri, “Speech analysis/synthesis based on a sinusoidal representation,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. ASSP-34, no. 4, pp. 744–754, Aug 1986. [3] Petr Tichavsk´y and Arye Nehorai, “Comparative study of four adaptive frequency trackers,” IEEE Transactions on Signal Processing, vol. 45, no. 6, pp. 1473–1484, Jun 1997. [4] Arye Nehorai and Boaz Porat, “Adaptive comb filtering for harmonic signal enhancement,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 34, no. 5, pp. 1124–1138, Oct. 1986. [5] Boualem Boashash, “Estimating and interpreting the instantaneous frequency of a signal—Part 1: Fundamentals,” Proceedings of the IEEE, vol. 80, no. 4, pp. 520– 538, Apr 1992. [6] J Tabrikian, S Dubnov, and Y Dickalov, “Maximum a-posteriori probability pitch tracking in noisy environments using harmonic model,” IEEE TRANSACTIONS ON SPEECH AND AUDIO PROCESSING, vol. 12, no. 1, pp. 76–87, Jan 2004. [7] Philip J. Parker and Brian D.O. Anderson, “Frequency tracking of nonsinusoidal periodic signals in noise,” Signal Processing, vol. 20, no. 2, pp. 127–152, Jun 1990. [8] Barbara F. La Scala, Robert R. Bitmead, and Barry G. Quinn, “An exteded Kalman filter frequency tracker for high-noise environments,” IEEE Transactions on Signal Processing, vol. 44, no. 2, pp. 431–434, Feb. 1996. [9] Sergio Bittanti and Sergio M. Savaresi, “On the parameterization and design of an extended Kalman filter frequency tracker,” IEEE Transactions on Automatic Control, vol. 45, no. 9, pp. 1718–1724, Sep 2000. [10] E. Fischler and B. Z. Bobrovsky, “Mean time to loose lock of phase tracking by particle filtering,” Signal Processing, vol. 86, pp. 3481–3485, 2006. [11] Corentin Dubois and Manuel Davy, “Joint detection and tracking of time-varying harmonic components: A flexible bayesian approach,” IEEE Transactions on Audio, Speech, and Language Processing, vol. 15, no. 4, pp. 1283–1295, May 2007.
[12] R. E. Kalman, “A new approach to linear filtering and prediction problems,” Transactions of the ASME - Journal of Basic Engineering, vol. 82, pp. 35–45, 1960. [13] S. J. Julier and J. K. Uhlmann, “Unscented filtering and nonlinear estimation,” in Processings of The IEEE, March 2004, vol. 92, pp. 401–422. [14] M. Norgaard, N. Poulsen, and O. Ravn, “New developments in state estimation for nonlinear systems,” Automatica, vol. 36, pp. 1627–1638, Nov 2000. [15] R. van der Merwe and E. Wan, “The square-root unscented kalman filter for state and parameter estimation,” in Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), May 2001, vol. 6, pp. 3461–3464.
