We present a multitaper algorithm to estimate the polarization of particle motion as a function of frequency from three-component seismic data. This algorithm is based on a singular value decomposition of a matrix of eigenspectra at a given frequency. The right complex eigenvector zˇ corresponding to the larget singular value of the matrix has the same direction as the dominant polarization of seismic motion at that frequency. The elements of the polarization vector zˇ specify the relative amplitudes and phases of motion measured along the recorded components within a chosen frequency band. The width of this frequency band is determined by the time-bandwidth product of the prolate spheroidal tapers used in the analysis. We manipulate the components of zˇ to determine the apparent azimuth and angle of incidence of seismic motion as a function of frequency. The orthogonality of the eigentapers allows one to calculate easily uncertainties in the estimated azimuth and angle of incidence. We apply this algorithm to data from the Anza Seismic Telemetered Array in the frequency band 0≤f≤30 Hz. The polarization is not always a smooth function of frequency and can exhibit sharp jumps, suggesting the existence of scattered modes within the crustal waveguide and/or receiver site resonances. ¿ American Geophysical Union 1987