General relativity (GR) lays the foundation for successfully explaining the current gravitational wave (GW) observations. We present a method for studying the orbital frequency evolution of GWs from binary black hole (BBH) systems based on their energy distribution on the time-frequency plane. The orbital frequency evolution of BBH systems is determined by the individual masses and spins of the component black holes and the governing gravity theory. If a beyond-GR theory of gravity governs the BBH orbital evolution for the same set of binary parameters, the time-frequency pixel energies will exhibit a different frequency evolution from what is predicted by GR. We develop a new consistency test to check whether GR explains the BBH orbital evolution. Through numerical simulation of beyond-GR theory of gravity, we demonstrate the efficiency of this new method in detecting any possible departure from GR in the framework of second-generation GW interferometers. Further, we discuss the utility of our method in probing missing physics in the GW waveform models. We apply our test to the GW190814 and GW190412 data from the LIGO-Livingston detector, assuming that the analyzing template waveform does not include higher-order modes. The lack of subdominant modes results in an incomplete representation of the GW signal, leading to systematic biases in the frequency evolution of the signal.