For a class of large linear input-output systems, we present a new model order reduction algorithm based on general orthogonal polynomials in the time domain. The main idea of the algorithm is first to expand the unknown state variables in the space spanned by orthogonal polynomials, then the coefficient terms of polynomial expansion are calculated by a recurrence formula. The basic procedure is to use the coefficient terms to generate a projection matrix. Many classic methods with orthogonal polynomials are special cases of the general approach. The proposed approach has a good computational efficiency and preserves the stability and passivity under certain condition. Numerical experiments are reported to verify the theoretical analysis.