1 /*
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
11 #include "common_audio/smoothing_filter.h"
12
13 #include <math.h>
14
15 #include <cmath>
16
17 #include "rtc_base/checks.h"
18 #include "rtc_base/time_utils.h"
19
20 namespace webrtc {
21
SmoothingFilterImpl(int init_time_ms)22 SmoothingFilterImpl::SmoothingFilterImpl(int init_time_ms)
23 : init_time_ms_(init_time_ms),
24 // Duing the initalization time, we use an increasing alpha. Specifically,
25 // alpha(n) = exp(-powf(init_factor_, n)),
26 // where |init_factor_| is chosen such that
27 // alpha(init_time_ms_) = exp(-1.0f / init_time_ms_),
28 init_factor_(init_time_ms_ == 0
29 ? 0.0f
30 : powf(init_time_ms_, -1.0f / init_time_ms_)),
31 // |init_const_| is to a factor to help the calculation during
32 // initialization phase.
33 init_const_(init_time_ms_ == 0
34 ? 0.0f
35 : init_time_ms_ -
36 powf(init_time_ms_, 1.0f - 1.0f / init_time_ms_)) {
37 UpdateAlpha(init_time_ms_);
38 }
39
40 SmoothingFilterImpl::~SmoothingFilterImpl() = default;
41
AddSample(float sample)42 void SmoothingFilterImpl::AddSample(float sample) {
43 const int64_t now_ms = rtc::TimeMillis();
44
45 if (!init_end_time_ms_) {
46 // 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;
49 init_end_time_ms_ = now_ms + init_time_ms_;
50 last_state_time_ms_ = now_ms;
51 return;
52 }
53
54 ExtrapolateLastSample(now_ms);
55 last_sample_ = sample;
56 }
57
GetAverage()58 absl::optional<float> SmoothingFilterImpl::GetAverage() {
59 if (!init_end_time_ms_) {
60 // |init_end_time_ms_| undefined since we have not received any sample.
61 return absl::nullopt;
62 }
63 ExtrapolateLastSample(rtc::TimeMillis());
64 return state_;
65 }
66
SetTimeConstantMs(int time_constant_ms)67 bool SmoothingFilterImpl::SetTimeConstantMs(int time_constant_ms) {
68 if (!init_end_time_ms_ || last_state_time_ms_ < *init_end_time_ms_) {
69 return false;
70 }
71 UpdateAlpha(time_constant_ms);
72 return true;
73 }
74
UpdateAlpha(int time_constant_ms)75 void SmoothingFilterImpl::UpdateAlpha(int time_constant_ms) {
76 alpha_ = time_constant_ms == 0 ? 0.0f : std::exp(-1.0f / time_constant_ms);
77 }
78
ExtrapolateLastSample(int64_t time_ms)79 void SmoothingFilterImpl::ExtrapolateLastSample(int64_t time_ms) {
80 RTC_DCHECK_GE(time_ms, last_state_time_ms_);
81 RTC_DCHECK(init_end_time_ms_);
82
83 float multiplier = 0.0f;
84
85 if (time_ms <= *init_end_time_ms_) {
86 // Current update is to be made during initialization phase.
87 // We update the state as if the |alpha| has been increased according
88 // alpha(n) = exp(-powf(init_factor_, n)),
89 // 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.
93 if (init_time_ms_ == 0) {
94 // This means |init_factor_| = 0.
95 multiplier = 0.0f;
96 } else if (init_time_ms_ == 1) {
97 // This means |init_factor_| = 1.
98 multiplier = std::exp(last_state_time_ms_ - time_ms);
99 } else {
100 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_);
104 }
105 } else {
106 if (last_state_time_ms_ < *init_end_time_ms_) {
107 // The latest state update was made during initialization phase.
108 // We first extrapolate to the initialization time.
109 ExtrapolateLastSample(*init_end_time_ms_);
110 // Then extrapolate the rest by the following.
111 }
112 multiplier = powf(alpha_, time_ms - last_state_time_ms_);
113 }
114
115 state_ = multiplier * state_ + (1.0f - multiplier) * last_sample_;
116 last_state_time_ms_ = time_ms;
117 }
118
119 } // namespace webrtc
120
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}
129 // Taking $\alpha_{n} = \exp(-\gamma^n)$, $\gamma$ denotes init\_factor\_, the
130 // multiplier becomes
131 // \begin{align}
132 // \prod_{i=m}^{n-1} \alpha_i
133 // &= \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}
139 // \end{align}
140 // We know $\gamma = T^{-\frac{1}{T}}$, where $T$ denotes init\_time\_ms\_. Then
141 // $1 - \gamma$ approaches zero when $T$ increases. This can cause numerical
142 // difficulties. We multiply $T$ (if $T > 0$) to both numerator and denominator
143 // in the fraction. See.
144 // \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}
148