• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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