解决方案
底线
$$(\theta_2-\theta_1) = 2\pi f(T_2-T_1)n -(\phi_2[n]-\phi_1[n]) \tag{1}$$
$f$ : 相同频率和固定相位偏移的两个音调的频率 (Hz)
$(\theta_2-\theta_1)$：被采样音调的弧度相位差
$T_1$ : 采样率$f_{s1}$以秒为单位的采样时钟 1 的周期
$T_2$ : 采样率$f_{s1}$以秒为单位的采样时钟 2 的周期
$\phi_1[n]$ :以弧度/样本为单位的$f_{s1}$采样色调的相位结果
$\phi_2[n]$ :以弧度/样本为单位的$f_{s2}$采样色调的相位结果
这显示了如何扩展以相同采样率采样的相同频率的两个音调之间的相位的任何标准方法（常见的相位检测器方法，包括乘法、相关等）可以扩展以处理两个采样率不同的情况.
先简单解释一下：
考虑等式 (1) 的指数频率形式：
$$e^{j(\theta_2-\theta_1)} = e^{j2\pi f(T_2-T_1)n}e^{-j(\phi_2[n]-\phi_1[n])} \tag {2}$$
术语$e^{j2\pi f(T_2-T_1)n}$是预测的两个音调之间的频率差异，这是由于以两种不同的采样率对单个音调进行采样（当在相同的归一化频率上观察两者时）规模）。
观察到的两种音调之间的频率差异为$e^{j(\phi_2[n]-\phi_1[n])} $。
这两项具有相同的频率，具有固定的相位偏移。该相位偏移是两个连续时间音调之间的实际相位差。通过共轭乘法，我们将两者相减，去除相位斜率和固定相位差结果。
推导
方法是仔细处理具有样本时间轴的所有单元。因此，频域以归一化频率为单位：当时间轴为秒时，周期/样本或弧度/样本对应于周期/秒或弧度/秒。因此，我们的采样率，无论它是以秒为单位的时间，将始终等于$1$周期/样本（或$2\pi$弧度/样本，如果在标准化弧度频率下工作）。
具有相同模拟频率的两个信号一旦在时域中分别以不同的速率被采样，将是各自具有不同归一化频率的两个信号。
这简化了问题，为我们提供了以下结果：
给定我们的原始信号为具有不同相位偏移的相同频率的归一化正弦音：
$$x_1(t) = \cos(2\pi ft + \theta_1) \tag{3}$$
$$x_1(t) = \cos(2\pi ft + \theta_2) \tag{4}$$
一次采样但仍以秒为单位的时间：
$$x_1(nT_1) = \cos(2\pi fn T_1 + \theta_1) \tag{5} $$
$$x_2(nT_2) = \cos(2\pi fn T_2 + \theta_2) \标签{6}$$
式（5）和式（6）以样本为单位表示的时间为：
$$x_1[n] = \cos(2\pi f T_1 n+ \theta_1) \tag{7}$$
$$x_2[n] = \cos(2\pi f T_2 n+ \theta_2) \tag{8} $$
转换为复指数形式，以便我们可以使用复共轭乘法轻松提取相位项，（对于单音，我们只需将输入信号拆分为正交分量；$\cos(\phi) \rightarrow [\cos(\ phi),\sin(\phi)]\rightarrow \cos(\phi)+j\sin(\phi) = e^{j\phi}$，这是使用希尔伯特变换描述为$h\{\} $ )
$$h\{x_1[n]\} =e^{-j(\phi_1[n])} = e^{2\pi f T_1 n+ \theta_1} = e^{2\pi f T_1 n}e ^{\theta_1} \tag{9}$$
$$h\{x_2[n]\} = e^{-j(\phi_2[n])} =e^{2\pi f T_2 n+ \theta_2} =e^{2\pi f T_2 n}e^{\theta_2} \tag{10}$$
复共轭乘法为我们提供了我们寻求的差分相位项及其与我们测量结果的关系：
$$e^{-j(\phi_2[n]-\phi_1[n])} = e^{2\pi f T_2 n}e^{\theta_2}e^{-2\pi f T_1 n}e ^{-\theta_1} \tag{11}$$
导致
$$e^{j(\theta_2-\theta_1)} = e^{j2\pi f(T_2-T_1)n}e^{-j(\phi_2[n]-\phi_1[n])} \tag {12}$$
注意$e^{-j(\phi_2[n]-\phi_1[n])}$表示对于单音将导致频率的测量值，该频率预计为$\omega = 2\pi f (T_2-T_1)n$，由 $e^{j2\pi f(T_2-T_1)n}$项给出。如果我们去除频率偏移（通过上面的乘法），结果就是原始信号的相位差。
取两边的自然对数显示以相位为单位（弧度）的结果：
$$(\theta_2-\theta_1) = 2\pi f(T_2-T_1)n-(\phi_2[n]-\phi_1[n]) \tag{13}$$
所以总而言之，$\phi_1[n]$ , $\phi_2[n]$来自我们给出的测量值$cos(\phi_1[n])$ , $cos(\phi_2[n])$并且我们建立了我们需要通过这些测量的希尔伯特变换的复共轭乘法来得到答案。
示范
我用下面的脚本演示了这一点，类似于 OP 的配置，结果绘制在下面，现在包括抽取，并针对 f = 25 MHz 和 f = 400 MHz（欠采样）进行了测试，结果相似。这显示了演示的每个步骤以上流程，可以进一步组合操作。实施中的希尔伯特变换将是任何将采样音延迟 90° 的选择方法（分数延迟全通滤波器是一个合理的选择）。
Len = 10000;
phase_diff_in = 45;
f=400e6; % Incoming signal frequency
D = 13
Fs1 = 157e6*D;
Fs2 = 121e6*D;
t1 = [0:Len-1]/Fs1; % Time vector channel 1
t2 = [0:Len-1]/Fs2; % Time vector channel 2
phi1 = 2*pi*f*t1;
phi2 = 2*pi*f*t2 + deg2rad(phase_diff_in);
sign1 = cos(phi1);
sign2 = cos(phi2);
% emulation of perfect Hilbert Transform for each tone:
c1_in = 2*(sign1 - 0.5*exp(j*phi1));
c2_in = 2*(sign2 - 0.5*exp(j*phi2));
% create expected phase slope to remove
n = [0:Len-1];
comp_in = exp(-j*2*pi*f*(1/Fs2-1/Fs1)*n);
% decimation
c1 = c1_in(1:D:end);
c2 = c2_in(1:D:end);
comp = comp_in(1:D:end);
pdout = c1.*conj(c2);
result = pdout.*comp;
%determine phase_diff
phase_out = rad2deg(unwrap(angle(result)));
mean_phase = mean(phase_out);
下面是两个测试用例的结果，OP 在他的示例中尝试 0°，然后是 45° 相移。
下面显示了频率为 $f$的输入信号副本的结果，在它们之间零度相位的情况下，由$f_{s1}$作为 sig1 和$f_{s2}$作为 sig2 采样。复共轭积 pdout 的实数是粗体红色正弦曲线，我们注意到它的相位偏移为零。
为了确认计算，下图将其直接与补偿项 $cos(2\pi f(T_2-T_1)) 的实数进行比较，以查看它们的频率是否与等式一致。
并重复$\theta_2-\theta_1 = 45°$
每个样本的原始相位数据的结果表明，每个样本单独具有极低的噪声（受数值精度限制，因此可以用很少的样本确定结果！）。这种性能将取决于希尔伯特变换的实际质量，以准确地将输入音延迟 90° 以创建正交副本。在噪声条件下，可以将结果平均到波形稳定性的范围内，以获得非常稳健的解决方案。
使用欠采样情况进行的扩展性能测试显示了出色的结果（f = 400e6）：
以 1 度的步长测试每个差异角：
10,000 个样本的 RMS 误差（注意垂直轴以 0.5e-11 为增量）
输入频率从 1e6 到 4000e6 以 1e6 为步长进行了极大扩展的频率扫描，具有 45 度相移，在每个频率处测量了 10,000 个点，结果表明，在所有频率（过采样和欠采样）下的相位确定结果都是一致的。这是 OP 配置的两个频率，包括 13 的抽取。（因此，该测试抽取后每个输入音调的采样率是 fs = 157e6 和 121e6，因此该图的最右边是频率被采样的音调为 4e9 明显欠采样。如图所示，RMS 误差与音调的频率成正比，但即使在这种极端条件下，误差仍然小于 5e-10 度。（8.7e-12弧度或-221 dB）。
实际限制
上述结果的准确性受限于对$f_{s1}$和$f_{s2}$给出的准确频率和相位关系的了解，以及对被采样音调的频率$f$的了解。
（如所写，该解决方案还假设两个采样时钟都在时间$t=0$开始，但是可以从等式（8）开始添加采样偏移量，结果相似；底线是两者之间的起始相位关系必须知道或测量采样时钟，因为它会引入额外的偏移）。
现实情况是，没有两个自由运行的时钟会保持完美同步。未锁定到公共参考的采样时钟之间的实际频率和相位差将不可避免地发生漂移（请参阅西格尔定律https://en.wikipedia.org/wiki/Segal%27s_law）。必须将其中一个时钟声明为我们的时间参考（我们的测量将与那个时钟的精度一致）。如果时钟不在物理上位于同一位置，则可以使用双向时间传输技术（参见https://tf.nist.gov/time/twoway.htm）来测量一个时钟与另一个时钟的对比。如果它们在物理上位于同一位置，那么简单的做法就是对一个时钟与另一个时钟进行采样。
下面我展示了这种方法如何从我们的解决方案的方程中完全消除一个采样时钟：（ 我尚未对此进行测试，因此可能包含数学错误）
考虑使用$f_{s2}=\frac{1}{T_2}$对$f_{s1}=\frac{1}{T_1}$进行抽样。这将最终通过使用$f_{s1}$作为公共参考完全从等式中删除$f_{s2}$ （我们基本上通过对$f_{ s1}进行采样来测量$f_{s2}$和$f_{s1}$ $和$f_{s2}$允许我们将$f_{s2}$的样本以$f_{s1}$计数为单位。）：
$f_{s1}$作为余弦：
$$x_{s1}(t) = cos(2\pi f_{s1}t) \tag{14}$$
在给定约束条件下使用$f_{s2}$进行采样时，它们都从 t=0 开始变为：
$$x_{s_1}(nT_2) = cos(2\pi f_{s1}nT_2) = cos(2\pi nT_2/T_1) \tag{15}$$
其中以样本为单位：
$$x_{s_1}[n] =cos(2\pi T_2/T_1 n) \tag{16}$$
导致以样本为单位的第三阶段测量，我们可以通过对$f_{s1}$和$f_{s2}$进行采样来获得（重要的是要同时完成$x_1(t)$和$x_2(t) $被抽样！）：
$$\phi_3[n] = 2\pi T_2/T_1 n \tag{17}$$
因此，如果我们不知道$T_2$但有$\phi_3$，我们可以代入上述等式得到：
$$T_2 = \frac{T_1 \phi_3[n]}{2\pi n} \tag{18}$$
代入(4)：
$$ \phi_2[n]- \phi_1[n] = 2\pi nf\bigg(\frac{T_1 \phi_3[n]}{2\pi n}-T_1\bigg) + (\theta_2-\theta_1) \标签{19} $$
对原始输入信号的期望相位差产生以下解决方案：
$$ \theta_2-\theta_1 = 2\pi f\bigg(\frac{T_1 \phi_3[n]}{2\pi n}-T_1\bigg) n - (\phi_2[n]- \phi_1[n] ) \标签{20} $$
在哪里
$f$ : 被采样的音调频率
$T_1$ : 采样频率为 $f_{s1}$的采样时钟 1 的周期
$\phi_1[n]$ ：使用$f_{s1}$采样音调的结果，值为$cos(\phi_1[n]) $
$\phi_2[n]$ ：使用$f_{s2}$采样音调的结果，值为$cos(\phi_2[n]) $
$\phi_3[n]$ :用$f_{s2}$采样$f_{s1}$的结果，值为$cos(\phi_3[n]) $
因此只知道$T_1$即$1/f_{s1}$，我们可以直接从$x_1(t)$的样本中测量$f$ ，通过采样$ x_1(t)来测量$\phi_1[n]$ $与$f_{s1}$，通过采样$x_2(t)$与$f_{s_2}$测量$\ phi_2 [n]$并通过采样 $f_{s1}$测量$\phi_3[n]$ $f_{s2}$并从这些测量中解析$\theta_2-\theta_1$。
同样，如果您的应用程序的相位偏移不会改变，那么您可以使用结果的斜率测量$f_{s2}$误差，而无需使用$ f_{s2 }$ 对$f_{s1}$进行采样。
真正的结果将取决于$f_{s1}$的实际时钟精度，但我们已经从等式中完全删除了 $f_{s2}$ 。如果您可以将$f_{s1}$视为您对时间的真正参考，这意味着它对于您的测量精度和准确性来说足够准确，那么结果将是您被采样的两个波形的相位差。这意味着最终你需要一些东西来作为你共同的时间参考。
