diff options
Diffstat (limited to 'alc/effects')
-rw-r--r-- | alc/effects/convolution.cpp | 132 |
1 files changed, 92 insertions, 40 deletions
diff --git a/alc/effects/convolution.cpp b/alc/effects/convolution.cpp index f46422d4..c7a342dc 100644 --- a/alc/effects/convolution.cpp +++ b/alc/effects/convolution.cpp @@ -34,6 +34,7 @@ #include "core/fmt_traits.h" #include "core/mixer.h" #include "intrusive_ptr.h" +#include "pffft.h" #include "polyphase_resampler.h" #include "vector.h" @@ -45,6 +46,10 @@ namespace { * each segment has an FFT applied with a 256-sample buffer (the latter half * left silent) to get its frequency-domain response. The resulting response * has its positive/non-mirrored frequencies saved (129 bins) in each segment. + * Note that since the 0th and 129th bins are real for a real signal, their + * imaginary components are always 0 and can be dropped, allowing their real + * components to be combined so only 128 complex values are stored for the 129 + * bins. * * Input samples are similarly broken up into 128-sample segments, with an FFT * applied to each new incoming segment to get its 129 bins. A history of FFT'd @@ -52,8 +57,8 @@ namespace { * * To apply the reverberation, each impulse response segment is convolved with * its paired input segment (using complex multiplies, far cheaper than FIRs), - * accumulating into a 256-bin FFT buffer. The input history is then shifted to - * align with later impulse response segments for next time. + * accumulating into a 129-bin FFT buffer. The input history is then shifted to + * align with later impulse response segments for the next input segment. * * An inverse FFT is then applied to the accumulated FFT buffer to get a 256- * sample time-domain response for output, which is split in two halves. The @@ -183,6 +188,35 @@ void apply_fir(al::span<float> dst, const float *RESTRICT src, const float *REST #endif } + +template<typename T> +struct AlignedDeleter { }; + +template<typename T> +struct AlignedDeleter<T[]> { + static_assert(std::is_trivially_destructible_v<T>); + using type = T; + + void operator()(T *ptr) { al_free(ptr); } +}; +template<typename T> +using AlignedUPtr = std::unique_ptr<T,AlignedDeleter<T>>; + +template<typename T, size_t A> +auto MakeAlignedPtr(size_t count) -> AlignedUPtr<T> +{ + using Type = typename AlignedDeleter<T>::type; + void *ptr{al_calloc(A, sizeof(Type)*count)}; + return AlignedUPtr<T>{::new(ptr) Type[count]}; +} + + +struct PFFFTSetupDeleter { + void operator()(PFFFT_Setup *ptr) { pffft_destroy_setup(ptr); } +}; +using PFFFTSetupPtr = std::unique_ptr<PFFFT_Setup,PFFFTSetupDeleter>; + + struct ConvolutionState final : public EffectState { FmtChannels mChannels{}; AmbiLayout mAmbiLayout{}; @@ -194,7 +228,9 @@ struct ConvolutionState final : public EffectState { al::vector<std::array<float,ConvolveUpdateSamples>,16> mFilter; al::vector<std::array<float,ConvolveUpdateSamples*2>,16> mOutput; - alignas(16) std::array<complex_f,ConvolveUpdateSize> mFftBuffer{}; + PFFFTSetupPtr mFft{}; + alignas(16) std::array<float,ConvolveUpdateSize> mFftBuffer{}; + alignas(16) std::array<float,ConvolveUpdateSize> mFftWorkBuffer{}; size_t mCurrentSegment{0}; size_t mNumConvolveSegs{0}; @@ -208,7 +244,7 @@ struct ConvolutionState final : public EffectState { }; using ChannelDataArray = al::FlexArray<ChannelData>; std::unique_ptr<ChannelDataArray> mChans; - std::unique_ptr<complex_f[]> mComplexData; + AlignedUPtr<float[]> mComplexData; ConvolutionState() = default; @@ -253,13 +289,17 @@ void ConvolutionState::deviceUpdate(const DeviceBase *device, const BufferStorag using UhjDecoderType = UhjDecoder<512>; static constexpr auto DecoderPadding = UhjDecoderType::sInputPadding; - constexpr uint MaxConvolveAmbiOrder{1u}; + static constexpr uint MaxConvolveAmbiOrder{1u}; + + if(!mFft) + mFft = PFFFTSetupPtr{pffft_new_setup(ConvolveUpdateSize, PFFFT_REAL)}; mFifoPos = 0; mInput.fill(0.0f); decltype(mFilter){}.swap(mFilter); decltype(mOutput){}.swap(mOutput); - mFftBuffer.fill(complex_f{}); + mFftBuffer.fill(0.0f); + mFftWorkBuffer.fill(0.0f); mCurrentSegment = 0; mNumConvolveSegs = 0; @@ -275,7 +315,6 @@ void ConvolutionState::deviceUpdate(const DeviceBase *device, const BufferStorag mAmbiScaling = IsUHJ(mChannels) ? AmbiScaling::UHJ : buffer->mAmbiScaling; mAmbiOrder = minu(buffer->mAmbiOrder, MaxConvolveAmbiOrder); - constexpr size_t m{ConvolveUpdateSize/2 + 1}; const auto bytesPerSample = BytesFromFmt(buffer->mType); const auto realChannels = buffer->channelsFromFmt(); const auto numChannels = (mChannels == FmtUHJ2) ? 3u : ChannelsFromFmt(mChannels, mAmbiOrder); @@ -309,9 +348,9 @@ void ConvolutionState::deviceUpdate(const DeviceBase *device, const BufferStorag mNumConvolveSegs = (resampledCount+(ConvolveUpdateSamples-1)) / ConvolveUpdateSamples; mNumConvolveSegs = maxz(mNumConvolveSegs, 2) - 1; - const size_t complex_length{mNumConvolveSegs * m * (numChannels+1)}; - mComplexData = std::make_unique<complex_f[]>(complex_length); - std::fill_n(mComplexData.get(), complex_length, complex_f{}); + const size_t complex_length{mNumConvolveSegs * ConvolveUpdateSize * (numChannels+1)}; + mComplexData = MakeAlignedPtr<float[],16>(complex_length); + std::fill_n(mComplexData.get(), complex_length, 0.0f); /* Load the samples from the buffer. */ const size_t srclinelength{RoundUp(buffer->mSampleLen+DecoderPadding, 16)}; @@ -332,7 +371,10 @@ void ConvolutionState::deviceUpdate(const DeviceBase *device, const BufferStorag auto ressamples = std::make_unique<double[]>(buffer->mSampleLen + (resampler ? resampledCount : 0)); - complex_f *filteriter = mComplexData.get() + mNumConvolveSegs*m; + auto ffttmp = al::vector<float,16>(ConvolveUpdateSize); + auto fftbuffer = std::vector<std::complex<double>>(ConvolveUpdateSize); + + float *filteriter = mComplexData.get() + mNumConvolveSegs*ConvolveUpdateSize; for(size_t c{0};c < numChannels;++c) { /* Resample to match the device. */ @@ -353,18 +395,34 @@ void ConvolutionState::deviceUpdate(const DeviceBase *device, const BufferStorag std::transform(ressamples.get(), ressamples.get()+first_size, mFilter[c].rbegin(), [](const double d) noexcept -> float { return static_cast<float>(d); }); - auto fftbuffer = std::vector<std::complex<double>>(ConvolveUpdateSize); size_t done{first_size}; for(size_t s{0};s < mNumConvolveSegs;++s) { const size_t todo{minz(resampledCount-done, ConvolveUpdateSamples)}; + /* Apply a double-precision forward FFT for more precise frequency + * measurements. + */ auto iter = std::copy_n(&ressamples[done], todo, fftbuffer.begin()); done += todo; std::fill(iter, fftbuffer.end(), std::complex<double>{}); - forward_fft(al::span{fftbuffer}); - filteriter = std::copy_n(fftbuffer.cbegin(), m, filteriter); + + /* Convert to, and pack in, a float buffer for PFFFT. Note that the + * first bin stores the real component of the half-frequency bin in + * the imaginary component. + */ + for(size_t i{0};i < ConvolveUpdateSamples;++i) + { + ffttmp[i*2 ] = static_cast<float>(fftbuffer[i].real()); + ffttmp[i*2 + 1] = static_cast<float>((i == 0) ? + fftbuffer[ConvolveUpdateSamples+1].real() : fftbuffer[i].imag()); + } + /* Reorder backward to make it suitable for pffft_zconvolve and the + * subsequent pffft_transform(..., PFFFT_BACKWARD). + */ + pffft_zreorder(mFft.get(), ffttmp.data(), al::to_address(filteriter), PFFFT_BACKWARD); + filteriter += ConvolveUpdateSize; } } } @@ -554,7 +612,6 @@ void ConvolutionState::process(const size_t samplesToDo, if(mNumConvolveSegs < 1) UNLIKELY return; - constexpr size_t m{ConvolveUpdateSize/2 + 1}; size_t curseg{mCurrentSegment}; auto &chans = *mChans; @@ -591,54 +648,49 @@ void ConvolutionState::process(const size_t samplesToDo, * frequency bins to the FFT history. */ auto fftiter = std::copy_n(mInput.cbegin(), ConvolveUpdateSamples, mFftBuffer.begin()); - std::fill(fftiter, mFftBuffer.end(), complex_f{}); - forward_fft(al::span{mFftBuffer}); + std::fill(fftiter, mFftBuffer.end(), 0.0f); + pffft_transform(mFft.get(), mFftBuffer.data(), + mComplexData.get() + curseg*ConvolveUpdateSize, mFftWorkBuffer.data(), PFFFT_FORWARD); - std::copy_n(mFftBuffer.cbegin(), m, &mComplexData[curseg*m]); - - const complex_f *RESTRICT filter{mComplexData.get() + mNumConvolveSegs*m}; + const float *RESTRICT filter{mComplexData.get() + mNumConvolveSegs*ConvolveUpdateSize}; for(size_t c{0};c < chans.size();++c) { - std::fill_n(mFftBuffer.begin(), m, complex_f{}); + /* The iFFT'd response is scaled up by the number of bins, so apply + * the inverse to normalize the output. + */ + static constexpr float fftscale{1.0f / float{ConvolveUpdateSize}}; /* Convolve each input segment with its IR filter counterpart * (aligned in time). */ - const complex_f *RESTRICT input{&mComplexData[curseg*m]}; + mFftBuffer.fill(0.0f); + const float *RESTRICT input{&mComplexData[curseg*ConvolveUpdateSize]}; for(size_t s{curseg};s < mNumConvolveSegs;++s) { - for(size_t i{0};i < m;++i,++input,++filter) - mFftBuffer[i] += *input * *filter; + pffft_zconvolve_accumulate(mFft.get(), input, filter, mFftBuffer.data(), fftscale); + input += ConvolveUpdateSize; + filter += ConvolveUpdateSize; } input = mComplexData.get(); for(size_t s{0};s < curseg;++s) { - for(size_t i{0};i < m;++i,++input,++filter) - mFftBuffer[i] += *input * *filter; + pffft_zconvolve_accumulate(mFft.get(), input, filter, mFftBuffer.data(), fftscale); + input += ConvolveUpdateSize; + filter += ConvolveUpdateSize; } - /* Reconstruct the mirrored/negative frequencies to do a proper - * inverse FFT. - */ - for(size_t i{m};i < ConvolveUpdateSize;++i) - mFftBuffer[i] = std::conj(mFftBuffer[ConvolveUpdateSize-i]); - /* Apply iFFT to get the 256 (really 255) samples for output. The * 128 output samples are combined with the last output's 127 * second-half samples (and this output's second half is * subsequently saved for next time). */ - inverse_fft(al::span{mFftBuffer}); + pffft_transform(mFft.get(), mFftBuffer.data(), mFftBuffer.data(), + mFftWorkBuffer.data(), PFFFT_BACKWARD); - /* The iFFT'd response is scaled up by the number of bins, so apply - * the inverse to normalize the output. - */ for(size_t i{0};i < ConvolveUpdateSamples;++i) - mOutput[c][i] = - (mFftBuffer[i].real()+mOutput[c][ConvolveUpdateSamples+i]) * - (1.0f/float{ConvolveUpdateSize}); + mOutput[c][i] = mFftBuffer[i] + mOutput[c][ConvolveUpdateSamples+i]; for(size_t i{0};i < ConvolveUpdateSamples;++i) - mOutput[c][ConvolveUpdateSamples+i] = mFftBuffer[ConvolveUpdateSamples+i].real(); + mOutput[c][ConvolveUpdateSamples+i] = mFftBuffer[ConvolveUpdateSamples+i]; } /* Shift the input history. */ |