  For a large Hermitian matrix $A\in \mathbb{C}^{N\times N}$, it is often the
case that the only affordable operation is matrix-vector multiplication. In
such case, randomized method is a powerful way to estimate the spectral density
(or density of states) of $A$. However, randomized methods developed so far for
estimating spectral densities only extract information from different random
vectors independently, and the accuracy is therefore inherently limited to
$\mathcal{O}(1/\sqrt{N_{v}})$ where $N_{v}$ is the number of random vectors. In
this paper we demonstrate that the "$\mathcal{O}(1/\sqrt{N_{v}})$ barrier" can
be overcome by taking advantage of the correlated information of random vectors
when properly filtered by polynomials of $A$. Our method uses the fact that the
estimation of the spectral density essentially requires the computation of the
trace of a series of matrix functions that are numerically low rank. By
repeatedly applying $A$ to the same set of random vectors and taking different
linear combination of the results, we can sweep through the entire spectrum of
$A$ by building such low rank decomposition at different parts of the spectrum.
Under some assumptions, we demonstrate that a robust and efficient
implementation of such spectrum sweeping method can compute the spectral
density accurately with $\mathcal{O}(N^2)$ computational cost and
$\mathcal{O}(N)$ memory cost. Numerical results indicate that the new method
can significantly outperform existing randomized methods in terms of accuracy.
As an application, we demonstrate a way to accurately compute a trace of a
smooth matrix function, by carefully balancing the smoothness of the integrand
and the regularized density of states using a deconvolution procedure.
