In this paper, we address the problem of frequency estimation when multiple stationary nonsinusoidal resonances oscillate about a trend in nonuniformly sampled data when the number and shape of the resonances are unknown. To solve this problem we postulate a model that relates the resonances to the data and then apply Bayesian probability theory to derive the posterior probability for the number of resonances. The calculation is implemented using simulated annealing in a Markov chain Monte Carlo simulation to draw samples from this posterior distribution. From these samples, using Monte Carlo integration, we compute the posterior probability for the resonance frequencies given the model indicators as well as a number of other posterior distributions of interest. For a single sinusoidal resonance, the Bayesian sufficient statistic is numerically equal to the Lomb‐Scargle periodogram. For a nonsinusoidal resonance this statistic is a straightforward generalization of both the discrete Fourier transform and the Lomb‐Scargle periodogram. Finally, we illustrate the calculations using data taken from two different astrophysical sources. © 2003 American Institute of Physics