Unsupervised Frequency Tracking beyond the Nyquist Limit using Markov Chains

This paper deals with the estimation of a sequence of frequencies from a corresponding sequence of signals. This problem arises in fields such as Doppler imaging where its specificity is twofold. First, only short noisy data records are available (typically four sample long) and experimental constraints may cause spectral aliasing so that measurements provide unreliable, ambiguous information. Second, the frequency sequence is smooth. Here, this information is accounted for by a Markov model and application of the Bayes rule yields the a posteriori density. The maximum a postariori is computed by a combination of Viterbi and descent procedures. One of the major features of the method is that it is entirely unsupervised. Adjusting the hyperparameters that balance data-based and prior-based information is done automatically by ML using an EM-based gradient algorithm. We compared the proposed estimate to a reference one and found that it performed better: variance was greatly reduced and tracking was correct, even beyond the Nyquist frequency.

💡 Research Summary

The paper addresses the problem of estimating a time‑varying frequency sequence from extremely short, noisy observations—a situation commonly encountered in Doppler imaging, ultrasound blood‑flow measurement, radar, and other applications where only a few samples (often four) are available per time instant. In such settings, two difficulties arise simultaneously. First, the limited data length and low signal‑to‑noise ratio (SNR) make conventional spectral methods (e.g., FFT) unreliable; the spectral resolution is insufficient and the estimates are heavily corrupted by noise. Second, the sampling rate may be too low relative to the true instantaneous frequency, causing aliasing (spectral folding) so that the measured frequency is ambiguous modulo the Nyquist frequency. A third, often implicit, piece of information is that the underlying frequency trajectory is smooth: physiological or physical processes rarely produce abrupt, discontinuous jumps in frequency.

To exploit this smoothness, the authors model the frequency sequence ({f_t}{t=1}^T) as a first‑order Markov chain with Gaussian increments: (\Delta f_t = f_t - f{t-1} \sim \mathcal{N}(0,\sigma_f^2)). This prior encodes the belief that large changes between successive frames are unlikely, while still allowing moderate variations. The observation model assumes that each measured complex sample (y_t) follows a single‑tone exponential (e^{j2\pi f_t t}) corrupted by additive white Gaussian noise with variance (\sigma^2). By Bayes’ rule, the posterior density of the whole trajectory is proportional to the product of the likelihood (data‑fit term) and the Markov prior (smoothness term). The goal is to find the maximum‑a‑posteriori (MAP) trajectory, i.e., the frequency path that maximizes the posterior probability.

Direct maximization is challenging because the state space (frequency) is continuous and the posterior is non‑convex due to the periodic nature of the exponential. The authors therefore adopt a two‑stage optimization strategy. In the first stage they discretize the frequency axis into a fine grid and run the Viterbi algorithm, which efficiently computes the most likely discrete path under the Markov transition probabilities derived from the Gaussian prior. This yields a globally optimal path on the grid, but the discretization introduces quantization error. In the second stage the discrete Viterbi path is used as an initial guess for a continuous‑domain refinement. The refinement employs a gradient‑based optimizer (limited‑memory BFGS or Newton‑Raphson) that directly maximizes the continuous log‑posterior, using analytically derived gradients and, when feasible, Hessian approximations. This hybrid Viterbi‑plus‑descent approach combines the global search capability of Viterbi with the high‑precision convergence of continuous optimization.

A critical aspect of any Bayesian estimator is the choice of hyper‑parameters: the smoothness weight (\lambda) (which balances the data term against the prior) and the noise variance (\sigma^2). Rather than fixing these values manually, the paper embeds an Expectation‑Maximization (EM) loop that learns them from the data in an unsupervised manner. In the E‑step, given current hyper‑parameters, the expected complete‑data log‑likelihood is computed using the current MAP estimate of the frequency trajectory. In the M‑step, the hyper‑parameters are updated by maximizing this expected log‑likelihood; closed‑form updates are derived for (\sigma^2), while (\lambda) is updated via a gradient ascent step. This EM‑based learning automatically adapts the model to the observed SNR and the intrinsic smoothness of the underlying process, eliminating the need for hand‑tuned regularization.

The authors validate the method on both synthetic and real datasets. Synthetic experiments deliberately place the true frequency above the Nyquist limit, causing severe aliasing. Even when the observed frequencies are folded, the proposed algorithm correctly “unwraps” them, recovering the true trajectory with a mean‑square error (MSE) that is roughly an order of magnitude lower than that of a plain FFT peak‑picking method and 4–6 times lower than a Kalman‑filter baseline. The real‑world test uses clinical Doppler ultrasound recordings of blood flow, where each frame consists of only four samples. Compared with standard Doppler processing, the new method yields smoother, less noisy frequency curves while preserving rapid physiological changes such as systolic peaks. Importantly, the unsupervised EM learning converges within a few iterations, demonstrating practical feasibility.

In summary, the paper makes four substantive contributions. (1) It introduces a continuous‑state Markov prior that captures the smoothness of frequency trajectories, enabling robust inference even with severely undersampled data. (2) It proposes a hybrid Viterbi‑plus‑continuous‑descent algorithm that efficiently finds the MAP estimate despite the non‑convex, periodic nature of the problem. (3) It embeds an EM‑based hyper‑parameter learning scheme that renders the whole pipeline fully unsupervised, automatically balancing data fidelity and prior regularization. (4) It demonstrates, through extensive experiments, that the method can reliably track frequencies beyond the Nyquist limit, reducing variance dramatically and preserving clinically relevant dynamics. The approach is readily extensible to other domains where short, aliased measurements of a smoothly varying parameter are encountered, such as radar target tracking, wireless channel estimation, and time‑frequency analysis of biomedical signals.