1 | /* |
2 | Copyright 2018 Google Inc. All Rights Reserved. |
3 | |
4 | Licensed under the Apache License, Version 2.0 (the "License"); |
5 | you may not use this file except in compliance with the License. |
6 | You may obtain a copy of the License at |
7 | |
8 | http://www.apache.org/licenses/LICENSE-2.0 |
9 | |
10 | Unless required by applicable law or agreed to in writing, software |
11 | distributed under the License is distributed on an "AS-IS" BASIS, |
12 | WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
13 | See the License for the specific language governing permissions and |
14 | limitations under the License. |
15 | */ |
16 | |
17 | // Prevent Visual Studio from complaining about std::copy_n. |
18 | #if defined(_WIN32) |
19 | #define _SCL_SECURE_NO_WARNINGS |
20 | #endif |
21 | |
22 | #include "dsp/resampler.h" |
23 | |
24 | #include <functional> |
25 | #include <numeric> |
26 | |
27 | #include "base/constants_and_types.h" |
28 | #include "base/logging.h" |
29 | #include "base/misc_math.h" |
30 | #include "base/simd_utils.h" |
31 | |
32 | #include "dsp/utils.h" |
33 | |
34 | namespace vraudio { |
35 | |
36 | namespace { |
37 | |
38 | // (kTransitionBandwidthRatio / 2) * (sample_rate / cutoff_frequency) |
39 | // = filter_length. |
40 | // The value below was chosen empirically as a tradeoff between execution time |
41 | // and filter rolloff wrt. cutoff frequency. |
42 | const size_t kTransitionBandwidthRatio = 13; |
43 | |
44 | // Maximum number of channels internally based upon the maximum supported |
45 | // ambisonic order. |
46 | const size_t kMaxNumChannels = |
47 | (kMaxSupportedAmbisonicOrder + 1) * (kMaxSupportedAmbisonicOrder + 1); |
48 | |
49 | } // namespace |
50 | |
51 | Resampler::Resampler() |
52 | : up_rate_(0), |
53 | down_rate_(0), |
54 | time_modulo_up_rate_(0), |
55 | last_processed_sample_(0), |
56 | num_channels_(0), |
57 | coeffs_per_phase_(0), |
58 | transposed_filter_coeffs_(kNumMonoChannels, kMaxSupportedNumFrames), |
59 | temporary_filter_coeffs_(kNumMonoChannels, kMaxSupportedNumFrames), |
60 | state_(kMaxNumChannels, kMaxSupportedNumFrames) { |
61 | state_.Clear(); |
62 | } |
63 | |
64 | void Resampler::Process(const AudioBuffer& input, AudioBuffer* output) { |
65 | |
66 | // See "Digital Signal Processing, 4th Edition, Prolakis and Manolakis, |
67 | // Pearson, Chapter 11 (specificaly Figures 11.5.10 and 11.5.13). |
68 | DCHECK_EQ(input.num_channels(), num_channels_); |
69 | const size_t input_length = input.num_frames(); |
70 | DCHECK_GE(output->num_frames(), GetNextOutputLength(input_length)); |
71 | DCHECK_LE(output->num_frames(), GetMaxOutputLength(input_length)); |
72 | DCHECK_EQ(output->num_channels(), num_channels_); |
73 | output->Clear(); |
74 | |
75 | if (up_rate_ == down_rate_) { |
76 | *output = input; |
77 | return; |
78 | } |
79 | |
80 | // |input_sample| is the last processed sample. |
81 | size_t input_sample = last_processed_sample_; |
82 | // |output_sample| is the output. |
83 | size_t output_sample = 0; |
84 | |
85 | const auto& filter_coefficients = transposed_filter_coeffs_[0]; |
86 | |
87 | while (input_sample < input_length) { |
88 | size_t filter_index = time_modulo_up_rate_ * coeffs_per_phase_; |
89 | size_t offset_input_index = input_sample - coeffs_per_phase_ + 1; |
90 | const int offset = -static_cast<int>(offset_input_index); |
91 | |
92 | if (offset > 0) { |
93 | // We will need to draw data from the |state_| buffer. |
94 | const int state_num_frames = static_cast<int>(coeffs_per_phase_ - 1); |
95 | int state_index = state_num_frames - offset; |
96 | while (state_index < state_num_frames) { |
97 | for (size_t channel = 0; channel < num_channels_; ++channel) { |
98 | (*output)[channel][output_sample] += |
99 | state_[channel][state_index] * filter_coefficients[filter_index]; |
100 | } |
101 | state_index++; |
102 | filter_index++; |
103 | } |
104 | // Move along by |offset| samples up as far as |input|. |
105 | offset_input_index += offset; |
106 | } |
107 | |
108 | // We now move back to where |input_sample| "points". |
109 | while (offset_input_index <= input_sample) { |
110 | for (size_t channel = 0; channel < num_channels_; ++channel) { |
111 | (*output)[channel][output_sample] += |
112 | input[channel][offset_input_index] * |
113 | filter_coefficients[filter_index]; |
114 | } |
115 | offset_input_index++; |
116 | filter_index++; |
117 | } |
118 | output_sample++; |
119 | |
120 | time_modulo_up_rate_ += down_rate_; |
121 | // Advance the input pointer. |
122 | input_sample += time_modulo_up_rate_ / up_rate_; |
123 | // Decide which phase of the polyphase filter to use next. |
124 | time_modulo_up_rate_ %= up_rate_; |
125 | } |
126 | DCHECK_GE(input_sample, input_length); |
127 | last_processed_sample_ = input_sample - input_length; |
128 | |
129 | // Take care of the |state_| buffer. |
130 | const int samples_left_in_input = |
131 | static_cast<int>(coeffs_per_phase_) - 1 - static_cast<int>(input_length); |
132 | if (samples_left_in_input > 0) { |
133 | for (size_t channel = 0; channel < num_channels_; ++channel) { |
134 | // Copy end of the |state_| buffer to the beginning. |
135 | auto& state_channel = state_[channel]; |
136 | DCHECK_GE(static_cast<int>(state_channel.size()), samples_left_in_input); |
137 | std::copy_n(first: state_channel.end() - samples_left_in_input, |
138 | n: samples_left_in_input, result: state_channel.begin()); |
139 | // Then copy input to the end of the buffer. |
140 | std::copy_n(first: input[channel].begin(), n: input_length, |
141 | result: state_channel.end() - input_length); |
142 | } |
143 | } else { |
144 | for (size_t channel = 0; channel < num_channels_; ++channel) { |
145 | // Copy the last of the |input| samples into the |state_| buffer. |
146 | DCHECK_GT(coeffs_per_phase_, 0U); |
147 | DCHECK_GT(input[channel].size(), coeffs_per_phase_ - 1); |
148 | std::copy_n(first: input[channel].end() - (coeffs_per_phase_ - 1), |
149 | n: coeffs_per_phase_ - 1, result: state_[channel].begin()); |
150 | } |
151 | } |
152 | } |
153 | |
154 | size_t Resampler::GetMaxOutputLength(size_t input_length) const { |
155 | if (up_rate_ == down_rate_) { |
156 | return input_length; |
157 | } |
158 | DCHECK_GT(down_rate_, 0U); |
159 | // The + 1 takes care of the case where: |
160 | // (time_modulo_up_rate_ + up_rate_ * last_processed_sample_) < |
161 | // ((input_length * up_rate_) % down_rate_) |
162 | // The output length will be equal to the return value or the return value -1. |
163 | return (input_length * up_rate_) / down_rate_ + 1; |
164 | } |
165 | |
166 | size_t Resampler::GetNextOutputLength(size_t input_length) const { |
167 | if (up_rate_ == down_rate_) { |
168 | return input_length; |
169 | } |
170 | const size_t max_length = GetMaxOutputLength(input_length); |
171 | if ((time_modulo_up_rate_ + up_rate_ * last_processed_sample_) >= |
172 | ((input_length * up_rate_) % down_rate_)) { |
173 | return max_length - 1; |
174 | } |
175 | return max_length; |
176 | } |
177 | |
178 | void Resampler::SetRateAndNumChannels(int source_frequency, |
179 | int destination_frequency, |
180 | size_t num_channels) { |
181 | |
182 | // Convert sampling rates to be relatively prime. |
183 | DCHECK_GT(source_frequency, 0); |
184 | DCHECK_GT(destination_frequency, 0); |
185 | DCHECK_GT(num_channels, 0U); |
186 | const int greatest_common_divisor = |
187 | FindGcd(a: destination_frequency, b: source_frequency); |
188 | const size_t destination = |
189 | static_cast<size_t>(destination_frequency / greatest_common_divisor); |
190 | const size_t source = |
191 | static_cast<size_t>(source_frequency / greatest_common_divisor); |
192 | |
193 | // Obtain the size of the |state_| before |coeffs_per_phase_| is updated in |
194 | // |GenerateInterpolatingFilter()|. |
195 | const size_t old_state_size = |
196 | coeffs_per_phase_ > 0 ? coeffs_per_phase_ - 1 : 0; |
197 | if ((destination != up_rate_) || (source != down_rate_)) { |
198 | up_rate_ = destination; |
199 | down_rate_ = source; |
200 | if (up_rate_ == down_rate_) { |
201 | return; |
202 | } |
203 | // Create transposed multirate filters from sincs. |
204 | GenerateInterpolatingFilter(sample_rate: source_frequency); |
205 | // Reset the time variable as it may be longer than the new filter length if |
206 | // we switched from upsampling to downsampling via a call to SetRate(). |
207 | time_modulo_up_rate_ = 0; |
208 | } |
209 | |
210 | // Update the |state_| buffer. |
211 | if (num_channels_ != num_channels) { |
212 | num_channels_ = num_channels; |
213 | InitializeStateBuffer(old_state_num_frames: old_state_size); |
214 | } |
215 | } |
216 | |
217 | bool Resampler::AreSampleRatesSupported(int source, int destination) { |
218 | DCHECK_GT(source, 0); |
219 | DCHECK_GT(destination, 0); |
220 | // Determines whether sample rates are supported based upon whether our |
221 | // maximul filter lenhgth is big enough to hold the corresponding |
222 | // interpolation filter. |
223 | const int max_rate = |
224 | std::max(a: source, b: destination) / FindGcd(a: source, b: destination); |
225 | size_t filter_length = max_rate * kTransitionBandwidthRatio; |
226 | filter_length += filter_length % 2; |
227 | return filter_length <= kMaxSupportedNumFrames; |
228 | } |
229 | |
230 | void Resampler::ResetState() { |
231 | |
232 | time_modulo_up_rate_ = 0; |
233 | last_processed_sample_ = 0; |
234 | state_.Clear(); |
235 | } |
236 | |
237 | void Resampler::InitializeStateBuffer(size_t old_state_num_frames) { |
238 | // Update the |state_| buffer if it is null or if the number of coefficients |
239 | // per phase in the polyphase filter has changed. |
240 | if (up_rate_ == down_rate_ || num_channels_ == 0) { |
241 | return; |
242 | } |
243 | // If the |state_| buffer is to be kept. For example in the case of a change |
244 | // in either source or destination sampling rate, maintaining the old |state_| |
245 | // buffers contents allows a glitch free transition. |
246 | const size_t new_state_num_frames = |
247 | coeffs_per_phase_ > 0 ? coeffs_per_phase_ - 1 : 0; |
248 | if (old_state_num_frames != new_state_num_frames) { |
249 | const size_t min_size = |
250 | std::min(a: new_state_num_frames, b: old_state_num_frames); |
251 | const size_t max_size = |
252 | std::max(a: new_state_num_frames, b: old_state_num_frames); |
253 | for (size_t channel = 0; channel < num_channels_; ++channel) { |
254 | auto& state_channel = state_[channel]; |
255 | DCHECK_LT(state_channel.begin() + max_size, state_channel.end()); |
256 | std::fill(first: state_channel.begin() + min_size, |
257 | last: state_channel.begin() + max_size, value: 0.0f); |
258 | } |
259 | } |
260 | } |
261 | |
262 | void Resampler::GenerateInterpolatingFilter(int sample_rate) { |
263 | // See "Digital Signal Processing, 4th Edition, Prolakis and Manolakis, |
264 | // Pearson, Chapter 11 (specificaly Figures 11.5.10 and 11.5.13). |
265 | const size_t max_rate = std::max(a: up_rate_, b: down_rate_); |
266 | const float cutoff_frequency = |
267 | static_cast<float>(sample_rate) / static_cast<float>(2 * max_rate); |
268 | size_t filter_length = max_rate * kTransitionBandwidthRatio; |
269 | filter_length += filter_length % 2; |
270 | auto* filter_channel = &temporary_filter_coeffs_[0]; |
271 | filter_channel->Clear(); |
272 | GenerateSincFilter(cutoff_frequency, sample_rate: static_cast<float>(sample_rate), |
273 | filter_length, buffer: filter_channel); |
274 | |
275 | // Pad out the filter length so that it can be arranged in polyphase fashion. |
276 | const size_t transposed_length = |
277 | filter_length + max_rate - (filter_length % max_rate); |
278 | coeffs_per_phase_ = transposed_length / max_rate; |
279 | ArrangeFilterAsPolyphase(filter_length, filter: *filter_channel); |
280 | } |
281 | |
282 | void Resampler::ArrangeFilterAsPolyphase(size_t filter_length, |
283 | const AudioBuffer::Channel& filter) { |
284 | // Coefficients are transposed and flipped. |
285 | // Suppose |up_rate_| is 3, and the input number of coefficients is 10, |
286 | // h[0], ..., h[9]. |
287 | // Then the |transposed_filter_coeffs_| buffer will look like this: |
288 | // h[9], h[6], h[3], h[0], flipped phase 0 coefs. |
289 | // 0, h[7], h[4], h[1], flipped phase 1 coefs (zero-padded). |
290 | // 0, h[8], h[5], h[2], flipped phase 2 coefs (zero-padded). |
291 | transposed_filter_coeffs_.Clear(); |
292 | auto& transposed_coefficients_channel = transposed_filter_coeffs_[0]; |
293 | for (size_t i = 0; i < up_rate_; ++i) { |
294 | for (size_t j = 0; j < coeffs_per_phase_; ++j) { |
295 | if (j * up_rate_ + i < filter_length) { |
296 | const size_t coeff_index = |
297 | (coeffs_per_phase_ - 1 - j) + i * coeffs_per_phase_; |
298 | transposed_coefficients_channel[coeff_index] = filter[j * up_rate_ + i]; |
299 | } |
300 | } |
301 | } |
302 | } |
303 | |
304 | void Resampler::GenerateSincFilter(float cutoff_frequency, float sample_rate, |
305 | size_t filter_length, |
306 | AudioBuffer::Channel* buffer) { |
307 | |
308 | DCHECK_GT(sample_rate, 0.0f); |
309 | const float angular_cutoff_frequency = |
310 | kTwoPi * cutoff_frequency / sample_rate; |
311 | |
312 | const size_t half_filter_length = filter_length / 2; |
313 | GenerateHannWindow(full_window: true /* Full Window */, window_length: filter_length, buffer); |
314 | auto* buffer_channel = &buffer[0]; |
315 | |
316 | for (size_t i = 0; i < filter_length; ++i) { |
317 | if (i == half_filter_length) { |
318 | (*buffer_channel)[half_filter_length] *= angular_cutoff_frequency; |
319 | } else { |
320 | const float denominator = |
321 | static_cast<float>(i) - (static_cast<float>(filter_length) / 2.0f); |
322 | DCHECK_GT(std::abs(denominator), kEpsilonFloat); |
323 | (*buffer_channel)[i] *= |
324 | std::sin(x: angular_cutoff_frequency * denominator) / denominator; |
325 | } |
326 | } |
327 | // Normalize. |
328 | const float normalizing_factor = |
329 | static_cast<float>(up_rate_) / |
330 | std::accumulate(first: buffer_channel->begin(), last: buffer_channel->end(), init: 0.0f); |
331 | ScalarMultiply(length: filter_length, gain: normalizing_factor, input: buffer_channel->begin(), |
332 | output: buffer_channel->begin()); |
333 | } |
334 | |
335 | } // namespace vraudio |
336 | |