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/fft_manager.h" |
18 | |
19 | #include <algorithm> |
20 | |
21 | #include "base/constants_and_types.h" |
22 | #include "base/logging.h" |
23 | #include "base/misc_math.h" |
24 | #include "base/simd_utils.h" |
25 | |
26 | |
27 | // Prevent Visual Studio from complaining about std::copy_n. |
28 | #if defined(_WIN32) |
29 | #define _SCL_SECURE_NO_WARNINGS |
30 | #endif |
31 | |
32 | namespace vraudio { |
33 | |
34 | namespace { |
35 | |
36 | // for |fft_size_|s less than 2^14, the stack is used, on the reccomendation of |
37 | // the author of the pffft library. |
38 | const size_t kPffftMaxStackSize = 16384; |
39 | |
40 | } // namespace |
41 | |
42 | // The pffft implementation requires a minimum fft size of 32 samples. |
43 | const size_t FftManager::kMinFftSize = 32; |
44 | |
45 | FftManager::FftManager(size_t frames_per_buffer) |
46 | : fft_size_(std::max(a: NextPowTwo(input: frames_per_buffer) * 2, b: kMinFftSize)), |
47 | frames_per_buffer_(frames_per_buffer), |
48 | inverse_fft_scale_(1.0f / static_cast<float>(fft_size_)), |
49 | temp_zeropad_buffer_(kNumMonoChannels, fft_size_), |
50 | temp_freq_buffer_(kNumMonoChannels, fft_size_) { |
51 | DCHECK_GT(frames_per_buffer, 0U); |
52 | DCHECK_GE(fft_size_, kMinFftSize); |
53 | DCHECK(!(fft_size_ & (fft_size_ - 1))); // Ensure its a power of two. |
54 | // Suggested pffft initialization. |
55 | if (fft_size_ > kPffftMaxStackSize) { |
56 | // Allocate memory for work space factors etc, Size reccomended by pffft. |
57 | const size_t num_bytes = 2 * fft_size_ * sizeof(float); |
58 | pffft_workspace_ = |
59 | reinterpret_cast<float*>(pffft_aligned_malloc(nb_bytes: num_bytes)); |
60 | } |
61 | |
62 | fft_ = pffft_new_setup(N: static_cast<int>(fft_size_), transform: PFFFT_REAL); |
63 | |
64 | temp_zeropad_buffer_.Clear(); |
65 | } |
66 | |
67 | FftManager::~FftManager() { |
68 | pffft_destroy_setup(fft_); |
69 | if (pffft_workspace_ != nullptr) { |
70 | pffft_aligned_free(pffft_workspace_); |
71 | } |
72 | } |
73 | |
74 | void FftManager::FreqFromTimeDomain(const AudioBuffer::Channel& time_channel, |
75 | AudioBuffer::Channel* freq_channel) { |
76 | |
77 | |
78 | DCHECK(freq_channel); |
79 | DCHECK_EQ(freq_channel->size(), fft_size_); |
80 | DCHECK_LE(time_channel.size(), fft_size_); |
81 | |
82 | // Perform forward FFT transform. |
83 | if (time_channel.size() == fft_size_) { |
84 | pffft_transform(setup: fft_, input: time_channel.begin(), output: freq_channel->begin(), |
85 | work: pffft_workspace_, direction: PFFFT_FORWARD); |
86 | } else { |
87 | std::copy_n(first: time_channel.begin(), n: frames_per_buffer_, |
88 | result: temp_zeropad_buffer_[0].begin()); |
89 | pffft_transform(setup: fft_, input: temp_zeropad_buffer_[0].begin(), |
90 | output: freq_channel->begin(), work: pffft_workspace_, direction: PFFFT_FORWARD); |
91 | } |
92 | } |
93 | |
94 | void FftManager::TimeFromFreqDomain(const AudioBuffer::Channel& freq_channel, |
95 | AudioBuffer::Channel* time_channel) { |
96 | |
97 | |
98 | DCHECK(time_channel); |
99 | DCHECK_EQ(freq_channel.size(), fft_size_); |
100 | |
101 | // Perform reverse FFT transform. |
102 | const size_t time_channel_size = time_channel->size(); |
103 | if (time_channel_size == fft_size_) { |
104 | pffft_transform(setup: fft_, input: freq_channel.begin(), output: time_channel->begin(), |
105 | work: pffft_workspace_, direction: PFFFT_BACKWARD); |
106 | } else { |
107 | DCHECK_EQ(time_channel_size, frames_per_buffer_); |
108 | auto& temp_channel = temp_freq_buffer_[0]; |
109 | pffft_transform(setup: fft_, input: freq_channel.begin(), output: temp_channel.begin(), |
110 | work: pffft_workspace_, direction: PFFFT_BACKWARD); |
111 | std::copy_n(first: temp_channel.begin(), n: frames_per_buffer_, |
112 | result: time_channel->begin()); |
113 | } |
114 | } |
115 | |
116 | void FftManager::ApplyReverseFftScaling(AudioBuffer::Channel* time_channel) { |
117 | DCHECK(time_channel->size() == frames_per_buffer_ || |
118 | time_channel->size() == fft_size_); |
119 | // Normalization must be performed here as we normally do this as part of the |
120 | // convolution. |
121 | ScalarMultiply(length: time_channel->size(), gain: inverse_fft_scale_, |
122 | input: time_channel->begin(), output: time_channel->begin()); |
123 | } |
124 | |
125 | void FftManager::GetCanonicalFormatFreqBuffer(const AudioBuffer::Channel& input, |
126 | AudioBuffer::Channel* output) { |
127 | DCHECK_EQ(input.size(), fft_size_); |
128 | DCHECK_EQ(output->size(), fft_size_); |
129 | pffft_zreorder(setup: fft_, input: input.begin(), output: output->begin(), direction: PFFFT_FORWARD); |
130 | } |
131 | |
132 | void FftManager::GetPffftFormatFreqBuffer(const AudioBuffer::Channel& input, |
133 | AudioBuffer::Channel* output) { |
134 | DCHECK_EQ(input.size(), fft_size_); |
135 | DCHECK_EQ(output->size(), fft_size_); |
136 | pffft_zreorder(setup: fft_, input: input.begin(), output: output->begin(), direction: PFFFT_BACKWARD); |
137 | } |
138 | |
139 | void FftManager::MagnitudeFromCanonicalFreqBuffer( |
140 | const AudioBuffer::Channel& freq_channel, |
141 | AudioBuffer::Channel* magnitude_channel) { |
142 | DCHECK(magnitude_channel); |
143 | DCHECK_EQ(freq_channel.size(), fft_size_); |
144 | DCHECK_EQ(magnitude_channel->size(), frames_per_buffer_ + 1); |
145 | |
146 | (*magnitude_channel)[0] = std::abs(x: freq_channel[0]); |
147 | ApproxComplexMagnitude(length: frames_per_buffer_ - 1, input: freq_channel.begin() + 2, |
148 | output: magnitude_channel->begin() + 1); |
149 | (*magnitude_channel)[frames_per_buffer_] = std::abs(x: freq_channel[1]); |
150 | } |
151 | |
152 | void FftManager::CanonicalFreqBufferFromMagnitudeAndPhase( |
153 | const AudioBuffer::Channel& magnitude_channel, |
154 | const AudioBuffer::Channel& phase_channel, |
155 | AudioBuffer::Channel* canonical_freq_channel) { |
156 | DCHECK(canonical_freq_channel); |
157 | DCHECK_EQ(magnitude_channel.size(), frames_per_buffer_ + 1); |
158 | DCHECK_EQ(phase_channel.size(), frames_per_buffer_ + 1); |
159 | DCHECK_EQ(canonical_freq_channel->size(), fft_size_); |
160 | |
161 | (*canonical_freq_channel)[0] = magnitude_channel[0]; |
162 | (*canonical_freq_channel)[1] = -magnitude_channel[frames_per_buffer_]; |
163 | for (size_t i = 1, j = 2; i < frames_per_buffer_; ++i, j += 2) { |
164 | (*canonical_freq_channel)[j] = |
165 | magnitude_channel[i] * std::cos(x: phase_channel[i]); |
166 | (*canonical_freq_channel)[j + 1] = |
167 | magnitude_channel[i] * std::sin(x: phase_channel[i]); |
168 | } |
169 | } |
170 | |
171 | void FftManager::CanonicalFreqBufferFromMagnitudeAndSinCosPhase( |
172 | size_t phase_offset, const AudioBuffer::Channel& magnitude_channel, |
173 | const AudioBuffer::Channel& sin_phase_channel, |
174 | const AudioBuffer::Channel& cos_phase_channel, |
175 | AudioBuffer::Channel* canonical_freq_channel) { |
176 | static const size_t kSimdLength = 4; |
177 | DCHECK(canonical_freq_channel); |
178 | DCHECK_EQ(magnitude_channel.size(), frames_per_buffer_ + 1); |
179 | DCHECK_GE(sin_phase_channel.size() + phase_offset, frames_per_buffer_ + 1); |
180 | DCHECK_GE(cos_phase_channel.size() + phase_offset, frames_per_buffer_ + 1); |
181 | DCHECK_EQ(canonical_freq_channel->size(), fft_size_); |
182 | |
183 | (*canonical_freq_channel)[0] = magnitude_channel[0]; |
184 | (*canonical_freq_channel)[1] = -magnitude_channel[frames_per_buffer_]; |
185 | // Continue on till we can gaurantee alignment in our audio buffer. |
186 | for (size_t i = 1, j = 2; i <= kSimdLength; ++i, j += 2) { |
187 | (*canonical_freq_channel)[j] = |
188 | magnitude_channel[i] * cos_phase_channel[i + phase_offset]; |
189 | (*canonical_freq_channel)[j + 1] = |
190 | magnitude_channel[i] * sin_phase_channel[i + phase_offset]; |
191 | } |
192 | ComplexInterleavedFormatFromMagnitudeAndSinCosPhase( |
193 | length: 2 * (frames_per_buffer_ - kSimdLength), magnitude: &magnitude_channel[kSimdLength], |
194 | cos_phase: &cos_phase_channel[kSimdLength + phase_offset], |
195 | sin_phase: &sin_phase_channel[kSimdLength + phase_offset], |
196 | complex_interleaved_format_output: &(*canonical_freq_channel)[2 * kSimdLength]); |
197 | } |
198 | |
199 | void FftManager::FreqDomainConvolution(const AudioBuffer::Channel& input_a, |
200 | const AudioBuffer::Channel& input_b, |
201 | AudioBuffer::Channel* scaled_output) { |
202 | DCHECK_EQ(input_a.size(), fft_size_); |
203 | DCHECK_EQ(input_b.size(), fft_size_); |
204 | DCHECK_EQ(scaled_output->size(), fft_size_); |
205 | pffft_zconvolve_accumulate(setup: fft_, dft_a: input_a.begin(), dft_b: input_b.begin(), |
206 | dft_ab: scaled_output->begin(), scaling: inverse_fft_scale_); |
207 | } |
208 | |
209 | } // namespace vraudio |
210 | |