A state–space integer–order approximation of commensurate–order systems is obtained using a data–driven interpolation approach based on Loewner matrices. Precisely, given the values of the original fractional–order transfer function at a number of generalised frequencies, a descriptor–form state–space model matching these frequency response values is constructed from a suitable Loewner matrix pencil, as already suggested for the reduction of high–dimensional integer–order systems. Even if the stability of the resulting integer–order system cannot be guaranteed, such an approach is particularly suitable for approximating (infinite–dimensional) fractional–order systems because: (i) the order of the approximation is bounded by half the number of interpolation points, (ii) the procedure is more robust and simple than alternative approximation methods, and (iii) the procedure is fairly flexible and often leads to satisfactory results, as shown by some examples discussed at the end of the article. Clearly, the approximation depends on the location, spacing and number of the generalised interpolation frequencies but there is no particular reason to choose the interpolation frequencies on the imaginary axis, which is a natural choice in integer–order model reduction, since this axis does not correspond to the stability boundary of the original fractional–order system.