fft.h

#ifndef SIGNALSMITH_AUDIO_LINEAR_FFT_H
#define SIGNALSMITH_AUDIO_LINEAR_FFT_H

#include <complex>
#include <vector>
#include <cmath>

#if defined(__FAST_MATH__) && (__apple_build_version__ >= 16000000) && (__apple_build_version__ <= 16000099) && !defined(SIGNALSMITH_IGNORE_BROKEN_APPLECLANG)
#	error Apple Clang 16.0.0 generates incorrect SIMD for ARM. If you HAVE to use this version of Clang, turn off -ffast-math.
#endif

#ifndef M_PI
#	define M_PI 3.14159265358979323846
#endif

namespace signalsmith { namespace linear {

namespace _impl {
	// Helpers for complex arithmetic, ignoring the NaN/Inf edge-cases you get without `-ffast-math`
	template<class V>
	void complexMul(std::complex<V> *a, const std::complex<V> *b, const std::complex<V> *c, size_t size) {
		for (size_t i = 0; i < size; ++i) {
			auto bi = b[i], ci = c[i];
			a[i] = {bi.real()*ci.real() - bi.imag()*ci.imag(), bi.imag()*ci.real() + bi.real()*ci.imag()};
		}
	}
	template<class V>
	void complexMulConj(std::complex<V> *a, const std::complex<V> *b, const std::complex<V> *c, size_t size) {
		for (size_t i = 0; i < size; ++i) {
			auto bi = b[i], ci = c[i];
			a[i] = {bi.real()*ci.real() + bi.imag()*ci.imag(), bi.imag()*ci.real() - bi.real()*ci.imag()};
		}
	}
	template<class V>
	void complexMul(V *ar, V *ai, const V *br, const V *bi, const V *cr, const V *ci, size_t size) {
		for (size_t i = 0; i < size; ++i) {
			V rr = br[i]*cr[i] - bi[i]*ci[i];
			V ri = br[i]*ci[i] + bi[i]*cr[i];
			ar[i] = rr;
			ai[i] = ri;
		}
	}
	template<class V>
	void complexMulConj(V *ar, V *ai, const V *br, const V *bi, const V *cr, const V *ci, size_t size) {
		for (size_t i = 0; i < size; ++i) {
			V rr = cr[i]*br[i] + ci[i]*bi[i];
			V ri = cr[i]*bi[i] - ci[i]*br[i];
			ar[i] = rr;
			ai[i] = ri;
		}
	}

	// Input: aStride elements next to each other -> output with bStride
	template<size_t aStride, class V>
	void interleaveCopy(const V *a, V *b, size_t bStride) {
		for (size_t bi = 0; bi < bStride; ++bi) {
			const V *offsetA = a + bi*aStride;
			V *offsetB = b + bi;
			for (size_t ai = 0; ai < aStride; ++ai) {
				offsetB[ai*bStride] = offsetA[ai];
			}
		}
	}
	template<class V>
	void interleaveCopy(const V *a, V *b, size_t aStride, size_t bStride) {
		for (size_t bi = 0; bi < bStride; ++bi) {
			const V *offsetA = a + bi*aStride;
			V *offsetB = b + bi;
			for (size_t ai = 0; ai < aStride; ++ai) {
				offsetB[ai*bStride] = offsetA[ai];
			}
		}
	}
	template<size_t aStride, class V>
	void interleaveCopy(const V *aReal, const V *aImag, V *bReal, V *bImag, size_t bStride) {
		for (size_t bi = 0; bi < bStride; ++bi) {
			const V *offsetAr = aReal + bi*aStride;
			const V *offsetAi = aImag + bi*aStride;
			V *offsetBr = bReal + bi;
			V *offsetBi = bImag + bi;
			for (size_t ai = 0; ai < aStride; ++ai) {
				offsetBr[ai*bStride] = offsetAr[ai];
				offsetBi[ai*bStride] = offsetAi[ai];
			}
		}
	}
	template<class V>
	void interleaveCopy(const V *aReal, const V *aImag, V *bReal, V *bImag, size_t aStride, size_t bStride) {
		for (size_t bi = 0; bi < bStride; ++bi) {
			const V *offsetAr = aReal + bi*aStride;
			const V *offsetAi = aImag + bi*aStride;
			V *offsetBr = bReal + bi;
			V *offsetBi = bImag + bi;
			for (size_t ai = 0; ai < aStride; ++ai) {
				offsetBr[ai*bStride] = offsetAr[ai];
				offsetBi[ai*bStride] = offsetAi[ai];
			}
		}
	}
}

/// Fairly simple and very portable power-of-2 FFT
template<typename Sample>
struct SimpleFFT {
	using Complex = std::complex<Sample>;
	
	SimpleFFT(size_t size=0) {
		resize(size);
	}
	
	void resize(size_t size) {
		twiddles.resize(size*3/4);
		for (size_t i = 0; i < size*3/4; ++i) {
			Sample twiddlePhase = -2*M_PI*i/size;
			twiddles[i] = std::polar(Sample(1), twiddlePhase);
		}
		working.resize(size);
	}
	
	void fft(const Complex *time, Complex *freq) {
		size_t size = working.size();
		if (size <= 1) {
			*freq = *time;
			return;
		}
		fftPass<false>(size, 1, time, freq, working.data());
	}

	void ifft(const Complex *freq, Complex *time) {
		size_t size = working.size();
		if (size <= 1) {
			*time = *freq;
			return;
		}
		fftPass<true>(size, 1, freq, time, working.data());
	}

	void fft(const Sample *inR, const Sample *inI, Sample *outR, Sample *outI) {
		size_t size = working.size();
		if (size <= 1) {
			*outR = *inR;
			*outI = *inI;
			return;
		}
		Sample *workingR = (Sample *)working.data(), *workingI = workingR + size;
		fftPass<false>(size, 1, inR, inI, outR, outI, workingR, workingI);
	}
	void ifft(const Sample *inR, const Sample *inI, Sample *outR, Sample *outI) {
		size_t size = working.size();
		if (size <= 1) {
			*outR = *inR;
			*outI = *inI;
			return;
		}
		Sample *workingR = (Sample *)working.data(), *workingI = workingR + size;
		fftPass<true>(size, 1, inR, inI, outR, outI, workingR, workingI);
	}
private:
	std::vector<Complex> twiddles;
	std::vector<Complex> working;
	
	template<bool conjB>
	static Complex mul(const Complex &a, const Complex &b) {
		return conjB ? Complex{
			a.real()*b.real() + a.imag()*b.imag(),
			a.imag()*b.real() - a.real()*b.imag()
		} : Complex{
			a.real()*b.real() - a.imag()*b.imag(),
			a.imag()*b.real() + a.real()*b.imag()
		};
	}

	// Calculate a [size]-point FFT, where each element is a block of [stride] values
	template<bool inverse>
	void fftPass(size_t size, size_t stride, const Complex *input, Complex *output, Complex *working) {
		if (size/4 > 1) {
			// Calculate four quarter-size FFTs
			fftPass<inverse>(size/4, stride*4, input, working, output);
			combine4<inverse>(size, stride, working, output);
		} else if (size == 4) {
			combine4<inverse>(4, stride, input, output);
		} else {
			// 2-point FFT
			for (size_t s = 0; s < stride; ++s) {
				Complex a = input[s];
				Complex b = input[s + stride];
				output[s] = a + b;
				output[s + stride] = a - b;
			}
		}
	}

	// Combine interleaved results into a single spectrum
	template<bool inverse>
	void combine4(size_t size, size_t stride, const Complex *input, Complex *output) const {
		auto twiddleStep = working.size()/size;
		for (size_t i = 0; i < size/4; ++i) {
			Complex twiddleB = twiddles[i*twiddleStep];
			Complex twiddleC = twiddles[i*2*twiddleStep];
			Complex twiddleD = twiddles[i*3*twiddleStep];
			
			const Complex *inputA = input + 4*i*stride;
			const Complex *inputB = input + (4*i + 1)*stride;
			const Complex *inputC = input + (4*i + 2)*stride;
			const Complex *inputD = input + (4*i + 3)*stride;
			Complex *outputA = output + i*stride;
			Complex *outputB = output + (i + size/4)*stride;
			Complex *outputC = output + (i + size/4*2)*stride;
			Complex *outputD = output + (i + size/4*3)*stride;
			for (size_t s = 0; s < stride; ++s) {
				Complex a = inputA[s];
				Complex b = mul<inverse>(inputB[s], twiddleB);
				Complex c = mul<inverse>(inputC[s], twiddleC);
				Complex d = mul<inverse>(inputD[s], twiddleD);
				Complex ac0 = a + c, ac1 = a - c;
				Complex bd0 = b + d, bd1 = inverse ? (b - d) : (d - b);
				Complex bd1i = {-bd1.imag(), bd1.real()};
				outputA[s] = ac0 + bd0;
				outputB[s] = ac1 + bd1i;
				outputC[s] = ac0 - bd0;
				outputD[s] = ac1 - bd1i;
			}
		}
	}

	// The same thing, but translated for split-complex input/output
	template<bool inverse>
	void fftPass(size_t size, size_t stride, const Sample *inputR, const Sample *inputI, Sample *outputR, Sample *outputI, Sample *workingR, Sample *workingI) const {
		if (size/4 > 1) {
			// Calculate four quarter-size FFTs
			fftPass<inverse>(size/4, stride*4, inputR, inputI, workingR, workingI, outputR, outputI);
			combine4<inverse>(size, stride, workingR, workingI, outputR, outputI);
		} else if (size == 4) {
			combine4<inverse>(4, stride, inputR, inputI, outputR, outputI);
		} else {
			// 2-point FFT
			for (size_t s = 0; s < stride; ++s) {
				Sample ar = inputR[s], ai = inputI[s];
				Sample br = inputR[s + stride], bi = inputI[s + stride];
				outputR[s] = ar + br;
				outputI[s] = ai + bi;
				outputR[s + stride] = ar - br;
				outputI[s + stride] = ai - bi;
			}
		}
	}

	// Combine interleaved results into a single spectrum
	template<bool inverse>
	void combine4(size_t size, size_t stride, const Sample *inputR, const Sample *inputI, Sample *outputR, Sample *outputI) const {
		auto twiddleStep = working.size()/size;
		for (size_t i = 0; i < size/4; ++i) {
			Complex twiddleB = twiddles[i*twiddleStep];
			Complex twiddleC = twiddles[i*2*twiddleStep];
			Complex twiddleD = twiddles[i*3*twiddleStep];
			
			const Sample *inputAr = inputR + 4*i*stride, *inputAi = inputI + 4*i*stride;
			const Sample *inputBr = inputR + (4*i + 1)*stride, *inputBi = inputI + (4*i + 1)*stride;
			const Sample *inputCr = inputR + (4*i + 2)*stride, *inputCi = inputI + (4*i + 2)*stride;
			const Sample *inputDr = inputR + (4*i + 3)*stride, *inputDi = inputI + (4*i + 3)*stride;
			Sample *outputAr = outputR + i*stride, *outputAi = outputI + i*stride;
			Sample *outputBr = outputR + (i + size/4)*stride, *outputBi = outputI + (i + size/4)*stride;
			Sample *outputCr = outputR + (i + size/4*2)*stride, *outputCi = outputI + (i + size/4*2)*stride;
			Sample *outputDr = outputR + (i + size/4*3)*stride, *outputDi = outputI + (i + size/4*3)*stride;
			for (size_t s = 0; s < stride; ++s) {
				Complex a = {inputAr[s], inputAi[s]};
				Complex b = mul<inverse>({inputBr[s], inputBi[s]}, twiddleB);
				Complex c = mul<inverse>({inputCr[s], inputCi[s]}, twiddleC);
				Complex d = mul<inverse>({inputDr[s], inputDi[s]}, twiddleD);
				Complex ac0 = a + c, ac1 = a - c;
				Complex bd0 = b + d, bd1 = inverse ? (b - d) : (d - b);
				Complex bd1i = {-bd1.imag(), bd1.real()};
				outputAr[s] = ac0.real() + bd0.real();
				outputAi[s] = ac0.imag() + bd0.imag();
				outputBr[s] = ac1.real() + bd1i.real();
				outputBi[s] = ac1.imag() + bd1i.imag();
				outputCr[s] = ac0.real() - bd0.real();
				outputCi[s] = ac0.imag() - bd0.imag();
				outputDr[s] = ac1.real() - bd1i.real();
				outputDi[s] = ac1.imag() - bd1i.imag();
			}
		}
	}
};

/// A power-of-2 only FFT, specialised with platform-specific fast implementations where available
template<typename Sample>
struct Pow2FFT {
	static constexpr bool prefersSplit = true; // whether this FFT implementation is faster when given split-complex inputs
	using Complex = std::complex<Sample>;
	
	Pow2FFT(size_t size=0) {
		resize(size);
	}
	// Allow move, but not copy
	Pow2FFT(const Pow2FFT &other) = delete;
	Pow2FFT(Pow2FFT &&other) : tmp(std::move(other.tmp)), simpleFFT(std::move(other.simpleFFT)) {}
	
	void resize(size_t size) {
		simpleFFT.resize(size);
		tmp.resize(size);
	}
	
	void fft(const Complex *time, Complex *freq) {
		simpleFFT.fft(time, freq);
	}
	void fft(const Sample *inR, const Sample *inI, Sample *outR, Sample *outI) {
		simpleFFT.fft(inR, inI, outR, outI);
	}

	void ifft(const Complex *freq, Complex *time) {
		simpleFFT.ifft(freq, time);
	}
	void ifft(const Sample *inR, const Sample *inI, Sample *outR, Sample *outI) {
		simpleFFT.ifft(inR, inI, outR, outI);
	}

private:
	std::vector<Complex> tmp;
	SimpleFFT<Sample> simpleFFT;
};

/// An FFT which can handle multiples of 3 and 5, and can be computed in chunks
template<typename Sample, bool splitComputation=false>
struct SplitFFT {
	using Complex = std::complex<Sample>;
	static constexpr bool prefersSplit = Pow2FFT<Sample>::prefersSplit;
	
	static constexpr size_t maxSplit = splitComputation ? 4 : 1;
	static constexpr size_t minInnerSize = 32;
	
	static size_t fastSizeAbove(size_t size) {
		size_t pow2 = 1;
		while (pow2 < 16 && pow2 < size) pow2 *= 2;
		while (pow2*8 < size) pow2 *= 2;
		size_t multiple = (size + pow2 - 1)/pow2; // will be 1-8
		if (multiple == 7) ++multiple;
		return multiple*pow2;
	}
	
	SplitFFT(size_t size=0) {
		resize(size);
	}
	
	void resize(size_t size) {
		innerSize = 1;
		outerSize = size;

		dftTmp.resize(0);
		dftTwists.resize(0);
		plan.resize(0);
		if (!size) return;
		
		// Inner size = largest power of 2 such that either the inner size >= minInnerSize, or we have the target number of splits
		while (!(outerSize&1) && (outerSize > maxSplit || innerSize < minInnerSize)) {
			innerSize *= 2;
			outerSize /= 2;
		}
		tmpFreq.resize(size);
		innerFFT.resize(innerSize);
		
		outerTwiddles.resize(innerSize*(outerSize - 1));
		outerTwiddlesR.resize(innerSize*(outerSize - 1));
		outerTwiddlesI.resize(innerSize*(outerSize - 1));
		for (size_t i = 0; i < innerSize; ++i) {
			for (size_t s = 1; s < outerSize; ++s) {
				Sample twiddlePhase = Sample(-2*M_PI*i/innerSize*s/outerSize);
				outerTwiddles[i + (s - 1)*innerSize] = std::polar(Sample(1), twiddlePhase);
			}
		}
		for (size_t i = 0; i < outerTwiddles.size(); ++i) {
			outerTwiddlesR[i] = outerTwiddles[i].real();
			outerTwiddlesI[i] = outerTwiddles[i].imag();
		}

		StepType interleaveStep = StepType::interleaveOrderN;
		StepType finalStep = StepType::finalOrderN;
		if (outerSize == 2) {
			interleaveStep = StepType::interleaveOrder2;
			finalStep = StepType::finalOrder2;
		}
		if (outerSize == 3) {
			interleaveStep = StepType::interleaveOrder3;
			finalStep = StepType::finalOrder3;
		}
		if (outerSize == 4) {
			interleaveStep = StepType::interleaveOrder4;
			finalStep = StepType::finalOrder4;
		}
		if (outerSize == 5) {
			interleaveStep = StepType::interleaveOrder5;
			finalStep = StepType::finalOrder5;
		}
		
		if (outerSize <= 1) {
			if (size > 0) plan.push_back(Step{StepType::passthrough, 0});
		} else {
			plan.push_back({interleaveStep, 0});
			plan.push_back({StepType::firstFFT, 0});
			for (size_t s = 1; s < outerSize; ++s) {
				plan.push_back({StepType::middleFFT, s*innerSize});
			}
			plan.push_back({StepType::twiddles, 0});
			plan.push_back({finalStep, 0});

			if (finalStep == StepType::finalOrderN) {
				dftTmp.resize(outerSize);
				dftTwists.resize(outerSize);
				for (size_t s = 0; s < outerSize; ++s) {
					Sample dftPhase = Sample(-2*M_PI*s/outerSize);
					dftTwists[s] = std::polar(Sample(1), dftPhase);
				}
			}
		}
	}
	
	size_t size() const {
		return innerSize*outerSize;
	}
	size_t steps() const {
		return plan.size();
	}
	
	void fft(const Complex *time, Complex *freq) {
		for (auto &step : plan) {
			fftStep<false>(step, time, freq);
		}
	}
	void fft(size_t step, const Complex *time, Complex *freq) {
		fftStep<false>(plan[step], time, freq);
	}
	void fft(const Sample *inR, const Sample *inI, Sample *outR, Sample *outI) {
		for (auto &step : plan) {
			fftStep<false>(step, inR, inI, outR, outI);
		}
	}
	void fft(size_t step, const Sample *inR, const Sample *inI, Sample *outR, Sample *outI) {
		fftStep<false>(plan[step], inR, inI, outR, outI);
	}
	
	void ifft(const Complex *freq, Complex *time) {
		for (auto &step : plan) {
			fftStep<true>(step, freq, time);
		}
	}
	void ifft(size_t step, const Complex *freq, Complex *time) {
		fftStep<true>(plan[step], freq, time);
	}
	void ifft(const Sample *inR, const Sample *inI, Sample *outR, Sample *outI) {
		for (auto &step : plan) {
			fftStep<true>(step, inR, inI, outR, outI);
		}
	}
	void ifft(size_t step, const Sample *inR, const Sample *inI, Sample *outR, Sample *outI) {
		fftStep<true>(plan[step], inR, inI, outR, outI);
	}
private:
	using InnerFFT = Pow2FFT<Sample>;
	InnerFFT innerFFT;

	size_t innerSize, outerSize;
	std::vector<Complex> tmpFreq;
	std::vector<Complex> outerTwiddles;
	std::vector<Sample> outerTwiddlesR, outerTwiddlesI;
	std::vector<Complex> dftTwists, dftTmp;

	enum class StepType {
		passthrough,
		interleaveOrder2, interleaveOrder3, interleaveOrder4, interleaveOrder5, interleaveOrderN,
		firstFFT, middleFFT,
		twiddles,
		finalOrder2, finalOrder3, finalOrder4, finalOrder5, finalOrderN
	};
	struct Step {
		StepType type;
		size_t offset;
	};
	std::vector<Step> plan;
	
	template<bool inverse>
	void fftStep(Step step, const Complex *time, Complex *freq) {
		switch (step.type) {
			case (StepType::passthrough): {
				if (inverse) {
					innerFFT.ifft(time, freq);
				} else {
					innerFFT.fft(time, freq);
				}
				break;
			}
			case (StepType::interleaveOrder2): {
				_impl::interleaveCopy<2>(time, tmpFreq.data(), innerSize);
				break;
			}
			case (StepType::interleaveOrder3): {
				_impl::interleaveCopy<3>(time, tmpFreq.data(), innerSize);
				break;
			}
			case (StepType::interleaveOrder4): {
				_impl::interleaveCopy<4>(time, tmpFreq.data(), innerSize);
				break;
			}
			case (StepType::interleaveOrder5): {
				_impl::interleaveCopy<5>(time, tmpFreq.data(), innerSize);
				break;
			}
			case (StepType::interleaveOrderN): {
				_impl::interleaveCopy(time, tmpFreq.data(), outerSize, innerSize);
				break;
			}
			case (StepType::firstFFT): {
				if (inverse) {
					innerFFT.ifft(tmpFreq.data(), freq);
				} else {
					innerFFT.fft(tmpFreq.data(), freq);
				}
				break;
			}
			case (StepType::middleFFT): {
				Complex *offsetOut = freq + step.offset;
				if (inverse) {
					innerFFT.ifft(tmpFreq.data() + step.offset, offsetOut);
				} else {
					innerFFT.fft(tmpFreq.data() + step.offset, offsetOut);
				}
				break;
			}
			case (StepType::twiddles): {
				if (inverse) {
					_impl::complexMulConj(freq + innerSize, freq + innerSize, outerTwiddles.data(), innerSize*(outerSize - 1));
				} else {
					_impl::complexMul(freq + innerSize, freq + innerSize, outerTwiddles.data(), innerSize*(outerSize - 1));
				}
				break;
			}
			case StepType::finalOrder2:
				finalPass2(freq);
				break;
			case StepType::finalOrder3:
				finalPass3<inverse>(freq);
				break;
			case StepType::finalOrder4:
				finalPass4<inverse>(freq);
				break;
			case StepType::finalOrder5:
				finalPass5<inverse>(freq);
				break;
			case StepType::finalOrderN:
				finalPassN<inverse>(freq);
				break;
		}
	}
	template<bool inverse>
	void fftStep(Step step, const Sample *inR, const Sample *inI, Sample *outR, Sample *outI) {
		Sample *tmpR = (Sample *)tmpFreq.data(), *tmpI = tmpR + tmpFreq.size();
		switch (step.type) {
			case (StepType::passthrough): {
				if (inverse) {
					innerFFT.ifft(inR, inI, outR, outI);
				} else {
					innerFFT.fft(inR, inI, outR, outI);
				}
				break;
			}
			case (StepType::interleaveOrder2): {
				_impl::interleaveCopy<2>(inR, tmpR, innerSize);
				_impl::interleaveCopy<2>(inI, tmpI, innerSize);
				break;
			}
			case (StepType::interleaveOrder3): {
				_impl::interleaveCopy<3>(inR, tmpR, innerSize);
				_impl::interleaveCopy<3>(inI, tmpI, innerSize);
				break;
			}
			case (StepType::interleaveOrder4): {
				_impl::interleaveCopy<4>(inR, tmpR, innerSize);
				_impl::interleaveCopy<4>(inI, tmpI, innerSize);
				break;
			}
			case (StepType::interleaveOrder5): {
				_impl::interleaveCopy<5>(inR, tmpR, innerSize);
				_impl::interleaveCopy<5>(inI, tmpI, innerSize);
				break;
			}
			case (StepType::interleaveOrderN): {
				_impl::interleaveCopy(inR, inI, tmpR, tmpI, outerSize, innerSize);
				break;
			}
			case (StepType::firstFFT): {
				if (inverse) {
					innerFFT.ifft(tmpR, tmpI, outR, outI);
				} else {
					innerFFT.fft(tmpR, tmpI, outR, outI);
				}
				break;
			}
			case (StepType::middleFFT): {
				size_t offset = step.offset;
				Sample *offsetOutR = outR + offset;
				Sample *offsetOutI = outI + offset;
				if (inverse) {
					innerFFT.ifft(tmpR + offset, tmpI + offset, offsetOutR, offsetOutI);
				} else {
					innerFFT.fft(tmpR + offset, tmpI + offset, offsetOutR, offsetOutI);
				}
				break;
			}
			case(StepType::twiddles): {
				auto *twiddlesR = outerTwiddlesR.data();
				auto *twiddlesI = outerTwiddlesI.data();
				if (inverse) {
					_impl::complexMulConj(outR + innerSize, outI + innerSize, outR + innerSize, outI + innerSize, twiddlesR, twiddlesI, innerSize*(outerSize - 1));
				} else {
					_impl::complexMul(outR + innerSize, outI + innerSize, outR + innerSize, outI + innerSize, twiddlesR, twiddlesI, innerSize*(outerSize - 1));
				}
				break;
			}
			case StepType::finalOrder2:
				finalPass2(outR, outI);
				break;
			case StepType::finalOrder3:
				finalPass3<inverse>(outR, outI);
				break;
			case StepType::finalOrder4:
				finalPass4<inverse>(outR, outI);
				break;
			case StepType::finalOrder5:
				finalPass5<inverse>(outR, outI);
				break;
			case StepType::finalOrderN:
				finalPassN<inverse>(outR, outI);
				break;
		}
	}
	
	void finalPass2(Complex *f0) {
		auto *f1 = f0 + innerSize;
		for (size_t i = 0; i < innerSize; ++i) {
			Complex a = f0[i], b = f1[i];
			f0[i] = a + b;
			f1[i] = a - b;
		}
	}
	void finalPass2(Sample *f0r, Sample *f0i) {
		auto *f1r = f0r + innerSize;
		auto *f1i = f0i + innerSize;
		for (size_t i = 0; i < innerSize; ++i) {
			Sample ar = f0r[i], ai = f0i[i];
			Sample br = f1r[i], bi = f1i[i];
			f0r[i] = ar + br;
			f0i[i] = ai + bi;
			f1r[i] = ar - br;
			f1i[i] = ai - bi;
		}
	}
	template<bool inverse>
	void finalPass3(Complex *f0) {
		auto *f1 = f0 + innerSize;
		auto *f2 = f0 + innerSize*2;
		const Complex tw1{Sample(-0.5), Sample(-std::sqrt(0.75)*(inverse ? -1 : 1))};
		for (size_t i = 0; i < innerSize; ++i) {
			Complex a = f0[i], b = f1[i], c = f2[i];
			Complex bc0 = b + c, bc1 = b - c;
			f0[i] = a + bc0;
			f1[i] = {
				a.real() + bc0.real()*tw1.real() - bc1.imag()*tw1.imag(),
				a.imag() + bc0.imag()*tw1.real() + bc1.real()*tw1.imag()
			};
			f2[i] = {
				a.real() + bc0.real()*tw1.real() + bc1.imag()*tw1.imag(),
				a.imag() + bc0.imag()*tw1.real() - bc1.real()*tw1.imag()
			};
		}
	}
	template<bool inverse>
	void finalPass3(Sample *f0r, Sample *f0i) {
		auto *f1r = f0r + innerSize;
		auto *f1i = f0i + innerSize;
		auto *f2r = f0r + innerSize*2;
		auto *f2i = f0i + innerSize*2;
		const Sample tw1r = -0.5, tw1i = -std::sqrt(0.75)*(inverse ? -1 : 1);
		
		for (size_t i = 0; i < innerSize; ++i) {
			Sample ar = f0r[i], ai = f0i[i], br = f1r[i], bi = f1i[i], cr = f2r[i], ci = f2i[i];

			f0r[i] = ar + br + cr;
			f0i[i] = ai + bi + ci;
			f1r[i] = ar + br*tw1r - bi*tw1i + cr*tw1r + ci*tw1i;
			f1i[i] = ai + bi*tw1r + br*tw1i - cr*tw1i + ci*tw1r;
			f2r[i] = ar + br*tw1r + bi*tw1i + cr*tw1r - ci*tw1i;
			f2i[i] = ai + bi*tw1r - br*tw1i + cr*tw1i + ci*tw1r;
		}
	}
	template<bool inverse>
	void finalPass4(Complex *f0) {
		auto *f1 = f0 + innerSize;
		auto *f2 = f0 + innerSize*2;
		auto *f3 = f0 + innerSize*3;
		for (size_t i = 0; i < innerSize; ++i) {
			Complex a = f0[i], b = f1[i], c = f2[i], d = f3[i];
			
			Complex ac0 = a + c, ac1 = a - c;
			Complex bd0 = b + d, bd1 = inverse ? (b - d) : (d - b);
			Complex bd1i = {-bd1.imag(), bd1.real()};
			f0[i] = ac0 + bd0;
			f1[i] = ac1 + bd1i;
			f2[i] = ac0 - bd0;
			f3[i] = ac1 - bd1i;
		}
	}
	template<bool inverse>
	void finalPass4(Sample *f0r, Sample *f0i) {
		auto *f1r = f0r + innerSize;
		auto *f1i = f0i + innerSize;
		auto *f2r = f0r + innerSize*2;
		auto *f2i = f0i + innerSize*2;
		auto *f3r = f0r + innerSize*3;
		auto *f3i = f0i + innerSize*3;
		for (size_t i = 0; i < innerSize; ++i) {
			Sample ar = f0r[i], ai = f0i[i], br = f1r[i], bi = f1i[i], cr = f2r[i], ci = f2i[i], dr = f3r[i], di = f3i[i];
			
			Sample ac0r = ar + cr, ac0i = ai + ci;
			Sample ac1r = ar - cr, ac1i = ai - ci;
			Sample bd0r = br + dr, bd0i = bi + di;
			Sample bd1r = br - dr, bd1i = bi - di;
			
			f0r[i] = ac0r + bd0r;
			f0i[i] = ac0i + bd0i;
			f1r[i] = inverse ? (ac1r - bd1i) : (ac1r + bd1i);
			f1i[i] = inverse ? (ac1i + bd1r) : (ac1i - bd1r);
			f2r[i] = ac0r - bd0r;
			f2i[i] = ac0i - bd0i;
			f3r[i] = inverse ? (ac1r + bd1i) : (ac1r - bd1i);
			f3i[i] = inverse ? (ac1i - bd1r) : (ac1i + bd1r);
		}
	}
	template<bool inverse>
	void finalPass5(Complex *f0) {
		auto *f1 = f0 + innerSize;
		auto *f2 = f0 + innerSize*2;
		auto *f3 = f0 + innerSize*3;
		auto *f4 = f0 + innerSize*4;
		const Sample tw1r = 0.30901699437494745;
		const Sample tw1i = -0.9510565162951535*(inverse ? -1 : 1);
		const Sample tw2r = -0.8090169943749473;
		const Sample tw2i = -0.5877852522924732*(inverse ? -1 : 1);
		for (size_t i = 0; i < innerSize; ++i) {
			Complex a = f0[i], b = f1[i], c = f2[i], d = f3[i], e = f4[i];

			Complex be0 = b + e, be1 = {e.imag() - b.imag(), b.real() - e.real()}; // (b - e)*i
			Complex cd0 = c + d, cd1 = {d.imag() - c.imag(), c.real() - d.real()};
			
			Complex bcde01 = be0*tw1r + cd0*tw2r;
			Complex bcde02 = be0*tw2r + cd0*tw1r;
			Complex bcde11 = be1*tw1i + cd1*tw2i;
			Complex bcde12 = be1*tw2i - cd1*tw1i;

			f0[i] = a + be0 + cd0;
			f1[i] = a + bcde01 + bcde11;
			f2[i] = a + bcde02 + bcde12;
			f3[i] = a + bcde02 - bcde12;
			f4[i] = a + bcde01 - bcde11;
		}
	}
	template<bool inverse>
	void finalPass5(Sample *f0r, Sample *f0i) {
		auto *f1r = f0r + innerSize;
		auto *f1i = f0i + innerSize;
		auto *f2r = f0r + innerSize*2;
		auto *f2i = f0i + innerSize*2;
		auto *f3r = f0r + innerSize*3;
		auto *f3i = f0i + innerSize*3;
		auto *f4r = f0r + innerSize*4;
		auto *f4i = f0i + innerSize*4;
		
		const Sample tw1r = 0.30901699437494745;
		const Sample tw1i = -0.9510565162951535*(inverse ? -1 : 1);
		const Sample tw2r = -0.8090169943749473;
		const Sample tw2i = -0.5877852522924732*(inverse ? -1 : 1);
		for (size_t i = 0; i < innerSize; ++i) {
			Sample ar = f0r[i], ai = f0i[i], br = f1r[i], bi = f1i[i], cr = f2r[i], ci = f2i[i], dr = f3r[i], di = f3i[i], er = f4r[i], ei = f4i[i];

			Sample be0r = br + er, be0i = bi + ei;
			Sample be1r = ei - bi, be1i = br - er;
			Sample cd0r = cr + dr, cd0i = ci + di;
			Sample cd1r = di - ci, cd1i = cr - dr;

			Sample bcde01r = be0r*tw1r + cd0r*tw2r, bcde01i = be0i*tw1r + cd0i*tw2r;
			Sample bcde02r = be0r*tw2r + cd0r*tw1r, bcde02i = be0i*tw2r + cd0i*tw1r;
			Sample bcde11r = be1r*tw1i + cd1r*tw2i, bcde11i = be1i*tw1i + cd1i*tw2i;
			Sample bcde12r = be1r*tw2i - cd1r*tw1i, bcde12i = be1i*tw2i - cd1i*tw1i;

			f0r[i] = ar + be0r + cd0r;
			f0i[i] = ai + be0i + cd0i;
			f1r[i] = ar + bcde01r + bcde11r;
			f1i[i] = ai + bcde01i + bcde11i;
			f2r[i] = ar + bcde02r + bcde12r;
			f2i[i] = ai + bcde02i + bcde12i;
			f3r[i] = ar + bcde02r - bcde12r;
			f3i[i] = ai + bcde02i - bcde12i;
			f4r[i] = ar + bcde01r - bcde11r;
			f4i[i] = ai + bcde01i - bcde11i;
		}
	}

	template<bool inverse>
	void finalPassN(Complex *f0) {
		for (size_t i = 0; i < innerSize; ++i) {
			Complex *offsetFreq = f0 + i;
			Complex sum = 0;
			for (size_t i2 = 0; i2 < outerSize; ++i2) {
				sum += (dftTmp[i2] = offsetFreq[i2*innerSize]);
			}
			offsetFreq[0] = sum;
			
			for (size_t f = 1; f < outerSize; ++f) {
				Complex sum = dftTmp[0];

				for (size_t i2 = 1; i2 < outerSize; ++i2) {
					size_t twistIndex = (i2*f)%outerSize;
					Complex twist = inverse ? std::conj(dftTwists[twistIndex]) : dftTwists[twistIndex];
					sum += Complex{
						dftTmp[i2].real()*twist.real() - dftTmp[i2].imag()*twist.imag(),
						dftTmp[i2].imag()*twist.real() + dftTmp[i2].real()*twist.imag()
					};
				}

				offsetFreq[f*innerSize] = sum;
			}
		}
	}
	template<bool inverse>
	void finalPassN(Sample *f0r, Sample *f0i) {
		Sample *tmpR = (Sample *)dftTmp.data(), *tmpI = tmpR + outerSize;
		
		for (size_t i = 0; i < innerSize; ++i) {
			Sample *offsetR = f0r + i;
			Sample *offsetI = f0i + i;
			Sample sumR = 0, sumI = 0;
			for (size_t i2 = 0; i2 < outerSize; ++i2) {
				sumR += (tmpR[i2] = offsetR[i2*innerSize]);
				sumI += (tmpI[i2] = offsetI[i2*innerSize]);
			}
			offsetR[0] = sumR;
			offsetI[0] = sumI;
			
			for (size_t f = 1; f < outerSize; ++f) {
				Sample sumR = *tmpR, sumI = *tmpI;

				for (size_t i2 = 1; i2 < outerSize; ++i2) {
					size_t twistIndex = (i2*f)%outerSize;
					Complex twist = inverse ? std::conj(dftTwists[twistIndex]) : dftTwists[twistIndex];
					sumR += tmpR[i2]*twist.real() - tmpI[i2]*twist.imag();
					sumI += tmpI[i2]*twist.real() + tmpR[i2]*twist.imag();
				}

				offsetR[f*innerSize] = sumR;
				offsetI[f*innerSize] = sumI;
			}
		}
	}
};

template<typename Sample, bool splitComputation=false>
using FFT = SplitFFT<Sample, splitComputation>;

// Wraps a complex FFT into a real one
template<typename Sample, class ComplexFFT=Pow2FFT<Sample>>
struct SimpleRealFFT {
	using Complex = std::complex<Sample>;
	
	static constexpr bool prefersSplit = ComplexFFT::prefersSplit;

	SimpleRealFFT(size_t size=0) {
		resize(size);
	}
	
	void resize(size_t size) {
		complexFft.resize(size);
		tmpTime.resize(size);
		tmpFreq.resize(size);
	}
	
	void fft(const Sample *time, Complex *freq) {
		for (size_t i = 0; i < tmpTime.size(); ++i) {
			tmpTime[i] = time[i];
		}
		complexFft.fft(tmpTime.data(), tmpFreq.data());
		for (size_t i = 0; i < tmpFreq.size()/2; ++i) {
			freq[i] = tmpFreq[i];
		}
		freq[0] = {
			tmpFreq[0].real(),
			tmpFreq[tmpFreq.size()/2].real()
		};
	}
	void fft(const Sample *inR, Sample *outR, Sample *outI) {
		Sample *tmpFreqR = (Sample *)tmpFreq.data(), *tmpFreqI = tmpFreqR + tmpFreq.size();
		for (size_t i = 0; i < tmpTime.size()/2; ++i) {
			tmpTime[i] = 0;
		}
		complexFft.fft(inR, (const Sample *)tmpTime.data(), tmpFreqR, tmpFreqI);
		for (size_t i = 0; i < tmpTime.size()/2; ++i) {
			outR[i] = tmpFreqR[i];
			outI[i] = tmpFreqI[i];
		}
		outI[0] = tmpFreqR[tmpFreq.size()/2];
	}

	void ifft(const Complex *freq, Sample *time) {
		tmpFreq[0] = freq[0].real();
		tmpFreq[tmpFreq.size()/2] = freq[0].imag();
		for (size_t i = 1; i < tmpFreq.size()/2; ++i) {
			tmpFreq[i] = freq[i];
			tmpFreq[tmpFreq.size() - i] = std::conj(freq[i]);
		}
		complexFft.ifft(tmpFreq.data(), tmpTime.data());
		for (size_t i = 0; i < tmpTime.size(); ++i) {
			time[i] = tmpTime[i].real();
		}
	}
	void ifft(const Sample *inR, const Sample *inI, Sample *outR) {
		Sample *tmpFreqR = (Sample *)tmpFreq.data(), *tmpFreqI = tmpFreqR + tmpFreq.size();
		tmpFreqR[0] = inR[0];
		tmpFreqR[tmpFreq.size()/2] = inI[0];
		tmpFreqI[0] = 0;
		tmpFreqI[tmpFreq.size()/2] = 0;
		for (size_t i = 1; i < tmpFreq.size()/2; ++i) {
			tmpFreqR[i] = inR[i];
			tmpFreqI[i] = inI[i];
			tmpFreqR[tmpFreq.size() - i] = inR[i];
			tmpFreqI[tmpFreq.size() - i] = -inI[i];
		}
		complexFft.ifft(tmpFreqR, tmpFreqI, outR, (Sample *)tmpTime.data());
	}

private:
	ComplexFFT complexFft;
	std::vector<Complex> tmpTime, tmpFreq;
};

/// A default power-of-2 FFT, specialised with platform-specific fast implementations where available
template<typename Sample>
struct Pow2RealFFT : public SimpleRealFFT<Sample> {
	static constexpr bool prefersSplit = SimpleRealFFT<Sample>::prefersSplit;
	
	using SimpleRealFFT<Sample>::SimpleRealFFT;

	// Prevent copying, since it might be a problem for specialisations
	Pow2RealFFT(const Pow2RealFFT &other) = delete;
	// Pass move-constructor through, just to be explicit about it
	Pow2RealFFT(Pow2RealFFT &&other) : SimpleRealFFT<Sample>(std::move(other)) {}
};

/// A Real FFT which can handle multiples of 3 and 5, and can be computed in chunks
template<typename Sample, bool splitComputation=false, bool halfBinShift=false>
struct RealFFT {
	using Complex = std::complex<Sample>;
	static constexpr bool prefersSplit = SplitFFT<Sample, splitComputation>::prefersSplit;

	static size_t fastSizeAbove(size_t size) {
		return ComplexFFT::fastSizeAbove((size + 1)/2)*2;
	}
	
	RealFFT(size_t size=0) {
		resize(size);
	}
	
	void resize(size_t size) {
		size_t hSize = size/2;
		complexFft.resize(hSize);
		tmpFreq.resize(hSize);
		tmpTime.resize(hSize);
		
		twiddles.resize(hSize/2 + 1);
		
		if (!halfBinShift) {
			for (size_t i = 0; i < twiddles.size(); ++i) {
				Sample rotPhase = i*(-2*M_PI/size) - M_PI/2; // bake rotation by (-i) into twiddles
				twiddles[i] = std::polar(Sample(1), rotPhase);
			}
		} else {
			for (size_t i = 0; i < twiddles.size(); ++i) {
				Sample rotPhase = (i + 0.5)*(-2*M_PI/size) - M_PI/2;
				twiddles[i] = std::polar(Sample(1), rotPhase);
			}
	
			halfBinTwists.resize(hSize);
			for (size_t i = 0; i < hSize; ++i) {
				Sample twistPhase = -2*M_PI*i/size;
				halfBinTwists[i] = std::polar(Sample(1), twistPhase);
			}
		}
	}
	
	size_t size() const {
		return complexFft.size()*2;
	}
	size_t steps() const {
		return complexFft.steps() + (splitComputation ? 3 : 2);
	}
	
	void fft(const Sample *time, Complex *freq) {
		for (size_t s = 0; s < steps(); ++s) {
			fft(s, time, freq);
		}
	}
	void fft(size_t step, const Sample *time, Complex *freq) {
		if (complexPrefersSplit) {
			size_t hSize = complexFft.size();
			Sample *tmpTimeR = (Sample *)tmpTime.data(), *tmpTimeI = tmpTimeR + hSize;
			Sample *tmpFreqR = (Sample *)tmpFreq.data(), *tmpFreqI = tmpFreqR + hSize;
			if (step-- == 0) {
				size_t hSize = complexFft.size();
				if (halfBinShift) {
					for (size_t i = 0; i < hSize; ++i) {
						Sample tr = time[2*i], ti = time[2*i + 1];
						Complex twist = halfBinTwists[i];
						tmpTimeR[i] = tr*twist.real() - ti*twist.imag();
						tmpTimeI[i] = ti*twist.real() + tr*twist.imag();
					}
				} else {
					for (size_t i = 0; i < hSize; ++i) {
						tmpTimeR[i] = time[2*i];
						tmpTimeI[i] = time[2*i + 1];
					}
				}
			} else if (step < complexFft.steps()) {
				complexFft.fft(step, tmpTimeR, tmpTimeI, tmpFreqR, tmpFreqI);
			} else {
				if (!halfBinShift) {
					Sample bin0r = tmpFreqR[0], bin0i = tmpFreqI[0];
					freq[0] = {bin0r + bin0i, bin0r - bin0i};
				}

				size_t startI = halfBinShift ? 0 : 1;
				size_t endI = hSize/2 + 1;
				if (splitComputation) { // Do this last twiddle in two halves
					if (step == complexFft.steps()) {
						endI = (startI + endI)/2;
					} else {
						startI = (startI + endI)/2;
					}
				}
				for (size_t i = startI; i < endI; ++i) {
					size_t conjI = halfBinShift ? (hSize - 1 - i) : (hSize - i);
					Complex twiddle = twiddles[i];

					Sample oddR = (tmpFreqR[i] + tmpFreqR[conjI])*Sample(0.5);
					Sample oddI = (tmpFreqI[i] - tmpFreqI[conjI])*Sample(0.5);
					Sample evenIR = (tmpFreqR[i] - tmpFreqR[conjI])*Sample(0.5);
					Sample evenII = (tmpFreqI[i] + tmpFreqI[conjI])*Sample(0.5);
					Sample evenRotMinusIR = evenIR*twiddle.real() - evenII*twiddle.imag();
					Sample evenRotMinusII = evenII*twiddle.real() + evenIR*twiddle.imag();

					freq[i] = {oddR + evenRotMinusIR, oddI + evenRotMinusII};
					freq[conjI] = {oddR - evenRotMinusIR, evenRotMinusII - oddI};
				}
			}
		} else {
			bool canUseTime = !halfBinShift && !(size_t(time)%alignof(Complex));
			if (step-- == 0) {
				size_t hSize = complexFft.size();
				if (halfBinShift) {
					for (size_t i = 0; i < hSize; ++i) {
						Sample tr = time[2*i], ti = time[2*i + 1];
						Complex twist = halfBinTwists[i];
						tmpTime[i] = {
							tr*twist.real() - ti*twist.imag(),
							ti*twist.real() + tr*twist.imag()
						};
					}
				} else if (!canUseTime) {
					std::memcpy(tmpTime.data(), time, sizeof(Complex)*hSize);
				}
			} else if (step < complexFft.steps()) {
				complexFft.fft(step, canUseTime ? (const Complex *)time : tmpTime.data(), tmpFreq.data());
			} else {
				if (!halfBinShift) {
					Complex bin0 = tmpFreq[0];
					freq[0] = { // pack DC & Nyquist together
						bin0.real() + bin0.imag(),
						bin0.real() - bin0.imag()
					};
				}
				
				size_t hSize = complexFft.size();
				size_t startI = halfBinShift ? 0 : 1;
				size_t endI = hSize/2 + 1;
				if (splitComputation) { // Do this last twiddle in two halves
					if (step == complexFft.steps()) {
						endI = (startI + endI)/2;
					} else {
						startI = (startI + endI)/2;
					}
				}
				for (size_t i = startI; i < endI; ++i) {
					size_t conjI = halfBinShift ? (hSize - 1 - i) : (hSize - i);
					Complex twiddle = twiddles[i];

					Complex odd = (tmpFreq[i] + std::conj(tmpFreq[conjI]))*Sample(0.5);
					Complex evenI = (tmpFreq[i] - std::conj(tmpFreq[conjI]))*Sample(0.5);
					Complex evenRotMinusI = { // twiddle includes a factor of -i
						evenI.real()*twiddle.real() - evenI.imag()*twiddle.imag(),
						evenI.imag()*twiddle.real() + evenI.real()*twiddle.imag()
					};

					freq[i] = odd + evenRotMinusI;
					freq[conjI] = {odd.real() - evenRotMinusI.real(), evenRotMinusI.imag() - odd.imag()};
				}
			}
		}
	}
	void fft(const Sample *inR, Sample *outR, Sample *outI) {
		for (size_t s = 0; s < steps(); ++s) {
			fft(s, inR, outR, outI);
		}
	}
	void fft(size_t step, const Sample *inR, Sample *outR, Sample *outI) {
		size_t hSize = complexFft.size();
		Sample *tmpTimeR = (Sample *)tmpTime.data(), *tmpTimeI = tmpTimeR + hSize;
		Sample *tmpFreqR = (Sample *)tmpFreq.data(), *tmpFreqI = tmpFreqR + hSize;
		if (step-- == 0) {
			size_t hSize = complexFft.size();
			if (halfBinShift) {
				for (size_t i = 0; i < hSize; ++i) {
					Sample tr = inR[2*i], ti = inR[2*i + 1];
					Complex twist = halfBinTwists[i];
					tmpTimeR[i] = tr*twist.real() - ti*twist.imag();
					tmpTimeI[i] = ti*twist.real() + tr*twist.imag();
				}
			} else {
				for (size_t i = 0; i < hSize; ++i) {
					tmpTimeR[i] = inR[2*i];
					tmpTimeI[i] = inR[2*i + 1];
				}
			}
		} else if (step < complexFft.steps()) {
			complexFft.fft(step, tmpTimeR, tmpTimeI, tmpFreqR, tmpFreqI);
		} else {
			if (!halfBinShift) {
				Sample bin0r = tmpFreqR[0], bin0i = tmpFreqI[0];
				outR[0] = bin0r + bin0i;
				outI[0] = bin0r - bin0i;
			}

			size_t startI = halfBinShift ? 0 : 1;
			size_t endI = hSize/2 + 1;
			if (splitComputation) { // Do this last twiddle in two halves
				if (step == complexFft.steps()) {
					endI = (startI + endI)/2;
				} else {
					startI = (startI + endI)/2;
				}
			}
			for (size_t i = startI; i < endI; ++i) {
				size_t conjI = halfBinShift ? (hSize - 1 - i) : (hSize - i);
				Complex twiddle = twiddles[i];

				Sample oddR = (tmpFreqR[i] + tmpFreqR[conjI])*Sample(0.5);
				Sample oddI = (tmpFreqI[i] - tmpFreqI[conjI])*Sample(0.5);
				Sample evenIR = (tmpFreqR[i] - tmpFreqR[conjI])*Sample(0.5);
				Sample evenII = (tmpFreqI[i] + tmpFreqI[conjI])*Sample(0.5);
				Sample evenRotMinusIR = evenIR*twiddle.real() - evenII*twiddle.imag();
				Sample evenRotMinusII = evenII*twiddle.real() + evenIR*twiddle.imag();

				outR[i] = oddR + evenRotMinusIR;
				outI[i] = oddI + evenRotMinusII;
				outR[conjI] = oddR - evenRotMinusIR;
				outI[conjI] = evenRotMinusII - oddI;
			}
		}
	}
	
	void ifft(const Complex *freq, Sample *time) {
		for (size_t s = 0; s < steps(); ++s) {
			ifft(s, freq, time);
		}
	}
	void ifft(size_t step, const Complex *freq, Sample *time) {
		if (complexPrefersSplit) {
			size_t hSize = complexFft.size();
			Sample *tmpTimeR = (Sample *)tmpTime.data(), *tmpTimeI = tmpTimeR + hSize;
			Sample *tmpFreqR = (Sample *)tmpFreq.data(), *tmpFreqI = tmpFreqR + hSize;

			bool splitFirst = splitComputation && (step-- == 0);
			if (splitFirst || step-- == 0) {
				Complex bin0 = freq[0];
				if (!halfBinShift) {
					tmpFreqR[0] = bin0.real() + bin0.imag();
					tmpFreqI[0] = bin0.real() - bin0.imag();
				}
				size_t startI = halfBinShift ? 0 : 1;
				size_t endI = hSize/2 + 1;
				if (splitComputation) { // Do this first twiddle in two halves
					if (splitFirst) {
						endI = (startI + endI)/2;
					} else {
						startI = (startI + endI)/2;
					}
				}
				for (size_t i = startI; i < endI; ++i) {
					size_t conjI = halfBinShift ? (hSize - 1 - i) : (hSize - i);
					Complex twiddle = twiddles[i];

					Complex odd = freq[i] + std::conj(freq[conjI]);
					Complex evenRotMinusI = freq[i] - std::conj(freq[conjI]);
					Complex evenI = { // Conjugate twiddle
						evenRotMinusI.real()*twiddle.real() + evenRotMinusI.imag()*twiddle.imag(),
						evenRotMinusI.imag()*twiddle.real() - evenRotMinusI.real()*twiddle.imag()
					};

					tmpFreqR[i] = odd.real() + evenI.real();
					tmpFreqI[i] = odd.imag() + evenI.imag();
					tmpFreqR[conjI] = odd.real() - evenI.real();
					tmpFreqI[conjI] = evenI.imag() - odd.imag();
				}
			} else if (step < complexFft.steps()) {
				complexFft.ifft(step, tmpFreqR, tmpFreqI, tmpTimeR, tmpTimeI);
			} else {
				size_t hSize = complexFft.size();
				if (halfBinShift) {
					for (size_t i = 0; i < hSize; ++i) {
						Sample tr = tmpTimeR[i], ti = tmpTimeI[i];
						Complex twist = halfBinTwists[i];
						time[2*i] = 	tr*twist.real() + ti*twist.imag();
						time[2*i + 1] = ti*twist.real() - tr*twist.imag();
					}
				} else {
					for (size_t i = 0; i < hSize; ++i) {
						time[2*i] = tmpTimeR[i];
						time[2*i + 1] = tmpTimeI[i];
					}
				}
			}
		} else {
			bool canUseTime = !halfBinShift && !(size_t(time)%alignof(Complex));
			bool splitFirst = splitComputation && (step-- == 0);
			if (splitFirst || step-- == 0) {
				Complex bin0 = freq[0];
				if (!halfBinShift) {
					tmpFreq[0] = {
						bin0.real() + bin0.imag(),
						bin0.real() - bin0.imag()
					};
				}
				size_t hSize = complexFft.size();
				size_t startI = halfBinShift ? 0 : 1;
				size_t endI = hSize/2 + 1;
				if (splitComputation) { // Do this first twiddle in two halves
					if (splitFirst) {
						endI = (startI + endI)/2;
					} else {
						startI = (startI + endI)/2;
					}
				}
				for (size_t i = startI; i < endI; ++i) {
					size_t conjI = halfBinShift ? (hSize - 1 - i) : (hSize - i);
					Complex twiddle = twiddles[i];

					Complex odd = freq[i] + std::conj(freq[conjI]);
					Complex evenRotMinusI = freq[i] - std::conj(freq[conjI]);
					Complex evenI = { // Conjugate twiddle
						evenRotMinusI.real()*twiddle.real() + evenRotMinusI.imag()*twiddle.imag(),
						evenRotMinusI.imag()*twiddle.real() - evenRotMinusI.real()*twiddle.imag()
					};

					tmpFreq[i] = odd + evenI;
					tmpFreq[conjI] = {odd.real() - evenI.real(), evenI.imag() - odd.imag()};
				}
			} else if (step < complexFft.steps()) {
				// Can't just use time as (Complex *), since it might not be aligned properly
				complexFft.ifft(step, tmpFreq.data(), canUseTime ? (Complex *)time : tmpTime.data());
			} else {
				size_t hSize = complexFft.size();
				if (halfBinShift) {
					for (size_t i = 0; i < hSize; ++i) {
						Complex t = tmpTime[i];
						Complex twist = halfBinTwists[i];
						time[2*i] = 	t.real()*twist.real() + t.imag()*twist.imag();
						time[2*i + 1] = t.imag()*twist.real() - t.real()*twist.imag();
					}
				} else if (!canUseTime) {
					std::memcpy(time, tmpTime.data(), sizeof(Complex)*hSize);
				}
			}
		}
	}
	void ifft(const Sample *inR, const Sample *inI, Sample *outR) {
		for (size_t s = 0; s < steps(); ++s) {
			ifft(s, inR, inI, outR);
		}
	}
	void ifft(size_t step, const Sample *inR, const Sample *inI, Sample *outR) {
		size_t hSize = complexFft.size();
		Sample *tmpTimeR = (Sample *)tmpTime.data(), *tmpTimeI = tmpTimeR + hSize;
		Sample *tmpFreqR = (Sample *)tmpFreq.data(), *tmpFreqI = tmpFreqR + hSize;

		bool splitFirst = splitComputation && (step-- == 0);
		if (splitFirst || step-- == 0) {
			Sample bin0r = inR[0], bin0i = inI[0];
			if (!halfBinShift) {
				tmpFreqR[0] = bin0r + bin0i;
				tmpFreqI[0] = bin0r - bin0i;
			}
			size_t startI = halfBinShift ? 0 : 1;
			size_t endI = hSize/2 + 1;
			if (splitComputation) { // Do this first twiddle in two halves
				if (splitFirst) {
					endI = (startI + endI)/2;
				} else {
					startI = (startI + endI)/2;
				}
			}
			for (size_t i = startI; i < endI; ++i) {
				size_t conjI = halfBinShift ? (hSize - 1 - i) : (hSize - i);
				Complex twiddle = twiddles[i];
				Sample fir = inR[i], fii = inI[i];
				Sample fcir = inR[conjI], fcii = inI[conjI];

				Complex odd = {fir + fcir, fii - fcii};
				Complex evenRotMinusI = {fir - fcir, fii + fcii};
				Complex evenI = { // Conjugate twiddle
					evenRotMinusI.real()*twiddle.real() + evenRotMinusI.imag()*twiddle.imag(),
					evenRotMinusI.imag()*twiddle.real() - evenRotMinusI.real()*twiddle.imag()
				};

				tmpFreqR[i] = odd.real() + evenI.real();
				tmpFreqI[i] = odd.imag() + evenI.imag();
				tmpFreqR[conjI] = odd.real() - evenI.real();
				tmpFreqI[conjI] = evenI.imag() - odd.imag();
			}
		} else if (step < complexFft.steps()) {
			// Can't just use time as (Complex *), since it might not be aligned properly
			complexFft.ifft(step, tmpFreqR, tmpFreqI, tmpTimeR, tmpTimeI);
		} else {
			if (halfBinShift) {
				for (size_t i = 0; i < hSize; ++i) {
					Sample tr = tmpTimeR[i], ti = tmpTimeI[i];
					Complex twist = halfBinTwists[i];
					outR[2*i] = 	tr*twist.real() + ti*twist.imag();
					outR[2*i + 1] = ti*twist.real() - tr*twist.imag();
				}
			} else {
				for (size_t i = 0; i < hSize; ++i) {
					outR[2*i] = tmpTimeR[i];
					outR[2*i + 1] = tmpTimeI[i];
				}
			}
		}
	}
private:
	using ComplexFFT = SplitFFT<Sample, splitComputation>;
	ComplexFFT complexFft;

	static constexpr bool complexPrefersSplit = ComplexFFT::prefersSplit;
	std::vector<Complex> tmpFreq, tmpTime;
	std::vector<Complex> twiddles, halfBinTwists;
};

template<typename Sample, bool splitComputation=false>
using ModifiedRealFFT = RealFFT<Sample, splitComputation, true>;

}} // namespace

// Override `Pow2FFT` / `Pow2RealFFT` templates with faster implementations
#if defined(SIGNALSMITH_USE_PFFFT) || defined(SIGNALSMITH_USE_PFFFT_DOUBLE)
#	if defined(SIGNALSMITH_USE_PFFFT)
#		include "./platform/fft-pffft.h"
#	endif
#	if defined(SIGNALSMITH_USE_PFFFT_DOUBLE)
#		include "./platform/fft-pffft-double.h"
#	endif
#elif defined(SIGNALSMITH_USE_ACCELERATE)
#	include "./platform/fft-accelerate.h"
#elif defined(SIGNALSMITH_USE_IPP)
#	include "./platform/fft-ipp.h"
#endif

#endif // include guard