alc/effects/convolution.cpp

   1
   2 #include "config.h"
   3
   4 #include <algorithm>
   5 #include <array>
   6 #include <complex>
   7 #include <cstddef>
   8 #include <functional>
   9 #include <iterator>
  10 #include <memory>
  11 #include <stdint.h>
  12 #include <utility>
  13
  14 #ifdef HAVE_SSE_INTRINSICS
  15 #include <xmmintrin.h>
  16 #elif defined(HAVE_NEON)
  17 #include <arm_neon.h>
  18 #endif
  19
  20 #include "albyte.h"
  21 #include "alcomplex.h"
  22 #include "almalloc.h"
  23 #include "alnumbers.h"
  24 #include "alnumeric.h"
  25 #include "alspan.h"
  26 #include "base.h"
  27 #include "core/ambidefs.h"
  28 #include "core/bufferline.h"
  29 #include "core/buffer_storage.h"
  30 #include "core/context.h"
  31 #include "core/devformat.h"
  32 #include "core/device.h"
  33 #include "core/effectslot.h"
  34 #include "core/filters/splitter.h"
  35 #include "core/fmt_traits.h"
  36 #include "core/mixer.h"
  37 #include "intrusive_ptr.h"
  38 #include "polyphase_resampler.h"
  39 #include "vector.h"
  40
  41
  42 namespace {
  43
  44 /* Convolution reverb is implemented using a segmented overlap-add method. The
  45  * impulse response is broken up into multiple segments of 128 samples, and
  46  * each segment has an FFT applied with a 256-sample buffer (the latter half
  47  * left silent) to get its frequency-domain response. The resulting response
  48  * has its positive/non-mirrored frequencies saved (129 bins) in each segment.
  49  *
  50  * Input samples are similarly broken up into 128-sample segments, with an FFT
  51  * applied to each new incoming segment to get its 129 bins. A history of FFT'd
  52  * input segments is maintained, equal to the length of the impulse response.
  53  *
  54  * To apply the reverberation, each impulse response segment is convolved with
  55  * its paired input segment (using complex multiplies, far cheaper than FIRs),
  56  * accumulating into a 256-bin FFT buffer. The input history is then shifted to
  57  * align with later impulse response segments for next time.
  58  *
  59  * An inverse FFT is then applied to the accumulated FFT buffer to get a 256-
  60  * sample time-domain response for output, which is split in two halves. The
  61  * first half is the 128-sample output, and the second half is a 128-sample
  62  * (really, 127) delayed extension, which gets added to the output next time.
  63  * Convolving two time-domain responses of lengths N and M results in a time-
  64  * domain signal of length N+M-1, and this holds true regardless of the
  65  * convolution being applied in the frequency domain, so these "overflow"
  66  * samples need to be accounted for.
  67  *
  68  * To avoid a delay with gathering enough input samples to apply an FFT with,
  69  * the first segment is applied directly in the time-domain as the samples come
  70  * in. Once enough have been retrieved, the FFT is applied on the input and
  71  * it's paired with the remaining (FFT'd) filter segments for processing.
  72  */
  73
  74
  75 void LoadSamples(float *RESTRICT dst, const al::byte *src, const size_t srcstep, FmtType srctype,
  76     const size_t samples) noexcept
  77 {
  78 #define HANDLE_FMT(T)  case T: al::LoadSampleArray<T>(dst, src, srcstep, samples); break
  79     switch(srctype)
  80     {
  81     HANDLE_FMT(FmtUByte);
  82     HANDLE_FMT(FmtShort);
  83     HANDLE_FMT(FmtFloat);
  84     HANDLE_FMT(FmtDouble);
  85     HANDLE_FMT(FmtMulaw);
  86     HANDLE_FMT(FmtAlaw);
  87     /* FIXME: Handle ADPCM decoding here. */
  88     case FmtIMA4:
  89     case FmtMSADPCM:
  90         std::fill_n(dst, samples, 0.0f);
  91         break;
  92     }
  93 #undef HANDLE_FMT
  94 }
  95
  96
  97 inline auto& GetAmbiScales(AmbiScaling scaletype) noexcept
  98 {
  99     switch(scaletype)
 100     {
 101     case AmbiScaling::FuMa: return AmbiScale::FromFuMa();
 102     case AmbiScaling::SN3D: return AmbiScale::FromSN3D();
 103     case AmbiScaling::UHJ: return AmbiScale::FromUHJ();
 104     case AmbiScaling::N3D: break;
 105     }
 106     return AmbiScale::FromN3D();
 107 }
 108
 109 inline auto& GetAmbiLayout(AmbiLayout layouttype) noexcept
 110 {
 111     if(layouttype == AmbiLayout::FuMa) return AmbiIndex::FromFuMa();
 112     return AmbiIndex::FromACN();
 113 }
 114
 115 inline auto& GetAmbi2DLayout(AmbiLayout layouttype) noexcept
 116 {
 117     if(layouttype == AmbiLayout::FuMa) return AmbiIndex::FromFuMa2D();
 118     return AmbiIndex::FromACN2D();
 119 }
 120
 121
 122 struct ChanMap {
 123     Channel channel;
 124     float angle;
 125     float elevation;
 126 };
 127
 128 constexpr float Deg2Rad(float x) noexcept
 129 { return static_cast<float>(al::numbers::pi / 180.0 * x); }
 130
 131
 132 using complex_f = std::complex<float>;
 133
 134 constexpr size_t ConvolveUpdateSize{256};
 135 constexpr size_t ConvolveUpdateSamples{ConvolveUpdateSize / 2};
 136
 137
 138 void apply_fir(al::span<float> dst, const float *RESTRICT src, const float *RESTRICT filter)
 139 {
 140 #ifdef HAVE_SSE_INTRINSICS
 141     for(float &output : dst)
 142     {
 143         __m128 r4{_mm_setzero_ps()};
 144         for(size_t j{0};j < ConvolveUpdateSamples;j+=4)
 145         {
 146             const __m128 coeffs{_mm_load_ps(&filter[j])};
 147             const __m128 s{_mm_loadu_ps(&src[j])};
 148
 149             r4 = _mm_add_ps(r4, _mm_mul_ps(s, coeffs));
 150         }
 151         r4 = _mm_add_ps(r4, _mm_shuffle_ps(r4, r4, _MM_SHUFFLE(0, 1, 2, 3)));
 152         r4 = _mm_add_ps(r4, _mm_movehl_ps(r4, r4));
 153         output = _mm_cvtss_f32(r4);
 154
 155         ++src;
 156     }
 157
 158 #elif defined(HAVE_NEON)
 159
 160     for(float &output : dst)
 161     {
 162         float32x4_t r4{vdupq_n_f32(0.0f)};
 163         for(size_t j{0};j < ConvolveUpdateSamples;j+=4)
 164             r4 = vmlaq_f32(r4, vld1q_f32(&src[j]), vld1q_f32(&filter[j]));
 165         r4 = vaddq_f32(r4, vrev64q_f32(r4));
 166         output = vget_lane_f32(vadd_f32(vget_low_f32(r4), vget_high_f32(r4)), 0);
 167
 168         ++src;
 169     }
 170
 171 #else
 172
 173     for(float &output : dst)
 174     {
 175         float ret{0.0f};
 176         for(size_t j{0};j < ConvolveUpdateSamples;++j)
 177             ret += src[j] * filter[j];
 178         output = ret;
 179         ++src;
 180     }
 181 #endif
 182 }
 183
 184 struct ConvolutionState final : public EffectState {
 185     FmtChannels mChannels{};
 186     AmbiLayout mAmbiLayout{};
 187     AmbiScaling mAmbiScaling{};
 188     uint mAmbiOrder{};
 189
 190     size_t mFifoPos{0};
 191     std::array<float,ConvolveUpdateSamples*2> mInput{};
 192     al::vector<std::array<float,ConvolveUpdateSamples>,16> mFilter;
 193     al::vector<std::array<float,ConvolveUpdateSamples*2>,16> mOutput;
 194
 195     alignas(16) std::array<complex_f,ConvolveUpdateSize> mFftBuffer{};
 196
 197     size_t mCurrentSegment{0};
 198     size_t mNumConvolveSegs{0};
 199
 200     struct ChannelData {
 201         alignas(16) FloatBufferLine mBuffer{};
 202         float mHfScale{}, mLfScale{};
 203         BandSplitter mFilter{};
 204         float Current[MAX_OUTPUT_CHANNELS]{};
 205         float Target[MAX_OUTPUT_CHANNELS]{};
 206     };
 207     using ChannelDataArray = al::FlexArray<ChannelData>;
 208     std::unique_ptr<ChannelDataArray> mChans;
 209     std::unique_ptr<complex_f[]> mComplexData;
 210
 211
 212     ConvolutionState() = default;
 213     ~ConvolutionState() override = default;
 214
 215     void NormalMix(const al::span<FloatBufferLine> samplesOut, const size_t samplesToDo);
 216     void UpsampleMix(const al::span<FloatBufferLine> samplesOut, const size_t samplesToDo);
 217     void (ConvolutionState::*mMix)(const al::span<FloatBufferLine>,const size_t)
 218     {&ConvolutionState::NormalMix};
 219
 220     void deviceUpdate(const DeviceBase *device, const BufferStorage *buffer) override;
 221     void update(const ContextBase *context, const EffectSlot *slot, const EffectProps *props,
 222         const EffectTarget target) override;
 223     void process(const size_t samplesToDo, const al::span<const FloatBufferLine> samplesIn,
 224         const al::span<FloatBufferLine> samplesOut) override;
 225
 226     DEF_NEWDEL(ConvolutionState)
 227 };
 228
 229 void ConvolutionState::NormalMix(const al::span<FloatBufferLine> samplesOut,
 230     const size_t samplesToDo)
 231 {
 232     for(auto &chan : *mChans)
 233         MixSamples({chan.mBuffer.data(), samplesToDo}, samplesOut, chan.Current, chan.Target,
 234             samplesToDo, 0);
 235 }
 236
 237 void ConvolutionState::UpsampleMix(const al::span<FloatBufferLine> samplesOut,
 238     const size_t samplesToDo)
 239 {
 240     for(auto &chan : *mChans)
 241     {
 242         const al::span<float> src{chan.mBuffer.data(), samplesToDo};
 243         chan.mFilter.processScale(src, chan.mHfScale, chan.mLfScale);
 244         MixSamples(src, samplesOut, chan.Current, chan.Target, samplesToDo, 0);
 245     }
 246 }
 247
 248
 249 void ConvolutionState::deviceUpdate(const DeviceBase *device, const BufferStorage *buffer)
 250 {
 251     using UhjDecoderType = UhjDecoder<512>;
 252     static constexpr auto DecoderPadding = UhjDecoderType::sInputPadding;
 253
 254     constexpr uint MaxConvolveAmbiOrder{1u};
 255
 256     mFifoPos = 0;
 257     mInput.fill(0.0f);
 258     decltype(mFilter){}.swap(mFilter);
 259     decltype(mOutput){}.swap(mOutput);
 260     mFftBuffer.fill(complex_f{});
 261
 262     mCurrentSegment = 0;
 263     mNumConvolveSegs = 0;
 264
 265     mChans = nullptr;
 266     mComplexData = nullptr;
 267
 268     /* An empty buffer doesn't need a convolution filter. */
 269     if(!buffer || buffer->mSampleLen < 1) return;
 270
 271     mChannels = buffer->mChannels;
 272     mAmbiLayout = IsUHJ(mChannels) ? AmbiLayout::FuMa : buffer->mAmbiLayout;
 273     mAmbiScaling = IsUHJ(mChannels) ? AmbiScaling::UHJ : buffer->mAmbiScaling;
 274     mAmbiOrder = minu(buffer->mAmbiOrder, MaxConvolveAmbiOrder);
 275
 276     constexpr size_t m{ConvolveUpdateSize/2 + 1};
 277     const auto bytesPerSample = BytesFromFmt(buffer->mType);
 278     const auto realChannels = buffer->channelsFromFmt();
 279     const auto numChannels = (mChannels == FmtUHJ2) ? 3u : ChannelsFromFmt(mChannels, mAmbiOrder);
 280
 281     mChans = ChannelDataArray::Create(numChannels);
 282
 283     /* The impulse response needs to have the same sample rate as the input and
 284      * output. The bsinc24 resampler is decent, but there is high-frequency
 285      * attenuation that some people may be able to pick up on. Since this is
 286      * called very infrequently, go ahead and use the polyphase resampler.
 287      */
 288     PPhaseResampler resampler;
 289     if(device->Frequency != buffer->mSampleRate)
 290         resampler.init(buffer->mSampleRate, device->Frequency);
 291     const auto resampledCount = static_cast<uint>(
 292         (uint64_t{buffer->mSampleLen}*device->Frequency+(buffer->mSampleRate-1)) /
 293         buffer->mSampleRate);
 294
 295     const BandSplitter splitter{device->mXOverFreq / static_cast<float>(device->Frequency)};
 296     for(auto &e : *mChans)
 297         e.mFilter = splitter;
 298
 299     mFilter.resize(numChannels, {});
 300     mOutput.resize(numChannels, {});
 301
 302     /* Calculate the number of segments needed to hold the impulse response and
 303      * the input history (rounded up), and allocate them. Exclude one segment
 304      * which gets applied as a time-domain FIR filter. Make sure at least one
 305      * segment is allocated to simplify handling.
 306      */
 307     mNumConvolveSegs = (resampledCount+(ConvolveUpdateSamples-1)) / ConvolveUpdateSamples;
 308     mNumConvolveSegs = maxz(mNumConvolveSegs, 2) - 1;
 309
 310     const size_t complex_length{mNumConvolveSegs * m * (numChannels+1)};
 311     mComplexData = std::make_unique<complex_f[]>(complex_length);
 312     std::fill_n(mComplexData.get(), complex_length, complex_f{});
 313
 314     /* Load the samples from the buffer. */
 315     const size_t srclinelength{RoundUp(buffer->mSampleLen+DecoderPadding, 16)};
 316     auto srcsamples = std::make_unique<float[]>(srclinelength * numChannels);
 317     std::fill_n(srcsamples.get(), srclinelength * numChannels, 0.0f);
 318     for(size_t c{0};c < numChannels && c < realChannels;++c)
 319         LoadSamples(srcsamples.get() + srclinelength*c, buffer->mData.data() + bytesPerSample*c,
 320             realChannels, buffer->mType, buffer->mSampleLen);
 321
 322     if(IsUHJ(mChannels))
 323     {
 324         auto decoder = std::make_unique<UhjDecoderType>();
 325         std::array<float*,4> samples{};
 326         for(size_t c{0};c < numChannels;++c)
 327             samples[c] = srcsamples.get() + srclinelength*c;
 328         decoder->decode({samples.data(), numChannels}, buffer->mSampleLen, buffer->mSampleLen);
 329     }
 330
 331     auto ressamples = std::make_unique<double[]>(buffer->mSampleLen +
 332         (resampler ? resampledCount : 0));
 333     complex_f *filteriter = mComplexData.get() + mNumConvolveSegs*m;
 334     for(size_t c{0};c < numChannels;++c)
 335     {
 336         /* Resample to match the device. */
 337         if(resampler)
 338         {
 339             std::copy_n(srcsamples.get() + srclinelength*c, buffer->mSampleLen,
 340                 ressamples.get() + resampledCount);
 341             resampler.process(buffer->mSampleLen, ressamples.get()+resampledCount,
 342                 resampledCount, ressamples.get());
 343         }
 344         else
 345             std::copy_n(srcsamples.get() + srclinelength*c, buffer->mSampleLen, ressamples.get());
 346
 347         /* Store the first segment's samples in reverse in the time-domain, to
 348          * apply as a FIR filter.
 349          */
 350         const size_t first_size{minz(resampledCount, ConvolveUpdateSamples)};
 351         std::transform(ressamples.get(), ressamples.get()+first_size, mFilter[c].rbegin(),
 352             [](const double d) noexcept -> float { return static_cast<float>(d); });
 353
 354         auto fftbuffer = std::vector<std::complex<double>>(ConvolveUpdateSize);
 355         size_t done{first_size};
 356         for(size_t s{0};s < mNumConvolveSegs;++s)
 357         {
 358             const size_t todo{minz(resampledCount-done, ConvolveUpdateSamples)};
 359
 360             auto iter = std::copy_n(&ressamples[done], todo, fftbuffer.begin());
 361             done += todo;
 362             std::fill(iter, fftbuffer.end(), std::complex<double>{});
 363
 364             forward_fft(al::as_span(fftbuffer));
 365             filteriter = std::copy_n(fftbuffer.cbegin(), m, filteriter);
 366         }
 367     }
 368 }
 369
 370
 371 void ConvolutionState::update(const ContextBase *context, const EffectSlot *slot,
 372     const EffectProps* /*props*/, const EffectTarget target)
 373 {
 374     /* NOTE: Stereo and Rear are slightly different from normal mixing (as
 375      * defined in alu.cpp). These are 45 degrees from center, rather than the
 376      * 30 degrees used there.
 377      *
 378      * TODO: LFE is not mixed to output. This will require each buffer channel
 379      * to have its own output target since the main mixing buffer won't have an
 380      * LFE channel (due to being B-Format).
 381      */
 382     static constexpr ChanMap MonoMap[1]{
 383         { FrontCenter, 0.0f, 0.0f }
 384     }, StereoMap[2]{
 385         { FrontLeft,  Deg2Rad(-45.0f), Deg2Rad(0.0f) },
 386         { FrontRight, Deg2Rad( 45.0f), Deg2Rad(0.0f) }
 387     }, RearMap[2]{
 388         { BackLeft,  Deg2Rad(-135.0f), Deg2Rad(0.0f) },
 389         { BackRight, Deg2Rad( 135.0f), Deg2Rad(0.0f) }
 390     }, QuadMap[4]{
 391         { FrontLeft,  Deg2Rad( -45.0f), Deg2Rad(0.0f) },
 392         { FrontRight, Deg2Rad(  45.0f), Deg2Rad(0.0f) },
 393         { BackLeft,   Deg2Rad(-135.0f), Deg2Rad(0.0f) },
 394         { BackRight,  Deg2Rad( 135.0f), Deg2Rad(0.0f) }
 395     }, X51Map[6]{
 396         { FrontLeft,   Deg2Rad( -30.0f), Deg2Rad(0.0f) },
 397         { FrontRight,  Deg2Rad(  30.0f), Deg2Rad(0.0f) },
 398         { FrontCenter, Deg2Rad(   0.0f), Deg2Rad(0.0f) },
 399         { LFE, 0.0f, 0.0f },
 400         { SideLeft,    Deg2Rad(-110.0f), Deg2Rad(0.0f) },
 401         { SideRight,   Deg2Rad( 110.0f), Deg2Rad(0.0f) }
 402     }, X61Map[7]{
 403         { FrontLeft,   Deg2Rad(-30.0f), Deg2Rad(0.0f) },
 404         { FrontRight,  Deg2Rad( 30.0f), Deg2Rad(0.0f) },
 405         { FrontCenter, Deg2Rad(  0.0f), Deg2Rad(0.0f) },
 406         { LFE, 0.0f, 0.0f },
 407         { BackCenter,  Deg2Rad(180.0f), Deg2Rad(0.0f) },
 408         { SideLeft,    Deg2Rad(-90.0f), Deg2Rad(0.0f) },
 409         { SideRight,   Deg2Rad( 90.0f), Deg2Rad(0.0f) }
 410     }, X71Map[8]{
 411         { FrontLeft,   Deg2Rad( -30.0f), Deg2Rad(0.0f) },
 412         { FrontRight,  Deg2Rad(  30.0f), Deg2Rad(0.0f) },
 413         { FrontCenter, Deg2Rad(   0.0f), Deg2Rad(0.0f) },
 414         { LFE, 0.0f, 0.0f },
 415         { BackLeft,    Deg2Rad(-150.0f), Deg2Rad(0.0f) },
 416         { BackRight,   Deg2Rad( 150.0f), Deg2Rad(0.0f) },
 417         { SideLeft,    Deg2Rad( -90.0f), Deg2Rad(0.0f) },
 418         { SideRight,   Deg2Rad(  90.0f), Deg2Rad(0.0f) }
 419     };
 420
 421     if(mNumConvolveSegs < 1) UNLIKELY
 422         return;
 423
 424     mMix = &ConvolutionState::NormalMix;
 425
 426     for(auto &chan : *mChans)
 427         std::fill(std::begin(chan.Target), std::end(chan.Target), 0.0f);
 428     const float gain{slot->Gain};
 429     if(IsAmbisonic(mChannels))
 430     {
 431         DeviceBase *device{context->mDevice};
 432         if(mChannels == FmtUHJ2 && !device->mUhjEncoder)
 433         {
 434             mMix = &ConvolutionState::UpsampleMix;
 435             (*mChans)[0].mHfScale = 1.0f;
 436             (*mChans)[0].mLfScale = DecoderBase::sWLFScale;
 437             (*mChans)[1].mHfScale = 1.0f;
 438             (*mChans)[1].mLfScale = DecoderBase::sXYLFScale;
 439             (*mChans)[2].mHfScale = 1.0f;
 440             (*mChans)[2].mLfScale = DecoderBase::sXYLFScale;
 441         }
 442         else if(device->mAmbiOrder > mAmbiOrder)
 443         {
 444             mMix = &ConvolutionState::UpsampleMix;
 445             const auto scales = AmbiScale::GetHFOrderScales(mAmbiOrder, device->mAmbiOrder,
 446                 device->m2DMixing);
 447             (*mChans)[0].mHfScale = scales[0];
 448             (*mChans)[0].mLfScale = 1.0f;
 449             for(size_t i{1};i < mChans->size();++i)
 450             {
 451                 (*mChans)[i].mHfScale = scales[1];
 452                 (*mChans)[i].mLfScale = 1.0f;
 453             }
 454         }
 455         mOutTarget = target.Main->Buffer;
 456
 457         auto&& scales = GetAmbiScales(mAmbiScaling);
 458         const uint8_t *index_map{Is2DAmbisonic(mChannels) ?
 459             GetAmbi2DLayout(mAmbiLayout).data() :
 460             GetAmbiLayout(mAmbiLayout).data()};
 461
 462         std::array<float,MaxAmbiChannels> coeffs{};
 463         for(size_t c{0u};c < mChans->size();++c)
 464         {
 465             const size_t acn{index_map[c]};
 466             coeffs[acn] = scales[acn];
 467             ComputePanGains(target.Main, coeffs.data(), gain, (*mChans)[c].Target);
 468             coeffs[acn] = 0.0f;
 469         }
 470     }
 471     else
 472     {
 473         DeviceBase *device{context->mDevice};
 474         al::span<const ChanMap> chanmap{};
 475         switch(mChannels)
 476         {
 477         case FmtMono: chanmap = MonoMap; break;
 478         case FmtSuperStereo:
 479         case FmtStereo: chanmap = StereoMap; break;
 480         case FmtRear: chanmap = RearMap; break;
 481         case FmtQuad: chanmap = QuadMap; break;
 482         case FmtX51: chanmap = X51Map; break;
 483         case FmtX61: chanmap = X61Map; break;
 484         case FmtX71: chanmap = X71Map; break;
 485         case FmtBFormat2D:
 486         case FmtBFormat3D:
 487         case FmtUHJ2:
 488         case FmtUHJ3:
 489         case FmtUHJ4:
 490             break;
 491         }
 492
 493         mOutTarget = target.Main->Buffer;
 494         if(device->mRenderMode == RenderMode::Pairwise)
 495         {
 496             auto ScaleAzimuthFront = [](float azimuth, float scale) -> float
 497             {
 498                 constexpr float half_pi{al::numbers::pi_v<float>*0.5f};
 499                 const float abs_azi{std::fabs(azimuth)};
 500                 if(!(abs_azi >= half_pi))
 501                     return std::copysign(minf(abs_azi*scale, half_pi), azimuth);
 502                 return azimuth;
 503             };
 504
 505             for(size_t i{0};i < chanmap.size();++i)
 506             {
 507                 if(chanmap[i].channel == LFE) continue;
 508                 const auto coeffs = CalcAngleCoeffs(ScaleAzimuthFront(chanmap[i].angle, 2.0f),
 509                     chanmap[i].elevation, 0.0f);
 510                 ComputePanGains(target.Main, coeffs.data(), gain, (*mChans)[i].Target);
 511             }
 512         }
 513         else for(size_t i{0};i < chanmap.size();++i)
 514         {
 515             if(chanmap[i].channel == LFE) continue;
 516             const auto coeffs = CalcAngleCoeffs(chanmap[i].angle, chanmap[i].elevation, 0.0f);
 517             ComputePanGains(target.Main, coeffs.data(), gain, (*mChans)[i].Target);
 518         }
 519     }
 520 }
 521
 522 void ConvolutionState::process(const size_t samplesToDo,
 523     const al::span<const FloatBufferLine> samplesIn, const al::span<FloatBufferLine> samplesOut)
 524 {
 525     if(mNumConvolveSegs < 1) UNLIKELY
 526         return;
 527
 528     constexpr size_t m{ConvolveUpdateSize/2 + 1};
 529     size_t curseg{mCurrentSegment};
 530     auto &chans = *mChans;
 531
 532     for(size_t base{0u};base < samplesToDo;)
 533     {
 534         const size_t todo{minz(ConvolveUpdateSamples-mFifoPos, samplesToDo-base)};
 535
 536         std::copy_n(samplesIn[0].begin() + base, todo,
 537             mInput.begin()+ConvolveUpdateSamples+mFifoPos);
 538
 539         /* Apply the FIR for the newly retrieved input samples, and combine it
 540          * with the inverse FFT'd output samples.
 541          */
 542         for(size_t c{0};c < chans.size();++c)
 543         {
 544             auto buf_iter = chans[c].mBuffer.begin() + base;
 545             apply_fir({buf_iter, todo}, mInput.data()+1 + mFifoPos, mFilter[c].data());
 546
 547             auto fifo_iter = mOutput[c].begin() + mFifoPos;
 548             std::transform(fifo_iter, fifo_iter+todo, buf_iter, buf_iter, std::plus<>{});
 549         }
 550
 551         mFifoPos += todo;
 552         base += todo;
 553
 554         /* Check whether the input buffer is filled with new samples. */
 555         if(mFifoPos < ConvolveUpdateSamples) break;
 556         mFifoPos = 0;
 557
 558         /* Move the newest input to the front for the next iteration's history. */
 559         std::copy(mInput.cbegin()+ConvolveUpdateSamples, mInput.cend(), mInput.begin());
 560
 561         /* Calculate the frequency domain response and add the relevant
 562          * frequency bins to the FFT history.
 563          */
 564         auto fftiter = std::copy_n(mInput.cbegin(), ConvolveUpdateSamples, mFftBuffer.begin());
 565         std::fill(fftiter, mFftBuffer.end(), complex_f{});
 566         forward_fft(al::as_span(mFftBuffer));
 567
 568         std::copy_n(mFftBuffer.cbegin(), m, &mComplexData[curseg*m]);
 569
 570         const complex_f *RESTRICT filter{mComplexData.get() + mNumConvolveSegs*m};
 571         for(size_t c{0};c < chans.size();++c)
 572         {
 573             std::fill_n(mFftBuffer.begin(), m, complex_f{});
 574
 575             /* Convolve each input segment with its IR filter counterpart
 576              * (aligned in time).
 577              */
 578             const complex_f *RESTRICT input{&mComplexData[curseg*m]};
 579             for(size_t s{curseg};s < mNumConvolveSegs;++s)
 580             {
 581                 for(size_t i{0};i < m;++i,++input,++filter)
 582                     mFftBuffer[i] += *input * *filter;
 583             }
 584             input = mComplexData.get();
 585             for(size_t s{0};s < curseg;++s)
 586             {
 587                 for(size_t i{0};i < m;++i,++input,++filter)
 588                     mFftBuffer[i] += *input * *filter;
 589             }
 590
 591             /* Reconstruct the mirrored/negative frequencies to do a proper
 592              * inverse FFT.
 593              */
 594             for(size_t i{m};i < ConvolveUpdateSize;++i)
 595                 mFftBuffer[i] = std::conj(mFftBuffer[ConvolveUpdateSize-i]);
 596
 597             /* Apply iFFT to get the 256 (really 255) samples for output. The
 598              * 128 output samples are combined with the last output's 127
 599              * second-half samples (and this output's second half is
 600              * subsequently saved for next time).
 601              */
 602             inverse_fft(al::as_span(mFftBuffer));
 603
 604             /* The iFFT'd response is scaled up by the number of bins, so apply
 605              * the inverse to normalize the output.
 606              */
 607             for(size_t i{0};i < ConvolveUpdateSamples;++i)
 608                 mOutput[c][i] =
 609                     (mFftBuffer[i].real()+mOutput[c][ConvolveUpdateSamples+i]) *
 610                     (1.0f/float{ConvolveUpdateSize});
 611             for(size_t i{0};i < ConvolveUpdateSamples;++i)
 612                 mOutput[c][ConvolveUpdateSamples+i] = mFftBuffer[ConvolveUpdateSamples+i].real();
 613         }
 614
 615         /* Shift the input history. */
 616         curseg = curseg ? (curseg-1) : (mNumConvolveSegs-1);
 617     }
 618     mCurrentSegment = curseg;
 619
 620     /* Finally, mix to the output. */
 621     (this->*mMix)(samplesOut, samplesToDo);
 622 }
 623
 624
 625 struct ConvolutionStateFactory final : public EffectStateFactory {
 626     al::intrusive_ptr<EffectState> create() override
 627     { return al::intrusive_ptr<EffectState>{new ConvolutionState{}}; }
 628 };
 629
 630 } // namespace
 631
 632 EffectStateFactory *ConvolutionStateFactory_getFactory()
 633 {
 634     static ConvolutionStateFactory ConvolutionFactory{};
 635     return &ConvolutionFactory;
 636 }