webRTC音频NetEq之音频包插入缓冲抖动BUFF处理过程_音频加抖动算法

    从远端接收的音频帧，经过解头部RTP后，会首先插入到抖动buff，然后统计延迟信息，绘制延迟直方图，根据直方图计算抖动延时的参数，后续dsp的处理根据这个参数以及其他参数，来决策何种策略处理音频数据。这部分根据webrtc源码详细讲解如何插入抖动buff以及统计延迟直方图。
    在webrtc中，NetEQ插入音频到抖动buff的函数为InsertPacketInternal，传入参数为音频数据的头部信息(包括时间戳或序列号等)和净荷数据。在此处会进行一个时间戳的转换，将外部时间戳转换为内部时间戳，外部时间戳即为RTP携带的时间戳，表示RTP报文发送的时钟频率，单位为样本数而非真正的时间单位秒等。在语音中，通常等于PCM语音的采样率，RTP携带Opus编码数据包时，时钟频率为固定的48kHz，但采样率可以有很多值；在视频中，无论是何种视频编码，外部时间戳（时钟频率）都设置为固定的90kHz。内部时间戳为WebRTC使用的时间戳。如接收到的音频格式采样率是16000HZ，内部按照48000Hz处理，则需要将时间戳转换到48000HZ下的时间戳单位处理。函数如下：

uint32_t TimestampScaler::ToInternal(uint32_t external_timestamp,
                                     uint8_t rtp_payload_type) {//external_timestamp为rtp包打包的时间戳，rtp_payload_type包类型
  const DecoderDatabase::DecoderInfo* info =
      decoder_database_.GetDecoderInfo(rtp_payload_type);
  if (!info) {
    // Payload type is unknown. Do not scale.
    return external_timestamp;
  }
  if (!(info->IsComfortNoise() || info->IsDtmf())) {//时间戳转换的前提是非舒适噪声和dtm包
    // Do not change the timestamp scaling settings for DTMF or CNG.
    numerator_ = info->SampleRateHz();
    if (info->GetFormat().clockrate_hz == 0) {
      // If the clockrate is invalid (i.e. with an old-style external codec)
      // we cannot do any timestamp scaling.
      denominator_ = numerator_;
    } else {
      denominator_ = info->GetFormat().clockrate_hz;
    }
  }
  if (numerator_ != denominator_) {//如果内外部处理时钟不同，则进行事件戳转换
    // We have a scale factor != 1.
    if (!first_packet_received_) {
      external_ref_ = external_timestamp;
      internal_ref_ = external_timestamp;
      first_packet_received_ = true;
    }
    const int64_t external_diff = int64_t{external_timestamp} - external_ref_;//先确定外部进入的包事件戳增量，如果采样率是16000hz，
    //一个包长度10ms=160增量，如果本端处理时钟48000hz，则本段时间戳增量=480，=160*(48000/16000)
    RTC_DCHECK_GT(denominator_, 0);
    external_ref_ = external_timestamp;
    internal_ref_ += (external_diff * numerator_) / denominator_;//转换后的事件戳
    return internal_ref_;
  } else {
    // No scaling.
    return external_timestamp;
  }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

    如果收到的是第一个包或包的源SSRC更改了，则需要初始化NetEq，清空packet_buffer_和dtmf_buffer_，更新sync_buffer_的end_timestamp_

    packet_buffer_->Flush();
    dtmf_buffer_->Flush();

    // Update audio buffer timestamp.
    sync_buffer_->IncreaseEndTimestamp(main_timestamp - timestamp_);//更新第一个包的时间戳到sync_buffer_的end_timestamp_

    // Update codecs.
    timestamp_ = main_timestamp;//记录timestamp_为第一个到达音频包的事件戳
1
2
3
4
5
6
7
8

处理RED类型的包以及DTMF类型的包后，将包插入到packet_buffer_之前先保存当前rtp的包类型，然后插入到packet_buffer_，代码流程如下：

int PacketBuffer::InsertPacketList(
    PacketList* packet_list,//接收的音频包，即待插入的音频包
    const DecoderDatabase& decoder_database,
    absl::optional<uint8_t>* current_rtp_payload_type,
    absl::optional<uint8_t>* current_cng_rtp_payload_type,
    StatisticsCalculator* stats) {
  RTC_DCHECK(stats);
  bool flushed = false;
  //确定当前rtp类型
  for (auto& packet : *packet_list) {
    if (decoder_database.IsComfortNoise(packet.payload_type)) {
      if (*current_cng_rtp_payload_type &&
          **current_cng_rtp_payload_type != packet.payload_type) {
        // 新的舒适噪声类型
        *current_rtp_payload_type = absl::nullopt;
        Flush();
        flushed = true;
      }
      *current_cng_rtp_payload_type = packet.payload_type;
    } else if (!decoder_database.IsDtmf(packet.payload_type)) {
      if ((*current_rtp_payload_type &&
           **current_rtp_payload_type != packet.payload_type) ||
          (*current_cng_rtp_payload_type &&
           !EqualSampleRates(packet.payload_type,
                             **current_cng_rtp_payload_type,
                             decoder_database))) {
        *current_cng_rtp_payload_type = absl::nullopt;
        Flush();
        flushed = true;
      }
      *current_rtp_payload_type = packet.payload_type;
    }
    //插入packet_buffer_
    int return_val = InsertPacket(std::move(packet), stats);
    if (return_val == kFlushed) {
      // The buffer flushed, but this is not an error. We can still continue.
      flushed = true;
    } else if (return_val != kOK) {
      // An error occurred. Delete remaining packets in list and return.
      packet_list->clear();
      return return_val;
    }
  }
  packet_list->clear();
  return flushed ? kFlushed : kOK;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

具体插入packet_buffer_的函数如下：

int PacketBuffer::InsertPacket(Packet&& packet, StatisticsCalculator* stats) {
  if (packet.empty()) {
    RTC_LOG(LS_WARNING) << "InsertPacket invalid packet";
    return kInvalidPacket;
  }

  RTC_DCHECK_GE(packet.priority.codec_level, 0);
  RTC_DCHECK_GE(packet.priority.red_level, 0);

  int return_val = kOK;

  packet.waiting_time = tick_timer_->GetNewStopwatch();//获取包入buff的时间

  if (buffer_.size() >= max_number_of_packets_) {//如果插入包的数量超过了buff的最大包数，则丢掉所有包
    // Buffer is full. Flush it.
    Flush();
    stats->FlushedPacketBuffer();
    RTC_LOG(LS_WARNING) << "Packet buffer flushed";
    return_val = kFlushed;
  }

  // 获取一个迭代器，指向缓冲区中应该插入新包的位置，从列表反向查找
//从后面搜索列表，因为最可能的情况是新包应该在列表的末尾
  PacketList::reverse_iterator rit = std::find_if(
      buffer_.rbegin(), buffer_.rend(), NewTimestampIsLarger(packet));//反向查找当前队列，包应该插入到队列中的包时间戳或者序列号大的位置

  //如果找到了位置，新包会插入到迭代器rit的右边，如果与新包时间戳相同，优先级比新入的包高，不插入
  if (rit != buffer_.rend() && packet.timestamp == rit->timestamp) {
    LogPacketDiscarded(packet.priority.codec_level, stats);
    return return_val;
  }
  
  PacketList::iterator it = rit.base();//如果没有找到位置，新包会插入到迭代器it的左边，如果与新包时间戳，优先级比新入的包低，移除掉，插入新收到的包
  if (it != buffer_.end() && packet.timestamp == it->timestamp) {//
    LogPacketDiscarded(it->priority.codec_level, stats);
    it = buffer_.erase(it);
  }
  buffer_.insert(it, std::move(packet));  // 根据时间戳或序列号的顺序将包插入到适当的位置

  return return_val;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

新包插入packe_buff后，如果不是事件戳乱序包则更新延迟信息。由DelayManager累统计包的延时，绘制延迟直方图，更新函数如下，解析标注在代码中。计算音频包的时间长度，即一个音频包的长度是10ms或者20ms或者更长等：
                    时间长度packet_len(ms)= 1000*（与上一包的index时间戳差/与上一包序列号差）/采样率；
IAT直方图
    IAT直方图统计2000ms内延迟统计概率，直方图划分为 2000ms/20ms = 100 个槽，编号index范围0~99，每个槽记录的是延迟为 index * 20ms 的概率，比如相对延迟50ms，则对应第2个槽，在第二个槽上增加概率值，直方图横坐标为index，如下图所示。

更新延迟信息代码如下：

int DelayManager::Update(uint16_t sequence_number,
                         uint32_t timestamp,
                         int sample_rate_hz) {
  if (sample_rate_hz <= 0) {
    return -1;
  }

  if (!first_packet_received_) {
    // 第一个包到达后，开始定时，并获取序列号和时间戳，准备接收下一个包
    packet_iat_stopwatch_ = tick_timer_->GetNewStopwatch();
    last_seq_no_ = sequence_number;
    last_timestamp_ = timestamp;
    first_packet_received_ = true;
    return 0;
  }

  int packet_len_ms;
  if (!IsNewerTimestamp(timestamp, last_timestamp_) ||
      !IsNewerSequenceNumber(sequence_number, last_seq_no_)) {
    // 如果事件戳或序列号乱序，则保存上一次的包长信息.
    packet_len_ms = packet_len_ms_;
  } else {
    // 根据时间戳来计算包长度
    int64_t packet_len_samp =
        static_cast<uint32_t>(timestamp - last_timestamp_) /
        static_cast<uint16_t>(sequence_number - last_seq_no_);
    packet_len_ms =
        rtc::saturated_cast<int>(1000 * packet_len_samp / sample_rate_hz);
  }

  bool reordered = false;
  if (packet_len_ms > 0) {//只有存在包长才会更新统计信息，计算包间到达时间（IAT），然后添加到包间到达时间直方图（IAT直方图）
    // Inter-arrival time (IAT) in integer "packet times" (rounding down). This
    // is the value added to the inter-arrival time histogram.
    int iat_ms = packet_iat_stopwatch_->ElapsedMs();
    int iat_packets = iat_ms / packet_len_ms;
    // Check for discontinuous packet sequence and re-ordering.
    if (IsNewerSequenceNumber(sequence_number, last_seq_no_ + 1)) {
      // Compensate for gap in the sequence numbers. Reduce IAT with the
      // expected extra time due to lost packets.
      int packet_offset =
          static_cast<uint16_t>(sequence_number - last_seq_no_ - 1);
      iat_packets -= packet_offset;
      iat_ms -= packet_offset * packet_len_ms;
    } else if (!IsNewerSequenceNumber(sequence_number, last_seq_no_)) {
      int packet_offset =
          static_cast<uint16_t>(last_seq_no_ + 1 - sequence_number);
      iat_packets += packet_offset;
      iat_ms += packet_offset * packet_len_ms;
      reordered = true;
    }

    int iat_delay = iat_ms - packet_len_ms;
    int relative_delay;
    if (reordered) {
      relative_delay = std::max(iat_delay, 0);
    } else {
      UpdateDelayHistory(iat_delay, timestamp, sample_rate_hz);
      relative_delay = CalculateRelativePacketArrivalDelay();
    }
    statistics_->RelativePacketArrivalDelay(relative_delay);

    switch (histogram_mode_) {
      case RELATIVE_ARRIVAL_DELAY: {
        const int index = relative_delay / kBucketSizeMs;
        if (index < histogram_->NumBuckets()) {
          // Maximum delay to register is 2000 ms.
          histogram_->Add(index);
        }
        break;
      }
      case INTER_ARRIVAL_TIME: {
        // Saturate IAT between 0 and maximum value.
        iat_packets =
            std::max(std::min(iat_packets, histogram_->NumBuckets() - 1), 0);
        histogram_->Add(iat_packets);
        break;
      }
    }
    // Calculate new |target_level_| based on updated statistics.
    target_level_ = CalculateTargetLevel(iat_packets, reordered);

    LimitTargetLevel();
  }  // End if (packet_len_ms > 0).

  if (enable_rtx_handling_ && reordered &&
      num_reordered_packets_ < kMaxReorderedPackets) {
    ++num_reordered_packets_;
    return 0;
  }
  num_reordered_packets_ = 0;
  // Prepare for next packet arrival.
  packet_iat_stopwatch_ = tick_timer_->GetNewStopwatch();
  last_seq_no_ = sequence_number;
  last_timestamp_ = timestamp;
  return 0;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

那么如何划分延迟信息呢？
    如：当一个长度len=20ms的新包到来时，统计与上一个包到达的时间差iat_ms，则延迟时间iat_delay=iat_ms-包长len，理想状态是iat_delay=0，即包长度len=时间差iat_ms。但是实际网络延迟或抖动会导致iat_delay不为0。将iat_delay记录到最大2000ms的历史延迟队列delay_history_中，如果到达的新包与队列第一个包的时间戳之差超过2000ms，则移除最早入队列的iat_delay，代码如下：

//存储包的延时信息到历史延迟队列
void DelayManager::UpdateDelayHistory(int iat_delay_ms,
                                      uint32_t timestamp,
                                      int sample_rate_hz) {
  PacketDelay delay;
  delay.iat_delay_ms = iat_delay_ms;
  delay.timestamp = timestamp;
  delay_history_.push_back(delay);
  while (timestamp - delay_history_.front().timestamp >
         static_cast<uint32_t>(kMaxHistoryMs * sample_rate_hz / 1000)) {
    delay_history_.pop_front();
  }
}
1
2
3
4
5
6
7
8
9
10
11
12
13

计算delay_history_队列中相对延迟relative_delay，因为该队列中记录的是包与包之间的延迟，延迟包括正延迟(iat_delay>0)和负延迟(iat_delay<0)，相对延迟relative_delay将队列中所有延迟时间依次叠加，如果叠加值不小于0，小于0的按0计算。

int DelayManager::CalculateRelativePacketArrivalDelay() const {
  int relative_delay = 0;
  for (const PacketDelay& delay : delay_history_) {
    relative_delay += delay.iat_delay_ms;
    relative_delay = std::max(relative_delay, 0);
  }
  return relative_delay;
}
1
2
3
4
5
6
7
8

直方图概率值采用定点Q15计算，在这里摘录下别人对定点运算的解释：
知识：Q格式DSP处理浮点数据转换成定点运算
许多DSP都是定点DSP，处理定点数据会相当快，但是处理浮点数据就会非常慢。可以利用Q格式进行浮点数据到定点的转化，节约CPU时间。实际应用中，浮点运算大都时候都是既有整数部分，也有小数部分的。所以要选择一个适当的定标格式才能更好的处理运算。

Q格式表示为：Qm.n，表示数据用m比特表示整数部分，n比特表示小数部分，共需要m+n+1位来表示这个数据，多余的一位用作符合位。假设小数点在n位的左边（从右向左数），从而确定小数的精度

例如Q15表示小数部分有15位，一个short型数据，占2个字节，最高位是符号位，后面15位是小数位，就假设小数点在第15位左边，表示的范围是：-1<X<0.9999695 。

浮点数据转化为Q15，将数据乘以2^{15；Q15数据转化为浮点数据，将数据除以2}15。

例如：假设数据存储空间为2个字节，0.333×2^{15=10911=0x2A9F，0.333的所有运算就可以用0x2A9F表示，同理10911×2}(-15)=0.332977294921875，可以看出浮点数据通过Q格式转化后是有误差的。

例：两个小数相乘，0.333*0.414=0.137862

0.333*2^{15=10911=0x2A9F，0.414*2}15=13565=0x34FD

short a = 0x2A9F;

short b = 0x34FD;

short c = a * b >> 15;　　// 两个Q15格式的数据相乘后为Q30格式数据，因此为了得到Q15的数据结果需要右移15位

这样c的结果是0x11A4=0001000110100100，这个数据同样是Q15格式的，它的小数点假设在第15位左边，即为0.001000110100100=0.1378173828125…和实际结果0.137862差距不大。或者0x11A4 / 2^15 = 0.1378173828125
Q格式的运算
　　1> 定点加减法：须转换成相同的Q格式才能加减

2> 定点乘法：不同Q格式的数据相乘，相当于Q值相加，即Q15数据乘以Q10数据后的结果是Q25格式的数据

3> 定点除法：不同Q格式的数据相除，相当于Q值相减

4> 定点左移：左移相当于Q值增加

5> 定点右移：右移相当于Q减少
如何统计到IAT直方图？
    首先确定相对延迟iat_delay的索引index= iat_delay / 20，20为直方图宽度20ms，然后在直方图中槽中找到index，在保证必须所有槽所对应的概率和为1的前提下，在对应index上添加概率。
    首先对原始直方图中的所有概率bucket(Q30)都与遗忘因子(Q15)相乘，并统计所有概率和vector_sum，为什么会乘遗忘因子，个人理解是，所有的概率和为1，如果再添加上去，打破了1的平衡，因此要将所有index对应的概率乘以一个小于1的因子，这样就可以叠加概率而维持概率总和为1的平衡了，叠加后的概率计算：
                bucket = bucket遗忘因子forget_factor_+（1-遗忘因子forget_factor_）
                vector_sum =vector_sum +（1-遗忘因子forget_factor_）
如果概率总和vector_sum !=1，则需要进行补偿使其为1，算法如下：
如果vector_sum>0，则需要减去一个纠正因子correction，
如果vector_sum<0，则需要加上一个纠正因子correction。
纠正因子correction =min((vector_sum -1），index对应概率的十六分之一)
遗忘因子forget_factor_不是固定不变的，而是变化的，每次接收到一个包都需要更新forget_factor_，并逐渐收敛于0.996（初始设定），收敛方程：
         forget_factor_=0.996（1- 2/延迟信息更新次数 ），随着延迟信息更新次数的逐渐增加，forget_factor_逐渐收敛于设定的0.996。

代码分析如下：

void Histogram::Add(int value) {
  RTC_DCHECK(value >= 0);
  RTC_DCHECK(value < static_cast<int>(buckets_.size()));
  int vector_sum = 0; 
  // 对原始直方图中的所有概率bucket(Q15)都与遗忘因子(Q15)相乘，并统计所有概率和vector_sum
  for (int& bucket : buckets_) {
    bucket = (static_cast<int64_t>(bucket) * forget_factor_) >> 15;
    vector_sum += bucket;
  }

  // 在对应index上增加概率值，forget_factor_是Q15格式，而buckets_值是Q30格式，所以在此还必须左翼<<15来转化成Q30格式
  buckets_[value] += (32768 - forget_factor_) << 15;
  vector_sum += (32768 - forget_factor_) << 15;  // 将（1-forget_factor_）叠加到概率总和上，vector_sum也是Q30格式

  // vector_sum应该为1，如果不为1 就需要补偿
  vector_sum -= 1 << 30; 
  if (vector_sum != 0) {
    // 更改前面一段的bucket值
    int flip_sign = vector_sum > 0 ? -1 : 1;
    for (int& bucket : buckets_) {
      int correction = flip_sign * std::min(std::abs(vector_sum), bucket >> 4);
      bucket += correction;
      vector_sum += correction;
      if (std::abs(vector_sum) == 0) {
        break;
      }
    }
  }
  RTC_DCHECK(vector_sum == 0);  // Verify that the above is correct.

  ++add_count_;

 
  if (start_forget_weight_) {
    if (forget_factor_ != base_forget_factor_) { // 更新forget_factor_，随着add_count_越来越大，forget_factor_逐渐接近base_forget_factor_=0.996
      int old_forget_factor = forget_factor_;
      int forget_factor =
          (1 << 15) * (1 - start_forget_weight_.value() / (add_count_ + 1));
      forget_factor_ =
          std::max(0, std::min(base_forget_factor_, forget_factor));
      RTC_DCHECK_GE((1 << 15) - forget_factor_,
                    ((1 << 15) - old_forget_factor) * forget_factor_ >> 15);
    }
  } else {
    forget_factor_ += (base_forget_factor_ - forget_factor_ + 3) >> 2;
  }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

至此，当前包的插入buff并延迟统计直方图统计结束，最后，还需要从直方图中计算出最终的网络延迟target_level，它能反映网络延迟的情况，target_level的计算有理由后面音频DSP处理决策，因此比较重要。
计算思路是：所有槽的值加起来是 1。 接下来遍历直方图，依次累加每个槽的概率， 直到累计值 >= 0.97, 这意味着，当前 index 代表的延时，可以覆盖 97% 的数据包， 我们将找到的 index 记为 bucket_index ，target_level >=1，且初始值为1，Q8格式，则
       target_level = target_level + （bucket_index * kBucketSizeMs） / packet_len_ms_    其中kBucketSizeMs 为直方图宽度为20ms。
计算bucket_index代码如下：

int Histogram::Quantile(int probability) {//probability是查找概率之和统计上限，这里是0.97，Q30格式

  int inverse_probability = (1 << 30) - probability;
  size_t index = 0;        // Start from the beginning of |buckets_|.
  int sum = 1 << 30;       // Assign to 1 in Q30.
  sum -= buckets_[index];
   // 思路：从index  0到99开始累加概率值，知道概率值总和大于等于probability，则返回指定index
  while ((sum > inverse_probability) && (index < buckets_.size() - 1)) {
 
    ++index;
    sum -= buckets_[index];
  }
  return static_cast<int>(index);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14

到此，音频NetEq模块之插入BUFF流程就结束了，每收到一个音频包都会执行这个流程。