80 lines
2.3 KiB
C#
80 lines
2.3 KiB
C#
namespace Monitor.Backend.Analysis;
|
|
|
|
/// <summary>
|
|
/// 估计两端的时钟偏移和线性漂移：单向时延 x(t) = 偏移 + 漂移·t + 真实的时延变化。
|
|
/// 最小二乘拟合出直线，残差再减去残差的低分位（时延“地板”），就是去除偏移和漂移之后的抖动。
|
|
/// </summary>
|
|
public static class ClockDriftFit
|
|
{
|
|
public const int MinPoints = 8;
|
|
|
|
/// <summary>拟合 y = a + b·x，x 先居中以保证数值稳定。样本不足或 x 全部相同时返回 false。</summary>
|
|
public static bool TryFit(ReadOnlySpan<double> x, ReadOnlySpan<double> y, out double intercept, out double slope, out double xMean)
|
|
{
|
|
intercept = 0;
|
|
slope = 0;
|
|
xMean = 0;
|
|
int n = x.Length;
|
|
if (n < 2)
|
|
{
|
|
return false;
|
|
}
|
|
|
|
double sumX = 0;
|
|
double sumY = 0;
|
|
for (int i = 0; i < n; i++)
|
|
{
|
|
sumX += x[i];
|
|
sumY += y[i];
|
|
}
|
|
|
|
xMean = sumX / n;
|
|
double yMean = sumY / n;
|
|
double sxx = 0;
|
|
double sxy = 0;
|
|
for (int i = 0; i < n; i++)
|
|
{
|
|
double dx = x[i] - xMean;
|
|
sxx += dx * dx;
|
|
sxy += dx * (y[i] - yMean);
|
|
}
|
|
|
|
if (sxx <= 0)
|
|
{
|
|
return false;
|
|
}
|
|
|
|
slope = sxy / sxx;
|
|
intercept = yMean - (slope * xMean);
|
|
return true;
|
|
}
|
|
|
|
/// <summary>
|
|
/// 去除偏移和漂移后的时延变化（非负）。样本少于 <see cref="MinPoints"/> 时返回 null（不产生观测值）。
|
|
/// </summary>
|
|
public static double[]? Jitter(ReadOnlySpan<double> timesUs, ReadOnlySpan<double> delaysUs)
|
|
{
|
|
int n = timesUs.Length;
|
|
if (n < MinPoints || !TryFit(timesUs, delaysUs, out double intercept, out double slope, out double xMean))
|
|
{
|
|
return null;
|
|
}
|
|
|
|
double[] residual = new double[n];
|
|
for (int i = 0; i < n; i++)
|
|
{
|
|
residual[i] = delaysUs[i] - (intercept + (slope * (timesUs[i] - xMean)));
|
|
}
|
|
|
|
double[] sorted = (double[])residual.Clone();
|
|
Array.Sort(sorted);
|
|
double floor = sorted[(int)(0.05 * (n - 1))];
|
|
for (int i = 0; i < n; i++)
|
|
{
|
|
residual[i] = Math.Max(0, residual[i] - floor);
|
|
}
|
|
|
|
return residual;
|
|
}
|
|
}
|