The multibit reference signal is typically not sinusoidal, can have offset and substantial noise, but is expected to contain one dominant frequency.

Against offset, you'd practically always start your processing with a high-pass filter.

If you know your dominant frequency to be sufficiently less than Nyquist, a fixed low-pass filter would also be a cheap way to increase your SNR before even doing anything complicated.

What you're looking for is the answer to the question

What's $f_d$ of the dominant frequency in the signal?

and that's a parametric estimation problem. So, I'd recommend using a parametric estimator.

I'm going to go ahead and recommend something that will seem quite like overkill (due to the scary complex algorithms mentioned therein), but it's really quite straightforward in the sizes of problem you deal with.

Do ROOT-MUSIC with a signal subspace size of $m=2$ and a noise subspace size of $n-m=1$ or $2$.

ESPRIT roughly works like the following:

1. You estimate a autocovariance matrix (and that you do by taking a vector $X$ of $n$ samples, and multiplying it with its transpose conjugate, so that you get a $n\times n$ matrix $\hat R_{XX}$. You typically accumulate and average a few of these to get an averaged $\bar R_{XX}$. Not that $n$ is "really small".
2. you do an Eigenvalue Decomposition. Because $\bar R_{XX}$ is real and symmetrical by construction, that becomes really trivial.
3. You take the two largest Eigenvalues (of your 3 to 5 eigenvalues...), and the matching eigenvector. These column vector $S$ spans your signal space! Yay! (Pick two, because your sine is real, and that means it's composed of two inherently linked complex sinusoids of $\pm f_d$).
4. You, however, consider the other eigenvectors. These span the noise space matrix $G$ (which has dimension 1 or 2 , depending on your choice of $n-m$), and by how eigenvalue decompositions of hermitian matrices like $\bar R_{XX}$ happen, that space is orthogonal to the signal space. If you calculated the "pure" sinusoid of the right frequencies $f_d$, their projection into the noise space would be zero. You exploit that. Let $z=e^{j2\pi ft}$, and define $a=\left[ z^{-0}, z^{-1}\right]$. Then, finding the zeros ("roots") of $$a^* GG^* a$$ (which is a very boring 1st or second degree polynomial) gives you your frequencies.

The multibit reference signal is typically not sinusoidal, can have offset and substantial noise, but is expected to contain one dominant frequency.

Against offset, you'd practically always start your processing with a high-pass filter.

If you know your dominant frequency to be sufficiently less than Nyquist, a fixed low-pass filter would also be a cheap way to increase your SNR before even doing anything complicated.

What you're looking for is the answer to the question

What's $f_d$ of the dominant frequency in the signal?

and that's a parametric estimation problem. So, I'd recommend using a parametric estimator.

I'm going to go ahead and recommend something that will seem quite like overkill (due to the scary complex algorithms mentioned therein), but it's really quite straightforward in the sizes of problem you deal with.

Do ROOT-MUSIC with a signal subspace size of $m=2$ and a noise subspace size of $n-m=1$ or $2$.

ROOT-MUSIC roughly works like the following:

1. You estimate a autocovariance matrix (and that you do by taking a vector $X$ of $n$ samples, and multiplying it with its transpose conjugate, so that you get a $n\times n$ matrix $\hat R_{XX}$. You typically accumulate and average a few of these to get an averaged $\bar R_{XX}$. Not that $n$ is "really small".
2. you do an Eigenvalue Decomposition. Because $\bar R_{XX}$ is real and symmetrical by construction, that becomes really trivial.
3. You take the two largest Eigenvalues (of your 3 to 5 eigenvalues...), and the matching eigenvector. These column vector $S$ spans your signal space! Yay! (Pick two, because your sine is real, and that means it's composed of two inherently linked complex sinusoids of $\pm f_d$).
4. You, however, consider the other eigenvectors. These span the noise space matrix $G$ (which has dimension 1 or 2 , depending on your choice of $n-m$), and by how eigenvalue decompositions of hermitian matrices like $\bar R_{XX}$ happen, that space is orthogonal to the signal space. If you calculated the "pure" sinusoid of the right frequencies $f_d$, their projection into the noise space would be zero. You exploit that. Let $z=e^{j2\pi ft}$, and define $a=\left[ z^{-0}, z^{-1}\right]$. Then, finding the zeros ("roots") of $$a^* GG^* a$$ (which is a very boring 1st or second degree polynomial) gives you your frequencies.