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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
|
/*
* Copyright (c) 2016 The WebRTC project authors. All Rights Reserved.
*
* Use of this source code is governed by a BSD-style license
* that can be found in the LICENSE file in the root of the source
* tree. An additional intellectual property rights grant can be found
* in the file PATENTS. All contributing project authors may
* be found in the AUTHORS file in the root of the source tree.
*/
#include "common_audio/smoothing_filter.h"
#include <math.h>
#include <cmath>
#include "rtc_base/checks.h"
#include "rtc_base/time_utils.h"
namespace webrtc {
SmoothingFilterImpl::SmoothingFilterImpl(int init_time_ms)
: init_time_ms_(init_time_ms),
// Duing the initalization time, we use an increasing alpha. Specifically,
// alpha(n) = exp(-powf(init_factor_, n)),
// where |init_factor_| is chosen such that
// alpha(init_time_ms_) = exp(-1.0f / init_time_ms_),
init_factor_(init_time_ms_ == 0
? 0.0f
: powf(init_time_ms_, -1.0f / init_time_ms_)),
// |init_const_| is to a factor to help the calculation during
// initialization phase.
init_const_(init_time_ms_ == 0
? 0.0f
: init_time_ms_ -
powf(init_time_ms_, 1.0f - 1.0f / init_time_ms_)) {
UpdateAlpha(init_time_ms_);
}
SmoothingFilterImpl::~SmoothingFilterImpl() = default;
void SmoothingFilterImpl::AddSample(float sample) {
const int64_t now_ms = rtc::TimeMillis();
if (!init_end_time_ms_) {
// This is equivalent to assuming the filter has been receiving the same
// value as the first sample since time -infinity.
state_ = last_sample_ = sample;
init_end_time_ms_ = now_ms + init_time_ms_;
last_state_time_ms_ = now_ms;
return;
}
ExtrapolateLastSample(now_ms);
last_sample_ = sample;
}
absl::optional<float> SmoothingFilterImpl::GetAverage() {
if (!init_end_time_ms_) {
// |init_end_time_ms_| undefined since we have not received any sample.
return absl::nullopt;
}
ExtrapolateLastSample(rtc::TimeMillis());
return state_;
}
bool SmoothingFilterImpl::SetTimeConstantMs(int time_constant_ms) {
if (!init_end_time_ms_ || last_state_time_ms_ < *init_end_time_ms_) {
return false;
}
UpdateAlpha(time_constant_ms);
return true;
}
void SmoothingFilterImpl::UpdateAlpha(int time_constant_ms) {
alpha_ = time_constant_ms == 0 ? 0.0f : std::exp(-1.0f / time_constant_ms);
}
void SmoothingFilterImpl::ExtrapolateLastSample(int64_t time_ms) {
RTC_DCHECK_GE(time_ms, last_state_time_ms_);
RTC_DCHECK(init_end_time_ms_);
float multiplier = 0.0f;
if (time_ms <= *init_end_time_ms_) {
// Current update is to be made during initialization phase.
// We update the state as if the |alpha| has been increased according
// alpha(n) = exp(-powf(init_factor_, n)),
// where n is the time (in millisecond) since the first sample received.
// With algebraic derivation as shown in the Appendix, we can find that the
// state can be updated in a similar manner as if alpha is a constant,
// except for a different multiplier.
if (init_time_ms_ == 0) {
// This means |init_factor_| = 0.
multiplier = 0.0f;
} else if (init_time_ms_ == 1) {
// This means |init_factor_| = 1.
multiplier = std::exp(last_state_time_ms_ - time_ms);
} else {
multiplier = std::exp(
-(powf(init_factor_, last_state_time_ms_ - *init_end_time_ms_) -
powf(init_factor_, time_ms - *init_end_time_ms_)) /
init_const_);
}
} else {
if (last_state_time_ms_ < *init_end_time_ms_) {
// The latest state update was made during initialization phase.
// We first extrapolate to the initialization time.
ExtrapolateLastSample(*init_end_time_ms_);
// Then extrapolate the rest by the following.
}
multiplier = powf(alpha_, time_ms - last_state_time_ms_);
}
state_ = multiplier * state_ + (1.0f - multiplier) * last_sample_;
last_state_time_ms_ = time_ms;
}
} // namespace webrtc
// Appendix: derivation of extrapolation during initialization phase.
// (LaTeX syntax)
// Assuming
// \begin{align}
// y(n) &= \alpha_{n-1} y(n-1) + \left(1 - \alpha_{n-1}\right) x(m) \\*
// &= \left(\prod_{i=m}^{n-1} \alpha_i\right) y(m) +
// \left(1 - \prod_{i=m}^{n-1} \alpha_i \right) x(m)
// \end{align}
// Taking $\alpha_{n} = \exp(-\gamma^n)$, $\gamma$ denotes init\_factor\_, the
// multiplier becomes
// \begin{align}
// \prod_{i=m}^{n-1} \alpha_i
// &= \exp\left(-\sum_{i=m}^{n-1} \gamma^i \right) \\*
// &= \begin{cases}
// \exp\left(-\frac{\gamma^m - \gamma^n}{1 - \gamma} \right)
// & \gamma \neq 1 \\*
// m-n & \gamma = 1
// \end{cases}
// \end{align}
// We know $\gamma = T^{-\frac{1}{T}}$, where $T$ denotes init\_time\_ms\_. Then
// $1 - \gamma$ approaches zero when $T$ increases. This can cause numerical
// difficulties. We multiply $T$ (if $T > 0$) to both numerator and denominator
// in the fraction. See.
// \begin{align}
// \frac{\gamma^m - \gamma^n}{1 - \gamma}
// &= \frac{T^\frac{T-m}{T} - T^\frac{T-n}{T}}{T - T^{1-\frac{1}{T}}}
// \end{align}
|
