Computationally scalable and accurate estimation, prediction, and simulation of wireless communication channels is critical to the development of more adaptive transceiver algorithms. Previously, the application of autoregressive moving average (ARMA) modeling to fading processes has been complicated by ill-conditioning and nonlinear parameter estimation. This correspondence presents a numerically stable and accurate method to synthesize ARMA rational approximations of correlated Rayleigh fading processes from more complex higher order representations. Here, the problem is decomposed into autoregressive (AR) model matching followed by linear system identification. Performance is compared to that of AR, inverse discrete Fourier transform, and sum of sinusoids techniques. Also, for the first time, the finite-precision performance of different methods is compared.