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/spectral_reverb.h"
23
24#include <algorithm>
25#include <cstdlib>
26#include <numeric>
27
28#include "base/constants_and_types.h"
29#include "base/logging.h"
30#include "base/simd_utils.h"
31#include "dsp/spectral_reverb_constants_and_tables.h"
32#include "dsp/utils.h"
33
34namespace vraudio {
35
36namespace {
37
38// FFT length and length of time domain data processed.
39const size_t kFftSize = 4096;
40
41// Length of magnitude and phase buffers.
42const size_t kMagnitudeLength = kFftSize / 2 + 1;
43
44// Number of overlaps per kFftSize time domain chunk of input.
45const size_t kNumOverlap = 4;
46
47// Number of samples per overlap.
48const size_t kOverlapLength = kFftSize / kNumOverlap;
49
50// Three quarters FFT size, used for copying chunks of input.
51const size_t kThreeQuarterFftSize = kFftSize - kOverlapLength;
52
53// Number of channels needed to overlap add output data.
54const size_t kNumQuadChannels = 4;
55
56// Number of buffers of delay applies to the magnitude spectrum.
57const size_t kMagnitudeDelay = 3;
58
59// Length of a buffer of noise used to provide random phase.
60const size_t kNoiseLength = 16384;
61
62// Length of the noise buffer from which a phase buffer can start.
63const size_t kAvailableNoiselength = kNoiseLength - kMagnitudeLength;
64
65// Returns a random integer in the range provided. Used to index into the random
66// buffer for phase values.
67inline size_t GetRandomIntegerInRange(size_t min, size_t max) {
68 DCHECK_GE(max, min);
69 return min + static_cast<size_t>(std::rand()) % (max - min);
70}
71
72} // namespace
73
74SpectralReverb::SpectralReverb(int sample_rate, size_t frames_per_buffer)
75 : sample_rate_(sample_rate),
76 frames_per_buffer_(frames_per_buffer),
77 magnitude_delay_index_(0),
78 overlap_add_index_(0),
79 fft_manager_(kFftSize / 2),
80 sin_cos_random_phase_buffer_(kNumStereoChannels, kNoiseLength),
81 unscaled_window_(kNumMonoChannels, kFftSize),
82 window_(kNumMonoChannels, kFftSize),
83 feedback_(kNumMonoChannels, kMagnitudeLength),
84 magnitude_compensation_(kNumMonoChannels, kMagnitudeLength),
85 magnitude_delay_(kMagnitudeDelay, kMagnitudeLength),
86 fft_size_input_(kNumMonoChannels, kFftSize),
87 input_circular_buffer_(kFftSize + frames_per_buffer_, frames_per_buffer_,
88 kOverlapLength),
89 output_circular_buffers_(kNumStereoChannels),
90 out_time_buffer_(kNumQuadChannels, kFftSize),
91 temp_freq_buffer_(kNumStereoChannels, kFftSize),
92 scaled_magnitude_buffer_(kNumMonoChannels, kMagnitudeLength),
93 temp_magnitude_buffer_(kNumMonoChannels, kMagnitudeLength),
94 temp_phase_buffer_(kNumStereoChannels, kMagnitudeLength),
95 output_accumulator_(kNumStereoChannels),
96 is_gain_near_zero_(false),
97 is_feedback_near_zero_(false) {
98 DCHECK_GT(sample_rate, 0);
99 DCHECK_GT(frames_per_buffer_, 0U);
100
101 // Seed std::rand, used for phase selection.
102 std::srand(seed: 1);
103 GenerateRandomPhaseBuffer();
104 GenerateAnalysisWindow();
105 InitializeCircularBuffersAndAccumulators();
106 fft_size_input_.Clear();
107 magnitude_compensation_.Clear();
108}
109
110void SpectralReverb::SetGain(float gain) {
111 DCHECK_GE(gain, 0.0f);
112 ScalarMultiply(length: window_.num_frames(), gain, input: &unscaled_window_[0][0],
113 output: &window_[0][0]);
114 // If the gain is less than -60dB we will bypass all processing.
115 is_gain_near_zero_ = gain <= kNegative60dbInAmplitude;
116 // If we are to bypass processing we clear the circular buffer so that we
117 // don't process old input when the spectral reverb is restarted.
118 if (is_gain_near_zero_ || is_feedback_near_zero_) {
119 input_circular_buffer_.Clear();
120 }
121}
122
123void SpectralReverb::SetRt60PerOctaveBand(const float* rt60_values) {
124 DCHECK(rt60_values);
125 const float sample_rate_float = static_cast<float>(sample_rate_);
126 // Fill the entries in the feedback channel with the feedback values up to
127 // that frequency. The feedback channels entries are per frequency in the same
128 // way as the magnitude vectors. Also fill the magnitude compensation vector
129 // such that the peak value at each frequency for any given reverberation time
130 // will be |kDefaultReverbGain|.
131 AudioBuffer::Channel* feedback_channel = &feedback_[0];
132 feedback_channel->Clear();
133 AudioBuffer::Channel* magnitude_compensation_channel =
134 &magnitude_compensation_[0];
135 magnitude_compensation_channel->Clear();
136 const float frequency_step = sample_rate_float / static_cast<float>(kFftSize);
137 int index = GetFeedbackIndexFromRt60(reverberation_time: rt60_values[0], sample_rate: sample_rate_float);
138 float current_feedback =
139 index == kInvalidIndex ? 0.0f : kSpectralReverbFeedback[index];
140 float magnitude_compensation_value =
141 index == kInvalidIndex ? 0.0f : kMagnitudeCompensation[index];
142 const size_t max_frequency_bin =
143 std::min(a: static_cast<size_t>(
144 (kOctaveBandCentres[kNumReverbOctaveBands - 1] * kSqrtTwo) /
145 frequency_step),
146 b: feedback_channel->size());
147 // The upper edge of the octave band is sqrt(2) * centre_frequency :
148 // https://en.wikipedia.org/wiki/Octave_band#Octave_Bands
149 float upper_octave_band_edge = kOctaveBandCentres[0] * kSqrtTwo;
150 for (size_t i = 0, j = 0; i < max_frequency_bin; ++i) {
151 const float current_freq = static_cast<float>(i) * frequency_step;
152 if (current_freq > upper_octave_band_edge) {
153 ++j;
154 DCHECK_LT(j, kNumReverbOctaveBands);
155 upper_octave_band_edge = kOctaveBandCentres[j] * kSqrtTwo;
156 index = GetFeedbackIndexFromRt60(reverberation_time: rt60_values[j], sample_rate: sample_rate_float);
157 current_feedback =
158 index == kInvalidIndex ? 0.0f : kSpectralReverbFeedback[index];
159 magnitude_compensation_value =
160 index == kInvalidIndex ? 0.0f : kMagnitudeCompensation[index];
161 }
162 (*feedback_channel)[i] = current_feedback;
163 (*magnitude_compensation_channel)[i] = magnitude_compensation_value;
164 }
165 // If the sum of all feedback values is below the minimum feedback value, it
166 // is safe to assume we can bypass the spectral reverb entirely.
167 is_feedback_near_zero_ =
168 std::accumulate(first: feedback_channel->begin(), last: feedback_channel->end(),
169 init: 0.0f) < kSpectralReverbFeedback[0];
170 // If we are to bypass processing we clear the circular buffer so that we
171 // don't process old input when the spectral reverb is restarted.
172 if (is_gain_near_zero_ || is_feedback_near_zero_) {
173 input_circular_buffer_.Clear();
174 }
175}
176
177void SpectralReverb::Process(const AudioBuffer::Channel& input,
178 AudioBuffer::Channel* left_out,
179 AudioBuffer::Channel* right_out) {
180 DCHECK_EQ(input.size(), left_out->size());
181 DCHECK_EQ(input.size(), right_out->size());
182 DCHECK_EQ(input.size(), frames_per_buffer_);
183
184
185 if (is_gain_near_zero_ || is_feedback_near_zero_) {
186 left_out->Clear();
187 right_out->Clear();
188 return;
189 }
190
191 // Here we insert |frames_per_buffer_| samples on each function call. Then,
192 // once there are |kOverlapLength| samples in the input circular buffer we
193 // retrieve |kOverlapLength| samples from it and process |kFftSize| samples of
194 // input at a time, sliding along by |kOverlapLength| samples. We then place
195 // |kOverlapLength| samples into the output buffer and we will extract
196 // |frames_per_buffer_| samples from the output buffer on each function call.
197 input_circular_buffer_.InsertBuffer(input);
198 while (input_circular_buffer_.GetOccupancy() >= kOverlapLength) {
199 std::copy_n(first: &fft_size_input_[0][kOverlapLength], n: kThreeQuarterFftSize,
200 result: &fft_size_input_[0][0]);
201 input_circular_buffer_.RetrieveBufferWithOffset(offset: kThreeQuarterFftSize,
202 output: &fft_size_input_[0]);
203 fft_manager_.FreqFromTimeDomain(time_channel: fft_size_input_[0], freq_channel: &temp_freq_buffer_[0]);
204 fft_manager_.GetCanonicalFormatFreqBuffer(input: temp_freq_buffer_[0],
205 output: &temp_freq_buffer_[1]);
206 fft_manager_.MagnitudeFromCanonicalFreqBuffer(freq_channel: temp_freq_buffer_[1],
207 magnitude_channel: &scaled_magnitude_buffer_[0]);
208 // Apply the magnitude compensation to the input magnitude spectrum before
209 // feedback is applied.
210 MultiplyPointwise(length: kMagnitudeLength, input_a: magnitude_compensation_[0].begin(),
211 input_b: scaled_magnitude_buffer_[0].begin(),
212 output: scaled_magnitude_buffer_[0].begin());
213 // Generate time domain reverb blocks.
214 GetNextReverbBlock(delay_index: magnitude_delay_index_, left_channel: &out_time_buffer_[0],
215 right_channel: &out_time_buffer_[1]);
216 magnitude_delay_index_ = (magnitude_delay_index_ + 1) % kMagnitudeDelay;
217 GetNextReverbBlock(delay_index: magnitude_delay_index_, left_channel: &out_time_buffer_[2],
218 right_channel: &out_time_buffer_[3]);
219
220 // Combine the reverb blocks for both left and right output.
221 AddPointwise(length: kFftSize, input_a: out_time_buffer_[0].begin(),
222 input_b: out_time_buffer_[2].begin(), output: out_time_buffer_[0].begin());
223 AddPointwise(length: kFftSize, input_a: out_time_buffer_[1].begin(),
224 input_b: out_time_buffer_[3].begin(), output: out_time_buffer_[1].begin());
225
226 // Window the left and right output (While applying inverse FFT scaling).
227 MultiplyPointwise(length: kFftSize, input_a: out_time_buffer_[0].begin(), input_b: window_[0].begin(),
228 output: out_time_buffer_[0].begin());
229 MultiplyPointwise(length: kFftSize, input_a: out_time_buffer_[1].begin(), input_b: window_[0].begin(),
230 output: out_time_buffer_[1].begin());
231
232 // Next perform the addition and the submission into the output buffer.
233 AccumulateOverlap(channel_index: 0 /*channel*/, buffer: out_time_buffer_[0]);
234 AccumulateOverlap(channel_index: 1 /*channel*/, buffer: out_time_buffer_[1]);
235 overlap_add_index_ = (overlap_add_index_ + 1) % kNumOverlap;
236 }
237 output_circular_buffers_[0]->RetrieveBuffer(output: left_out);
238 output_circular_buffers_[1]->RetrieveBuffer(output: right_out);
239}
240
241void SpectralReverb::AccumulateOverlap(size_t channel_index,
242 const AudioBuffer::Channel& buffer) {
243 // Use a modulo indexed multi channel audio buffer with each channel of length
244 // |kOverlapLength| to perform an overlap add.
245 for (size_t i = 0, index = overlap_add_index_; i < kNumOverlap;
246 ++i, index = (index + 1) % kNumOverlap) {
247 float* accumulator_start_point =
248 output_accumulator_[channel_index][index].begin();
249 AddPointwise(length: kOverlapLength, input_a: buffer.begin() + i * kOverlapLength,
250 input_b: accumulator_start_point, output: accumulator_start_point);
251 }
252 output_circular_buffers_[channel_index]->InsertBuffer(
253 input: output_accumulator_[channel_index][overlap_add_index_]);
254 output_accumulator_[channel_index][overlap_add_index_].Clear();
255}
256
257void SpectralReverb::GenerateAnalysisWindow() {
258 // Genarate a pseudo tukey window from three overlapping hann windows, scaled
259 // by the inverse fft scale.
260 AudioBuffer::Channel* window_channel = &window_[0];
261 // Use the unscaled window buffer as temporary storage.
262 GenerateHannWindow(full_window: true /* full */, window_length: kMagnitudeLength, buffer: &unscaled_window_[0]);
263 float* hann_window = &unscaled_window_[0][0];
264 // Scale the hann window such that the sum of three will have unity peak.
265 const float kThreeQuarters = 0.75f;
266 ScalarMultiply(length: kMagnitudeLength, gain: kThreeQuarters, input: hann_window, output: hann_window);
267 for (size_t offset = 0; offset < kThreeQuarterFftSize;
268 offset += kOverlapLength) {
269 float* tripple_hann_window = &(*window_channel)[offset];
270 AddPointwise(length: kMagnitudeLength, input_a: hann_window, input_b: tripple_hann_window,
271 output: tripple_hann_window);
272 }
273 fft_manager_.ApplyReverseFftScaling(time_channel: window_channel);
274 unscaled_window_[0] = *window_channel;
275}
276
277void SpectralReverb::GenerateRandomPhaseBuffer() {
278 AudioBuffer::Channel* sin_phase_channel = &sin_cos_random_phase_buffer_[0];
279 AudioBuffer::Channel* cos_phase_channel = &sin_cos_random_phase_buffer_[1];
280 // Initially use the sin channel to store the random data before taking sine
281 // and cosine.
282 GenerateUniformNoise(/*min=*/0.0f, /*max=*/kPi, /*seed=*/1U,
283 noise_channel: sin_phase_channel);
284 for (size_t i = 0; i < sin_cos_random_phase_buffer_.num_frames(); ++i) {
285
286 (*cos_phase_channel)[i] = std::cos(x: (*sin_phase_channel)[i]);
287 (*sin_phase_channel)[i] = std::sin(x: (*sin_phase_channel)[i]);
288 }
289}
290
291void SpectralReverb::GetNextReverbBlock(size_t delay_index,
292 AudioBuffer::Channel* left_channel,
293 AudioBuffer::Channel* right_channel) {
294 DCHECK(left_channel);
295 DCHECK(right_channel);
296
297 // Generate reverb magnitude by combining the delayed magnitude with the
298 // current magnitude.
299 AudioBuffer::Channel* temp_magnitude_channel = &temp_magnitude_buffer_[0];
300 *temp_magnitude_channel = scaled_magnitude_buffer_[0];
301 MultiplyAndAccumulatePointwise(
302 length: kMagnitudeLength, input_a: magnitude_delay_[delay_index].begin(),
303 input_b: feedback_[0].begin(), accumulator: temp_magnitude_channel->begin());
304 // Reinsert this new reverb magnitude into the delay buffer.
305 magnitude_delay_[delay_index] = *temp_magnitude_channel;
306
307 for (size_t i = 0; i < kNumStereoChannels; ++i) {
308 // Extract a random phase buffer.
309 const size_t random_offset =
310 GetRandomIntegerInRange(min: 0, max: kAvailableNoiselength);
311 // We gaurantee an aligned offset as when SSE is used we need it.
312 const size_t phase_offset = FindNextAlignedArrayIndex(
313 length: random_offset, type_size_bytes: sizeof(float), memory_alignment_bytes: kMemoryAlignmentBytes);
314 // Convert from magnitude and phase to a time domain output.
315 fft_manager_.CanonicalFreqBufferFromMagnitudeAndSinCosPhase(
316 phase_offset, magnitude_channel: (*temp_magnitude_channel),
317 sin_phase_channel: sin_cos_random_phase_buffer_[0], cos_phase_channel: sin_cos_random_phase_buffer_[1],
318 canonical_freq_channel: &temp_freq_buffer_[0]);
319
320 AudioBuffer::Channel* out_channel = i == 0 ? left_channel : right_channel;
321 fft_manager_.GetPffftFormatFreqBuffer(input: temp_freq_buffer_[0],
322 output: &temp_freq_buffer_[1]);
323 fft_manager_.TimeFromFreqDomain(freq_channel: temp_freq_buffer_[1], time_channel: out_channel);
324 }
325}
326
327void SpectralReverb::InitializeCircularBuffersAndAccumulators() {
328 AudioBuffer zeros(kNumMonoChannels, kOverlapLength);
329 zeros.Clear();
330 for (size_t channel = 0; channel < kNumStereoChannels; ++channel) {
331 // Prefill the |output_circular_buffers_| with |kOverlapLength| /
332 // |frames_per_buffer_| calls worth of zeros.
333 output_circular_buffers_[channel].reset(
334 p: new CircularBuffer(kOverlapLength + frames_per_buffer_, kOverlapLength,
335 frames_per_buffer_));
336 // Due to differences in the |frames_per_buffer_| used for input and output
337 // and |kOverlapLength| used for processing, a certain number of buffers of
338 // zeros must be inserted into the output buffers such that enough input can
339 // build up to process |kOverlapLength| worth, and enough output will build
340 // up to return |frames_per_buffer_| worth.
341 const size_t zeroed_buffers_of_output = kOverlapLength / frames_per_buffer_;
342 for (size_t i = 0; i < zeroed_buffers_of_output; ++i) {
343 output_circular_buffers_[channel]->InsertBuffer(input: zeros[0]);
344 }
345 output_accumulator_[channel] = AudioBuffer(kNumOverlap, kOverlapLength);
346 output_accumulator_[channel].Clear();
347 }
348}
349
350} // namespace vraudio
351

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