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 also 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 \\ \\ &= \tfrac{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 filter with $2L+1$ taps (of which $L$ interlaced 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 & |n| > 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 & |n| > 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 ripple and transition bandwidth around DC of the Hilbert transformer FIR. A good value of $\beta$ is 6 or 7, to get you about -63 dB or -72 dB of ripple around DC.
This Hilbert transformer is not causal but this implementation 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 (and gets doubled):
$$\begin{align} x_\mathrm{a}[n] &= x[n] + j \hat{x}[n] \\ \\ &= A \, 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{A e^{j(\omega_0 n + \theta)}}{A 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] = \tfrac{\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.
So I'm gonna rob a Wikipedia graphic and make lotsa assumptions that I think are prudent. We'll say that $q[n]$ is additive Nyquist-bandlimited white Gaussian noise with zero mean and $\overline{q^2[n]}$ variance.
So it's added:
$$\begin{align} x[n] &= A \cos(\omega_0 n + \theta) + q[n] \qquad\qquad\qquad A>0, \quad 0 < \omega_0 < \pi \\ \\ &= \tfrac{A}{2}\left(e^{j(\omega_0 n + \theta)} + e^{-j(\omega_0 n + \theta)}\right) + q[n] \\ \end{align}$$
The Hilbert transform
$$ \hat{x}[n] = A \sin(\omega_0 n + \theta) + \hat{q}[n] $$
where $\hat{q}[n] = \mathscr{H}\big\{ q[n] \big\}$ and the analytic signal is
$$\begin{align} x_\mathrm{a}[n] &= x[n] + j \hat{x}[n] \\ \\ &= A \, e^{j(\omega_0 n + \theta)} + q[n] + j \hat{q}[n] \\ \end{align}$$
While $q[n]$ and $\hat{q}[n]$ are certainly coupled, they're both still Gaussian and white and have virtually the same energy or variance because the Hilbert transformer has 0 dB gain for all frequencies except DC and Nyquist.
$$ \overline{\hat{q}^2[n]} = \overline{q^2[n]}$$
So, I think it's gonna look a bit like this figure
and, assuming a decent S/N ratio,
$$\begin{align} x_\mathrm{a}[n] &= x[n] + j \hat{x}[n] \\ \\ &= A \, e^{j(\omega_0 n + \theta)} + q[n] + j \hat{q}[n] \\ \\ &= A \, e^{j(\omega_0 n + \theta)} \left(1 + \tfrac{1}{A}e^{-j(\omega_0 n + \theta)}\big(q[n] + j \hat{q}[n]\big) \right) \\ \\ &= A \, e^{j(\omega_0 n + \theta)} \big(1 + \epsilon[n] \big) \\ \\ &\approx A \, e^{j(\omega_0 n + \theta)} e^{\epsilon[n]} \\ \\ &= A \, e^{j(\omega_0 n + \theta) + \epsilon[n]} \\ \end{align} $$
$\epsilon[n]$ is complex with the magnitude square having mean
$$\begin{align} \overline{|\epsilon[n]|^2} &= \frac{1}{A^2}\big(\overline{q^2[n]}+\overline{\hat{q}^2[n]}\big) \\ \\ &= \frac{2\overline{q^2[n]}}{A^2} \\ \end{align}$$
And the variance of either the real part or the imaginary part:
$$\begin{align} \overline{\big(\Re e \{\epsilon[n]\}\big)^2} &= \tfrac12 \overline{|\epsilon[n]|^2} \\ \\ \\\overline{\big(\Im m \{\epsilon[n]\}\big)^2} &= \tfrac12 \overline{\big|e[n]\big|^2} \\ \\ &= \frac{\overline{q^2[n]}}{A^2} \\ \end{align}$$
It is the imaginary part, $\Im m \{\epsilon[n]\}$, that is the additive error to the increasing phase of the sinusoid.
Still not quite done with this ...