webrtc源码分析(6)- jitter delay计算详解

1. 前言

本文主要介绍webrtc jitter buffer中的对于视频帧抖动的计算,关于jitter buffer如何处理乱序组帧的可以参考WebRTC视频JitterBuffer详解,关于处理的抖动后,如何保证视频和音频的同步的可以参考WebRTC音视频同步详解

webrtc版本:M91

2. 正文

2.1 jitter buffer的思想

视频帧从发送端发出后到接收端会经历如下过程:

时刻线: ---------收到包----------开始解码---------完成解码/开始渲染----------完成渲染
耗时:  *网络传输延时*|****延时*****|****解码耗时*********|********渲染耗时*****|

其中网络耗时,解码耗时间和渲染耗时都一定会有,而延时则是为了平滑播放而人为增添的一个耗时;

什么是平滑播放? webrtc的平滑播放是指接收端的帧能够趋向于原始采集源的帧间隔进行播放,由于网络抖动的存在,导致每个帧来的或早或慢,早没问题,但是慢就会导致前面的帧要播放的长而当前帧播放的晚,所以收到一帧后不立刻解码,而是添加一个延时(wait_ms)后再开始解码,这样就能确保下一个帧来的慢也能到;

如下图所示,有三个帧A,B,C,接收到A帧之后不马上进行解码渲染,而是做了一个延时,如此一来,哪怕B帧会因为抖动的关系,导致到达时刻是一个区间,也能保证在A的播放时刻内到达

这个添加的延时的大小非常重要,大了就导致端到端的延迟过高,小了就无法起到有效缓冲的作用,而帧到达的抖动本身是由于传输时间的抖动导致的,所以最合适精确的值是帧传输时间的抖动,这个传输时间的抖动会因为帧的大小以及网络情况的变化而变化,所以要根据这些因素进行调整,简单的思路是对于一段时间内的帧传输时间抖动做加权平均,但这不够精细,也无法考虑到网络路由中发送队列引入的延迟等等;webrtc的思路是这样的,统计一段时间内的最大帧大小(MaxFrameSize)和平均帧大小(AvgFrameSize),计算当前的网络传输速率(C),那就可以通过以下公式得到最保险的抖动大小:

\[jitter = \frac{MaxFrameSize - AvgFrameSize}{C} \]

由于观测过程会存在误差,导致帧间抖动以及网络传输速率不准,引入了卡尔曼滤波去修正以获得一个更准确的值

2.2 卡尔曼滤波

jitter的计算思想如2.1所述是很简单的,卡尔曼滤波更多的是起一个修正的作用,因为我们观测到的网络抖动受噪声影响并不准确,卡尔曼滤波在实际代码中涉及的比较广,所以要看懂jitter的计算;卡尔曼滤波是绕不过的点,卡尔曼滤波的思想是很朴素的,基于贝叶斯概率,贝叶斯概率就是说A和B有关联,知道A发生了就能知道实际上B发生的概率是多少,而表征一个概率的手段就是均值和方差,如果能知道这个系统的均值和方差就知道了这个系统的概率分布,那A发生的时候就能推算出B发生的概率,但这个均值和方差不是上帝告诉我们的,那怎么得到呢?是先随便假设一个,然后通过不断反馈和迭代得到的,现在假设当前的系统的一个概率分布,然后现在A发生了,那么基于当前的概率分布预测到了B发生的概率是50%,但实际上观察到了B发生概率是60%,说明我们当前的概率分布模型有问题(也就是均值和方差有问题),那就让当前的概率模型朝着这个方向去收敛,经过不断迭代之后慢慢的,这个系统的概率分布模型就准确了;要完全理解,请参考这篇笔记--从贝叶斯到卡尔曼和其中的课程视频,只有理解了其思想和本质,到了理解jitter的算法部分,只要理解了其状态方程和观测方程,剩下的就是很简单的套公式了。

对于如下预测方程和观测方程

\[预测方程: X_{k} = F(X_{k-1}) + Q_{K} \hspace{4cm}\\ \]

\[观测方程: Y_k = H(X_{K}) + R_{k} \hspace{4cm}\\ \]

\[Q_{K}为预测误差, Q_{K}服从正态分布N(0 ,Q) \hspace{4cm}\\\ \]

\[R_{k}为观测误差, R_{k}服从正态分布N(0 ,R) \hspace{4cm}\\\ \]

卡尔曼滤波公式如下:

\[先验期望: \vec{u_{k}^{-}} = F \cdot \vec{u_{k-1}^{+}} \hspace{4cm} (2.2.1)\\ \]

\[先验方差: {\Sigma_{k}^{-}}=F {\Sigma_{k-1}^{+}}F^T+Q \hspace{2.5cm} (2.2.2)\\ \]

\[卡尔曼增益: k=\frac{H{\Sigma_{k}^{-}}}{H\Sigma_{k}^{-}H^T + R} \hspace{3cm} (2.2.3)\\ \]

\[后验期望: \vec{u_{k}^{+}} = \vec{u_{k}^{-}} + k*(\vec{y_{k}} - H\vec{u_{k}^{-}}) \hspace{1.1cm} (2.2.4)\\ \]

\[后验方差: {\Sigma_{k}^{+}}= {(I-kH)\Sigma_{k}^{-}} \hspace{3cm} (2.2.5)\\ \]

2.3 jitter滤波算法分析

2.3.1 算法阐述

本节必须在理解透了2.2卡尔曼滤波的思想和本质,才能继续;

首先需要对系统建出两个概率模型-状态方程和观测方程;

帧的传输延时抖动和三个因素有关:帧大小,网络传输速率,网络排队延迟,建模如下:

状态方程:

\[theta(i) = theta(i-1) + u(i-1) \hspace{3cm} (2.3.1) \]

theta(i): 一个二维矩阵\(\left[\begin{array}{c} \frac{1}{C(i)} \\m(i) \end{array} \right]\)\(C(i)\)是网络传输速率,\(m(i)\)是网络排队延迟,网络路由发送队列堆积引入的延迟,这两个变化都很混沌和很多已知未知的因素有关,所以用了直接加噪声这种最简单暴力的建模

u(i-1): 预估误差,高斯噪声矩阵, \(P(u)\sim (0, Q)\), 此处Q是个矩阵

即: 认为theta(i)是与theta(i -1)呈斜率为1的线性关系,这种建模误差由u(i -1) 去弥补


观测方程:

\[d(i) = H * theta(i) + v(i) \hspace{3cm} (2.3.2) \]

d(i): 帧间传输延迟

H: 观测系数变换矩阵\(\left[\begin{array}{c} dL(i) & 1 \end{array} \right]\)\(dL(i)\)是当上一帧和当前帧的大小差

v(i): 观测误差,\(P(v)\sim (0, R)\), 此处\(R\)是一个数

其中,将公式2.3.1代入2.3.2中,得以下方程, 看起来更直观

\[\begin{align*} d(i) &= \left[\begin{array}{} dL(i) & 1 \end{array} \right] * \left[\begin{array}{} \frac{1}{C(i)} \\ m(i)\end{array} \right] + v(i) \\ &= \frac{dL(i)}{C(i)} + m(i) + v(i) \hspace{3cm} (2.3.3) \end{align*} \]

即: 帧间传输延迟的观测值等于 帧变化大小/传输速率 + 网络变化延迟 + 观测噪声


有了两个概率模型之后,参照2.2中的卡尔曼滤波公式就可以得到以下方程了:

先验期望: \(theta\_hat^{-}(i) = theta\_hat(i-1) \hspace{4cm} (2.3.4)\)

先验方差: $theta_cov^{-}(i) = theta_cov(i-1) + Q(i) \hspace{2.5cm} (2.3.5) $

卡尔曼增益: \(k = \frac{H * theta\_cov^{-}(i)}{H * theta\_cov^{-}(i) * H^T + R} \hspace{6.0cm} (2.3.6)\)

后验期望: \(theta\_hat(i) = theta\_hat^{-}(i) + k * (d(i) - H*theta\_hat^{-}(i)) \hspace{2cm} (2.3.7)\)

后验方差: $theta_cov(i) = (I - kH) * theta_cov^{-}(i) \hspace{2.8cm} (2.3.8) $

这里为了贴近google相关文档, 变量起的名字有点四不像, 后缀hat表示期望,cov表示协方差,对于先验的期望和方差,右上角添加了负号

2.3.2 源码分析

2.3.2.1 计算帧间传输抖动观测值d(i)

jitter这一块的计算在webrtc里很集中,在解码线程调用GetNextFrame( )时,会从FramesToDecode队列取出最新要解码的帧,这时就根据该帧的大小和网络传输时间计算抖动d(i)

EncodedFrame* FrameBuffer::GetNextFrame() {
......

  if (!superframe_delayed_by_retransmission) {
    // 正常情况下,非重传
    int64_t frame_delay;
	// 计算当前帧和上一帧在传输时间上的抖动
    if (inter_frame_delay_.CalculateDelay(first_frame->Timestamp(),
                                          &frame_delay, receive_time_ms)) {
      // 通过当前帧大小和传输时间上的抖动更新预估器
      jitter_estimator_.UpdateEstimate(frame_delay, superframe_size);
    }
	
    float rtt_mult = protection_mode_ == kProtectionNackFEC ? 0.0 : 1.0;
    absl::optional<float> rtt_mult_add_cap_ms = absl::nullopt;
    if (rtt_mult_settings_.has_value()) {
      rtt_mult = rtt_mult_settings_->rtt_mult_setting;
      rtt_mult_add_cap_ms = rtt_mult_settings_->rtt_mult_add_cap_ms;
    }
    // 获取当前帧的用于平滑网络延迟,并设置到延迟器中用于后续解码播放
    timing_->SetJitterDelay(
        jitter_estimator_.GetJitterEstimate(rtt_mult, rtt_mult_add_cap_ms));
    //获取之后马上就会做解码,now_ms就是decode_time
    timing_->UpdateCurrentDelay(render_time_ms, now_ms);
  } else {
    if (RttMultExperiment::RttMultEnabled() || add_rtt_to_playout_delay_)
      jitter_estimator_.FrameNacked();
  }

......
}

帧间传输抖动d(i)的计算很简单,每一帧都有rtp时间与接收时间,可以通过计算两帧传输时间,然后相减就可以得到传输抖动,如下图所示:

具体的计算过程如下:

/**
 * @description: 计算帧间传输延迟
 * @param {timestamp} 到达帧的rtp timestamp
 * @param {delay} 计算出的结果
 * @param {CurrentWallClock} 当前时钟
 * @return {} 如果有计算出帧间delay,返回true
 */
bool VCMInterFrameDelay::CalculateDelay(uint32_t timestamp,
                                        int64_t* delay,
                                        int64_t currentWallClock) {
  if (_prevWallClock == 0) {
    // First set of data, initialization, wait for next frame.
    _prevWallClock = currentWallClock;
    _prevTimestamp = timestamp;
    *delay = 0;
    return true;
  }

  int32_t prevWrapArounds = _wrapArounds;
  // 检查是否回环,调整_wrapArounds
  CheckForWrapArounds(timestamp);

  // This will be -1 for backward wrap arounds and +1 for forward wrap arounds.
  // 看是前向回环还是后向回环
  int32_t wrapAroundsSincePrev = _wrapArounds - prevWrapArounds;

  // Account for reordering in jitter variance estimate in the future?
  // Note that this also captures incomplete frames which are grabbed for
  // decoding after a later frame has been complete, i.e. real packet losses.
  // 如果没有发生时间戳回环而且时间戳小于上一个帧,那么当前帧是一个过时的帧,
  if ((wrapAroundsSincePrev == 0 && timestamp < _prevTimestamp) ||
      wrapAroundsSincePrev < 0) {
    *delay = 0;
    return false;
  }

  // Compute the compensated timestamp difference and convert it to ms and round
  // it to closest integer.


  // divide 90 是因为rtp timestamp中的以90khz为单位, ->1ms采样90hz
  // 要想将timestamp转化成ms,就需要除90
  // 根据wrapAroundsSincePrev发生回环数,增大1<<32
  // +0.5 上取整
  _dTS = static_cast<int64_t>(
      (timestamp + wrapAroundsSincePrev * (static_cast<int64_t>(1) << 32) -
       _prevTimestamp) /
          90.0 +
      0.5);
  

  // frameDelay is the difference of dT and dTS -- i.e. the difference of the
  // wall clock time difference and the timestamp difference between two
  // following frames.
  // 这个地方计算两帧的delay是分开来计算的
  //(到达时间差) dt =  currentWallClock - _prevWallClock
  //(发送时间差) dts = currentTimestamp - prevTimeStamp
  *delay = static_cast<int64_t>(currentWallClock - _prevWallClock - _dTS);

  _prevTimestamp = timestamp;
  _prevWallClock = currentWallClock;

  return true;
}

2.3.2.2 计算deviation

inter_frame_delay_.CalculateDelay( )计算除这个抖动值之后,就开始通过函数VCMJitterEstimator::UpdateEstimate( )开始准备卡尔曼滤波。

// Updates the estimates with the new measurements.
void VCMJitterEstimator::UpdateEstimate(int64_t frameDelayMS,
                                        uint32_t frameSizeBytes,
                                        bool incompleteFrame /* = false */) {
  if (frameSizeBytes == 0) {
    return;
  }
  int deltaFS = frameSizeBytes - _prevFrameSize;
  // 启动阶段,先收集几个包,计算一个平均的frameSize
  if (_fsCount < kFsAccuStartupSamples) {
    _fsSum += frameSizeBytes;
    _fsCount++;
  } else if (_fsCount == kFsAccuStartupSamples) {
    // Give the frame size filter.
    _avgFrameSize = static_cast<double>(_fsSum) / static_cast<double>(_fsCount);
    _fsCount++;
  }


  if (!incompleteFrame || frameSizeBytes > _avgFrameSize) {
    // 当frameSize大于_avgFrameSize时,_avgFrameSize平滑趋近于frameSizeBytes
    double avgFrameSize = _phi * _avgFrameSize + (1 - _phi) * frameSizeBytes;
    if (frameSizeBytes < _avgFrameSize + 2 * sqrt(_varFrameSize)) {
      // Only update the average frame size if this sample wasn't a key frame.
      _avgFrameSize = avgFrameSize;
    }
    // Update the variance anyway since we want to capture cases where we only
    // get key frames.
    // 更新framesize方差
    _varFrameSize = VCM_MAX(
        _phi * _varFrameSize + (1 - _phi) * (frameSizeBytes - avgFrameSize) *
                                   (frameSizeBytes - avgFrameSize),
        1.0);
  }

  // Update max frameSize estimate.
  // 在关键帧到达后maxFrameSize会突然增大,这里使用了_psi不断的去收敛_maxFrameSize
  // 到普通P帧大小
  _maxFrameSize =
      VCM_MAX(_psi * _maxFrameSize, static_cast<double>(frameSizeBytes));

  if (_prevFrameSize == 0) {
    _prevFrameSize = frameSizeBytes;
    return;
  }
  _prevFrameSize = frameSizeBytes;

  // Cap frameDelayMS based on the current time deviation noise.
  // A*sqrt(_varNoise)这样的形式在整个文件中出现数次,A实际上是标准正态分布N~(0, 1)下的
  // 分位点,当扩大sqrt(_varNoise)倍时就是 N~(0, _varNoise)下的分位点
  int64_t max_time_deviation_ms =
      static_cast<int64_t>(time_deviation_upper_bound_ * sqrt(_varNoise) + 0.5);
  // 将frameDelayMs 控制在[-max_time_deviation_ms, max_time_deviation_ms]范围内
  frameDelayMS = std::max(std::min(frameDelayMS, max_time_deviation_ms),
                          -max_time_deviation_ms);


  // 利用先验期望计算deviation: d(i) - H * theta_hatˉ(i)
  double deviation = DeviationFromExpectedDelay(frameDelayMS, deltaFS);

  {
    
     RTC_LOG(LS_WARNING)<<"deviation" << deviation
                        <<", _numStdDevDelayOutlier" << _numStdDevDelayOutlier
                        <<", sqrt(_varNoise)"<<sqrt(_varNoise)
                        <<", _numStdDevDelayOutlier * sqrt(_varNoise)" << _numStdDevDelayOutlier * sqrt(_varNoise);
  }


  // Only update the Kalman filter if the sample is not considered an extreme
  // outlier. Even if it is an extreme outlier from a delay point of view, if
  // the frame size also is large the deviation is probably due to an incorrect
  // line slope.
  // 不使用异常点(outlier)更新卡尔曼滤波器
  // 这里异常点的判断依据是:观测噪声值正常会分布于区间[-range, range]
  // range 作为分位点被设置成了: _numStdDevDelayOutlier * sqrt(_varNoise) 
  // _numStdDevDelayOutlier: 标准正态分布中的分位点, 但被设置成了15,这个值怎么会这么大..
  // 但还有一点,当前帧的frame size很大的时候带来的观测噪声很大,有可能是
  // 因为数据量小的时候测得传输速率1/C(即_theta[0])不准导致的,那这就不是一个异常点

  if (fabs(deviation) < _numStdDevDelayOutlier * sqrt(_varNoise) ||
      frameSizeBytes >
          _avgFrameSize + _numStdDevFrameSizeOutlier * sqrt(_varFrameSize)) {
    // Update the variance of the deviation from the line given by the Kalman
    // filter.
    // 利用deviation更新观测噪声的概率分布
    EstimateRandomJitter(deviation, incompleteFrame);
    // Prevent updating with frames which have been congested by a large frame,
    // and therefore arrives almost at the same time as that frame.
    // This can occur when we receive a large frame (key frame) which has been
    // delayed. The next frame is of normal size (delta frame), and thus deltaFS
    // will be << 0. This removes all frame samples which arrives after a key
    // frame.
    if ((!incompleteFrame || deviation >= 0.0) &&
        static_cast<double>(deltaFS) > -0.25 * _maxFrameSize) {
      // Update the Kalman filter with the new data
      KalmanEstimateChannel(frameDelayMS, deltaFS);
    }
  } else {
    // 出现了异常值,
    // 异常值可能是因为某些意外引起的,不希望引起指数滤波器的太大变化
    // 所以使用异常值边界_numStdDevDelayOutlier去更新_varNoise
    int nStdDev =
        (deviation >= 0) ? _numStdDevDelayOutlier : -_numStdDevDelayOutlier;
    EstimateRandomJitter(nStdDev * sqrt(_varNoise), incompleteFrame);
  }

  // Post process the total estimated jitter
  if (_startupCount >= kStartupDelaySamples) {
    PostProcessEstimate();
  } else {
    _startupCount++;
  }
}

首先会统计平均帧大小的_avgFrameSize和最大的帧大小_maxFrameSize。

然后开始计算deviation,deviation是公式2.3.7中的\(d(i) - H*theta\_hat^{-}(i)\)。最开始的时候,我曾把它直接理解成观测误差v(i),因为它很像公式2.3.2的变形,也就是观测值-H*实际值=观测噪声,看起来很合理,但可惜\(theta\_hat^{-}(i)\)不是一个实际值,它是所有实际值的期望,所以这里就把它称为deviation,虽然不是观测误差,但是deviation和观测误差是相近的。

deviation的具体计算过程如下, 有一个隐藏点,由于\(theta\_hat^{-}(i) = theta\_hat(i-1)\),所以在计算过程中直接把\(theta\_hat(i-1)\)拿来用了

/**
 * @description: 利用先验期望计算deviation,计算公式为: v(i) = d(i) - H *theta_hat(i)
 * v(i): 观测值和先验
 * d(i): frameDelayMs 观测值
 * theta_hat(i): theta是一个二维向量[1/C(i)  m(i)]^T, 
 * 1/C(i)是传输速率倒数,m(i-1)是网络时延,theta_hat(i-1) 表示i-1时刻该向量期望,也就是后验值
 * H(i): [dL(i)  1]^T, dL(i)为第i帧和第i-1帧的数据量之差
 * @param {frameDelayMS} 此帧距上帧到达的时间间隔
 * @param {deltaFSBytes} 此帧和上帧大小相差
 * @return {观测噪声值}
 */
double VCMJitterEstimator::DeviationFromExpectedDelay(
    int64_t frameDelayMS,
    int32_t deltaFSBytes) const {
  // d(i) - H(i) *theta_hatˉ(i) =d(i) - _theta[0] * deltaFSBytes + _theta[1] * 1
  // 注意! 此处使用的_theta[0]和 _theta[1] 实际上是theta_hat(i-1)
  // 但由于theta_hatˉ(i) = theta_hat(i-1),所以直接拿来用了
  return frameDelayMS - (_theta[0] * deltaFSBytes + _theta[1]);
}

在计算完deviation之后,会根据这个deviation去判断当前采样点在概率分布上是否为一个超过一个异常边界的异常值,如果不是异常值,会使用它对观测噪声概率分布进行更新,但是这里设置的,如果是,则直接用异常边界去更新观测噪声概率分布。但这个异常边界的值被设置为15,这个分位点在标准正态分布中是非常非常大的,根据标准正态分布表,设置到3.5已经是万分之一的概率了,至今想不明白这个点。不知道是否哪里理解错了,希望知道有缘人告知一下呀。

一般的卡尔曼滤波的观测噪声的概率分布都是固定的,但这里是动态变化的,更新噪声概率分布的具体过程如下

// Estimates the random jitter by calculating the variance of the sample
// distance from the line given by theta.
/**
 * @description: 根据deviation,更新观测噪声均值和噪声方差
 * @param {d_dT} deviation
 * @param {incompleteFrame} 
 * @return {*}
 */
void VCMJitterEstimator::EstimateRandomJitter(double d_dT,
                                              bool incompleteFrame) {
  uint64_t now = clock_->TimeInMicroseconds();
  if (_lastUpdateT != -1) {
    fps_counter_.AddSample(now - _lastUpdateT);
  }
  _lastUpdateT = now;

  if (_alphaCount == 0) {
    assert(false);
    return;
  }
  double alpha =
      static_cast<double>(_alphaCount - 1) / static_cast<double>(_alphaCount);
  _alphaCount++;
  if (_alphaCount > _alphaCountMax)
    _alphaCount = _alphaCountMax;

  // In order to avoid a low frame rate stream to react slower to changes,
  // scale the alpha weight relative a 30 fps stream.
  // alpha会在下面用于指数移动加权滤波器的权重系数, 用于对观测噪声进行更新
  // 此处会通过帧率去修改alpha的值,以30帧为目标,当帧率低于30帧的时候,此时
  // 来帧太慢,导致观测噪声更改的太慢,会通过30.0/帧率对alpha进行pow操作,让alpha
  // 变得更小,从而让观测噪声变化更倾向于当前
  double fps = GetFrameRate();
  if (fps > 0.0) {
    double rate_scale = 30.0 / fps;
    // At startup, there can be a lot of noise in the fps estimate.
    // Interpolate rate_scale linearly, from 1.0 at sample #1, to 30.0 / fps
    // at sample #kStartupDelaySamples.
    if (_alphaCount < kStartupDelaySamples) {
      rate_scale =
          (_alphaCount * rate_scale + (kStartupDelaySamples - _alphaCount)) /
          kStartupDelaySamples;
    }
    // 根据帧率调整alpha 大小使得EWMA中观测噪声调整不因为帧率低而慢
    alpha = pow(alpha, rate_scale);
  }

  // 和传统卡尔曼滤波不同的是,此处观测方程中噪声R不是固定的,
  // 而是会根据deviation去跟新
  // 更新的方式是通过指数移动加权滤波(EWMA)对其进行更新,EWMA如下滤波如下
  // EWMA(i) = λY(i)+ ( 1-λ)EWMA(i-1) 
  // EWMA(i): i时刻的估计值
  // Y(i): i时刻的测量值
  // λ: 权重系数,[0,1],越大说明历史趋势越重要,越小说明即时性更重要
  // 在这里alpha就是λ

  // EWMA对噪声期望更新
  double avgNoise = alpha * _avgNoise + (1 - alpha) * d_dT;
  // EWMA对噪声协方差更新
  double varNoise =
      alpha * _varNoise + (1 - alpha) * (d_dT - _avgNoise) * (d_dT - _avgNoise);

  if (!incompleteFrame || varNoise > _varNoise) {
    _avgNoise = avgNoise;
    _varNoise = varNoise;
  }
  if (_varNoise < 1.0) {
    // The variance should never be zero, since we might get stuck and consider
    // all samples as outliers.
    _varNoise = 1.0;
  }
}

2.3.2.3 卡尔曼滤波

更新完观测噪声的概率分布后,就开始卡尔曼滤波计算后验期望和方差,卡尔曼滤波的计算过程基本就是公式\(2.3.4\)到公式\(2.3.8\), 只要把符号对好,并不复杂

// Updates Kalman estimate of the channel.
// The caller is expected to sanity check the inputs.
void VCMJitterEstimator::KalmanEstimateChannel(int64_t frameDelayMS,
                                               int32_t deltaFSBytes) {
  double Mh[2];         // theta_covˉ * H
  double hMh_sigma;     // H*theta_covˉH' + R
  double kalmanGain[2]; //卡尔曼增益k
  double measureRes;
  double t00, t01;

  // Kalman filtering

  // Prediction,先验协方差预测
  // theta_covˉ(i) = theta_cov(i-1) + Q(i)
  _thetaCov[0][0] += _Qcov[0][0];
  _thetaCov[0][1] += _Qcov[0][1];
  _thetaCov[1][0] += _Qcov[1][0];
  _thetaCov[1][1] += _Qcov[1][1];

  // 卡尔曼增益k的更新,
  // k = H * theta_covˉ(i) / H * theta_covˉ(i) * H' + R
  Mh[0] = _thetaCov[0][0] * deltaFSBytes + _thetaCov[0][1];
  Mh[1] = _thetaCov[1][0] * deltaFSBytes + _thetaCov[1][1];
  // sigma weights measurements with a small deltaFS as noisy and
  // measurements with large deltaFS as good
  if (_maxFrameSize < 1.0) {
    return;
  }
  //sigma 观测噪声方差的计算使用指数滤波的方式从噪声协方差获得
  // deltaFs越大,噪声就越小
  double sigma = (300.0 * exp(-fabs(static_cast<double>(deltaFSBytes)) /
                              (1e0 * _maxFrameSize)) +
                  1) *
                 sqrt(_varNoise);
  if (sigma < 1.0) {
    sigma = 1.0;
  }
  hMh_sigma = deltaFSBytes * Mh[0] + Mh[1] + sigma;
  if ((hMh_sigma < 1e-9 && hMh_sigma >= 0) ||
      (hMh_sigma > -1e-9 && hMh_sigma <= 0)) {
    assert(false);
    return;
  }
  kalmanGain[0] = Mh[0] / hMh_sigma;
  kalmanGain[1] = Mh[1] / hMh_sigma;

  // 后验期望的计算
  // theta_hat(i) = theta_hatˉ(i) + k*(d(i) - H *theta_hatˉ(i))
  measureRes = frameDelayMS - (deltaFSBytes * _theta[0] + _theta[1]);//d(i) - H *theta_hatˉ(i)
  _theta[0] += kalmanGain[0] * measureRes;
  _theta[1] += kalmanGain[1] * measureRes;

  if (_theta[0] < _thetaLow) {
    _theta[0] = _thetaLow;
  }

  // 后验方差的计算
  // theta_cov(i) = (I -kH) * theta_covˉ(i)
  t00 = _thetaCov[0][0];
  t01 = _thetaCov[0][1];
  _thetaCov[0][0] = (1 - kalmanGain[0] * deltaFSBytes) * t00 -
                    kalmanGain[0] * _thetaCov[1][0];
  _thetaCov[0][1] = (1 - kalmanGain[0] * deltaFSBytes) * t01 -
                    kalmanGain[0] * _thetaCov[1][1];
  _thetaCov[1][0] = _thetaCov[1][0] * (1 - kalmanGain[1]) -
                    kalmanGain[1] * deltaFSBytes * t00;
  _thetaCov[1][1] = _thetaCov[1][1] * (1 - kalmanGain[1]) -
                    kalmanGain[1] * deltaFSBytes * t01;

  // Covariance matrix, must be positive semi-definite.
  assert(_thetaCov[0][0] + _thetaCov[1][1] >= 0 &&
         _thetaCov[0][0] * _thetaCov[1][1] -
                 _thetaCov[0][1] * _thetaCov[1][0] >=
             0 &&
         _thetaCov[0][0] >= 0);
}

2.3.2.4 计算抖动

卡尔曼滤波处理完之后,此时就得到了一个相对准确的传输速率C(i)网络排队延迟m(i), 则可以通过这两个值取计算抖动了,在VCMJitterEstimator::UpdateEstimate( )的最后调用了PostProcessEstimate(),计算了jitter

void VCMJitterEstimator::PostProcessEstimate() {
  _filterJitterEstimate = CalculateEstimate();
}

// Calculates the current jitter estimate from the filtered estimates.
/**
 * @description: 计算当前的jitter
 * @param {*}
 * @return {jitter}值
 */
double VCMJitterEstimator::CalculateEstimate() {
  // 这里jitter的计算是一个保守的jitter
  // 使用_maxFrameSize - _avgFrameSize预估下一个帧最大可能增长作为deltaFs
  //计算其在当前预估的传输速率C下可能呆的最大抖动,
  double ret = _theta[0] * (_maxFrameSize - _avgFrameSize) + NoiseThreshold();

  // A very low estimate (or negative) is neglected.
  if (ret < 1.0) {
    if (_prevEstimate <= 0.01) {
      ret = 1.0;
    } else {
      ret = _prevEstimate;
    }
  }
  if (ret > 10000.0) {  // Sanity
    ret = 10000.0;
  }
  _prevEstimate = ret;
  return ret;
}

在估计过程中添加NoiseThreshold( )让我费解,在这里,jitter的计算只是一种估计,基于当前的噪声概率分布添加噪声无可厚非,但是这个噪声阈值的计算过程不好理解,_noiseStdDevs * sqrt(_varNoise)为1%噪声上分位 ,而_noiseStdDevOffset的含义不太清楚,其值为30ms,下面的if分支做了判断,仅做推测其含义为,噪声不大的时候使用1.0作为噪声,大于30ms的时候,使用原值

/**
 * @description: 计算噪声阈值
 * @param {*}
 * @return {*}
 */
double VCMJitterEstimator::NoiseThreshold() const {
  // 而是使用一个基础噪声系数 * 噪声标准差 -  差值
  double noiseThreshold = _noiseStdDevs * sqrt(_varNoise) - _noiseStdDevOffset;
  if (noiseThreshold < 1.0) {
    // 噪声不大的时候使用1.0,大于30ms的时候,使用原值
    noiseThreshold = 1.0;
  }
  return noiseThreshold;
}

至此,预测的抖动计算完毕

2.3.2.5 外部获取jitter

再回到GetNextFrame( )这个函数中,外部会通过VCMJitterEstimator::GetJitterEstimate( )去获取jitter,获取的jitter delay会传递到VCMTiming,用它来计算一个帧解码时间,下面是外部获取jitter delay的具体实现:

// Returns the current filtered estimate if available,
// otherwise tries to calculate an estimate.
int VCMJitterEstimator::GetJitterEstimate(
    double rttMultiplier,
    absl::optional<double> rttMultAddCapMs) {
  // 计算jitter
  double jitterMS = CalculateEstimate() + OPERATING_SYSTEM_JITTER;
  uint64_t now = clock_->TimeInMicroseconds();

  // 此处注意! nack在jitter模块一分钟才过期清零,影响下面jitter的变大
  // 可以优化
  if (now - _latestNackTimestamp > kNackCountTimeoutMs * 1000)
    _nackCount = 0;

  // 暂时不知道什么情况会触发
  if (_filterJitterEstimate > jitterMS)
    jitterMS = _filterJitterEstimate;
  
  if (_nackCount >= _nackLimit) {
    // 丢包太多的时候,将jitterMs增大
    // 不太理解为什么添加额外处理丢包延时,按道理来说,卡尔曼滤波建模的
    // 传输速率的计算是deltaFs/delta 会考虑到平均丢包时间
    if (rttMultAddCapMs.has_value()) {
      jitterMS +=
          std::min(_rttFilter.RttMs() * rttMultiplier, rttMultAddCapMs.value());
    } else {
      jitterMS += _rttFilter.RttMs() * rttMultiplier;
    }
  }

  // 低帧率的时候允许降低jitter
  if (enable_reduced_delay_) {
    static const double kJitterScaleLowThreshold = 5.0;
    static const double kJitterScaleHighThreshold = 10.0;
    double fps = GetFrameRate();
    // Ignore jitter for very low fps streams.
    // 低帧率的时候,网络jitter引起的抖动已经不明显了
    // 这个时候的主要问题是卡顿,可以适当降低jitter了
    if (fps < kJitterScaleLowThreshold) {
      // 只有5帧时,每一帧都超过了200ms,网络抖动没有太大意义
      if (fps == 0.0) {
        return rtc::checked_cast<int>(std::max(0.0, jitterMS) + 0.5);
      }
      return 0;
    }

    // Semi-low frame rate; scale by factor linearly interpolated from 0.0 at
    // kJitterScaleLowThreshold to 1.0 at kJitterScaleHighThreshold.
    if (fps < kJitterScaleHighThreshold) {
      // 处于semid_low区间时,根据帧率调整
      jitterMS =
          (1.0 / (kJitterScaleHighThreshold - kJitterScaleLowThreshold)) *
          (fps - kJitterScaleLowThreshold) * jitterMS;
    }
  }

  return rtc::checked_cast<int>(std::max(0.0, jitterMS) + 0.5);
}

3. REF

1. A Google Congestion Control Algorithm for Real-Time Communication draft-alvestrand-rmcat-congestion-03

2. WebRtc Video Receiver(八)-基于Kalman filter模型的JitterDelay原理分析

3.从贝叶斯到卡尔曼滤波

4.WebRTC视频JitterBuffer详解

posted @ 2021-07-12 10:26  woder  阅读(6781)  评论(6编辑  收藏  举报