A recursive algorithm which maps a general class of nonlinear rational model, defined as the ratio of two polynomial functions, into the frequency domain is derived using the harmonic expansion method. The new algorithm provides, for the first time, a direct analytic map from the time domain rational model parameters to the higher-order frequency response functions. Complex nonlinear time domain behaviours can be analysed and interpreted in the frequency domain, and simulated examples are included to illustrate the concepts involved.