Data from gravitational wave detectors are recorded as time series that include contributions from myriad noise sources in addition to any gravitational wave signals. When regularly sampled data are available, such as for ground based and future space based interferometers, analyses are typically performed in the frequency domain, where stationary (time invariant) noise processes can be modeled very efficiently. In reality, detector noise is not stationary due to a combination of short duration noise transients and longer duration drifts in the power spectrum. This non-stationarity produces correlations across samples at different frequencies, obviating the main advantage of a frequency domain analysis. Here an alternative time-frequency approach to gravitational wave data analysis is proposed that uses discrete, orthogonal wavelet wavepackets. The time domain data is mapped onto a uniform grid of time-frequency pixels. For locally stationary noise - that is, noise with an adiabatically varying spectrum - the time-frequency pixels are uncorrelated, which greatly simplifies the calculation of quantities such as the likelihood. Moreover, the gravitational wave signals from binary systems can be compactly represented as a collection of lines in time-frequency space, resulting in a computational cost for computing waveforms and likelihoods that scales as the square root of the number of time samples, as opposed to the linear scaling for time or frequency based analyses. Key to this approach is having fast methods for computing binary signals directly in the wavelet domain. Multiple fast transform methods are developed in detail.