A model reduction method for stable delay systems under L-2 optimality is introduced in this paper. The reduced models may take the form of either a stable finite dimensional system or a delay system with a reduced-order finite dimensional part. The two cases are studied under a unified framework. Contrary to many existing algorithms, the formulation is entirely given in terms of state-space matrices. The Routh canonical realization is used to parametrize the stable reduced models. The parameters which form the reduced order system are obtained by minimizing the weighted L-2 norm of the approximation error through an unconstrained gradient-based minimization procedure. Numerical examples are used to illustrate the effectiveness of the proposed method.