The viscoelastic relaxation spectrum provides deep insights into the complex behavior of polymers. The spectrum is not directly measurable and must be recovered from oscillatory shear or relaxation stress data. The paper deals with the problem of recovery of the relaxation spectrum of linear viscoelastic materials from discrete-time noise-corrupted measurements of relaxation modulus obtained in the stress relaxation test. A class of robust algorithms of approximation of the continuous spectrum of relaxation frequencies by finite series of orthonormal functions is proposed. A quadratic identification index, which refers to the measured relaxation modulus, is adopted. Since the problem of relaxation spectrum identification is an ill-posed inverse problem, Tikhonov regularization combined with generalized cross-validation is used to guarantee the stability of the scheme. It is proved that the accuracy of the spectrum approximation depends both on measurement noises and the regularization parameter and on the proper selection of the basis functions. The series expansions using the Laguerre, Legendre, Hermite and Chebyshev functions were studied in this paper as examples. The numerical realization of the scheme by the singular value decomposition technique is discussed and the resulting computer algorithm is outlined. Numerical calculations on model data and relaxation spectrum of polydisperse polymer are presented. Analytical analysis and numerical studies proved that by choosing an appropriate model through selection of orthonormal basis functions from the proposed class of models and using a developed algorithm of least-square regularized identification, it is possible to determine the relaxation spectrum model for a wide class of viscoelastic materials. The model is smoothed and robust on measurement noises; small model approximation errors are obtained. The identification scheme can be easily implemented in available computing environments.