blob: eaaf3a003377e5efe3ecf222033e65cd53e2a1e8 [file] [log] [blame]
minyue35483572016-09-21 06:13:081/*
2 * Copyright (c) 2016 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
Mirko Bonadei92ea95e2017-09-15 04:47:3111#include "common_audio/smoothing_filter.h"
minyue35483572016-09-21 06:13:0812
Mirko Bonadeib00dec52019-03-25 08:31:0613#include <math.h>
Jonas Olssona4d87372019-07-05 17:08:3314
minyue301fc4a2016-12-13 14:52:5615#include <cmath>
16
Yves Gerey988cc082018-10-23 10:03:0117#include "rtc_base/checks.h"
Steve Anton10542f22019-01-11 17:11:0018#include "rtc_base/time_utils.h"
michaelt92aef172017-04-18 07:11:4819
minyue35483572016-09-21 06:13:0820namespace webrtc {
21
michaelt92aef172017-04-18 07:11:4822SmoothingFilterImpl::SmoothingFilterImpl(int init_time_ms)
minyue7667db42016-12-28 10:57:5023 : init_time_ms_(init_time_ms),
minyue301fc4a2016-12-13 14:52:5624 // Duing the initalization time, we use an increasing alpha. Specifically,
minyue7667db42016-12-28 10:57:5025 // alpha(n) = exp(-powf(init_factor_, n)),
Artem Titov96315752021-07-26 10:15:2926 // where `init_factor_` is chosen such that
minyue301fc4a2016-12-13 14:52:5627 // alpha(init_time_ms_) = exp(-1.0f / init_time_ms_),
michaelt92aef172017-04-18 07:11:4828 init_factor_(init_time_ms_ == 0
29 ? 0.0f
30 : powf(init_time_ms_, -1.0f / init_time_ms_)),
Artem Titov96315752021-07-26 10:15:2931 // `init_const_` is to a factor to help the calculation during
minyue301fc4a2016-12-13 14:52:5632 // initialization phase.
minyue7667db42016-12-28 10:57:5033 init_const_(init_time_ms_ == 0
34 ? 0.0f
35 : init_time_ms_ -
michaelt92aef172017-04-18 07:11:4836 powf(init_time_ms_, 1.0f - 1.0f / init_time_ms_)) {
minyue301fc4a2016-12-13 14:52:5637 UpdateAlpha(init_time_ms_);
38}
39
40SmoothingFilterImpl::~SmoothingFilterImpl() = default;
minyue35483572016-09-21 06:13:0841
42void SmoothingFilterImpl::AddSample(float sample) {
michaelt92aef172017-04-18 07:11:4843 const int64_t now_ms = rtc::TimeMillis();
minyue35483572016-09-21 06:13:0844
minyue7667db42016-12-28 10:57:5045 if (!init_end_time_ms_) {
minyue301fc4a2016-12-13 14:52:5646 // This is equivalent to assuming the filter has been receiving the same
47 // value as the first sample since time -infinity.
48 state_ = last_sample_ = sample;
Oskar Sundbom5e1a7492017-11-16 09:57:1949 init_end_time_ms_ = now_ms + init_time_ms_;
minyue301fc4a2016-12-13 14:52:5650 last_state_time_ms_ = now_ms;
minyue35483572016-09-21 06:13:0851 return;
52 }
53
minyue301fc4a2016-12-13 14:52:5654 ExtrapolateLastSample(now_ms);
55 last_sample_ = sample;
56}
57
Danil Chapovalov196100e2018-06-21 08:17:2458absl::optional<float> SmoothingFilterImpl::GetAverage() {
minyue7667db42016-12-28 10:57:5059 if (!init_end_time_ms_) {
Artem Titov96315752021-07-26 10:15:2960 // `init_end_time_ms_` undefined since we have not received any sample.
Danil Chapovalov196100e2018-06-21 08:17:2461 return absl::nullopt;
minyue7667db42016-12-28 10:57:5062 }
michaelt92aef172017-04-18 07:11:4863 ExtrapolateLastSample(rtc::TimeMillis());
Oskar Sundbom5e1a7492017-11-16 09:57:1964 return state_;
minyue301fc4a2016-12-13 14:52:5665}
66
67bool SmoothingFilterImpl::SetTimeConstantMs(int time_constant_ms) {
minyue7667db42016-12-28 10:57:5068 if (!init_end_time_ms_ || last_state_time_ms_ < *init_end_time_ms_) {
minyue301fc4a2016-12-13 14:52:5669 return false;
70 }
71 UpdateAlpha(time_constant_ms);
72 return true;
73}
74
75void SmoothingFilterImpl::UpdateAlpha(int time_constant_ms) {
Mirko Bonadeib00dec52019-03-25 08:31:0676 alpha_ = time_constant_ms == 0 ? 0.0f : std::exp(-1.0f / time_constant_ms);
minyue301fc4a2016-12-13 14:52:5677}
78
79void SmoothingFilterImpl::ExtrapolateLastSample(int64_t time_ms) {
80 RTC_DCHECK_GE(time_ms, last_state_time_ms_);
minyue7667db42016-12-28 10:57:5081 RTC_DCHECK(init_end_time_ms_);
minyue301fc4a2016-12-13 14:52:5682
83 float multiplier = 0.0f;
minyue7667db42016-12-28 10:57:5084
85 if (time_ms <= *init_end_time_ms_) {
minyue301fc4a2016-12-13 14:52:5686 // Current update is to be made during initialization phase.
Artem Titov96315752021-07-26 10:15:2987 // We update the state as if the `alpha` has been increased according
minyue7667db42016-12-28 10:57:5088 // alpha(n) = exp(-powf(init_factor_, n)),
minyue301fc4a2016-12-13 14:52:5689 // where n is the time (in millisecond) since the first sample received.
90 // With algebraic derivation as shown in the Appendix, we can find that the
91 // state can be updated in a similar manner as if alpha is a constant,
92 // except for a different multiplier.
minyue7667db42016-12-28 10:57:5093 if (init_time_ms_ == 0) {
Artem Titov96315752021-07-26 10:15:2994 // This means `init_factor_` = 0.
minyue7667db42016-12-28 10:57:5095 multiplier = 0.0f;
96 } else if (init_time_ms_ == 1) {
Artem Titov96315752021-07-26 10:15:2997 // This means `init_factor_` = 1.
Mirko Bonadeib00dec52019-03-25 08:31:0698 multiplier = std::exp(last_state_time_ms_ - time_ms);
minyue7667db42016-12-28 10:57:5099 } else {
Mirko Bonadeib00dec52019-03-25 08:31:06100 multiplier = std::exp(
101 -(powf(init_factor_, last_state_time_ms_ - *init_end_time_ms_) -
102 powf(init_factor_, time_ms - *init_end_time_ms_)) /
103 init_const_);
minyue7667db42016-12-28 10:57:50104 }
minyue301fc4a2016-12-13 14:52:56105 } else {
minyue7667db42016-12-28 10:57:50106 if (last_state_time_ms_ < *init_end_time_ms_) {
minyue301fc4a2016-12-13 14:52:56107 // The latest state update was made during initialization phase.
108 // We first extrapolate to the initialization time.
minyue7667db42016-12-28 10:57:50109 ExtrapolateLastSample(*init_end_time_ms_);
minyue301fc4a2016-12-13 14:52:56110 // Then extrapolate the rest by the following.
minyue35483572016-09-21 06:13:08111 }
minyue7667db42016-12-28 10:57:50112 multiplier = powf(alpha_, time_ms - last_state_time_ms_);
minyue35483572016-09-21 06:13:08113 }
114
minyue301fc4a2016-12-13 14:52:56115 state_ = multiplier * state_ + (1.0f - multiplier) * last_sample_;
116 last_state_time_ms_ = time_ms;
michaelt2fedf9c2016-11-28 10:34:18117}
118
minyue35483572016-09-21 06:13:08119} // namespace webrtc
minyue301fc4a2016-12-13 14:52:56120
121// Appendix: derivation of extrapolation during initialization phase.
122// (LaTeX syntax)
123// Assuming
124// \begin{align}
125// y(n) &= \alpha_{n-1} y(n-1) + \left(1 - \alpha_{n-1}\right) x(m) \\*
126// &= \left(\prod_{i=m}^{n-1} \alpha_i\right) y(m) +
127// \left(1 - \prod_{i=m}^{n-1} \alpha_i \right) x(m)
128// \end{align}
minyue7667db42016-12-28 10:57:50129// Taking $\alpha_{n} = \exp(-\gamma^n)$, $\gamma$ denotes init\_factor\_, the
minyue301fc4a2016-12-13 14:52:56130// multiplier becomes
131// \begin{align}
132// \prod_{i=m}^{n-1} \alpha_i
minyue7667db42016-12-28 10:57:50133// &= \exp\left(-\sum_{i=m}^{n-1} \gamma^i \right) \\*
134// &= \begin{cases}
135// \exp\left(-\frac{\gamma^m - \gamma^n}{1 - \gamma} \right)
136// & \gamma \neq 1 \\*
137// m-n & \gamma = 1
138// \end{cases}
minyue301fc4a2016-12-13 14:52:56139// \end{align}
minyue7667db42016-12-28 10:57:50140// We know $\gamma = T^{-\frac{1}{T}}$, where $T$ denotes init\_time\_ms\_. Then
minyue301fc4a2016-12-13 14:52:56141// $1 - \gamma$ approaches zero when $T$ increases. This can cause numerical
minyue7667db42016-12-28 10:57:50142// difficulties. We multiply $T$ (if $T > 0$) to both numerator and denominator
143// in the fraction. See.
minyue301fc4a2016-12-13 14:52:56144// \begin{align}
145// \frac{\gamma^m - \gamma^n}{1 - \gamma}
146// &= \frac{T^\frac{T-m}{T} - T^\frac{T-n}{T}}{T - T^{1-\frac{1}{T}}}
147// \end{align}