A multipoint Routh gamma-delta canonical continued-fraction expansion for the transfer function of a linear system is derived. Based on this general form, a multifrequency Routh approximant to the system is derived by selecting the expansion points on the imaginary axis and truncating the resulting continued-fraction expansion. A connection between the stability preservation property of the multifrequency Routh approximants and the expansion points on the imaginary axis is established. Thus the multifrequency Routh approximation procedure is flexible in deriving stable reduced-order models while fitting the frequency responses and retaining the time moments and/or Markov parameters of the impulse response of the system.