We investigate a method to incorporate signal models that allow an additional frequency harmonic in searches for gravitational waves from spinning neutron stars. We assume that emission is given by the general triaxial non-aligned model of Jones, whose waveform under certain conditions reduces to that of a biaxial precessing star, or a simple rigidly rotating triaxial aligned star. The triaxial non-aligned and biaxial models can produce emission at both the star's rotation frequency (f) and 2f, whilst the latter only emits at 2f. We have studied parameter estimation for signal models using both a set of physical source parameters and a set of waveform parameters that remove a degeneracy. We have assessed the signal detection efficiency, and used Bayesian model selection to investigate how well we can distinguish between the three models. We found that for signal-to-noise ratios (SNRs) ≳6, there is no significant loss in efficiency if performing a search for a signal at f and 2f when the source is only producing emission at 2f. However, for sources with emission at both f and 2f, signals could be missed by a search only at 2f. We also find that for a triaxial aligned source, the correct model is always favoured, but for a triaxial non-aligned source it can be hard to distinguish between the triaxial non-aligned model and the biaxial model, even at high SNR. Finally, we apply the method to a selection of known pulsars using data from the LIGO fifth science run. We give the first upper limits on gravitational wave amplitude at both f and 2f and apply the model selection criteria on real data.