minyue | 3548357 | 2016-09-20 23:13:08 -0700 | [diff] [blame] | 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 | |
Mirko Bonadei | 92ea95e | 2017-09-15 06:47:31 +0200 | [diff] [blame] | 11 | #include "common_audio/smoothing_filter.h" |
minyue | 3548357 | 2016-09-20 23:13:08 -0700 | [diff] [blame] | 12 | |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 13 | #include <cmath> |
| 14 | |
Mirko Bonadei | 92ea95e | 2017-09-15 06:47:31 +0200 | [diff] [blame] | 15 | #include "rtc_base/timeutils.h" |
michaelt | 92aef17 | 2017-04-18 00:11:48 -0700 | [diff] [blame] | 16 | |
minyue | 3548357 | 2016-09-20 23:13:08 -0700 | [diff] [blame] | 17 | namespace webrtc { |
| 18 | |
michaelt | 92aef17 | 2017-04-18 00:11:48 -0700 | [diff] [blame] | 19 | SmoothingFilterImpl::SmoothingFilterImpl(int init_time_ms) |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 20 | : init_time_ms_(init_time_ms), |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 21 | // Duing the initalization time, we use an increasing alpha. Specifically, |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 22 | // alpha(n) = exp(-powf(init_factor_, n)), |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 23 | // where |init_factor_| is chosen such that |
| 24 | // alpha(init_time_ms_) = exp(-1.0f / init_time_ms_), |
michaelt | 92aef17 | 2017-04-18 00:11:48 -0700 | [diff] [blame] | 25 | init_factor_(init_time_ms_ == 0 |
| 26 | ? 0.0f |
| 27 | : powf(init_time_ms_, -1.0f / init_time_ms_)), |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 28 | // |init_const_| is to a factor to help the calculation during |
| 29 | // initialization phase. |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 30 | init_const_(init_time_ms_ == 0 |
| 31 | ? 0.0f |
| 32 | : init_time_ms_ - |
michaelt | 92aef17 | 2017-04-18 00:11:48 -0700 | [diff] [blame] | 33 | powf(init_time_ms_, 1.0f - 1.0f / init_time_ms_)) { |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 34 | UpdateAlpha(init_time_ms_); |
| 35 | } |
| 36 | |
| 37 | SmoothingFilterImpl::~SmoothingFilterImpl() = default; |
minyue | 3548357 | 2016-09-20 23:13:08 -0700 | [diff] [blame] | 38 | |
| 39 | void SmoothingFilterImpl::AddSample(float sample) { |
michaelt | 92aef17 | 2017-04-18 00:11:48 -0700 | [diff] [blame] | 40 | const int64_t now_ms = rtc::TimeMillis(); |
minyue | 3548357 | 2016-09-20 23:13:08 -0700 | [diff] [blame] | 41 | |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 42 | if (!init_end_time_ms_) { |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 43 | // This is equivalent to assuming the filter has been receiving the same |
| 44 | // value as the first sample since time -infinity. |
| 45 | state_ = last_sample_ = sample; |
Oskar Sundbom | 5e1a749 | 2017-11-16 10:57:19 +0100 | [diff] [blame] | 46 | init_end_time_ms_ = now_ms + init_time_ms_; |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 47 | last_state_time_ms_ = now_ms; |
minyue | 3548357 | 2016-09-20 23:13:08 -0700 | [diff] [blame] | 48 | return; |
| 49 | } |
| 50 | |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 51 | ExtrapolateLastSample(now_ms); |
| 52 | last_sample_ = sample; |
| 53 | } |
| 54 | |
| 55 | rtc::Optional<float> SmoothingFilterImpl::GetAverage() { |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 56 | if (!init_end_time_ms_) { |
| 57 | // |init_end_time_ms_| undefined since we have not received any sample. |
Oskar Sundbom | 5e1a749 | 2017-11-16 10:57:19 +0100 | [diff] [blame] | 58 | return rtc::nullopt; |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 59 | } |
michaelt | 92aef17 | 2017-04-18 00:11:48 -0700 | [diff] [blame] | 60 | ExtrapolateLastSample(rtc::TimeMillis()); |
Oskar Sundbom | 5e1a749 | 2017-11-16 10:57:19 +0100 | [diff] [blame] | 61 | return state_; |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 62 | } |
| 63 | |
| 64 | bool SmoothingFilterImpl::SetTimeConstantMs(int time_constant_ms) { |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 65 | if (!init_end_time_ms_ || last_state_time_ms_ < *init_end_time_ms_) { |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 66 | return false; |
| 67 | } |
| 68 | UpdateAlpha(time_constant_ms); |
| 69 | return true; |
| 70 | } |
| 71 | |
| 72 | void SmoothingFilterImpl::UpdateAlpha(int time_constant_ms) { |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 73 | alpha_ = time_constant_ms == 0 ? 0.0f : exp(-1.0f / time_constant_ms); |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 74 | } |
| 75 | |
| 76 | void SmoothingFilterImpl::ExtrapolateLastSample(int64_t time_ms) { |
| 77 | RTC_DCHECK_GE(time_ms, last_state_time_ms_); |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 78 | RTC_DCHECK(init_end_time_ms_); |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 79 | |
| 80 | float multiplier = 0.0f; |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 81 | |
| 82 | if (time_ms <= *init_end_time_ms_) { |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 83 | // Current update is to be made during initialization phase. |
| 84 | // We update the state as if the |alpha| has been increased according |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 85 | // alpha(n) = exp(-powf(init_factor_, n)), |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 86 | // where n is the time (in millisecond) since the first sample received. |
| 87 | // With algebraic derivation as shown in the Appendix, we can find that the |
| 88 | // state can be updated in a similar manner as if alpha is a constant, |
| 89 | // except for a different multiplier. |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 90 | if (init_time_ms_ == 0) { |
| 91 | // This means |init_factor_| = 0. |
| 92 | multiplier = 0.0f; |
| 93 | } else if (init_time_ms_ == 1) { |
| 94 | // This means |init_factor_| = 1. |
| 95 | multiplier = exp(last_state_time_ms_ - time_ms); |
| 96 | } else { |
| 97 | multiplier = |
| 98 | exp(-(powf(init_factor_, last_state_time_ms_ - *init_end_time_ms_) - |
| 99 | powf(init_factor_, time_ms - *init_end_time_ms_)) / |
| 100 | init_const_); |
| 101 | } |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 102 | } else { |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 103 | if (last_state_time_ms_ < *init_end_time_ms_) { |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 104 | // The latest state update was made during initialization phase. |
| 105 | // We first extrapolate to the initialization time. |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 106 | ExtrapolateLastSample(*init_end_time_ms_); |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 107 | // Then extrapolate the rest by the following. |
minyue | 3548357 | 2016-09-20 23:13:08 -0700 | [diff] [blame] | 108 | } |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 109 | multiplier = powf(alpha_, time_ms - last_state_time_ms_); |
minyue | 3548357 | 2016-09-20 23:13:08 -0700 | [diff] [blame] | 110 | } |
| 111 | |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 112 | state_ = multiplier * state_ + (1.0f - multiplier) * last_sample_; |
| 113 | last_state_time_ms_ = time_ms; |
michaelt | 2fedf9c | 2016-11-28 02:34:18 -0800 | [diff] [blame] | 114 | } |
| 115 | |
minyue | 3548357 | 2016-09-20 23:13:08 -0700 | [diff] [blame] | 116 | } // namespace webrtc |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 117 | |
| 118 | // Appendix: derivation of extrapolation during initialization phase. |
| 119 | // (LaTeX syntax) |
| 120 | // Assuming |
| 121 | // \begin{align} |
| 122 | // y(n) &= \alpha_{n-1} y(n-1) + \left(1 - \alpha_{n-1}\right) x(m) \\* |
| 123 | // &= \left(\prod_{i=m}^{n-1} \alpha_i\right) y(m) + |
| 124 | // \left(1 - \prod_{i=m}^{n-1} \alpha_i \right) x(m) |
| 125 | // \end{align} |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 126 | // Taking $\alpha_{n} = \exp(-\gamma^n)$, $\gamma$ denotes init\_factor\_, the |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 127 | // multiplier becomes |
| 128 | // \begin{align} |
| 129 | // \prod_{i=m}^{n-1} \alpha_i |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 130 | // &= \exp\left(-\sum_{i=m}^{n-1} \gamma^i \right) \\* |
| 131 | // &= \begin{cases} |
| 132 | // \exp\left(-\frac{\gamma^m - \gamma^n}{1 - \gamma} \right) |
| 133 | // & \gamma \neq 1 \\* |
| 134 | // m-n & \gamma = 1 |
| 135 | // \end{cases} |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 136 | // \end{align} |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 137 | // We know $\gamma = T^{-\frac{1}{T}}$, where $T$ denotes init\_time\_ms\_. Then |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 138 | // $1 - \gamma$ approaches zero when $T$ increases. This can cause numerical |
minyue | 7667db4 | 2016-12-28 02:57:50 -0800 | [diff] [blame] | 139 | // difficulties. We multiply $T$ (if $T > 0$) to both numerator and denominator |
| 140 | // in the fraction. See. |
minyue | 301fc4a | 2016-12-13 06:52:56 -0800 | [diff] [blame] | 141 | // \begin{align} |
| 142 | // \frac{\gamma^m - \gamma^n}{1 - \gamma} |
| 143 | // &= \frac{T^\frac{T-m}{T} - T^\frac{T-n}{T}}{T - T^{1-\frac{1}{T}}} |
| 144 | // \end{align} |