14 #ifdef HAVE_SSE_INTRINSICS
15 #include <xmmintrin.h>
16 #elif defined(HAVE_NEON)
21 #include "alcomplex.h"
23 #include "alnumbers.h"
24 #include "alnumeric.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"
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.
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.
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.
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.
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.
75 void LoadSamples(float *RESTRICT dst, const al::byte *src, const size_t srcstep, FmtType srctype,
76 const size_t samples) noexcept
78 #define HANDLE_FMT(T) case T: al::LoadSampleArray<T>(dst, src, srcstep, samples); break
84 HANDLE_FMT(FmtDouble);
87 /* FIXME: Handle ADPCM decoding here. */
90 std::fill_n(dst, samples, 0.0f);
97 inline auto& GetAmbiScales(AmbiScaling scaletype) noexcept
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;
106 return AmbiScale::FromN3D();
109 inline auto& GetAmbiLayout(AmbiLayout layouttype) noexcept
111 if(layouttype == AmbiLayout::FuMa) return AmbiIndex::FromFuMa();
112 return AmbiIndex::FromACN();
115 inline auto& GetAmbi2DLayout(AmbiLayout layouttype) noexcept
117 if(layouttype == AmbiLayout::FuMa) return AmbiIndex::FromFuMa2D();
118 return AmbiIndex::FromACN2D();
128 constexpr float Deg2Rad(float x) noexcept
129 { return static_cast<float>(al::numbers::pi / 180.0 * x); }
132 using complex_f = std::complex<float>;
134 constexpr size_t ConvolveUpdateSize{256};
135 constexpr size_t ConvolveUpdateSamples{ConvolveUpdateSize / 2};
138 void apply_fir(al::span<float> dst, const float *RESTRICT src, const float *RESTRICT filter)
140 #ifdef HAVE_SSE_INTRINSICS
141 for(float &output : dst)
143 __m128 r4{_mm_setzero_ps()};
144 for(size_t j{0};j < ConvolveUpdateSamples;j+=4)
146 const __m128 coeffs{_mm_load_ps(&filter[j])};
147 const __m128 s{_mm_loadu_ps(&src[j])};
149 r4 = _mm_add_ps(r4, _mm_mul_ps(s, coeffs));
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);
158 #elif defined(HAVE_NEON)
160 for(float &output : dst)
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);
173 for(float &output : dst)
176 for(size_t j{0};j < ConvolveUpdateSamples;++j)
177 ret += src[j] * filter[j];
184 struct ConvolutionState final : public EffectState {
185 FmtChannels mChannels{};
186 AmbiLayout mAmbiLayout{};
187 AmbiScaling mAmbiScaling{};
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;
195 alignas(16) std::array<complex_f,ConvolveUpdateSize> mFftBuffer{};
197 size_t mCurrentSegment{0};
198 size_t mNumConvolveSegs{0};
201 alignas(16) FloatBufferLine mBuffer{};
202 float mHfScale{}, mLfScale{};
203 BandSplitter mFilter{};
204 float Current[MAX_OUTPUT_CHANNELS]{};
205 float Target[MAX_OUTPUT_CHANNELS]{};
207 using ChannelDataArray = al::FlexArray<ChannelData>;
208 std::unique_ptr<ChannelDataArray> mChans;
209 std::unique_ptr<complex_f[]> mComplexData;
212 ConvolutionState() = default;
213 ~ConvolutionState() override = default;
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};
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;
226 DEF_NEWDEL(ConvolutionState)
229 void ConvolutionState::NormalMix(const al::span<FloatBufferLine> samplesOut,
230 const size_t samplesToDo)
232 for(auto &chan : *mChans)
233 MixSamples({chan.mBuffer.data(), samplesToDo}, samplesOut, chan.Current, chan.Target,
237 void ConvolutionState::UpsampleMix(const al::span<FloatBufferLine> samplesOut,
238 const size_t samplesToDo)
240 for(auto &chan : *mChans)
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);
249 void ConvolutionState::deviceUpdate(const DeviceBase *device, const BufferStorage *buffer)
251 using UhjDecoderType = UhjDecoder<512>;
252 static constexpr auto DecoderPadding = UhjDecoderType::sInputPadding;
254 constexpr uint MaxConvolveAmbiOrder{1u};
258 decltype(mFilter){}.swap(mFilter);
259 decltype(mOutput){}.swap(mOutput);
260 mFftBuffer.fill(complex_f{});
263 mNumConvolveSegs = 0;
266 mComplexData = nullptr;
268 /* An empty buffer doesn't need a convolution filter. */
269 if(!buffer || buffer->mSampleLen < 1) return;
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);
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);
281 mChans = ChannelDataArray::Create(numChannels);
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.
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);
295 const BandSplitter splitter{device->mXOverFreq / static_cast<float>(device->Frequency)};
296 for(auto &e : *mChans)
297 e.mFilter = splitter;
299 mFilter.resize(numChannels, {});
300 mOutput.resize(numChannels, {});
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.
307 mNumConvolveSegs = (resampledCount+(ConvolveUpdateSamples-1)) / ConvolveUpdateSamples;
308 mNumConvolveSegs = maxz(mNumConvolveSegs, 2) - 1;
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{});
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);
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);
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)
336 /* Resample to match the device. */
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());
345 std::copy_n(srcsamples.get() + srclinelength*c, buffer->mSampleLen, ressamples.get());
347 /* Store the first segment's samples in reverse in the time-domain, to
348 * apply as a FIR filter.
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); });
354 auto fftbuffer = std::vector<std::complex<double>>(ConvolveUpdateSize);
355 size_t done{first_size};
356 for(size_t s{0};s < mNumConvolveSegs;++s)
358 const size_t todo{minz(resampledCount-done, ConvolveUpdateSamples)};
360 auto iter = std::copy_n(&ressamples[done], todo, fftbuffer.begin());
362 std::fill(iter, fftbuffer.end(), std::complex<double>{});
364 forward_fft(al::as_span(fftbuffer));
365 filteriter = std::copy_n(fftbuffer.cbegin(), m, filteriter);
371 void ConvolutionState::update(const ContextBase *context, const EffectSlot *slot,
372 const EffectProps* /*props*/, const EffectTarget target)
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.
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).
382 static constexpr ChanMap MonoMap[1]{
383 { FrontCenter, 0.0f, 0.0f }
385 { FrontLeft, Deg2Rad(-45.0f), Deg2Rad(0.0f) },
386 { FrontRight, Deg2Rad( 45.0f), Deg2Rad(0.0f) }
388 { BackLeft, Deg2Rad(-135.0f), Deg2Rad(0.0f) },
389 { BackRight, Deg2Rad( 135.0f), Deg2Rad(0.0f) }
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) }
396 { FrontLeft, Deg2Rad( -30.0f), Deg2Rad(0.0f) },
397 { FrontRight, Deg2Rad( 30.0f), Deg2Rad(0.0f) },
398 { FrontCenter, Deg2Rad( 0.0f), Deg2Rad(0.0f) },
400 { SideLeft, Deg2Rad(-110.0f), Deg2Rad(0.0f) },
401 { SideRight, Deg2Rad( 110.0f), Deg2Rad(0.0f) }
403 { FrontLeft, Deg2Rad(-30.0f), Deg2Rad(0.0f) },
404 { FrontRight, Deg2Rad( 30.0f), Deg2Rad(0.0f) },
405 { FrontCenter, Deg2Rad( 0.0f), Deg2Rad(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) }
411 { FrontLeft, Deg2Rad( -30.0f), Deg2Rad(0.0f) },
412 { FrontRight, Deg2Rad( 30.0f), Deg2Rad(0.0f) },
413 { FrontCenter, Deg2Rad( 0.0f), Deg2Rad(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) }
421 if(mNumConvolveSegs < 1) UNLIKELY
424 mMix = &ConvolutionState::NormalMix;
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))
431 DeviceBase *device{context->mDevice};
432 if(mChannels == FmtUHJ2 && !device->mUhjEncoder)
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;
442 else if(device->mAmbiOrder > mAmbiOrder)
444 mMix = &ConvolutionState::UpsampleMix;
445 const auto scales = AmbiScale::GetHFOrderScales(mAmbiOrder, device->mAmbiOrder,
447 (*mChans)[0].mHfScale = scales[0];
448 (*mChans)[0].mLfScale = 1.0f;
449 for(size_t i{1};i < mChans->size();++i)
451 (*mChans)[i].mHfScale = scales[1];
452 (*mChans)[i].mLfScale = 1.0f;
455 mOutTarget = target.Main->Buffer;
457 auto&& scales = GetAmbiScales(mAmbiScaling);
458 const uint8_t *index_map{Is2DAmbisonic(mChannels) ?
459 GetAmbi2DLayout(mAmbiLayout).data() :
460 GetAmbiLayout(mAmbiLayout).data()};
462 std::array<float,MaxAmbiChannels> coeffs{};
463 for(size_t c{0u};c < mChans->size();++c)
465 const size_t acn{index_map[c]};
466 coeffs[acn] = scales[acn];
467 ComputePanGains(target.Main, coeffs.data(), gain, (*mChans)[c].Target);
473 DeviceBase *device{context->mDevice};
474 al::span<const ChanMap> chanmap{};
477 case FmtMono: chanmap = MonoMap; break;
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;
493 mOutTarget = target.Main->Buffer;
494 if(device->mRenderMode == RenderMode::Pairwise)
496 auto ScaleAzimuthFront = [](float azimuth, float scale) -> float
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);
505 for(size_t i{0};i < chanmap.size();++i)
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);
513 else for(size_t i{0};i < chanmap.size();++i)
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);
522 void ConvolutionState::process(const size_t samplesToDo,
523 const al::span<const FloatBufferLine> samplesIn, const al::span<FloatBufferLine> samplesOut)
525 if(mNumConvolveSegs < 1) UNLIKELY
528 constexpr size_t m{ConvolveUpdateSize/2 + 1};
529 size_t curseg{mCurrentSegment};
530 auto &chans = *mChans;
532 for(size_t base{0u};base < samplesToDo;)
534 const size_t todo{minz(ConvolveUpdateSamples-mFifoPos, samplesToDo-base)};
536 std::copy_n(samplesIn[0].begin() + base, todo,
537 mInput.begin()+ConvolveUpdateSamples+mFifoPos);
539 /* Apply the FIR for the newly retrieved input samples, and combine it
540 * with the inverse FFT'd output samples.
542 for(size_t c{0};c < chans.size();++c)
544 auto buf_iter = chans[c].mBuffer.begin() + base;
545 apply_fir({buf_iter, todo}, mInput.data()+1 + mFifoPos, mFilter[c].data());
547 auto fifo_iter = mOutput[c].begin() + mFifoPos;
548 std::transform(fifo_iter, fifo_iter+todo, buf_iter, buf_iter, std::plus<>{});
554 /* Check whether the input buffer is filled with new samples. */
555 if(mFifoPos < ConvolveUpdateSamples) break;
558 /* Move the newest input to the front for the next iteration's history. */
559 std::copy(mInput.cbegin()+ConvolveUpdateSamples, mInput.cend(), mInput.begin());
561 /* Calculate the frequency domain response and add the relevant
562 * frequency bins to the FFT history.
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));
568 std::copy_n(mFftBuffer.cbegin(), m, &mComplexData[curseg*m]);
570 const complex_f *RESTRICT filter{mComplexData.get() + mNumConvolveSegs*m};
571 for(size_t c{0};c < chans.size();++c)
573 std::fill_n(mFftBuffer.begin(), m, complex_f{});
575 /* Convolve each input segment with its IR filter counterpart
578 const complex_f *RESTRICT input{&mComplexData[curseg*m]};
579 for(size_t s{curseg};s < mNumConvolveSegs;++s)
581 for(size_t i{0};i < m;++i,++input,++filter)
582 mFftBuffer[i] += *input * *filter;
584 input = mComplexData.get();
585 for(size_t s{0};s < curseg;++s)
587 for(size_t i{0};i < m;++i,++input,++filter)
588 mFftBuffer[i] += *input * *filter;
591 /* Reconstruct the mirrored/negative frequencies to do a proper
594 for(size_t i{m};i < ConvolveUpdateSize;++i)
595 mFftBuffer[i] = std::conj(mFftBuffer[ConvolveUpdateSize-i]);
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).
602 inverse_fft(al::as_span(mFftBuffer));
604 /* The iFFT'd response is scaled up by the number of bins, so apply
605 * the inverse to normalize the output.
607 for(size_t i{0};i < ConvolveUpdateSamples;++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();
615 /* Shift the input history. */
616 curseg = curseg ? (curseg-1) : (mNumConvolveSegs-1);
618 mCurrentSegment = curseg;
620 /* Finally, mix to the output. */
621 (this->*mMix)(samplesOut, samplesToDo);
625 struct ConvolutionStateFactory final : public EffectStateFactory {
626 al::intrusive_ptr<EffectState> create() override
627 { return al::intrusive_ptr<EffectState>{new ConvolutionState{}}; }
632 EffectStateFactory *ConvolutionStateFactory_getFactory()
634 static ConvolutionStateFactory ConvolutionFactory{};
635 return &ConvolutionFactory;