Evaluation of phase difference between two acquired signals (expected to be sine waves) with the same frequency
I really dislike the convention of notation in the question.  "$f$" is for frequency.  Functions having capital letters are for Fourier or Laplace transforms of the time-domain signal having the lower-case letter denoting it.

Let's say our reference signal is 

$$ x[n] \triangleq A_x \cos(\omega n + \theta_x) $$

Somehow we can get the Hilbert transform of that signal

$$ \hat{x}[n] \triangleq \mathscr{H} \big \{ x[n] \big \} = A_x \sin(\omega n + \theta_x) $$

If you're generating the reference signal $x[n]$, you should also be able to synchronously generate $\hat{x}[n]$.   If you're measuring $x[n]$, you have to make a hilbert transformer to create $\hat{x}[n]$ (which will be a special FIR filter) and that will introduce a known delay.  You must delay all other signals by the same amount to get this relative phase measurement right.  The general relation of the discrete Hilbert transform is:

$$ \hat{x}[n] = \sum\limits_{m=-\infty}^{\infty} \frac{1 - (-1)^m}{m \, \pi} x[n-m] $$

Of course, you will need to window that infinite sum to a finite sum and delay it to make it causal.  Note that all even $m$ terms in the above summation are zero.

$$\begin{align}
 \hat{x}[n] &amp;= \sum\limits_{m=-\infty}^{\infty} \frac{1 - (-1)^m}{m \, \pi} x[n-m] \\
 &amp;\approx \sum\limits_{m=-2M+1}^{2M-1} w[m] \frac{1 - (-1)^m}{m \, \pi} x[n-m] \\
 &amp;= \sum\limits_{m=1}^{2M-1} w[m] \frac{1 - (-1)^m}{m \, \pi} (x[n-m]-x[n+m]) \\
 &amp;= \sum\limits_{m=0}^{M-1} w[2m+1]\frac{2}{(2m+1) \pi} (x[n-2m-1]-x[n+2m+1]) \\
 &amp;= \sum\limits_{m=0}^{M-1} h[m] \, (x[n-2m-1]-x[n+2m+1]) \\
\end{align}$$

where $w[n]=w[-n]$ is a decent window (like a Kaiser) of width $4M-1$ and the hilbert coefficients are

$$ h[n] \triangleq \frac{2\,w[2n+1]}{(2n+1) \pi} $$

Since this symmetrical, acausal FIR filter looks $2M-1$ samples into the "future", it must be delayed by $2M-1$ samples to be causal and realizable in a real-time context.  So you can not realize in real-time a true Hilbert transformer, but you can realize a very close approximation with the corresponding delay of $2M-1$ samples.  You must delay all other signals by the same amount to make the following calculations valid.  This means you will really be working with $\hat{x}[n-2M+1]$ and you have to explicitly delay $x[n]$ and $y[n]$ by the same amount.

So, while in real-time reality you will be working with $x[n-2M+1]$, $\hat{x}[n-2M+1]$, and $y[n-2M+1]$, for clarity, I am not showing the signals below as so delayed.

Now the measured signal that has the same frequency $\omega$ but different phase is:

$$ y[n] = A_y \cos(\omega n + \theta_y) $$

and the number we are looking for is $\theta_y - \theta_x$ .

Consider these two products:

$$\begin{align}
y[n] \cdot x[n] &amp;= A_y A_x \cos(\omega n + \theta_y) \cos(\omega n + \theta_x) \\
   &amp;= \frac{A_y A_x}{2} \big( \cos(2\omega n + \theta_y+\theta_x) + \cos(\theta_y - \theta_x)  \big) \\
\\
y[n] \cdot \hat{x}[n] &amp;= A_y A_x \cos(\omega n + \theta_y) \sin(\omega n + \theta_x) \\
   &amp;= \frac{A_y A_x}{2} \big( \sin(2\omega n + \theta_y+\theta_x) - \sin(\theta_y - \theta_x)  \big) \\
\end{align}$$

Like you have done before, you want to compute the mean of these two products:

$$\begin{align}
 \overline{y[n] \, x[n]} &amp; = \frac{A_y A_x}{2} \cos(\theta_y - \theta_x) \\
\\
 \overline{y[n] \, \hat{x}[n]} &amp; = -\frac{A_y A_x}{2} \sin(\theta_y - \theta_x) \\
\end{align}$$

We get the mean by use of any decent low-pass filtering.  We don't even need to worry that the DC gain = 1 as you normally do with the mean, since the scaling will come out in the wash.
Then, using complex expression:

$$\begin{align}
 e^{j(\theta_y - \theta_x)} &amp;= \cos(\theta_y - \theta_x) + j \, \sin(\theta_y - \theta_x) \\
  &amp;= \frac{2}{A_y A_x} \left( \overline{y[n] \, x[n]} - j \, \overline{y[n] \, \hat{x}[n]} \right) \\
\end{align}$$

and

$$\begin{align}
\theta_y - \theta_x &amp;= \arg \bigg\{ \overline{y[n] \, x[n]} - j \,\overline{y[n] \, \hat{x}[n]} \bigg\} \\
 &amp;= -\operatorname{atan2}\left(\overline{y[n] \, \hat{x}[n]}, \ \overline{y[n] \, x[n]}\right) \\
\end{align}$$

To do the four-quadrant arctangent:

$$\begin{align}
\operatorname{atan2}(y,x) &amp;= \begin{cases}
\arctan(\frac y x) &amp;\text{if } x &gt; 0, \\
\arctan(\frac y x) + \pi &amp;\text{if } x &lt; 0 \text{ and } y \ge 0, \\
\arctan(\frac y x) - \pi &amp;\text{if } x &lt; 0 \text{ and } y &lt; 0, \\
+\frac{\pi}{2} &amp;\text{if } x = 0 \text{ and } y &gt; 0, \\
-\frac{\pi}{2} &amp;\text{if } x = 0 \text{ and } y &lt; 0, \\
\text{undefined} &amp;\text{if } x = 0 \text{ and } y = 0 \\
\end{cases} \\
\\
&amp;= \begin{cases}
\arctan\left(\frac{y}{x}\right) &amp;\text{if } x &gt; 0, \\
\frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &amp;\text{if } y &gt; 0, \\
-\frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &amp;\text{if } y &lt; 0, \\
\arctan\left(\frac{y}{x}\right) \pm \pi &amp;\text{if } x &lt; 0, \\
\text{undefined} &amp;\text{if } x = 0 \text{ and } y = 0 \\
\end{cases} \\
\end{align}$$