We have developed a new seismic array data processing method to produce slowness vector estimates and an objective measure of their uncertainties in the form of statistical confidence intervals. The slowness vector, which is typically transformed into bearing and velocity, is a key parameter used for identifying seismic phases and for event source location. Our method, multi-wavelet beamforming, is closely related to both time-domain and frequency-domain beamforming. The major advantage of multi-wavelet beamforming is that it produces multiple estimates of the slowness vector that are approximately statistically independent. First, a set of wavelet transforms is applied to the data in a manner analogous to the use of the windowed Fourier transform. Next, for each wavelet transform, we calculate semblance, a measure of signal coherence, for a range of possible slowness vectors. Then, the slowness vector estimate associated with that transform is the vector that produces the largest semblance value. The multiple slowness vector estimates can be treated as samples from a probability distribution, whose “center” we estimate using the mean, the median, and an M-estimator. Uncertainty intervals are calculated for these estimators by applying the jackknife statistical method. The intervals for the mean estimator appear to be true statistical confidence intervals, but the estimates can be biased by a directional noise field in low signal-to-noise circumstances. The median estimates are less biased by a directional noise field but sometimes underestimate the uncertainty. The M-estimator produces less-biased estimates while appearing to estimate correctly their uncertainty.