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 | #include "dsp/utils.h" |
18 | |
19 | #include <algorithm> |
20 | #include <cmath> |
21 | #include <limits> |
22 | #include <random> |
23 | #include <vector> |
24 | |
25 | #include "base/constants_and_types.h" |
26 | #include "base/logging.h" |
27 | #include "base/misc_math.h" |
28 | |
29 | #include "dsp/biquad_filter.h" |
30 | #include "dsp/filter_coefficient_generators.h" |
31 | |
32 | namespace { |
33 | |
34 | // The mean and standard deviation of the normal distribution for bandlimited |
35 | // Gaussian noise. |
36 | const float kMean = 0.0f; |
37 | const float kStandardDeviation = 1.0f; |
38 | |
39 | // Maximum group delay in seconds for each filter. In order to avoid audible |
40 | // distortion, the maximum phase shift of a re-combined stereo sequence should |
41 | // not exceed 5ms at high frequencies. That is why, maximum phase shift of |
42 | // each filter is set to 1/2 of that value. |
43 | const float kMaxGroupDelaySeconds = 0.0025f; |
44 | |
45 | // Phase modulation depth, chosen so that for a given max group delay filters |
46 | // provide the lowest cross-correlation coefficient. |
47 | const float kPhaseModulationDepth = 1.18f; |
48 | |
49 | // Constants used in the generation of uniform random number distributions. |
50 | // https://en.wikipedia.org/wiki/Linear_congruential_generator |
51 | const uint64 kMultiplier = 1664525L; |
52 | const uint64 kIncrement = 1013904223L; |
53 | const float kInt32ToFloat = |
54 | 1.0f / static_cast<float>(std::numeric_limits<uint32>::max()); |
55 | |
56 | } // namespace |
57 | |
58 | namespace vraudio { |
59 | |
60 | void GenerateGaussianNoise(float mean, float std_deviation, unsigned seed, |
61 | AudioBuffer::Channel* noise_channel) { |
62 | DCHECK(noise_channel); |
63 | // First generate uniform noise. |
64 | GenerateUniformNoise(min: 0.0f, max: 1.0f, seed, noise_channel); |
65 | const size_t length = noise_channel->size(); |
66 | |
67 | // Gaussian distribution with mean and standard deviation in pairs via the |
68 | // box-muller transform |
69 | // https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform. |
70 | for (size_t i = 0; i < length - 1; i += 2) { |
71 | const float part_one = std::sqrt(x: -2.0f * std::log(x: (*noise_channel)[i])); |
72 | const float part_two = kTwoPi * (*noise_channel)[i + 1]; |
73 | const float z0 = part_one * std::cos(x: part_two); |
74 | const float z1 = part_one * std::sin(x: part_two); |
75 | (*noise_channel)[i] = std_deviation * z0 + mean; |
76 | (*noise_channel)[i + 1] = std_deviation * z1 + mean; |
77 | } |
78 | // Handle the odd buffer length case cheaply. |
79 | if (length % 2 > 0) { |
80 | (*noise_channel)[length - 1] = (*noise_channel)[0]; |
81 | } |
82 | } |
83 | |
84 | void GenerateUniformNoise(float min, float max, unsigned seed, |
85 | AudioBuffer::Channel* noise_channel) { |
86 | // Simple random generator to avoid the use of std::uniform_real_distribution |
87 | // affected by https://gcc.gnu.org/bugzilla/show_bug.cgi?id=56202 |
88 | DCHECK(noise_channel); |
89 | DCHECK_LT(min, max); |
90 | const float scaled_conversion_factor = kInt32ToFloat * (max - min); |
91 | uint32 state = static_cast<uint32>(seed); |
92 | for (float& sample : *noise_channel) { |
93 | state = static_cast<uint32>(state * kMultiplier + kIncrement); |
94 | sample = min + static_cast<float>(state) * scaled_conversion_factor; |
95 | } |
96 | } |
97 | |
98 | void GenerateBandLimitedGaussianNoise(float center_frequency, int sampling_rate, |
99 | unsigned seed, |
100 | AudioBuffer* noise_buffer) { |
101 | |
102 | |
103 | DCHECK(noise_buffer); |
104 | DCHECK_GT(sampling_rate, 0); |
105 | DCHECK_LT(center_frequency, static_cast<float>(sampling_rate) / 2.0f); |
106 | const size_t num_frames = noise_buffer->num_frames(); |
107 | |
108 | BiquadCoefficients bandpass_coefficients = ComputeBandPassBiquadCoefficients( |
109 | sample_rate: sampling_rate, centre_frequency: center_frequency, /*bandwidth=*/1); |
110 | BiquadFilter bandpass_filter(bandpass_coefficients, num_frames); |
111 | |
112 | for (auto& channel : *noise_buffer) { |
113 | GenerateGaussianNoise(mean: kMean, std_deviation: kStandardDeviation, seed, noise_channel: &channel); |
114 | bandpass_filter.Filter(input_channel: channel, output_channel: &channel); |
115 | bandpass_filter.Clear(); |
116 | } |
117 | } |
118 | |
119 | std::unique_ptr<AudioBuffer> GenerateDecorrelationFilters(int sampling_rate) { |
120 | |
121 | const int kMaxGroupDelaySamples = static_cast<int>( |
122 | roundf(x: kMaxGroupDelaySeconds * static_cast<float>(sampling_rate))); |
123 | |
124 | // Filter coefficients according to: |
125 | // [1] F. Zotter, M. Frank, "Efficient Phantom Source Widening", Archives of |
126 | // Acoustics, Vol. 38, No. 1, pp. 27–37 (2013). |
127 | const float g0 = 1.0f - 0.25f * IntegerPow(base: kPhaseModulationDepth, exp: 2); |
128 | const float g1 = 0.5f * kPhaseModulationDepth - |
129 | 0.0625f * IntegerPow(base: kPhaseModulationDepth, exp: 3); |
130 | const float g2 = 0.1250f * IntegerPow(base: kPhaseModulationDepth, exp: 2); |
131 | std::vector<float> filter1_coefficients{g2, g1, g0, -g1, g2}; |
132 | std::vector<float> filter2_coefficients{g2, -g1, g0, g1, g2}; |
133 | |
134 | const size_t filter_length = |
135 | filter1_coefficients.size() * kMaxGroupDelaySamples; |
136 | std::unique_ptr<AudioBuffer> decorrelation_filters( |
137 | new AudioBuffer(kNumStereoChannels, filter_length)); |
138 | decorrelation_filters->Clear(); |
139 | |
140 | for (size_t coefficient = 0; coefficient < filter1_coefficients.size(); |
141 | ++coefficient) { |
142 | (*decorrelation_filters)[0][coefficient * kMaxGroupDelaySamples] = |
143 | filter1_coefficients[coefficient]; |
144 | (*decorrelation_filters)[1][coefficient * kMaxGroupDelaySamples] = |
145 | filter2_coefficients[coefficient]; |
146 | } |
147 | |
148 | return decorrelation_filters; |
149 | } |
150 | |
151 | size_t GetNumReverbOctaveBands(int sampling_rate) { |
152 | DCHECK_GT(sampling_rate, 0); |
153 | |
154 | const float max_band = |
155 | log2f(x: 0.5f * static_cast<float>(sampling_rate) / kLowestOctaveBandHz); |
156 | return std::min(a: kNumReverbOctaveBands, b: static_cast<size_t>(roundf(x: max_band))); |
157 | } |
158 | |
159 | size_t GetNumSamplesFromMilliseconds(float milliseconds, int sampling_rate) { |
160 | DCHECK_GE(milliseconds, 0.0f); |
161 | DCHECK_GT(sampling_rate, 0); |
162 | return static_cast<size_t>(milliseconds * kSecondsFromMilliseconds * |
163 | static_cast<float>(sampling_rate)); |
164 | } |
165 | |
166 | size_t CeilToMultipleOfFramesPerBuffer(size_t size, size_t frames_per_buffer) { |
167 | DCHECK_NE(frames_per_buffer, 0U); |
168 | const size_t remainder = size % frames_per_buffer; |
169 | return remainder == 0 ? std::max(a: size, b: frames_per_buffer) |
170 | : size + frames_per_buffer - remainder; |
171 | } |
172 | |
173 | void GenerateHannWindow(bool full_window, size_t window_length, |
174 | AudioBuffer::Channel* buffer) { |
175 | |
176 | DCHECK(buffer); |
177 | DCHECK_LE(window_length, buffer->size()); |
178 | const float full_window_scaling_factor = |
179 | kTwoPi / (static_cast<float>(window_length) - 1.0f); |
180 | const float half_window_scaling_factor = |
181 | kTwoPi / (2.0f * static_cast<float>(window_length) - 1.0f); |
182 | const float scaling_factor = |
183 | (full_window) ? full_window_scaling_factor : half_window_scaling_factor; |
184 | for (size_t i = 0; i < window_length; ++i) { |
185 | (*buffer)[i] = |
186 | 0.5f * (1.0f - std::cos(x: scaling_factor * static_cast<float>(i))); |
187 | } |
188 | } |
189 | |
190 | } // namespace vraudio |
191 | |