Phase difference measurement of a signal sampled with two different sampling frequencies
This is an untested time domain solution, but the math looks solid.

This will be impossible to implement unless you solve the reciever synchronization problem first.  That is either a hardware fix or a calibration operation.  Assume it is solved and your two signals are coming in as time aligned sequences.  Assume also your sampling rates (I don't like "sampling frequency" as a term) are high compared to the frequency of the incoming signal.  Also the incoming signal is assumed to be a pure real tone of a known frequency ($\omega$).

You will be searching for peak values on each sequence.  Use a parabolic fit with the center point being the max.  You will then need to find the signal values from the sequences at that time instance by interpolation.

Another problem you will have is normalizing the magnitude of one signal against the other.  This solution does that implicitly.

With that in mind, we can do the math with a continuous definition.

$$ x_1(t) = A_1 \cos( \omega t + \phi ) $$

$$
\begin{aligned}
x_2(t) &amp;= A_2 \cos( \omega ( t - d ) + \phi )  \\
&amp;= A_2 \left[ \cos( \omega t + \phi ) \cos( \omega d ) + \sin( \omega t + \phi ) \sin( \omega d )  \right] \\
\end{aligned}
$$

Divide the second signal by the first.

$$ \frac{x_2(t)}{x_1(t)} = \frac{A_2}{A_1} \left[ \cos( \omega d ) + \tan( \omega t + \phi ) \sin( \omega d ) \right] $$

At a peak of signal 1, $ \tan( \omega t + \phi ) = 0 $, since the $\sin$ does which is the derivative.

$$ (\frac{x_2}{x_1})_{peak1} = \frac{A_2}{A_1} \cos( \omega d ) $$

If your amplitudes are normalized, you have your answer from this.

By symmetry the reverse argument can be made.

$$ s = t - d $$

$$ x_2(s) = A_2 \cos( \omega s + \phi ) $$

$$
\begin{aligned}
x_1(s) &amp;= A_1 \cos( \omega ( s + d ) + \phi ) \\
&amp;= A_1 \left[ \cos( \omega s + \phi ) \cos( \omega d ) - \sin( \omega s + \phi ) \sin( \omega d )  \right] \\
\end{aligned}
$$

$$ \frac{x_1(s)}{x_2(s)} = \frac{A_1}{A_2} \left[ \cos( \omega d ) - \tan( \omega s + \phi ) \sin( \omega d ) \right] $$

$ \tan( \omega s + \phi ) $ will be zero at a peak of signal 2.

$$ (\frac{x_1}{x_2})_{peak2} = \frac{A_1}{A_2} \cos( \omega d ) $$

We can now combine the ratios from the two different peak locations:

$$ \cos( \omega d ) = \sqrt{ (\frac{x_2}{x_1})_{peak1} \cdot (\frac{x_1}{x_2})_{peak2} } = V $$

Since we are squaring, the sign of $V$ is lost.  It can be deduced from the individual parts.  

$$ d = \frac{ \cos^{-1} (V) }{ \omega } $$

Note, the solution is not unique as the shift can be greater than one cycle length.

This should take a lot fewer calculations than an FFT approach.





Had a "duh moment."

If you can measure the peaks directly, then $d$ is simply their distance apart plus possibly a number of whole cycles.

You can calibrate your system if you can move your source straight ahead for a baseline shift measurement, ideally zero.  Then a measurement the same distance on either side will give you an idea of the angle of the first cycle overlap if there is one.



For a third solution, there is a simple frequency domain approach that will work well in this situation if your signal is too noisy for the time domain one.  Select a duration which is a whole number of cycles, say $k$, of the signal you are receiving.  It will be roughly $M$ number of samples at your one rate and $N$ on the second.  Calculate the $k$th DFT bin for each using two sets of basis vectors (two different DFTs, technically, but one bin each).  You should be able to find an interval so the $M$ and $N$ are pretty good fits.  Read the phase difference from the difference of the two bin angles.  Translate that to a time shift using the frequency.  This should still take a significantly fewer number of calculations than doing FFTs (DFTs) on each signal.