A nonuniform, fast Fourier transform can be used to reduce the computational cost of the empirical characteristic function (ECF) by a factor of 100. This fast ECF calculation method is applied to a new, objective, and robust method for estimating the probability distribution of univariate data, which effectively modulates and filters the ECF of a dataset in a way that yields an optimal estimate of the (Fourier transformed) underlying distribution. This improvement in computational efficiency is leveraged to estimate probability densities from a large ensemble of atmospheric velocity increments (gradients), with the purpose of characterizing the statistical and fractal properties of the velocity field. It is shown that the distribution of velocity increments depends on location in an atmospheric model and that the increments are clearly not normally distributed. The estimated increment distributions exhibit self-similar and distinctly multifractal behavior, as shown by structure functions that exhibit power-law scaling with a non-linear dependence of the power-law exponent on the structure function order.