1 /*
2 * Copyright (c) 2022 The WebRTC project authors. All Rights Reserved.
3 *
4 * Use of this source code is governed by a BSD-style license
5 * that can be found in the LICENSE file in the root of the source
6 * tree. An additional intellectual property rights grant can be found
7 * in the file PATENTS. All contributing project authors may
8 * be found in the AUTHORS file in the root of the source tree.
9 */
10
11 #include "modules/video_coding/timing/frame_delay_variation_kalman_filter.h"
12
13 #include "api/units/data_size.h"
14 #include "api/units/time_delta.h"
15
16 namespace webrtc {
17
18 namespace {
19 // TODO(brandtr): The value below corresponds to 8 Gbps. Is that reasonable?
20 constexpr double kMaxBandwidth = 0.000001; // Unit: [1 / bytes per ms].
21 } // namespace
22
FrameDelayVariationKalmanFilter()23 FrameDelayVariationKalmanFilter::FrameDelayVariationKalmanFilter() {
24 // TODO(brandtr): Is there a factor 1000 missing here?
25 estimate_[0] = 1 / (512e3 / 8); // Unit: [1 / bytes per ms]
26 estimate_[1] = 0; // Unit: [ms]
27
28 // Initial estimate covariance.
29 estimate_cov_[0][0] = 1e-4; // Unit: [(1 / bytes per ms)^2]
30 estimate_cov_[1][1] = 1e2; // Unit: [ms^2]
31 estimate_cov_[0][1] = estimate_cov_[1][0] = 0;
32
33 // Process noise covariance.
34 process_noise_cov_diag_[0] = 2.5e-10; // Unit: [(1 / bytes per ms)^2]
35 process_noise_cov_diag_[1] = 1e-10; // Unit: [ms^2]
36 }
37
PredictAndUpdate(double frame_delay_variation_ms,double frame_size_variation_bytes,double max_frame_size_bytes,double var_noise)38 void FrameDelayVariationKalmanFilter::PredictAndUpdate(
39 double frame_delay_variation_ms,
40 double frame_size_variation_bytes,
41 double max_frame_size_bytes,
42 double var_noise) {
43 // Sanity checks.
44 if (max_frame_size_bytes < 1) {
45 return;
46 }
47 if (var_noise <= 0.0) {
48 return;
49 }
50
51 // This member function follows the data flow in
52 // https://en.wikipedia.org/wiki/Kalman_filter#Details.
53
54 // 1) Estimate prediction: `x = F*x`.
55 // For this model, there is no need to explicitly predict the estimate, since
56 // the state transition matrix is the identity.
57
58 // 2) Estimate covariance prediction: `P = F*P*F' + Q`.
59 // Again, since the state transition matrix is the identity, this update
60 // is performed by simply adding the process noise covariance.
61 estimate_cov_[0][0] += process_noise_cov_diag_[0];
62 estimate_cov_[1][1] += process_noise_cov_diag_[1];
63
64 // 3) Innovation: `y = z - H*x`.
65 // This is the part of the measurement that cannot be explained by the current
66 // estimate.
67 double innovation =
68 frame_delay_variation_ms -
69 GetFrameDelayVariationEstimateTotal(frame_size_variation_bytes);
70
71 // 4) Innovation variance: `s = H*P*H' + r`.
72 double estim_cov_times_obs[2];
73 estim_cov_times_obs[0] =
74 estimate_cov_[0][0] * frame_size_variation_bytes + estimate_cov_[0][1];
75 estim_cov_times_obs[1] =
76 estimate_cov_[1][0] * frame_size_variation_bytes + estimate_cov_[1][1];
77 double observation_noise_stddev =
78 (300.0 * exp(-fabs(frame_size_variation_bytes) /
79 (1e0 * max_frame_size_bytes)) +
80 1) *
81 sqrt(var_noise);
82 if (observation_noise_stddev < 1.0) {
83 observation_noise_stddev = 1.0;
84 }
85 // TODO(brandtr): Shouldn't we add observation_noise_stddev^2 here? Otherwise,
86 // the dimensional analysis fails.
87 double innovation_var = frame_size_variation_bytes * estim_cov_times_obs[0] +
88 estim_cov_times_obs[1] + observation_noise_stddev;
89 if ((innovation_var < 1e-9 && innovation_var >= 0) ||
90 (innovation_var > -1e-9 && innovation_var <= 0)) {
91 RTC_DCHECK_NOTREACHED();
92 return;
93 }
94
95 // 5) Optimal Kalman gain: `K = P*H'/s`.
96 // How much to trust the model vs. how much to trust the measurement.
97 double kalman_gain[2];
98 kalman_gain[0] = estim_cov_times_obs[0] / innovation_var;
99 kalman_gain[1] = estim_cov_times_obs[1] / innovation_var;
100
101 // 6) Estimate update: `x = x + K*y`.
102 // Optimally weight the new information in the innovation and add it to the
103 // old estimate.
104 estimate_[0] += kalman_gain[0] * innovation;
105 estimate_[1] += kalman_gain[1] * innovation;
106
107 // (This clamping is not part of the linear Kalman filter.)
108 if (estimate_[0] < kMaxBandwidth) {
109 estimate_[0] = kMaxBandwidth;
110 }
111
112 // 7) Estimate covariance update: `P = (I - K*H)*P`
113 double t00 = estimate_cov_[0][0];
114 double t01 = estimate_cov_[0][1];
115 estimate_cov_[0][0] =
116 (1 - kalman_gain[0] * frame_size_variation_bytes) * t00 -
117 kalman_gain[0] * estimate_cov_[1][0];
118 estimate_cov_[0][1] =
119 (1 - kalman_gain[0] * frame_size_variation_bytes) * t01 -
120 kalman_gain[0] * estimate_cov_[1][1];
121 estimate_cov_[1][0] = estimate_cov_[1][0] * (1 - kalman_gain[1]) -
122 kalman_gain[1] * frame_size_variation_bytes * t00;
123 estimate_cov_[1][1] = estimate_cov_[1][1] * (1 - kalman_gain[1]) -
124 kalman_gain[1] * frame_size_variation_bytes * t01;
125
126 // Covariance matrix, must be positive semi-definite.
127 RTC_DCHECK(estimate_cov_[0][0] + estimate_cov_[1][1] >= 0 &&
128 estimate_cov_[0][0] * estimate_cov_[1][1] -
129 estimate_cov_[0][1] * estimate_cov_[1][0] >=
130 0 &&
131 estimate_cov_[0][0] >= 0);
132 }
133
GetFrameDelayVariationEstimateSizeBased(double frame_size_variation_bytes) const134 double FrameDelayVariationKalmanFilter::GetFrameDelayVariationEstimateSizeBased(
135 double frame_size_variation_bytes) const {
136 // Unit: [1 / bytes per millisecond] * [bytes] = [milliseconds].
137 return estimate_[0] * frame_size_variation_bytes;
138 }
139
GetFrameDelayVariationEstimateTotal(double frame_size_variation_bytes) const140 double FrameDelayVariationKalmanFilter::GetFrameDelayVariationEstimateTotal(
141 double frame_size_variation_bytes) const {
142 double frame_transmission_delay_ms =
143 GetFrameDelayVariationEstimateSizeBased(frame_size_variation_bytes);
144 double link_queuing_delay_ms = estimate_[1];
145 return frame_transmission_delay_ms + link_queuing_delay_ms;
146 }
147
148 } // namespace webrtc
149