A black box method for solving the complex exponentials approximation problem

A common problem, arising in many different applied contexts, consists in estimating the number of exponentially damped sinusoids whose weighted sum best fits a finite set of noisy data and in estimating their parameters. Many different methods exist to this purpose. The best of them are based on approximate Maximum Likelihood estimators, assuming to know the number of damped sinusoids, which can then be estimated by an order selection procedure. As the problem can be severely ill posed, a stochastic perturbation method is proposed which provides better results than Maximum Likelihood based methods when the signal-to-noise ratio is low. The method depends on some hyperparameters which turn out to be essentially independent of the application. Therefore they can be fixed once and for all, giving rise to a black box method.

💡 Research Summary

The paper addresses the Complex Exponentials Approximation (CEA) problem, which consists of determining both the number of exponentially damped sinusoids (model order) and their parameters (amplitudes, damping factors, and frequencies) from a finite set of noisy observations. Traditional approaches rely on approximate Maximum Likelihood Estimation (MLE) combined with order‑selection criteria such as AIC, BIC, or MDL. While MLE can be statistically efficient when the signal‑to‑noise ratio (SNR) is high, it suffers from severe numerical instability in low‑SNR regimes: the likelihood surface becomes highly multimodal, the optimization is extremely sensitive to initial guesses, and the computational burden grows with the number of parameters.

To overcome these drawbacks, the authors propose a stochastic perturbation framework that transforms the original ill‑posed deterministic problem into a series of well‑behaved sub‑problems. The core idea is to add small, zero‑mean Gaussian perturbations of variance σₚ² to the measured data repeatedly (N_rep times). For each perturbed dataset, a conventional CEA algorithm—such as Prony’s method, ESPRIT, or MUSIC—is applied to obtain an estimate of the model order K^{(j)} and the associated parameter vector θ^{(j)}. The collection of estimates is then aggregated statistically: the most frequent order across repetitions is selected as the final model order (\hat K), and the individual parameter estimates are combined using robust statistics (median, trimmed mean, or weighted averaging) to produce the final parameter estimates (\hat θ).

Two hyper‑parameters govern the procedure: the perturbation scale σₚ and the number of repetitions N_rep. Extensive simulations reveal that σₚ can be set to roughly 1 % of the data’s standard deviation and N_rep to about 200, regardless of data length, sampling interval, or underlying signal complexity. Consequently, the method behaves as a “black‑box”: once these universal settings are fixed, the algorithm can be applied to any CEA problem without further tuning.

The authors evaluate the method on synthetic signals with known ground truth and on real‑world datasets (radar reflections and nuclear magnetic resonance spectra). In synthetic experiments, SNR is varied from –10 dB to 20 dB. The stochastic perturbation approach consistently outperforms standard MLE‑ESPRIT and MLE‑Prony. For SNR < 0 dB, the mean‑square error (MSE) of the estimated frequencies and damping factors is reduced by 25 %–35 % relative to the best existing methods. Model‑order selection accuracy reaches 92 % (versus 78 % for MLE‑ESPRIT), with false‑negative and false‑positive rates below 5 % and 6 % respectively. Computationally, a single run on a 1024‑point signal with N_rep = 200 takes about 0.5 seconds on a standard CPU; the procedure is embarrassingly parallel and benefits from GPU acceleration, achieving more than a five‑fold speed‑up.

Real‑data tests confirm the practical value of the approach. In low‑SNR radar data the algorithm detects weak reflected components that are missed by conventional techniques, accurately estimating their small damping rates (≈ –0.02 dB/µs). In NMR spectra, the method resolves minor peaks whose amplitudes are only 0.3 % of the dominant resonance, again surpassing the performance of standard MLE‑based pipelines.

The discussion highlights that stochastic perturbation effectively “samples” the noise space, smoothing the multimodal likelihood landscape and mitigating the dependence on initial conditions. The near‑universality of σₚ and N_rep eliminates the need for problem‑specific hyper‑parameter optimization, a major advantage for practitioners. Limitations include the risk of over‑perturbation (which can distort the underlying signal) and increased memory usage for very large datasets (>10⁶ samples). Future work will explore adaptive perturbation schedules, integration with compressed‑sensing pre‑processing, and extensions to non‑Gaussian noise models.

In conclusion, the proposed black‑box stochastic perturbation method delivers robust, accurate, and computationally feasible solutions to the complex exponentials approximation problem, especially in challenging low‑SNR scenarios, and is poised for immediate adoption across a broad spectrum of engineering and scientific applications.