A model reduction method for stable delay systems under L2 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 L2 norm of the approximation error through an unconstrained gradient-based minimization procedure. Numerical examples are used to illustrate the effectiveness of the proposed method.