3 The Tracking Algorithm

For a single sample $r[n]$, there is no reason to prefer one solution for $\theta[n]$ and $\phi[n]$ over the other possible solution. However, a sequence of solutions $\ldots, \theta[n-2]$, $\theta[n-1]$, $\theta[n]$, $\ldots$, and $\ldots, \phi[n-2]$, $\phi[n-1]$, $\phi[n]$, $\ldots$, can be chosen that has the bandwidth (or spectral density, if known) expected of the underlying modulated phase.

A two-state trellis is set up. The first state corresponds to $s=1$ in Eq. (3), and the second state corresponds to $s=-1$. The sequence of solution choices will be traced through the trellis using a Viterbi algorithm. Stored at each state at time $n$ are hypothesized phase solutions $(\hat \theta[n], \hat \phi[n])$ determined from Eq. (3), as well as the instantaneous frequencies $(\hat \theta'[n], \hat \phi'[n])$, calculated by finite differences with phases traced back to time $n-1$. In addition, each state at time $n$ stores a prediction $(\hat \theta_p'[n+1], \hat \phi_p'[n+1])$ of the instantaneous frequencies at time $n+1$.

The branch metric from state $i$ at time $n-1$ to state $j$ at time $n$ is the squared difference between the predicted instantaneous frequency $\theta_p'[n]$ from state $i$ and the hypothesized solution $\theta'[n]$ from state $j$, added to the similar squared difference for $\phi'[n]$. The Viterbi algorithm operates by computing the four branch metrics at each time step, storing an accumulated metric, and tracing backwards in the trellis to find the correct solution.

An $m$th order Levinson-Durbin linear predictor is used to compute the one-step instantaneous frequency prediction. It has the form

$$\hat \theta_p'[n+1] = \sum_{i=1}^m a_i \hat \theta'[n+1-i].$$

A similar linear predictor is used for $\hat \phi'[n+1]$. The Levinson-Durbin algorithm is the linear minimum-mean squared error (LMMSE) estimator for the instantaneous frequencies. This is a standard technique, and is closely related to a Kalman filter and the EM method for parameter estimation for Markov models [9]. The coefficients $a_i$ are determined as follows. Let $r_i = E(\theta'[n]\theta'[n-i])$, let $v = (r_1, \ldots, r_m)^T$, and let

$$R_m = \left[ \begin{array}{cccc} r_0 & r_1 & \cdots & r_{m-1... ...vdots \\ r_{m-1} & r_{m-2} & \cdots & r_0 \end{array}\right].$$

Then the coefficients are given by $a = R_{k-1}^{-1} v$. There is an efficient iterative technique to determine the $m$th order coefficients from the $(m-1)$th order coefficients, so that computing the inverse of the large autocorrelation matrix is not necessary [10]. An estimate of $r_i$ may be obtained by assuming a flat power spectrum density for $\theta'$ with a sharp cutoff at $B$ Hz. Taking the inverse Fourier transform gives $r_i =$   sinc$(2 i B T_s)$, where $T_s$ is the sampling rate. For example, if the bandwidth is 4 kHz. and $T_s = 1/132300$ sec., then for a fourth order linear predictor, $v = (0.993996, 0.976115, 0.946741, 0.906506)^T$ and $a = (3.96459, -5.92481, 3.95553, -0.995421)^T$. If more is known about the spectral characteristics of $\theta'$, then the coefficients may be determined from more accurately determined $r_i$ values.