In RTP to NTP estimator use linear regression instead of ad hoc filter

Make averaging test in NtpEstimator less sensitive.

TESTED=Locally patched into chrome and tested on 1st party software and in video_loopback. All produced parameters looked reasonable.

Bug: webrtc:9698
Change-Id: Idc5e80c657ef190dc95da1e27d1288ff9eddd139
Reviewed-on: https://webrtc-review.googlesource.com/c/110500
Commit-Queue: Ilya Nikolaevskiy <ilnik@webrtc.org>
Reviewed-by: Niels Moller <nisse@webrtc.org>
Reviewed-by: Danil Chapovalov <danilchap@webrtc.org>
Cr-Commit-Position: refs/heads/master@{#25603}

modules/rtp_rtcp/source/remote_ntp_time_estimator_unittest.cc

@@ -111,7 +111,17 @@ TEST_F(RemoteNtpTimeEstimatorTest, Estimate) {
 }
 
 TEST_F(RemoteNtpTimeEstimatorTest, AveragesErrorsOut) {
-  // Remote peer sends first 5 RTCP SR without errors.
+  // Remote peer sends first 10 RTCP SR without errors.
+  AdvanceTimeMilliseconds(1000);
+  SendRtcpSr();
+  AdvanceTimeMilliseconds(1000);
+  SendRtcpSr();
+  AdvanceTimeMilliseconds(1000);
+  SendRtcpSr();
+  AdvanceTimeMilliseconds(1000);
+  SendRtcpSr();
+  AdvanceTimeMilliseconds(1000);
+  SendRtcpSr();
   AdvanceTimeMilliseconds(1000);
   SendRtcpSr();
   AdvanceTimeMilliseconds(1000);
@@ -123,18 +133,17 @@ TEST_F(RemoteNtpTimeEstimatorTest, AveragesErrorsOut) {
   AdvanceTimeMilliseconds(1000);
   SendRtcpSr();
 
-  AdvanceTimeMilliseconds(15);
+  AdvanceTimeMilliseconds(150);
   uint32_t rtp_timestamp = GetRemoteTimestamp();
   int64_t capture_ntp_time_ms = local_clock_.CurrentNtpInMilliseconds();
-
   // Local peer gets enough RTCP SR to calculate the capture time.
   EXPECT_EQ(capture_ntp_time_ms, estimator_->Estimate(rtp_timestamp));
 
   // Remote sends corrupted RTCP SRs
   AdvanceTimeMilliseconds(1000);
-  SendRtcpSrInaccurately(10, 10);
+  SendRtcpSrInaccurately(/*ntp_error_ms=*/2, /*networking_delay_ms=*/-1);
   AdvanceTimeMilliseconds(1000);
-  SendRtcpSrInaccurately(-20, 5);
+  SendRtcpSrInaccurately(/*ntp_error_ms=*/-2, /*networking_delay_ms=*/1);
 
   // New RTP packet to estimate timestamp.
   AdvanceTimeMilliseconds(150);

system_wrappers/BUILD.gn

@@ -34,6 +34,7 @@ rtc_static_library("system_wrappers") {
   deps = [
     ":cpu_features_api",
     "..:webrtc_common",
+    "../api:array_view",
     "../modules:module_api_public",
     "../rtc_base:checks",
     "../rtc_base/synchronization:rw_lock_wrapper",

system_wrappers/include/rtp_to_ntp_estimator.h

@@ -42,23 +42,13 @@ class RtpToNtpEstimator {
 
   // Estimated parameters from RTP and NTP timestamp pairs in |measurements_|.
   struct Parameters {
-    // Implicit conversion from int because MovingMedianFilter returns 0
-    // internally if no samples are present. However, it should never happen as
-    // we don't ask smoothing_filter_ to return anything if there were no
-    // samples.
-    Parameters(const int& value) {  // NOLINT
-      RTC_NOTREACHED();
-    }
     Parameters() : frequency_khz(0.0), offset_ms(0.0) {}
 
+    Parameters(double frequency_khz, double offset_ms)
+        : frequency_khz(frequency_khz), offset_ms(offset_ms) {}
+
     double frequency_khz;
     double offset_ms;
-
-    // Needed to make it work inside MovingMedianFilter
-    bool operator<(const Parameters& other) const;
-    bool operator==(const Parameters& other) const;
-    bool operator<=(const Parameters& other) const;
-    bool operator!=(const Parameters& other) const;
   };
 
   // Updates measurements with RTP/NTP timestamp pair from a RTCP sender report.
@@ -70,7 +60,7 @@ class RtpToNtpEstimator {
 
   // Converts an RTP timestamp to the NTP domain in milliseconds.
   // Returns true on success, false otherwise.
-  bool Estimate(int64_t rtp_timestamp, int64_t* rtp_timestamp_ms) const;
+  bool Estimate(int64_t rtp_timestamp, int64_t* ntp_timestamp_ms) const;
 
   // Returns estimated rtp to ntp linear transform parameters.
   const absl::optional<Parameters> params() const;
@@ -82,8 +72,7 @@ class RtpToNtpEstimator {
 
   int consecutive_invalid_samples_;
   std::list<RtcpMeasurement> measurements_;
-  MovingMedianFilter<Parameters> smoothing_filter_;
+  absl::optional<Parameters> params_;
-  bool params_calculated_;
   mutable TimestampUnwrapper unwrapper_;
 };
 }  // namespace webrtc

system_wrappers/source/rtp_to_ntp_estimator.cc

@@ -11,36 +11,23 @@
 #include "system_wrappers/include/rtp_to_ntp_estimator.h"
 
 #include <stddef.h>
+#include <cmath>
+#include <vector>
 
+#include "api/array_view.h"
 #include "rtc_base/checks.h"
 #include "rtc_base/logging.h"
 
 namespace webrtc {
 namespace {
-// Number of RTCP SR reports to use to map between RTP and NTP.
+// Maximum number of RTCP SR reports to use to map between RTP and NTP.
-const size_t kNumRtcpReportsToUse = 2;
+const size_t kNumRtcpReportsToUse = 20;
-// Number of parameters samples used to smooth.
-const size_t kNumSamplesToSmooth = 20;
 // Don't allow NTP timestamps to jump more than 1 hour. Chosen arbitrary as big
 // enough to not affect normal use-cases. Yet it is smaller than RTP wrap-around
 // half-period (90khz RTP clock wrap-arounds every 13.25 hours). After half of
 // wrap-around period it is impossible to unwrap RTP timestamps correctly.
 const int kMaxAllowedRtcpNtpIntervalMs = 60 * 60 * 1000;
 
-// Calculates the RTP timestamp frequency from two pairs of NTP/RTP timestamps.
-bool CalculateFrequency(int64_t ntp_ms1,
-                        uint32_t rtp_timestamp1,
-                        int64_t ntp_ms2,
-                        uint32_t rtp_timestamp2,
-                        double* frequency_khz) {
-  if (ntp_ms1 <= ntp_ms2)
-    return false;
-
-  *frequency_khz = static_cast<double>(rtp_timestamp1 - rtp_timestamp2) /
-                   static_cast<double>(ntp_ms1 - ntp_ms2);
-  return true;
-}
-
 bool Contains(const std::list<RtpToNtpEstimator::RtcpMeasurement>& measurements,
               const RtpToNtpEstimator::RtcpMeasurement& other) {
   for (const auto& measurement : measurements) {
@@ -49,29 +36,47 @@ bool Contains(const std::list<RtpToNtpEstimator::RtcpMeasurement>& measurements,
   }
   return false;
 }
-}  // namespace
 
-bool RtpToNtpEstimator::Parameters::operator<(const Parameters& other) const {
+// Given x[] and y[] writes out such k and b that line y=k*x+b approximates
-  if (frequency_khz < other.frequency_khz - 1e-6) {
+// given points in the best way (Least Squares Method).
-    return true;
+bool LinearRegression(rtc::ArrayView<const double> x,
-  } else if (frequency_khz > other.frequency_khz + 1e-6) {
+                      rtc::ArrayView<const double> y,
+                      double* k,
+                      double* b) {
+  size_t n = x.size();
+  if (n < 2)
     return false;
-  } else {
+
-    return offset_ms < other.offset_ms - 1e-6;
+  if (y.size() != n)
+    return false;
+
+  double avg_x = 0;
+  double avg_y = 0;
+  for (size_t i = 0; i < n; ++i) {
+    avg_x += x[i];
+    avg_y += y[i];
   }
+  avg_x /= n;
+  avg_y /= n;
+
+  double variance_x = 0;
+  double covariance_xy = 0;
+  for (size_t i = 0; i < n; ++i) {
+    double normalized_x = x[i] - avg_x;
+    double normalized_y = y[i] - avg_y;
+    variance_x += normalized_x * normalized_x;
+    covariance_xy += normalized_x * normalized_y;
+  }
+
+  if (std::fabs(variance_x) < 1e-8)
+    return false;
+
+  *k = static_cast<double>(covariance_xy / variance_x);
+  *b = static_cast<double>(avg_y - (*k) * avg_x);
+  return true;
 }
 
-bool RtpToNtpEstimator::Parameters::operator==(const Parameters& other) const {
+}  // namespace
-  return !(other < *this || *this < other);
-}
-
-bool RtpToNtpEstimator::Parameters::operator!=(const Parameters& other) const {
-  return other < *this || *this < other;
-}
-
-bool RtpToNtpEstimator::Parameters::operator<=(const Parameters& other) const {
-  return !(other < *this);
-}
 
 RtpToNtpEstimator::RtcpMeasurement::RtcpMeasurement(uint32_t ntp_secs,
                                                     uint32_t ntp_frac,
@@ -88,31 +93,29 @@ bool RtpToNtpEstimator::RtcpMeasurement::IsEqual(
 }
 
 // Class for converting an RTP timestamp to the NTP domain.
-RtpToNtpEstimator::RtpToNtpEstimator()
+RtpToNtpEstimator::RtpToNtpEstimator() : consecutive_invalid_samples_(0) {}
-    : consecutive_invalid_samples_(0),
-      smoothing_filter_(kNumSamplesToSmooth),
-      params_calculated_(false) {}
 
 RtpToNtpEstimator::~RtpToNtpEstimator() {}
 
 void RtpToNtpEstimator::UpdateParameters() {
-  if (measurements_.size() != kNumRtcpReportsToUse)
+  if (measurements_.size() < 2)
     return;
 
-  Parameters params;
+  std::vector<double> x;
-  int64_t timestamp_new = measurements_.front().unwrapped_rtp_timestamp;
+  std::vector<double> y;
-  int64_t timestamp_old = measurements_.back().unwrapped_rtp_timestamp;
+  x.reserve(measurements_.size());
+  y.reserve(measurements_.size());
+  for (auto it = measurements_.begin(); it != measurements_.end(); ++it) {
+    x.push_back(it->unwrapped_rtp_timestamp);
+    y.push_back(it->ntp_time.ToMs());
+  }
+  double slope, offset;
 
-  int64_t ntp_ms_new = measurements_.front().ntp_time.ToMs();
+  if (!LinearRegression(x, y, &slope, &offset)) {
-  int64_t ntp_ms_old = measurements_.back().ntp_time.ToMs();
-
-  if (!CalculateFrequency(ntp_ms_new, timestamp_new, ntp_ms_old, timestamp_old,
-                          &params.frequency_khz)) {
     return;
   }
-  params.offset_ms = timestamp_new - params.frequency_khz * ntp_ms_new;
+
-  params_calculated_ = true;
+  params_.emplace(1 / slope, offset);
-  smoothing_filter_.Insert(params);
 }
 
 bool RtpToNtpEstimator::UpdateMeasurements(uint32_t ntp_secs,
@@ -159,8 +162,7 @@ bool RtpToNtpEstimator::UpdateMeasurements(uint32_t ntp_secs,
     RTC_LOG(LS_WARNING) << "Multiple consecutively invalid RTCP SR reports, "
                            "clearing measurements.";
     measurements_.clear();
-    smoothing_filter_.Reset();
+    params_ = absl::nullopt;
-    params_calculated_ = false;
   }
   consecutive_invalid_samples_ = 0;
 
@@ -177,35 +179,29 @@ bool RtpToNtpEstimator::UpdateMeasurements(uint32_t ntp_secs,
 }
 
 bool RtpToNtpEstimator::Estimate(int64_t rtp_timestamp,
-                                 int64_t* rtp_timestamp_ms) const {
+                                 int64_t* ntp_timestamp_ms) const {
-  if (!params_calculated_)
+  if (!params_)
     return false;
 
   int64_t rtp_timestamp_unwrapped = unwrapper_.Unwrap(rtp_timestamp);
 
-  Parameters params = smoothing_filter_.GetFilteredValue();
-
   // params_calculated_ should not be true unless ms params.frequency_khz has
   // been calculated to something non zero.
-  RTC_DCHECK_NE(params.frequency_khz, 0.0);
+  RTC_DCHECK_NE(params_->frequency_khz, 0.0);
   double rtp_ms =
-      (static_cast<double>(rtp_timestamp_unwrapped) - params.offset_ms) /
+      static_cast<double>(rtp_timestamp_unwrapped) / params_->frequency_khz +
-          params.frequency_khz +
+      params_->offset_ms + 0.5f;
-      0.5f;
 
   if (rtp_ms < 0)
     return false;
 
-  *rtp_timestamp_ms = rtp_ms;
+  *ntp_timestamp_ms = rtp_ms;
+
   return true;
 }
 
 const absl::optional<RtpToNtpEstimator::Parameters> RtpToNtpEstimator::params()
     const {
-  absl::optional<Parameters> res;
+  return params_;
-  if (params_calculated_) {
-    res.emplace(smoothing_filter_.GetFilteredValue());
-  }
-  return res;
 }
 }  // namespace webrtc

system_wrappers/source/rtp_to_ntp_estimator_unittest.cc

@@ -9,11 +9,13 @@
  */
 
 #include "system_wrappers/include/rtp_to_ntp_estimator.h"
+#include "rtc_base/random.h"
 #include "test/gtest.h"
 
 namespace webrtc {
 namespace {
 const uint32_t kOneMsInNtpFrac = 4294967;
+const uint32_t kOneHourInNtpSec = 60 * 60;
 const uint32_t kTimestampTicksPerMs = 90;
 }  // namespace
 
@@ -224,6 +226,22 @@ TEST(UpdateRtcpMeasurementTests, FailsForOldNtp) {
       estimator.UpdateMeasurements(ntp_sec, ntp_frac, timestamp, &new_sr));
 }
 
+TEST(UpdateRtcpMeasurementTests, FailsForTooNewNtp) {
+  RtpToNtpEstimator estimator;
+  uint32_t ntp_sec = 1;
+  uint32_t ntp_frac = 699925050;
+  uint32_t timestamp = 0x12345678;
+  bool new_sr;
+  EXPECT_TRUE(
+      estimator.UpdateMeasurements(ntp_sec, ntp_frac, timestamp, &new_sr));
+  EXPECT_TRUE(new_sr);
+  // Ntp time from far future, list not updated.
+  ntp_sec += kOneHourInNtpSec * 2;
+  timestamp += kTimestampTicksPerMs * 10;
+  EXPECT_FALSE(
+      estimator.UpdateMeasurements(ntp_sec, ntp_frac, timestamp, &new_sr));
+}
+
 TEST(UpdateRtcpMeasurementTests, FailsForEqualTimestamp) {
   RtpToNtpEstimator estimator;
   uint32_t ntp_sec = 0;
@@ -292,4 +310,37 @@ TEST(RtpToNtpTests, FailsForNoParameters) {
   EXPECT_FALSE(estimator.Estimate(timestamp, &timestamp_ms));
 }
 
+TEST(RtpToNtpTests, AveragesErrorOut) {
+  RtpToNtpEstimator estimator;
+  uint32_t ntp_sec = 1;
+  uint32_t ntp_frac = 90000000;  // More than 1 ms.
+  uint32_t timestamp = 0x12345678;
+  const int kNtpSecStep = 1;  // 1 second.
+  const int kRtpTicksPerMs = 90;
+  const int kRtpStep = kRtpTicksPerMs * 1000;
+  bool new_sr;
+  EXPECT_TRUE(
+      estimator.UpdateMeasurements(ntp_sec, ntp_frac, timestamp, &new_sr));
+  EXPECT_TRUE(new_sr);
+
+  Random rand(1123536L);
+  for (size_t i = 0; i < 1000; i++) {
+    // Advance both timestamps by exactly 1 second.
+    ntp_sec += kNtpSecStep;
+    timestamp += kRtpStep;
+    // Add upto 1ms of errors to NTP and RTP timestamps passed to estimator.
+    EXPECT_TRUE(estimator.UpdateMeasurements(
+        ntp_sec,
+        ntp_frac + rand.Rand(-static_cast<int>(kOneMsInNtpFrac),
+                             static_cast<int>(kOneMsInNtpFrac)),
+        timestamp + rand.Rand(-kRtpTicksPerMs, kRtpTicksPerMs), &new_sr));
+    EXPECT_TRUE(new_sr);
+
+    int64_t estimated_ntp_ms;
+    EXPECT_TRUE(estimator.Estimate(timestamp, &estimated_ntp_ms));
+    // Allow upto 2 ms of error.
+    EXPECT_NEAR(NtpTime(ntp_sec, ntp_frac).ToMs(), estimated_ntp_ms, 2);
+  }
+}
+
 };  // namespace webrtc