We consider Lyapunov exponents and Sacker–Sell spectrum for linear, nonautonomous retarded functional differential equations posed on an appropriate Hilbert space. A numerical method is proposed to approximate such quantities, based on the reduction to finite dimension of the evolution family associated to the system, to which a classic discrete QR method is then applied. The discretization of the evolution family is accomplished by a combination of collocation and generalized Fourier projection. A rigorous error analysis is developed to bound the difference between the computed stability spectra and the exact stability spectra. The efficacy of the results is illustrated with some numerical examples.