The boundary element method (BEM) is a well-established method and particularly well suited to treat wave propagation phenomena in unbounded domains. However, the occurrence of dense system matrices is prohibitive, limiting the classical BEM to small and mid-sized problems. In the present work we propose a Chebyshev interpolation based multi-level fast multipole method (FMM) to reduce memory and computational cost of the 3D elastodynamic boundary integral operators. We present two versions for the proposed algorithm: Firstly, the direct approximation of the tensorial elastodynamic displacement and traction kernels and secondly, a version using a representation of the fundamental solutions based on scalar Helmholtz kernels. The former offers easy extensibility to more complicated kernel functions, which arise for instance in poroelastic problems. The latter minimizes the number of moment-to-local (M2L) operations and, additionally, offers the possibility to exploit the rotational invariance of the scalar kernel to further reduce memory requirements. For both approaches a directional clustering scheme in combination with a plane wave modification of the kernel function is implemented to treat the high frequency case. In order to validate the proposed numerical schemes, the FMM approximation error is investigated for both the low and high frequency regime. Furthermore, convergence results are given for a Dirichlet as well as a mixed boundary value problem in Laplace domain. Finally, the applicability of the proposed FMM to transient problems treated with the Convolution Quadrature Method is investigated.