I wasn't going to answer this, since the question is stale. But I'm a little bit dissatisfied with Royi's answer and with the Kay algorithm as presented.
The Kay method is, at first glance, simply using the analytic signal, which is what comes out of MATLAB's hilbert() operator, and then deriving the phase increment out of that to get away from the phase wrapping problem. This phase increment is your instantaneous frequency, and if your sinusoid has mostly unchanging frequency (a true sinusoid has constant frequency), then the instantaneous frequency should be mostly DC and low-pass filtering will attenuate noise and increase the resolution and accuracy of the estimated frequency. This is a very old idea and even predates the Kay paper. I think the Rife and Boorstin paper that Kay cited might be the earliest reference. I dunno. It's in O&S or Rabiner and Schafer or Papoulis or some old text.
First, let's start with noise free input. Let's assume that the input is DC free (run it through a DC-blocking HPF if you need to).
$$\begin{align} x[n] &= A \cos(\omega_0 n + \theta) \qquad\qquad\qquad A>0, \quad 0 < \omega_0 < \pi \\ \\ &= \frac{A}{2}\left(e^{j(\omega_0 n + \theta)} + e^{-j(\omega_0 n + \theta)}\right) \\ \end{align}$$
We compute the Hilbert transform
$$\begin{align} \hat{x}[n] &= \mathscr{H}\big\{ x[n] \big\} \\ \\ &\approx \sum\limits_{i=-L}^{+L} \frac{1-(-1)^i}{\pi \ i} w[i] \ x[n-i] \end{align}$$
and then analytic signal:
$$ x_\mathrm{a}[n] \triangleq x[n] + j \hat{x}[n] $$
$L$ should be a large, positive, and odd integer. The Hilbert transformer is an FIR with $2L+1$ taps (of which $L-1$ taps have zero coefficients and need not be computed) and $w[n]$ is a good window function. If it were a Hamming Window, it would be:
$$ w[n] = \begin{cases} 0.54 \ + \ 0.46 \cdot \cos\left(\pi \frac{n}{L} \right) \qquad & |n| \le L \\ 0 & |t| > L \\ \end{cases}$$
If it were a Kaiser window it would be
$$ w[n] = \begin{cases} \frac{1}{I_0(\beta)} \, I_0\left(\beta \sqrt{1 - \left(\frac{n}{L}\right)^2 } \right) \qquad & |n| \le L \\ 0 & |t| > L \\ \end{cases}$$
where
$$ I_0(u) \triangleq \sum\limits_{k=0}^{\infty} \frac{(-1)^k \big( \tfrac{u}{2} \big)^{2k}}{(k!)^2} $$
is the zeroth-order Bessel function of the first kind. With the Kaiser window, knowing the length $L$ (like $L$ might be 31 samples, or 63 if you want it really good) you can choose $\beta$ to tradeoff between the stopband attenuation and transition bandwidth of the reconstruction brickwall filter. A good value of $\beta$ is 6 or 7, to get you about 63 dB or 72 dB of stopband attenuation.
This Hilbert transformer is not causal but it is FIR and a delay of $L$ samples will make it causal and "realizable" for real-time operation. So instead of $x_\mathrm{a}[n]$, you're really going to have $x_\mathrm{a}[n-L]$, but that delay should not be a problem if $L$ is reasonably small (not in the thousands). This means, when constructing the analytic function that, since $\hat{x}[n]$ is delayed by $L$ samples, so also must the input $x[n]$, to do this correctly. For simplicity, I will not depict that delay until maybe the very end.
Now after applying the Hilbert transform and computing the analytic function, only the positive frequency component survives:
$$\begin{align} x_\mathrm{a}[n] &= x[n] + j \hat{x}[n] \\ \\ &= \frac{A}{2} e^{j(\omega_0 n + \theta)} \\ \end{align}$$
To get the phase increment, we can compute the $\log(\cdot)$ (which involves the $\arg\{\cdot\}$) and subtract the exponents of two adjacent samples of $x_\mathrm{a}[n]$, but that will require phase unwrapping (otherwise you will get a spurious spikes of $2 \pi$ added to some phase increments. A better way to get the phase increment is simply to divide the two values and get the phase increment from that:
$$\begin{align} \frac{x_\mathrm{a}[n]}{x_\mathrm{a}[n-1]} &= \frac{\frac{A}{2} e^{j(\omega_0 n + \theta)}}{\frac{A}{2} e^{j(\omega_0 (n-1) + \theta)}} \\ \\ &= e^{j \omega_0} \end{align}$$
See how the amplitude $A$ and constant angle offset $\theta$ get killed off?
$$\begin{align} \omega_0[n] &= \arg \left\{ \frac{x_\mathrm{a}[n]}{x_\mathrm{a}[n-1]} \right\} \\ \\ &= \arg \left\{ \frac{x_\mathrm{a}[n] \, (x_\mathrm{a}[n-1])^*}{x_\mathrm{a}[n-1] \, (x_\mathrm{a}[n-1])^*} \right\} \qquad \qquad (\cdot)^* \text{ is complex conjugate}\\ \\ &= \arg \left\{ \frac{x_\mathrm{a}[n] \, (x_\mathrm{a}[n-1])^*}{\big| x_\mathrm{a}[n-1] \big|^2} \right\} \\ \\ &= \arg \Big\{ x_\mathrm{a}[n] \, (x_\mathrm{a}[n-1])^* \Big\} \\ \\ &= \arg \Big\{ (x[n] + j \hat{x}[n]) \, (x[n-1] - j \hat{x}[n-1]) \Big\} \\ \\ &= \arg \Big\{ (x[n]x[n-1] + \hat{x}[n]\hat{x}[n-1]) \, + \, j \big(\hat{x}[n]x[n-1] - x[n]\hat{x}[n-1] \big) \Big\} \\ \end{align}$$
Now I am showing the instantaneous frequency $\omega_0$ as a function of time $n$, in case the phase increment is not constant. And we know that $0 < \omega_0[n] < \pi$ so we have a single half-plane for the $\arg\{\cdot \}$ function. We'll use the second line.
$$ \arg \big\{ u+jv \big\} = \begin{cases} \arctan\left(\frac{v}{u}\right) &\text{if } u > 0, \\ \frac{\pi}{2} - \arctan\left(\frac{u}{v}\right) &\text{if } v > 0, \\ -\frac{\pi}{2} - \arctan\left(\frac{u}{v}\right) &\text{if } v < 0, \\ \arctan\left(\frac{v}{u}\right) \pm \pi &\text{if } u < 0, \\ \text{undefined} &\text{if } u = 0 \text{ and } v = 0 \end{cases} $$
So, I think, after computing the Hilbert transform $\hat{x}[n]$, we get for instantaneous frequency:
$$ \omega_0[n] = \frac{\pi}{2} - \arctan\left(\frac{x[n]x[n-1] + \hat{x}[n]\hat{x}[n-1]}{\hat{x}[n]x[n-1] - x[n]\hat{x}[n-1]}\right) $$
So that's how you get the instantaneous frequency of a single real sinusoid $x[n]$ in a noise-free environment. Note that all of this is really delayed by $L$ samples because the Hilbert transformer is delayed by that.
How it does with Gaussian noise added is getting dealt with later. Stay tuned.