In searching for gravitational waves from known pulsars if the gravitational wave signal phase is assumed to exactly follow the best fit phase evolution provided by electromagnetic observations then there is a single template for the phase evolution. However, if there is some uncertainty on the phase evolution from the electromagnetic observations, or if the gravitational wave phase evolution can deviate from the rotational phase, then searches that allow some variation in the parameters defining the phase evolution are required. In our analysis the phase parameters are included in Bayesian parameter estimation, which requires the likelihood of the data given the signal model to be evaluated. This requires recalculating the phase evolution across the entire time series of the data, which for the lengths of data being used becomes computationally intensive. We show that a reduced basis of the signal model can be calculated, and from this an empirical interpolant can allow the likelihood to be calculated far more efficiently, with only minor differences to the full likelihood calculation. This method is often referred to as a reduced order quadrature.