1/*
2Copyright 2018 Google Inc. All Rights Reserved.
3
4Licensed under the Apache License, Version 2.0 (the "License");
5you may not use this file except in compliance with the License.
6You may obtain a copy of the License at
7
8 http://www.apache.org/licenses/LICENSE-2.0
9
10Unless required by applicable law or agreed to in writing, software
11distributed under the License is distributed on an "AS-IS" BASIS,
12WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13See the License for the specific language governing permissions and
14limitations 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
34namespace vraudio {
35
36namespace {
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.
42const size_t kTransitionBandwidthRatio = 13;
43
44// Maximum number of channels internally based upon the maximum supported
45// ambisonic order.
46const size_t kMaxNumChannels =
47 (kMaxSupportedAmbisonicOrder + 1) * (kMaxSupportedAmbisonicOrder + 1);
48
49} // namespace
50
51Resampler::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
64void 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
154size_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
166size_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
178void 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
217bool 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
230void Resampler::ResetState() {
231
232 time_modulo_up_rate_ = 0;
233 last_processed_sample_ = 0;
234 state_.Clear();
235}
236
237void 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
262void 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
282void 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
304void 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

source code of qtmultimedia/src/3rdparty/resonance-audio/resonance_audio/dsp/resampler.cc