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/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 | |
34 | namespace vraudio { |
35 | |
36 | namespace { |
37 | |
38 | // FFT length and length of time domain data processed. |
39 | const size_t kFftSize = 4096; |
40 | |
41 | // Length of magnitude and phase buffers. |
42 | const size_t kMagnitudeLength = kFftSize / 2 + 1; |
43 | |
44 | // Number of overlaps per kFftSize time domain chunk of input. |
45 | const size_t kNumOverlap = 4; |
46 | |
47 | // Number of samples per overlap. |
48 | const size_t kOverlapLength = kFftSize / kNumOverlap; |
49 | |
50 | // Three quarters FFT size, used for copying chunks of input. |
51 | const size_t kThreeQuarterFftSize = kFftSize - kOverlapLength; |
52 | |
53 | // Number of channels needed to overlap add output data. |
54 | const size_t kNumQuadChannels = 4; |
55 | |
56 | // Number of buffers of delay applies to the magnitude spectrum. |
57 | const size_t kMagnitudeDelay = 3; |
58 | |
59 | // Length of a buffer of noise used to provide random phase. |
60 | const size_t kNoiseLength = 16384; |
61 | |
62 | // Length of the noise buffer from which a phase buffer can start. |
63 | const 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. |
67 | inline 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 | |
74 | SpectralReverb::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 | |
110 | void 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 | |
123 | void 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 | |
177 | void 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 | |
241 | void 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 | |
257 | void 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 | |
277 | void 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 | |
291 | void 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 | |
327 | void 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 | |