Since I don't really know what you're using software-wise, I just whipped up a really short demo in GNU Radio Companion. This uses a PLL in a live processing demo:
If you're inclined to play around with it, here is the flow graph ready to be run if you've got a running installation of GNU Radio. The flow graph looks like this:
The left part up to the conversion to complex (which actually is a bit superfluous) is just the generation of a noise cosine.
The "PLL Freq Det" Block does pretty much what you describe: it compares the phase of the signal it observes with an internally generated oscillator, and generates a "phase signal" that is used to tell the internal frequency reference "how fast" it should oscillate. That signal is just a phase difference to be compensated, and needs to be internally filtered (the loop filter does that) before being applied to the internal oscillator.
At the same time, that signal is actually also a phase per sample - i.e., a frequency. So, you get a frequency estimate out of this.
Now, this is really practical for things like radio receivers that need to update their estimate of the received signal's frequency while operating, but in your case, you can do all this much more accurate, since you can actually work on the full signal.
A very popular (albeit, seldom optimal) estimator for the spectrum of a signal (and you're looking for which frequencies are making up your signal, so that's a question about the spectrum of your signal) is in fact the (magnitude of the) discrete fourier transform (DFT), often implemented by an FFT.
If you, for example, have your data in a Python program as vector
signal, the following (untested) code should do interesting things to your signal:
from matplotlib import pyplot
## Replace with your own samples and rate
samples = …
sampling_rate = 32e3
## Do a fourier transform of the whole data set
fourier_transform = numpy.fft.fft(samples)
positive_freq_components = fourier_transform[0:len(fourier_transform) / 2]
positive_freq_components_mag_squared = numpy.abs(positive_freq_components)**2
frequencies = numpy.linspace(0, sampling_rate / 2, len(fourier_transform) / 2)
## Plot that stuff!
This should give you a plot of power over frequency. Your noise will probably be relatively white, i.e. it will not have a lot of concentration in frequency domain. The sine, on the other hand, is very periodic and thus has all its energy at one frequency.
Note that I said "seldom optimal" about the DFT as spectrum estimator: It gives you the same frequency resolution all over your bandwidth, and you don't care about most frequencies. The (usable) frequency resolution depends linearly on the observation length, and you can't just get arbitrary frequencies by looking for peaks – all frequencies represented are sampling rate / length of DFT.
Other, "cooler", methods can give you a spectrum estimate that you can evaluate at arbitrary frequencies (such as MUSIC¹), or directly an estimate for the "strongest" frequency (or the strongest N frequencies) (such as ESPRIT¹). However, using these often requires a few strong assumptions, so I won't go deeply into that: If you want to know more about parametric estimators for fundamental frequencies, feel encouraged to research and ask a new question – there's folks on here way more experienced with such estimation tasks than I am.