doi:10.1155/2010/467150 Research Article Multiharmonic Frequency Tracking Method Using The Sigma-Point Kalman Smoother

Several groups have proposed the state-space approach to tracking time-varying frequencies of multiharmonic quasiperiodic signals. The extended Kalman filter/smoother (EKF/EKS) is one of the common frequency tracking approaches seen in the literature. We introduce a multiharmonic frequency tracker based on the forward-backward statistical linearized Sigma-Point Kalman smoother (FBSL-SPKS) and compare its performance to that of the extended Kalman smoother (EKS). In all cases the FBSL-SPKS tracker outperformed the EKS tracker over a wide range of signal-to-noise (SNR) ratios. We also demonstrate its superior performance on real signals.


Introduction
Many natural signals contain nearly periodic rhythms with slowly varying morphologies. Example signals with this property include tremor, speech, electrocardiogram (ECG), and arterial blood pressure (ABP). In many applications the instantaneous frequency (IF) of these signals contains useful information for further analysis.
Many signal processing methods have been applied to the problem of multiharmonic frequency tracking in quasiperiodic signals. Especially, the pitch tracking in the speech signal analysis is one of the most common applications of multiharmonic frequency tracking. Pitch detection/tracking algorithms can be roughly categorized into three groups: time-domain methods such as zero-crossing, frequencydomain methods, and time-frequency-domain methods. All pitch tracking methods apply the frame-by-frame analysis due to the nature of human voice [1]. Recently Tabrikian et al. proposed the maximum a posteriori (MAP) probability pitch tracking method using harmonic model [2]. They implemented the MAP estimator by a dynamic programming procedure based on measurement collected over several frames. However, these frame-by-frame based algorithms are always not applicable especially when a local signal stationarity cannot be assumed. There are other methods that have been applied to track rhythmicity (harmonic components) in nonstationary quasiperiodic signals based on adaptive schemes [3]. The advantage of using these adaptive schemes is that one can track rhythmicity (frequencies) recursively as signal samples are acquired.
In this paper we use a Fourier series representation, which is shown in (1) Section 2.1, of multiharmonic quasiperiodic signals in which the amplitudes, phases, and frequencies are allowed to change slowly over time. The application of state space methods to continuously track the amplitudes, phases, and frequencies was pioneered by Parker and Anderson in [4] with many subsequent investigations [5][6][7][8][9]. Recently there have been several proposed methods based on particle filters [10,11] which are highly computationally intensive and hence practically intractable.
The Kalman filter (KF) recursively estimates the optimal state of a linear state space system driven by Gaussian noise by minimizing the MSE [12]. However, it cannot be applied directly to frequency tracking because our state space model has nonlinearity due to the relationship between frequencies and observed signals. There are many types of generalizations of the KF for the case of a nonlinear state space model. The extended Kalman filter (EKF) uses a local linear approximation of the model. The algorithm is relatively simple and faster than other generalizations of the KF because it relies 2 EURASIP Journal on Advances in Signal Processing on a first-order Taylor series approximation of the nonlinear system around the estimate of the current state. The Sigma-Point Kalman filter (SPKF) is another generalization to nonlinear state-space models, which includes the Unscented Kalman filter (UKF) [13], Central Difference Kalman filter (CDKF) [14], and their square-root variants [15]. Like the EKF, the SPKF approximates the state distribution by a Gaussian Random Variable. The SPKF uses a deterministic sampling approach to approximate the probability density of the state-error and noise covariances by a set of carefully chosen sample points known as sigma-points. These sigmapoints are chosen in such a way that they completely capture the mean and covariance of the corresponding densities. These sigma-points are then propagated through the true nonlinear system, with the posterior mean and covariance estimated using simple weighted averaging. This approach captures the posterior mean and covariance accurately to the 2nd order (3rd order is achieved for symmetric distributions) compared to EKF which only achieves 1storder accuracy. Another advantage of the SPKF over other Kalman generalizations is that it maintains the same order of computational complexity as the EKF.
The Kalman smoother (KS) is a noncausal version of the KF. Typically, smoothers can achieve better estimates than filters since they deal with more measurements with proper design. We proposed a tremor frequency tracking method utilizing the extended Kalman smoother (EKS) in [16,17]. However, we are unaware of any literature that investigates the estimation accuracy of smoothers in the multiharmonic frequency tracking application.
Forward-backward statistical linearized sigma-point Kalman smoother (FBSL-SPKS), which is recently proposed in [18], presents a new formulation for nonlinear smoothing using Sigma-Point Kalman filtering method. The derivation of the FBSL-SPKS is obtained by making use of the relationship between the SPKF and weighted statistical linear regression (WSLR). WSLR takes into account both the mean and covariance of the prior distribution to pseudolinearize the nonlinear dynamics. Therefore, it is more accurate than the first-order Taylor series-based linearization approach, which completely neglects the prior covariance at the point of linearization. In [18], the FBSL-SPKS is shown to obtain superior estimates than the EKS in general. To our best knowledge, however, the head-to-head performance comparison between the EKS and SPKS has not been made explicitly for the multiharmonic frequency tracking application. Especially, FBSL-SPKS has never been applied to any practical applications such as multiharmonic frequency tracking.
The first objective of our study was to implement two multiharmonic frequency trackers utilizing the EKS and FBSL-SPKS and demonstrate their feasibility of tracking the frequency of multiharmonic signals. The second objective was to compare the performance of the EKS and FBSL-SPKS trackers based on the Monte Carlo simulations and real biomedical signals. We used three performance metrics to quantify different aspects of the multiharmonic tracking performance. We only examined the smoothers since our work was focused on an offline analysis of prerecorded signals.

Methodology
We apply two nonlinear smoothing schemes using the EKF and SPKF approaches for multiharmonic frequency tracking problem. The EKF-based smoother, that is, the EKS, has many mathematically equivalent expressions. Here, we use a variant similar to that developed in [19] (see [20, page 374]). The nonlinear SPKF-based smoother was derived from the first principle in [18] and is referred as the FBSL-SPKS. The FBSL-SPKS is a fixed interval smoother, which uses two independent forward and backward filters for smoothing. The standard SPKF is used as a forward filter. The backward filter requires the inverse dynamics of the forward filter. While the EKS can easily invert the Taylor series based linearized dynamics, the SPKS requires a new approach to linearize the forward nonlinear dynamic model. There are two major variants of SPKS available in the literature which can solve this problem in a roundabout way. In [21], the inverse dynamic model was learned from the data by training a backward nonlinear predictor (e.g., neural network). The major disadvantages of this method are that it is application and data specific and requires a learning phase. Recently an Unscented Rauch-Tung-Striebel-(URTSS-) based smoother was proposed in [22], where a joint distribution of the current and future state is maintained in order to smoothen the current state. This method requires more computation due to doubling of the state dimension.
The FBSL-SPKS introduced a direct and straightforward formulation for forward-backward smoothing [18]. Instead of learning a backward dynamical model from the data, the proposed smoother (FBSL-SPKS) makes use of weighted statistical linear regression (WSLR) formulation of SPKF (see [18] for details). WSLR is a linearization technique that takes into account the uncertainty of the prior random variable when linearizing the nonlinear model. In this way, WSLR is more accurate in the statistical sense than the first-order Taylor series-based linearization employed by the EKF which only considers the mean of the prior densities while linearizing. By representing the forward nonlinear dynamics in terms of WSLR, a linear backward filter was derived from first principle in [18]. The forward and backward estimates were then statistically combined to obtain a smoothed estimate. This newly proposed FBSL-SPKS performed comparably with the smoothers presented in [21,22] but with higher computational efficiency and ease of implementation.

State Space
Model. We use 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 y n where n = 0, . . . , N represent discrete time.
Our state space model is based on the one proposed in [4] with some modifications. The measurement model is based on a Fourier series representation in which the amplitudes, phases, and frequencies are allowed to change slowly over time. It can be expressed as where m is the number of the harmonics assumed to be known, θ n the instantaneous angle, a k,n and b k,n the amplitudes of the kth harmonic sinusoidal components, y n the trend of y n , and v n is a white noise process with zeromean and variance r. The instantaneous angle θ is modeled as where f is the mean frequency, ξ n is the difference between the instantaneous frequency f n and the mean frequency f , φ n the accumulative sum of ξ n , and T s the sample interval. This is one of the major differences between our state-space model and the one proposed in [4]. This modification was necessary because the FBSL-SPKS requires the state variables to have zero mean. Since φ n is the accumulative sum of ξ n = f n − f , its mean is zero. This increases numerical stability and makes it easier to invert the model for the backward filter.
Each state-space variable was modeled as follows: where γ n is the fluctuating component in φ n , α an autoregressive (AR) process coefficient of γ n , and u ·,n mutually uncorrelated white noise processes. A value of α = 1 results in a random walk model of φ n and α = 0 results in a white noise model. The variance q of u ·,n determines how quickly the parameters are expected to change over time.
The state vector x n is defined as Then, the state-space model can be written as follows: where f (·) and h(·) are the linear state transition and nonlinear observation functions, respectively.

Forward Updates.
The filtered and predicted state estimates can be computed directly from the well-known EKF recursions, which can be found in [20]. In the recursions, the derivatives of the state transition function f n (x) and observation function h n (x) have to be computed as part of time-update and measurement-update equations, respectively. They can be expressed as follows.
(i) Derivative of f n (x) for time-update equations is (ii) Derivative of h n (x) for measurement-update equations is

EURASIP Journal on Advances in Signal Processing
. . .
The further detail of the EKF recursions can be found in [20].

Smoothing.
There are many mathematically equivalent expressions for the extended Kalman smoother (EKS). Here we use a variant similar to that developed in [19] (see [20, page 374]). The backward recursive update equations for the EKS start with initialization at time N such as where ψ is called the adjoint variable. The smoothed estimates can then be computed as follows.
(i) Backward-update equations are

SPKS Multiharmonic Frequency Tracker Recursions.
Our proposed FBSL-SPKS uses a forward-backward approach. A standard SPKF is run in the forward direction using the nonlinear model shown in (5) and (7). A backward filter then computes the estimates operating on the inverse dynamics of the forward filter. WSLR formulation as described below is used to pseudolinearize the nonlinear state-space model so that it is inverse can be computed. The forward and backward estimates are then optimally combined to generate the smoothed estimates. In order to better understand the equations of FBSL-SPKS, first we describe how SPKF performs an inherent linearization called WSLR, which considers both the mean and covariance of the prior random variable (RV) at the point of linearization.
Consider a prior RV x which is propagated through a nonlinear function g(x) to obtain a posterior RV y. Sigmapoints χ i , i = 0, 1, . . . , 2M are selected as the prior mean x plus and minus the columns of the square root of the prior covariance P x : where M is the RV dimension and γ is the composite scaling parameter. The sigma-points set χ completely captures the mean x and the covariance P x of the prior RV x: where w i is the normalized scaler weight for each sigmapoint. Each prior sigma-point is propagated through the nonlinearity to form the posterior sigma-points set γ i : The posterior statistics can then be calculated using weighted averaging of the posterior sigma-points, An alternate view is to consider the estimates arising from the sigma-point approach as a weighted statistical linearization of the nonlinear dynamics: where A and b are the statistical linearization parameters and can be determined by minimizing the expected mean square error which takes into account the uncertainty of the prior RV x. Defining J = E[ T W ] is the expected mean square error with sigma-point weighting matrix W: The true expectation can be replaced as a finite sample approximation: where the point wise linearization error is i = γ i − Aχ i − b. Now taking partial derivative on J with respect to b we obtain EURASIP Journal on Advances in Signal Processing the equation can be rewritten as follows: After cross multiplication and differentiation, (20) simplifies to Solving for b from (21) we get Substituting the value of b obtained in (22) into J and taking the partial derivative with respect to A we get Then, the equation can be rewritten as Cross multiplication and differentiation with respect to A on (24) provides where x = (x − x); solving for A from (25) we get where the prior mean (x) and covariance (P x ) are calculated in (13) from the prior sigma-points. Similarly, the posterior mean ( z) and covariances ( P z and P xz ) are calculated from the posterior sigma-points. The linearization error has zero mean and covariance P which is defined as follows Replacing P T xz = AP x from(26) From (28), P z = AP x A T + P , we observe that the covariance of the linearization error P is added when calculating the posterior covariance P z . The uncertainty feedback scheme is very important especially when there is severe nonlinearity over the uncertainty region of prior RV. First-order Taylor series-based linearization employed by EKF often diverges in highly nonlinear region as it only performs linearization around the mean of the RV but neglects this error term. In general, the WSLR technique is an optimal way of linearizing any nonlinear function in the minimum mean square error (MMSE) sense as this approach explicitly takes into account the prior RV statistics (e.g., mean and covariance).
To form the SPKF estimator, we consider the nonlinear state-space model: x n+1 = f n (x n , u n ), where x n ∈ R M is the state, y n ∈ R P is the observation at time index n, u n and v n are Gaussian distributed process and observation noises, f (·) is the nonlinear dynamic model and h(·) is the nonlinear observation model function. The process and observation noise has zero mean and covariances Q n and R n , respectively. The SPKF is then derived by recursively applying the sigma-point selection scheme shown above at every time index to these dynamic equations (see [13] for more details).
Alternatively, we may form the statistically linearized state-space using the WSLR technique: x n+1 = A f ,n x n + b f ,n + G f ,n u n + G f ,n f ,n , where A f ,n , A h,n , b f ,n , and b h,n are the statistical linearization parameters and f ,n , h,n are the linearization error with mean zero and covariance P f ,n and P h ,n . All the parameters can be obtained by applying (22) and (26) iteratively at each time index n. Deriving the KF using the linearized state-space shown in (30) also leads to SPKF (see [21]). This statistically linearized form allows to form the dynamics of the backward filter used in forward-backward smoothing approach. As the statistically linearized state space shown in (30) is different from the standard linear state space used by the Kalman filter, the detailed derivation of the FBSL-SPKS which is demonstrated in the next sections needs to be done from the first principle. The pseudocode for the FBSL-SPKS can now be specified as follows.

Forward Updates.
A standard SPKF is used as the forward filter. The task of the SPKF is to estimate x n at time index n given all past and current measurements. The SPKF recursions, which operates on the nonlinear statespace model defined in (29), are written below with WSLR. 6 EURASIP Journal on Advances in Signal Processing (i) Initialization: (ii) Calculation of sigma-points: where Λ = (L + λ)P a xn|n .

(vi) Weighted Statistical Linearization of h(·):
A h,n = P T xn+1|nyn+1|n P xn+1|n (vii) where (viii) Parameters: λ is the composite scaling parameter which is given by where w (c) i and w (m) i are the scaler sigma-point weights and they are defined as EURASIP Journal on Advances in Signal Processing 7 where α controls the size of the sigma-point distribution and should be within 0 ≤ α ≤ 1 to avoid sampling nonlocal points when the nonlinearities are strong [21]. β ≥ 0 is the weighting term which incorporates the higherorder moments of the prior distribution. As generally sigmapoints can effectively capture the first 2 moments (mean and covariance) of the distribution (for gaussian any symmetrical sigma-points set also capture the third-order moment, i.e., skewness), the parameter β also can be used to minimize the error of higher-order moments of the distribution due to sigma-point approximation effects. For Gaussian prior, β = 2 [13]. The parameter κ is used to make sure that the positive definiteness of the covariance matrices and the default choice of κ ≥ 0 should work for most of the cases. L is the dimension of the augmented state; Q and R are the process and observation noise covariances.

Backward
Updates. An information filter is used to estimate the states from the backward direction given all the present and future measurements. As the statistically linearized state-space is different from the standard linear statespace used by the Kalman filter, the time and measurement update equations had to be derived from the first principle [18]. The backward filter recursion which operates on the statistically linearized state-space shown in (30) is given as fpllows.
(i) Initializations: where S n|n = (P b xn|n ) −1 is the information matrix and z n|n = S n|n x b n|n is defined as the information state. The state estimate and error covariance matrix for the backward filter can be denoted as x b n|n and (P b xn|n ), respectively.

Synthetic Time-Variant Harmonic Signals.
We generated two sets of synthetic signals with time-variant harmonics whose sample rate was f s = 2 kHz, mean frequency f = 100 Hz, and duration 3 s using (1)-(3). The first set of synthetic signals contains the rhythmicity during the entire 3 seconds duration. The second set of synthetic signals contains the rhythmicity only during the first and last one seconds, 0-1 and 2-3 seconds. Between 1 and 2 seconds the signals are simply white Gaussian noise. The second set of synthetic signals mimics those signals whose rhythmicity is intermittent. Tables 1 and 2 list the user-specified parameters that we used for the results and examples in this paper.

Parameter Selection.
In [17] we demonstrated that the ratio (γ) between the measurement noise variance and process noise variance is the critical factor that affects the performance of the EKS frequency tracker. We used the same value for the ratio (γ = q f /r), which was 100. The other values such as q s and q y were selected based on empirical results obtained during the development of the multiharmonic frequency trackers. The user-specified parameters were chosen for the best performance of the EKS tracker [17]. We did not perform any additional tuning process for the SPKS tracker. Therefore,

Performance Measures.
There are two main issues that need to be addressed when comparing the performance of frequency trackers: accuracy and lock-on time. The accuracy quantifies how closely the tracker estimates the state. The lock-on time is a measure of how quickly the tracker can converge to the true state.
Depending on the application, the primary objective of frequency tracking may be accurate tracking of an instantaneous frequency or "signal denoising". When the rhythmicity in a given signal is intermittent, it is also important that the frequency tracker can regain its track of the intermittent instantaneous frequency as quickly as possible [10].
We used three metrics to compare the accuracy and speed of the EKS and FBSL-SPKS multiharmonic frequency trackers. The first metric is the normalized mean-squareerror (NMSE): where N is the signal duration. The second metric is normalized frequency meansquare-error (NFMSE): where f n is the instantaneous frequency (IF), f n is the estimated IF, and f is the mean IF. 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 poorer frequency tracking than a simple mean estimator and those of NFMSE 1 indicate accurate frequency tracking.
The third metric is the square-frequency-error (SFE(n)), which can be written as (48) When this metric is averaged over an ensemble of synthetic signals, it visualizes how rapidly the trackers lock on to the true frequency. In contrast to NMSE and NFMSE, SFE(n) is a function of time that shows the squared difference between the true IF and its estimate at a given time. For all of our results we calculated the NFMSE, NMSE, and SFE(n) over an ensemble of 300 synthetic signals. Figure 1 show the estimated multiharmonic frequencies using the EKS (a) and FBSL-SPKS (b) trackers on top of the spectrogram of a synthetic signal generated using (1)-(3) whose SNR was −3 dB. At 1.1 s the EKS tracker lost track of the true frequency because the estimated third harmonic started tracking the fourth harmonic of the signal. The same situation occurred toward to the end of the signal at 2.7 s. However, the FBSL-SPKS tracker never lost its track of the true IT during the entire signal duration. Two plots in Figure 2 show the spectrograms of estimation residuals using the EKS (a) and FBSL-SPKS (b) trackers. The residual spectrogram (a) in Figure 2 depicts some harmonic structures between 1.1-1.7 s and 2.7-3.0 s due to the estimation error of the EKS tracker. Figure 3(a) shows NMSE versus SNR of the EKS and FBSL-SPKS trackers. It demonstrates that the FBSL-SPKS tracker can estimate the true signal better than the EKS tracker over a wide range of SNR. Figure 3(b) depicts NFMSE versus SNR of two multiharmonic trackers. The FBSL-SPKS tracker outperformed the EKS tracker over the entire range of SNR. The performance difference is larger at low SNR values. This is probably due to a better approximation of the  state and error covariances with the sampling approach of the FBSL-SPKS as compared to the local linearization approach of the EKS. Figure 4 depicts the SFE(n) of the two multiharmonic trackers. It demonstrates that on average the FBSL-SPKS tracker can regain its track of the true IF faster than the EKS tracker.

Synthetic Signals. Two plots in
Plots in Figure 5 show the estimated instantaneous frequencies using the EKS and FBSL-SPKS trackers on top of the spectrogram of a synthetic signal whose rhythmicity is present only during 0-1 s and 2-3 s. The FBSL-SPKS tracker started tracking the true IF accurately at 2.5 s while the EKS frequency tracker failed to regain its track of the true IF. Plots in Figure 6 show the spectrograms of estimation residuals using the EKS Figure 6(a) and FBSL-SPKS Figure 6(b) trackers. The EKS tracker barely started tracking the true frequency toward the very end of the signal. However it took only 0.4 seconds for the FBSL-SPKS to start tracking the true frequency after the rhythmicity came back at 2 seconds.

Real Signal Examples.
We applied both trackers to two different types of real signals: a photosenor insect activity signal and an arterial blood pressure (ABP) signal. The photosensor insect activity signal has a clear harmonic structure, which carries important entomological information. The instantaneous frequency and the harmonic amplitudes help entomologists determine what kind of insects flew over the photosensor [23,24]. The ABP signal also has many harmonics by nature. Accurate tracking of the harmonics in the ABP signal is critical to check a patient's heart condition. However, the ABP signal can often be noisy due to signal drops and medical device interference. The following example will demonstrate that the FBSL-SPKS harmonic tracker is more robust to this type of noise than the EKS harmonic tracker.
The sampling frequency of the photosensor insect activity signal was 16 kHz and its duration was 10s. Figures 7(a) and 7(b) show the estimated harmonics using the EKS (a) and FBSL-SPKS (b) multiharmonic trackers on top of the spectrogram of a photosensor insect activity signal. Figures  8(a) and 8(b) are the spectrograms of estimation residuals using the EKS and FBSL-SPKS, respectively. The NMSE between the true and reconstructed bug signals using the FBSL-SPKS was 0.038 while that using the EKS was 0.104. The FBSL-SPKS tracker could track the harmonics during the entire time period while the EKS tracker lost its track between 2.3 s and 2.9 s, which is marked with two dark grey bars. The performance difference may not be apparent in Figures 8(a) and 8(b). However, the estimated harmonic frequencies between two grey bars in Figure 7(a) show that the slight error in fundamental frequency estimation results in the complete mismatch of higher harmonic frequency estimation. This result matches the simulation results shown in Figure 3.
The ABP signal was sampled at 500 Hz and its duration was 30 minutes. Figures 9(a) and 9(b) depict the estimated harmonics using the EKS (a) and FBSL-SPKS (b) multiharmonic trackers on top of the spectrogram of an arterial blood pressure (ABP) signal. Figures 10(a) and 10(b) are the spectrograms of estimation residuals using the EKS and FBSL-SPKS trackers, respectively. Figure 9 shows a typical example of signal drops at 25 minutes, which is common in ABP signals. While the EKS tracker could not regain its track of the right frequencies after this signal drop, the FBSL-SPKS tracker was able to regain its track. This result again demonstrated that the FBSL-SPKS harmonic tracker is more reliable than the EKS harmonic tracker.

Conclusion
We implemented the multiharmonic tracker using the recently proposed FBSL-SPKS technique and made the headto-head performance comparison between the FBSL-SPKS and EKS multiharmonic trackers based on synthetic and realworld signals. Using three difference performance metrics, we demonstrated that the FBSL-SPKS multiharmonic tracker is more accurate and robust to noise than the EKS multiharmonic tracker.