In this first of two papers, we present a new method for searching for oscillatory features in the primordial power spectrum. A wide variety of models predict these features in one of two different flavors: logarithmically spaced oscillations and linearly spaced oscillations. The proposed method treats the oscillations as perturbations on top of the scale-invariant power spectrum, allowing us to vary all cosmological parameters. This perturbative approach reduces the computational requirements for the search as the transfer functions and their derivatives can be precomputed. We show that the most significant degeneracy in the analysis is between the distance to last scattering and the overall amplitude at low frequencies. For models with logarithmic oscillations, this degeneracy leads to an uncertainty in the phase. For linear spaced oscillations, it affects the frequency of the oscillations. In this first of two papers, we test our code on simulated Planck-like data, and show we are able to recover fiducial input oscillations with an amplitude of a few times order 10^{-2}. We apply the code to WMAP9-year data and confirm the existence of two intriguing resonant frequencies for log spaced oscillations. For linear spaced oscillations we find a single resonance peak. We use numerical simulations to assess the significance of these features and conclude that the data do not provide compelling evidence for the existence of oscillatory features in the primordial spectrum.