| 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 | |