We develop and apply a new technique to determine array-averaged particle motions from three-component seismic array data. The method is based on multi-wavelets, which are an extension of multi-taper spectral methods, and is a hybrid of Fourier and time-domain methods of array processing. Particle motions are determined by a time-domain principal-component method. A complex singular value decomposition is used on wavelet transformed signals assembled into multiple matrices (one for each wavelet). The eigenvector of the largest singular value of each matrix is used to estimate the phase between individual signals. We determine the relative phase between components to estimate an average particle motion ellipse for the array. The estimation procedure is made more stable by the redundancy inherent in the multi-wavelets and by M-estimators applied to individual phase factors in the complex plane. The method is applied to data from three-component array experiments conducted at Piñon Flats, California, in 1990 and 1991. We find remarkable departures of P-wave particle motions from the pure longitudinal motion expected for an isotropic media. Anomalies as large as 40° are measured from some azimuths. The azimuthally varying particle-motion anomalies are frequency dependent, generally increasing in magnitude as frequency increases. Borehole measurements from sensors at 153 and 274 m depth below the array show a pattern indistinguishable from the surface sensors. The data are fit with a dipping, transversely isotropic medium with a symmetry plane having a strike of 70° and a dip of 30° to the northwest. We attribute our results to three superimposed effects: (1) an anisotropy of the near surface induced by preferential weathering of the grandiorite bedrock along joints, (2) a larger scale anisotropy induced by structural and intrinsic anisotropy related to the Santa Rosa mylonite, and (3) near-surface scattering.