diff options
author | Sven Gothel <[email protected]> | 2023-11-28 12:51:46 +0100 |
---|---|---|
committer | Sven Gothel <[email protected]> | 2023-11-28 12:51:46 +0100 |
commit | 1aaf4f070011490bcece50394b9b32dfa593fd9e (patch) | |
tree | 17d68284e401a35eea3d3a574d986d446a60763a /alc/effects/convolution.cpp | |
parent | 6e7cee4fa9a8af03f28ca26cd89f8357390dfc90 (diff) | |
parent | 571b546f35eead77ce109f8d4dd6c3de3199d573 (diff) |
Merge remote-tracking branch 'upstream/master'
Diffstat (limited to 'alc/effects/convolution.cpp')
-rw-r--r-- | alc/effects/convolution.cpp | 388 |
1 files changed, 229 insertions, 159 deletions
diff --git a/alc/effects/convolution.cpp b/alc/effects/convolution.cpp index 7f36c415..517e6b08 100644 --- a/alc/effects/convolution.cpp +++ b/alc/effects/convolution.cpp @@ -17,7 +17,6 @@ #include <arm_neon.h> #endif -#include "albyte.h" #include "alcomplex.h" #include "almalloc.h" #include "alnumbers.h" @@ -35,44 +34,49 @@ #include "core/fmt_traits.h" #include "core/mixer.h" #include "intrusive_ptr.h" +#include "pffft.h" #include "polyphase_resampler.h" #include "vector.h" namespace { -/* Convolution reverb is implemented using a segmented overlap-add method. The - * impulse response is broken up into multiple segments of 128 samples, and - * 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. +/* Convolution is implemented using a segmented overlap-add method. The impulse + * response is split into multiple segments of 128 samples, and 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 0- and half-frequency 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 - * input segments is maintained, equal to the length of the impulse response. + * Input samples are similarly broken up into 128-sample segments, with a 256- + * sample FFT applied to each new incoming segment to get its 129 bins. A + * history of FFT'd input segments is maintained, equal to the number of + * impulse response segments. * - * To apply the reverberation, each impulse response segment is convolved with + * To apply the convolution, 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 * first half is the 128-sample output, and the second half is a 128-sample * (really, 127) delayed extension, which gets added to the output next time. - * Convolving two time-domain responses of lengths N and M results in a time- - * domain signal of length N+M-1, and this holds true regardless of the - * convolution being applied in the frequency domain, so these "overflow" - * samples need to be accounted for. + * Convolving two time-domain responses of length N results in a time-domain + * signal of length N*2 - 1, and this holds true regardless of the convolution + * being applied in the frequency domain, so these "overflow" samples need to + * be accounted for. * - * To avoid a delay with gathering enough input samples to apply an FFT with, - * the first segment is applied directly in the time-domain as the samples come - * in. Once enough have been retrieved, the FFT is applied on the input and - * it's paired with the remaining (FFT'd) filter segments for processing. + * To avoid a delay with gathering enough input samples for the FFT, the first + * segment is applied directly in the time-domain as the samples come in. Once + * enough have been retrieved, the FFT is applied on the input and it's paired + * with the remaining (FFT'd) filter segments for processing. */ -void LoadSamples(float *RESTRICT dst, const al::byte *src, const size_t srcstep, FmtType srctype, +void LoadSamples(float *RESTRICT dst, const std::byte *src, const size_t srcstep, FmtType srctype, const size_t samples) noexcept { #define HANDLE_FMT(T) case T: al::LoadSampleArray<T>(dst, src, srcstep, samples); break @@ -80,6 +84,7 @@ void LoadSamples(float *RESTRICT dst, const al::byte *src, const size_t srcstep, { HANDLE_FMT(FmtUByte); HANDLE_FMT(FmtShort); + HANDLE_FMT(FmtInt); HANDLE_FMT(FmtFloat); HANDLE_FMT(FmtDouble); HANDLE_FMT(FmtMulaw); @@ -94,40 +99,43 @@ void LoadSamples(float *RESTRICT dst, const al::byte *src, const size_t srcstep, } -inline auto& GetAmbiScales(AmbiScaling scaletype) noexcept +constexpr auto GetAmbiScales(AmbiScaling scaletype) noexcept { switch(scaletype) { - case AmbiScaling::FuMa: return AmbiScale::FromFuMa(); - case AmbiScaling::SN3D: return AmbiScale::FromSN3D(); - case AmbiScaling::UHJ: return AmbiScale::FromUHJ(); + case AmbiScaling::FuMa: return al::span{AmbiScale::FromFuMa}; + case AmbiScaling::SN3D: return al::span{AmbiScale::FromSN3D}; + case AmbiScaling::UHJ: return al::span{AmbiScale::FromUHJ}; case AmbiScaling::N3D: break; } - return AmbiScale::FromN3D(); + return al::span{AmbiScale::FromN3D}; } -inline auto& GetAmbiLayout(AmbiLayout layouttype) noexcept +constexpr auto GetAmbiLayout(AmbiLayout layouttype) noexcept { - if(layouttype == AmbiLayout::FuMa) return AmbiIndex::FromFuMa(); - return AmbiIndex::FromACN(); + if(layouttype == AmbiLayout::FuMa) return al::span{AmbiIndex::FromFuMa}; + return al::span{AmbiIndex::FromACN}; } -inline auto& GetAmbi2DLayout(AmbiLayout layouttype) noexcept +constexpr auto GetAmbi2DLayout(AmbiLayout layouttype) noexcept { - if(layouttype == AmbiLayout::FuMa) return AmbiIndex::FromFuMa2D(); - return AmbiIndex::FromACN2D(); + if(layouttype == AmbiLayout::FuMa) return al::span{AmbiIndex::FromFuMa2D}; + return al::span{AmbiIndex::FromACN2D}; } -struct ChanMap { +constexpr float sin30{0.5f}; +constexpr float cos30{0.866025403785f}; +constexpr float sin45{al::numbers::sqrt2_v<float>*0.5f}; +constexpr float cos45{al::numbers::sqrt2_v<float>*0.5f}; +constexpr float sin110{ 0.939692620786f}; +constexpr float cos110{-0.342020143326f}; + +struct ChanPosMap { Channel channel; - float angle; - float elevation; + std::array<float,3> pos; }; -constexpr float Deg2Rad(float x) noexcept -{ return static_cast<float>(al::numbers::pi / 180.0 * x); } - using complex_f = std::complex<float>; @@ -181,6 +189,13 @@ void apply_fir(al::span<float> dst, const float *RESTRICT src, const float *REST #endif } + +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{}; @@ -188,11 +203,13 @@ struct ConvolutionState final : public EffectState { uint mAmbiOrder{}; size_t mFifoPos{0}; - std::array<float,ConvolveUpdateSamples*2> mInput{}; + alignas(16) std::array<float,ConvolveUpdateSamples*2> mInput{}; 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}; @@ -204,9 +221,8 @@ struct ConvolutionState final : public EffectState { float Current[MAX_OUTPUT_CHANNELS]{}; float Target[MAX_OUTPUT_CHANNELS]{}; }; - using ChannelDataArray = al::FlexArray<ChannelData>; - std::unique_ptr<ChannelDataArray> mChans; - std::unique_ptr<complex_f[]> mComplexData; + std::vector<ChannelData> mChans; + al::vector<float,16> mComplexData; ConvolutionState() = default; @@ -229,7 +245,7 @@ struct ConvolutionState final : public EffectState { void ConvolutionState::NormalMix(const al::span<FloatBufferLine> samplesOut, const size_t samplesToDo) { - for(auto &chan : *mChans) + for(auto &chan : mChans) MixSamples({chan.mBuffer.data(), samplesToDo}, samplesOut, chan.Current, chan.Target, samplesToDo, 0); } @@ -237,7 +253,7 @@ void ConvolutionState::NormalMix(const al::span<FloatBufferLine> samplesOut, void ConvolutionState::UpsampleMix(const al::span<FloatBufferLine> samplesOut, const size_t samplesToDo) { - for(auto &chan : *mChans) + for(auto &chan : mChans) { const al::span<float> src{chan.mBuffer.data(), samplesToDo}; chan.mFilter.processScale(src, chan.mHfScale, chan.mLfScale); @@ -251,19 +267,23 @@ 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; - mChans = nullptr; - mComplexData = nullptr; + decltype(mChans){}.swap(mChans); + decltype(mComplexData){}.swap(mComplexData); /* An empty buffer doesn't need a convolution filter. */ if(!buffer || buffer->mSampleLen < 1) return; @@ -273,12 +293,11 @@ 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); - mChans = ChannelDataArray::Create(numChannels); + mChans.resize(numChannels); /* The impulse response needs to have the same sample rate as the input and * output. The bsinc24 resampler is decent, but there is high-frequency @@ -293,7 +312,7 @@ void ConvolutionState::deviceUpdate(const DeviceBase *device, const BufferStorag buffer->mSampleRate); const BandSplitter splitter{device->mXOverFreq / static_cast<float>(device->Frequency)}; - for(auto &e : *mChans) + for(auto &e : mChans) e.mFilter = splitter; mFilter.resize(numChannels, {}); @@ -307,9 +326,8 @@ 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.resize(complex_length, 0.0f); /* Load the samples from the buffer. */ const size_t srclinelength{RoundUp(buffer->mSampleLen+DecoderPadding, 16)}; @@ -330,7 +348,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.data() + mNumConvolveSegs*ConvolveUpdateSize; for(size_t c{0};c < numChannels;++c) { /* Resample to match the device. */ @@ -351,71 +372,85 @@ 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}); - forward_fft(al::as_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. Also scale the FFT by its length so the + * iFFT'd output will be normalized. + */ + static constexpr float fftscale{1.0f / float{ConvolveUpdateSize}}; + for(size_t i{0};i < ConvolveUpdateSamples;++i) + { + ffttmp[i*2 ] = static_cast<float>(fftbuffer[i].real()) * fftscale; + ffttmp[i*2 + 1] = static_cast<float>((i == 0) ? + fftbuffer[ConvolveUpdateSamples].real() : fftbuffer[i].imag()) * fftscale; + } + /* 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; } } } void ConvolutionState::update(const ContextBase *context, const EffectSlot *slot, - const EffectProps* /*props*/, const EffectTarget target) + const EffectProps *props, const EffectTarget target) { - /* NOTE: Stereo and Rear are slightly different from normal mixing (as - * defined in alu.cpp). These are 45 degrees from center, rather than the - * 30 degrees used there. - * - * TODO: LFE is not mixed to output. This will require each buffer channel + /* TODO: LFE is not mixed to output. This will require each buffer channel * to have its own output target since the main mixing buffer won't have an * LFE channel (due to being B-Format). */ - static constexpr ChanMap MonoMap[1]{ - { FrontCenter, 0.0f, 0.0f } + static constexpr ChanPosMap MonoMap[1]{ + { FrontCenter, std::array{0.0f, 0.0f, -1.0f} } }, StereoMap[2]{ - { FrontLeft, Deg2Rad(-45.0f), Deg2Rad(0.0f) }, - { FrontRight, Deg2Rad( 45.0f), Deg2Rad(0.0f) } + { FrontLeft, std::array{-sin30, 0.0f, -cos30} }, + { FrontRight, std::array{ sin30, 0.0f, -cos30} }, }, RearMap[2]{ - { BackLeft, Deg2Rad(-135.0f), Deg2Rad(0.0f) }, - { BackRight, Deg2Rad( 135.0f), Deg2Rad(0.0f) } + { BackLeft, std::array{-sin30, 0.0f, cos30} }, + { BackRight, std::array{ sin30, 0.0f, cos30} }, }, QuadMap[4]{ - { FrontLeft, Deg2Rad( -45.0f), Deg2Rad(0.0f) }, - { FrontRight, Deg2Rad( 45.0f), Deg2Rad(0.0f) }, - { BackLeft, Deg2Rad(-135.0f), Deg2Rad(0.0f) }, - { BackRight, Deg2Rad( 135.0f), Deg2Rad(0.0f) } + { FrontLeft, std::array{-sin45, 0.0f, -cos45} }, + { FrontRight, std::array{ sin45, 0.0f, -cos45} }, + { BackLeft, std::array{-sin45, 0.0f, cos45} }, + { BackRight, std::array{ sin45, 0.0f, cos45} }, }, X51Map[6]{ - { FrontLeft, Deg2Rad( -30.0f), Deg2Rad(0.0f) }, - { FrontRight, Deg2Rad( 30.0f), Deg2Rad(0.0f) }, - { FrontCenter, Deg2Rad( 0.0f), Deg2Rad(0.0f) }, - { LFE, 0.0f, 0.0f }, - { SideLeft, Deg2Rad(-110.0f), Deg2Rad(0.0f) }, - { SideRight, Deg2Rad( 110.0f), Deg2Rad(0.0f) } + { FrontLeft, std::array{-sin30, 0.0f, -cos30} }, + { FrontRight, std::array{ sin30, 0.0f, -cos30} }, + { FrontCenter, std::array{ 0.0f, 0.0f, -1.0f} }, + { LFE, {} }, + { SideLeft, std::array{-sin110, 0.0f, -cos110} }, + { SideRight, std::array{ sin110, 0.0f, -cos110} }, }, X61Map[7]{ - { FrontLeft, Deg2Rad(-30.0f), Deg2Rad(0.0f) }, - { FrontRight, Deg2Rad( 30.0f), Deg2Rad(0.0f) }, - { FrontCenter, Deg2Rad( 0.0f), Deg2Rad(0.0f) }, - { LFE, 0.0f, 0.0f }, - { BackCenter, Deg2Rad(180.0f), Deg2Rad(0.0f) }, - { SideLeft, Deg2Rad(-90.0f), Deg2Rad(0.0f) }, - { SideRight, Deg2Rad( 90.0f), Deg2Rad(0.0f) } + { FrontLeft, std::array{-sin30, 0.0f, -cos30} }, + { FrontRight, std::array{ sin30, 0.0f, -cos30} }, + { FrontCenter, std::array{ 0.0f, 0.0f, -1.0f} }, + { LFE, {} }, + { BackCenter, std::array{ 0.0f, 0.0f, 1.0f} }, + { SideLeft, std::array{-1.0f, 0.0f, 0.0f} }, + { SideRight, std::array{ 1.0f, 0.0f, 0.0f} }, }, X71Map[8]{ - { FrontLeft, Deg2Rad( -30.0f), Deg2Rad(0.0f) }, - { FrontRight, Deg2Rad( 30.0f), Deg2Rad(0.0f) }, - { FrontCenter, Deg2Rad( 0.0f), Deg2Rad(0.0f) }, - { LFE, 0.0f, 0.0f }, - { BackLeft, Deg2Rad(-150.0f), Deg2Rad(0.0f) }, - { BackRight, Deg2Rad( 150.0f), Deg2Rad(0.0f) }, - { SideLeft, Deg2Rad( -90.0f), Deg2Rad(0.0f) }, - { SideRight, Deg2Rad( 90.0f), Deg2Rad(0.0f) } + { FrontLeft, std::array{-sin30, 0.0f, -cos30} }, + { FrontRight, std::array{ sin30, 0.0f, -cos30} }, + { FrontCenter, std::array{ 0.0f, 0.0f, -1.0f} }, + { LFE, {} }, + { BackLeft, std::array{-sin30, 0.0f, cos30} }, + { BackRight, std::array{ sin30, 0.0f, cos30} }, + { SideLeft, std::array{ -1.0f, 0.0f, 0.0f} }, + { SideRight, std::array{ 1.0f, 0.0f, 0.0f} }, }; if(mNumConvolveSegs < 1) UNLIKELY @@ -423,7 +458,7 @@ void ConvolutionState::update(const ContextBase *context, const EffectSlot *slot mMix = &ConvolutionState::NormalMix; - for(auto &chan : *mChans) + for(auto &chan : mChans) std::fill(std::begin(chan.Target), std::end(chan.Target), 0.0f); const float gain{slot->Gain}; if(IsAmbisonic(mChannels)) @@ -432,46 +467,66 @@ void ConvolutionState::update(const ContextBase *context, const EffectSlot *slot if(mChannels == FmtUHJ2 && !device->mUhjEncoder) { mMix = &ConvolutionState::UpsampleMix; - (*mChans)[0].mHfScale = 1.0f; - (*mChans)[0].mLfScale = DecoderBase::sWLFScale; - (*mChans)[1].mHfScale = 1.0f; - (*mChans)[1].mLfScale = DecoderBase::sXYLFScale; - (*mChans)[2].mHfScale = 1.0f; - (*mChans)[2].mLfScale = DecoderBase::sXYLFScale; + mChans[0].mHfScale = 1.0f; + mChans[0].mLfScale = DecoderBase::sWLFScale; + mChans[1].mHfScale = 1.0f; + mChans[1].mLfScale = DecoderBase::sXYLFScale; + mChans[2].mHfScale = 1.0f; + mChans[2].mLfScale = DecoderBase::sXYLFScale; } else if(device->mAmbiOrder > mAmbiOrder) { mMix = &ConvolutionState::UpsampleMix; const auto scales = AmbiScale::GetHFOrderScales(mAmbiOrder, device->mAmbiOrder, device->m2DMixing); - (*mChans)[0].mHfScale = scales[0]; - (*mChans)[0].mLfScale = 1.0f; - for(size_t i{1};i < mChans->size();++i) + mChans[0].mHfScale = scales[0]; + mChans[0].mLfScale = 1.0f; + for(size_t i{1};i < mChans.size();++i) { - (*mChans)[i].mHfScale = scales[1]; - (*mChans)[i].mLfScale = 1.0f; + mChans[i].mHfScale = scales[1]; + mChans[i].mLfScale = 1.0f; } } mOutTarget = target.Main->Buffer; - auto&& scales = GetAmbiScales(mAmbiScaling); + alu::Vector N{props->Convolution.OrientAt[0], props->Convolution.OrientAt[1], + props->Convolution.OrientAt[2], 0.0f}; + N.normalize(); + alu::Vector V{props->Convolution.OrientUp[0], props->Convolution.OrientUp[1], + props->Convolution.OrientUp[2], 0.0f}; + V.normalize(); + /* Build and normalize right-vector */ + alu::Vector U{N.cross_product(V)}; + U.normalize(); + + const float mixmatrix[4][4]{ + {1.0f, 0.0f, 0.0f, 0.0f}, + {0.0f, U[0], -U[1], U[2]}, + {0.0f, -V[0], V[1], -V[2]}, + {0.0f, -N[0], N[1], -N[2]}, + }; + + const auto scales = GetAmbiScales(mAmbiScaling); const uint8_t *index_map{Is2DAmbisonic(mChannels) ? GetAmbi2DLayout(mAmbiLayout).data() : GetAmbiLayout(mAmbiLayout).data()}; std::array<float,MaxAmbiChannels> coeffs{}; - for(size_t c{0u};c < mChans->size();++c) + for(size_t c{0u};c < mChans.size();++c) { const size_t acn{index_map[c]}; - coeffs[acn] = scales[acn]; - ComputePanGains(target.Main, coeffs.data(), gain, (*mChans)[c].Target); - coeffs[acn] = 0.0f; + const float scale{scales[acn]}; + + for(size_t x{0};x < 4;++x) + coeffs[x] = mixmatrix[acn][x] * scale; + + ComputePanGains(target.Main, coeffs, gain, mChans[c].Target); } } else { DeviceBase *device{context->mDevice}; - al::span<const ChanMap> chanmap{}; + al::span<const ChanPosMap> chanmap{}; switch(mChannels) { case FmtMono: chanmap = MonoMap; break; @@ -493,28 +548,55 @@ void ConvolutionState::update(const ContextBase *context, const EffectSlot *slot mOutTarget = target.Main->Buffer; if(device->mRenderMode == RenderMode::Pairwise) { - auto ScaleAzimuthFront = [](float azimuth, float scale) -> float + /* Scales the azimuth of the given vector by 3 if it's in front. + * Effectively scales +/-30 degrees to +/-90 degrees, leaving > +90 + * and < -90 alone. + */ + auto ScaleAzimuthFront = [](std::array<float,3> pos) -> std::array<float,3> { - constexpr float half_pi{al::numbers::pi_v<float>*0.5f}; - const float abs_azi{std::fabs(azimuth)}; - if(!(abs_azi >= half_pi)) - return std::copysign(minf(abs_azi*scale, half_pi), azimuth); - return azimuth; + if(pos[2] < 0.0f) + { + /* Normalize the length of the x,z components for a 2D + * vector of the azimuth angle. Negate Z since {0,0,-1} is + * angle 0. + */ + const float len2d{std::sqrt(pos[0]*pos[0] + pos[2]*pos[2])}; + float x{pos[0] / len2d}; + float z{-pos[2] / len2d}; + + /* Z > cos(pi/6) = -30 < azimuth < 30 degrees. */ + if(z > cos30) + { + /* Triple the angle represented by x,z. */ + x = x*3.0f - x*x*x*4.0f; + z = z*z*z*4.0f - z*3.0f; + + /* Scale the vector back to fit in 3D. */ + pos[0] = x * len2d; + pos[2] = -z * len2d; + } + else + { + /* If azimuth >= 30 degrees, clamp to 90 degrees. */ + pos[0] = std::copysign(len2d, pos[0]); + pos[2] = 0.0f; + } + } + return pos; }; for(size_t i{0};i < chanmap.size();++i) { if(chanmap[i].channel == LFE) continue; - const auto coeffs = CalcAngleCoeffs(ScaleAzimuthFront(chanmap[i].angle, 2.0f), - chanmap[i].elevation, 0.0f); - ComputePanGains(target.Main, coeffs.data(), gain, (*mChans)[i].Target); + const auto coeffs = CalcDirectionCoeffs(ScaleAzimuthFront(chanmap[i].pos), 0.0f); + ComputePanGains(target.Main, coeffs, gain, mChans[i].Target); } } else for(size_t i{0};i < chanmap.size();++i) { if(chanmap[i].channel == LFE) continue; - const auto coeffs = CalcAngleCoeffs(chanmap[i].angle, chanmap[i].elevation, 0.0f); - ComputePanGains(target.Main, coeffs.data(), gain, (*mChans)[i].Target); + const auto coeffs = CalcDirectionCoeffs(chanmap[i].pos, 0.0f); + ComputePanGains(target.Main, coeffs, gain, mChans[i].Target); } } } @@ -525,9 +607,7 @@ 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; for(size_t base{0u};base < samplesToDo;) { @@ -539,9 +619,9 @@ void ConvolutionState::process(const size_t samplesToDo, /* Apply the FIR for the newly retrieved input samples, and combine it * with the inverse FFT'd output samples. */ - for(size_t c{0};c < chans.size();++c) + for(size_t c{0};c < mChans.size();++c) { - auto buf_iter = chans[c].mBuffer.begin() + base; + auto buf_iter = mChans[c].mBuffer.begin() + base; apply_fir({buf_iter, todo}, mInput.data()+1 + mFifoPos, mFilter[c].data()); auto fifo_iter = mOutput[c].begin() + mFifoPos; @@ -557,59 +637,49 @@ void ConvolutionState::process(const size_t samplesToDo, /* Move the newest input to the front for the next iteration's history. */ std::copy(mInput.cbegin()+ConvolveUpdateSamples, mInput.cend(), mInput.begin()); + std::fill(mInput.begin()+ConvolveUpdateSamples, mInput.end(), 0.0f); - /* Calculate the frequency domain response and add the relevant + /* Calculate the frequency-domain response and add the relevant * 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::as_span(mFftBuffer)); + pffft_transform(mFft.get(), mInput.data(), mComplexData.data() + curseg*ConvolveUpdateSize, + mFftWorkBuffer.data(), PFFFT_FORWARD); - std::copy_n(mFftBuffer.cbegin(), m, &mComplexData[curseg*m]); - - const complex_f *RESTRICT filter{mComplexData.get() + mNumConvolveSegs*m}; - for(size_t c{0};c < chans.size();++c) + const float *filter{mComplexData.data() + mNumConvolveSegs*ConvolveUpdateSize}; + for(size_t c{0};c < mChans.size();++c) { - std::fill_n(mFftBuffer.begin(), m, complex_f{}); - /* 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 *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()); + input += ConvolveUpdateSize; + filter += ConvolveUpdateSize; } - input = mComplexData.get(); + input = mComplexData.data(); 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()); + 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::as_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. - */ + /* The filter was attenuated, so the response is already scaled. */ 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. */ |