🛠️🐜 Antkeeper superbuild with dependencies included https://antkeeper.com
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

612 lines
22 KiB

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