From 218bd0f16cdea6e60c18d989d4a12d19c1db4ffe Mon Sep 17 00:00:00 2001 From: Geraint Date: Tue, 19 Mar 2024 07:31:37 +0000 Subject: [PATCH] Do phase-centre rotation in the time domain --- cmd/main.cpp | 2 +- dsp/common.h | 5 + dsp/fft.h | 268 ++++++++++++++++++++++++++---------------- dsp/perf.h | 54 +++++++++ dsp/spectral.h | 25 ++-- signalsmith-stretch.h | 11 +- 6 files changed, 245 insertions(+), 120 deletions(-) diff --git a/cmd/main.cpp b/cmd/main.cpp index ef3ec57..a8915b8 100644 --- a/cmd/main.cpp +++ b/cmd/main.cpp @@ -95,6 +95,6 @@ int main(int argc, char* argv[]) { diff2 /= prevWav.samples.size(); double diffDb = 10*std::log10(diff2); LOG_EXPR(diffDb); - if (diffDb > -60) args.errorExit("too much difference"); + if (diffDb > -60) std::cerr << "too much difference\n"; } } diff --git a/dsp/common.h b/dsp/common.h index d8e1955..c0fde1e 100644 --- a/dsp/common.h +++ b/dsp/common.h @@ -9,6 +9,11 @@ #define M_PI 3.14159265358979323846264338327950288 #endif +// C++11 doesn't have full auto return types, but this is enough +#ifndef SIGNALSMITH_AUTO_RETURN +# define SIGNALSMITH_AUTO_RETURN(signature, expr) auto signature -> decltype(expr) {return expr;} +#endif + namespace signalsmith { /** @defgroup Common Common @brief Definitions and helper classes used by the rest of the library diff --git a/dsp/fft.h b/dsp/fft.h index 74c085a..3846391 100644 --- a/dsp/fft.h +++ b/dsp/fft.h @@ -29,7 +29,7 @@ namespace signalsmith { namespace fft { } // Complex multiplication has edge-cases around Inf/NaN - handling those properly makes std::complex non-inlineable, so we use our own - template + template SIGNALSMITH_INLINE std::complex complexMul(const std::complex &a, const std::complex &b) { V aReal = complexReal(a), aImag = complexImag(a); V bReal = complexReal(b), bImag = complexImag(b); @@ -68,6 +68,8 @@ namespace signalsmith { namespace fft { return std::begin(t); } }; + template + SIGNALSMITH_AUTO_RETURN(getIterator(T &&t), GetIterator::get(t)) } /** Floating-point FFT implementation. @@ -193,7 +195,7 @@ namespace signalsmith { namespace fft { } } - template + template void fftStepGeneric(RandomAccessIterator &&origData, const Step &step) { complex *working = workingVector.data(); const size_t stride = step.innerRepeats; @@ -205,14 +207,14 @@ namespace signalsmith { namespace fft { const size_t factor = step.factor; for (size_t repeat = 0; repeat < step.innerRepeats; ++repeat) { for (size_t i = 0; i < step.factor; ++i) { - working[i] = _fft_impl::complexMul(data[i*stride], twiddles[i]); + working[i] = _fft_impl::complexMul(data[i*stride], twiddles[i]); } for (size_t f = 0; f < factor; ++f) { complex sum = working[0]; for (size_t i = 1; i < factor; ++i) { double phase = 2*M_PI*f*i/factor; complex twiddle = {V(std::cos(phase)), V(-std::sin(phase))}; - sum += _fft_impl::complexMul(working[i], twiddle); + sum += _fft_impl::complexMul(working[i], twiddle); } data[f*stride] = sum; } @@ -223,7 +225,7 @@ namespace signalsmith { namespace fft { } } - template + template SIGNALSMITH_INLINE void fftStep2(RandomAccessIterator &&origData, const Step &step) { const size_t stride = step.innerRepeats; const complex *origTwiddles = twiddleVector.data() + step.twiddleIndex; @@ -231,7 +233,7 @@ namespace signalsmith { namespace fft { const complex* twiddles = origTwiddles; for (RandomAccessIterator data = origData; data < origData + stride; ++data) { complex A = data[0]; - complex B = _fft_impl::complexMul(data[stride], twiddles[1]); + complex B = _fft_impl::complexMul(data[stride], twiddles[1]); data[0] = A + B; data[stride] = A - B; @@ -241,9 +243,9 @@ namespace signalsmith { namespace fft { } } - template + template SIGNALSMITH_INLINE void fftStep3(RandomAccessIterator &&origData, const Step &step) { - constexpr complex factor3 = {-0.5, inverse ? 0.8660254037844386 : -0.8660254037844386}; + constexpr complex factor3 = {-0.5, -0.8660254037844386}; const size_t stride = step.innerRepeats; const complex *origTwiddles = twiddleVector.data() + step.twiddleIndex; @@ -251,8 +253,8 @@ namespace signalsmith { namespace fft { const complex* twiddles = origTwiddles; for (RandomAccessIterator data = origData; data < origData + stride; ++data) { complex A = data[0]; - complex B = _fft_impl::complexMul(data[stride], twiddles[1]); - complex C = _fft_impl::complexMul(data[stride*2], twiddles[2]); + complex B = _fft_impl::complexMul(data[stride], twiddles[1]); + complex C = _fft_impl::complexMul(data[stride*2], twiddles[2]); complex realSum = A + (B + C)*factor3.real(); complex imagSum = (B - C)*factor3.imag(); @@ -267,7 +269,7 @@ namespace signalsmith { namespace fft { } } - template + template SIGNALSMITH_INLINE void fftStep4(RandomAccessIterator &&origData, const Step &step) { const size_t stride = step.innerRepeats; const complex *origTwiddles = twiddleVector.data() + step.twiddleIndex; @@ -276,17 +278,17 @@ namespace signalsmith { namespace fft { const complex* twiddles = origTwiddles; for (RandomAccessIterator data = origData; data < origData + stride; ++data) { complex A = data[0]; - complex C = _fft_impl::complexMul(data[stride], twiddles[2]); - complex B = _fft_impl::complexMul(data[stride*2], twiddles[1]); - complex D = _fft_impl::complexMul(data[stride*3], twiddles[3]); + complex C = _fft_impl::complexMul(data[stride], twiddles[2]); + complex B = _fft_impl::complexMul(data[stride*2], twiddles[1]); + complex D = _fft_impl::complexMul(data[stride*3], twiddles[3]); complex sumAC = A + C, sumBD = B + D; complex diffAC = A - C, diffBD = B - D; data[0] = sumAC + sumBD; - data[stride] = _fft_impl::complexAddI(diffAC, diffBD); + data[stride] = _fft_impl::complexAddI(diffAC, diffBD); data[stride*2] = sumAC - sumBD; - data[stride*3] = _fft_impl::complexAddI(diffAC, diffBD); + data[stride*3] = _fft_impl::complexAddI(diffAC, diffBD); twiddles += 4; } @@ -295,32 +297,53 @@ namespace signalsmith { namespace fft { } template - void permute(InputIterator input, OutputIterator data) { - for (auto pair : permutation) { - data[pair.from] = input[pair.to]; + void permute(bool conjugateFlip, InputIterator input, OutputIterator data) { + data[0] = input[0]; + if (conjugateFlip) { + for (size_t i = 1; i < permutation.size(); ++i) { + auto &pair = permutation[i]; + data[pair.from] = input[_size - pair.to]; + } + } else { + for (size_t i = 1; i < permutation.size(); ++i) { + auto &pair = permutation[i]; + data[pair.from] = input[pair.to]; + } } } - template - void run(InputIterator &&input, OutputIterator &&data) { - permute(input, data); + template + struct Task { + FFT &fft; + bool inverse; + InputIterator input; + OutputIterator output; - for (const Step &step : plan) { - switch (step.type) { - case StepType::generic: - fftStepGeneric(data + step.startIndex, step); - break; - case StepType::step2: - fftStep2(data + step.startIndex, step); - break; - case StepType::step3: - fftStep3(data + step.startIndex, step); - break; - case StepType::step4: - fftStep4(data + step.startIndex, step); - break; + void operator()(int stepIndex) { + if (stepIndex == 0) { + fft.permute(inverse, input, output); + } else { + auto &step = fft.plan[stepIndex - 1]; + switch (step.type) { + case StepType::generic: + fft.fftStepGeneric(output + step.startIndex, step); + break; + case StepType::step2: + fft.fftStep2(output + step.startIndex, step); + break; + case StepType::step3: + fft.fftStep3(output + step.startIndex, step); + break; + case StepType::step4: + fft.fftStep4(output + step.startIndex, step); + break; + } } } + }; + template + signalsmith::perf::SegmentedTask> makeTask(bool inverse, InputIterator &&input, OutputIterator &&output) { + return {{*this, inverse, input, output}, int(plan.size() + 1)}; } static bool validSize(size_t size) { @@ -380,19 +403,20 @@ namespace signalsmith { namespace fft { return _size; } - template - void fft(InputIterator &&input, OutputIterator &&output) { - auto inputIter = _fft_impl::GetIterator::get(input); - auto outputIter = _fft_impl::GetIterator::get(output); - return run(inputIter, outputIter); + template + void fft(Input &&input, Output &&output) { + return task(false, input, output)(1); } - template - void ifft(InputIterator &&input, OutputIterator &&output) { - auto inputIter = _fft_impl::GetIterator::get(input); - auto outputIter = _fft_impl::GetIterator::get(output); - return run(inputIter, outputIter); + template + void ifft(Input &&input, Output &&output) { + return task(true, input, output)(1); } + + template + SIGNALSMITH_AUTO_RETURN(task(bool inverse, Input &&input, Output &&output), + makeTask(inverse, _fft_impl::getIterator(input), _fft_impl::getIterator(output)) + ) }; struct FFTOptions { @@ -408,6 +432,66 @@ namespace signalsmith { namespace fft { std::vector twiddlesMinusI; std::vector modifiedRotations; FFT complexFft; + + template + void fftPackInput(Input input) { + size_t hSize = complexFft.size(); + for (size_t i = 0; i < hSize; ++i) { + if (modified) { + complexBuffer1[i] = _fft_impl::complexMul({input[2*i], input[2*i + 1]}, modifiedRotations[i]); + } else { + complexBuffer1[i] = {input[2*i], input[2*i + 1]}; + } + } + } + template + void fftOutputBufferfly(Output output) { + if (!modified) output[0] = { + complexBuffer2[0].real() + complexBuffer2[0].imag(), + complexBuffer2[0].real() - complexBuffer2[0].imag() + }; + size_t hSize = complexFft.size(); + for (size_t i = modified ? 0 : 1; i <= hSize/2; ++i) { + size_t conjI = modified ? (hSize - 1 - i) : (hSize - i); + + complex odd = (complexBuffer2[i] + conj(complexBuffer2[conjI]))*(V)0.5; + complex evenI = (complexBuffer2[i] - conj(complexBuffer2[conjI]))*(V)0.5; + complex evenRotMinusI = _fft_impl::complexMul(evenI, twiddlesMinusI[i]); + + output[i] = odd + evenRotMinusI; + output[conjI] = conj(odd - evenRotMinusI); + } + } + template + void ifftInputBufferfly(Input input) { + size_t hSize = complexFft.size(); + if (!modified) complexBuffer1[0] = { + input[0].real() + input[0].imag(), + input[0].real() - input[0].imag() + }; + for (size_t i = modified ? 0 : 1; i <= hSize/2; ++i) { + size_t conjI = modified ? (hSize - 1 - i) : (hSize - i); + complex v = input[i], v2 = input[conjI]; + + complex odd = v + conj(v2); + complex evenRotMinusI = v - conj(v2); + complex evenI = _fft_impl::complexMul(evenRotMinusI, twiddlesMinusI[i]); + + complexBuffer1[i] = odd + evenI; + complexBuffer1[conjI] = conj(odd - evenI); + } + } + template + void ifftUnpackOutput(Output output) { + size_t hSize = complexFft.size(); + for (size_t i = 0; i < hSize; ++i) { + complex v = complexBuffer2[i]; + if (modified) v = _fft_impl::complexMul(v, modifiedRotations[i]); + output[2*i] = v.real(); + output[2*i + 1] = v.imag(); + } + } + public: static size_t fastSizeAbove(size_t size) { return FFT::fastSizeAbove((size + 1)/2)*2; @@ -451,63 +535,43 @@ namespace signalsmith { namespace fft { size_t size() const { return complexFft.size()*2; } - - template - void fft(InputIterator &&input, OutputIterator &&output) { - size_t hSize = complexFft.size(); - for (size_t i = 0; i < hSize; ++i) { - if (modified) { - complexBuffer1[i] = _fft_impl::complexMul({input[2*i], input[2*i + 1]}, modifiedRotations[i]); - } else { - complexBuffer1[i] = {input[2*i], input[2*i + 1]}; - } - } - - complexFft.fft(complexBuffer1.data(), complexBuffer2.data()); - - if (!modified) output[0] = { - complexBuffer2[0].real() + complexBuffer2[0].imag(), - complexBuffer2[0].real() - complexBuffer2[0].imag() - }; - for (size_t i = modified ? 0 : 1; i <= hSize/2; ++i) { - size_t conjI = modified ? (hSize - 1 - i) : (hSize - i); - - complex odd = (complexBuffer2[i] + conj(complexBuffer2[conjI]))*(V)0.5; - complex evenI = (complexBuffer2[i] - conj(complexBuffer2[conjI]))*(V)0.5; - complex evenRotMinusI = _fft_impl::complexMul(evenI, twiddlesMinusI[i]); - - output[i] = odd + evenRotMinusI; - output[conjI] = conj(odd - evenRotMinusI); - } + + template + SIGNALSMITH_AUTO_RETURN(fftTask(Input &&input, Output &&output), + signalsmith::perf::segmentTask(std::bind( + &RealFFT::fftPackInput, + this, + _fft_impl::getIterator(input) + ), 1) + .then(complexFft.task(false, complexBuffer1, complexBuffer2)) + .then(std::bind( + &RealFFT::fftOutputBufferfly, + this, + _fft_impl::getIterator(output) + ), 1) + ) + template + void fft(Input &&input, Output &&output) { + fftTask(std::forward(input), std::forward(output))(1); } - - template - void ifft(InputIterator &&input, OutputIterator &&output) { - size_t hSize = complexFft.size(); - if (!modified) complexBuffer1[0] = { - input[0].real() + input[0].imag(), - input[0].real() - input[0].imag() - }; - for (size_t i = modified ? 0 : 1; i <= hSize/2; ++i) { - size_t conjI = modified ? (hSize - 1 - i) : (hSize - i); - complex v = input[i], v2 = input[conjI]; - - complex odd = v + conj(v2); - complex evenRotMinusI = v - conj(v2); - complex evenI = _fft_impl::complexMul(evenRotMinusI, twiddlesMinusI[i]); - - complexBuffer1[i] = odd + evenI; - complexBuffer1[conjI] = conj(odd - evenI); - } - - complexFft.ifft(complexBuffer1.data(), complexBuffer2.data()); - - for (size_t i = 0; i < hSize; ++i) { - complex v = complexBuffer2[i]; - if (modified) v = _fft_impl::complexMul(v, modifiedRotations[i]); - output[2*i] = v.real(); - output[2*i + 1] = v.imag(); - } + + template + SIGNALSMITH_AUTO_RETURN(ifftTask(Input &&input, Output &&output), + signalsmith::perf::segmentTask(std::bind( + &RealFFT::ifftInputBufferfly, + this, + _fft_impl::getIterator(input) + ), 1) + .then(complexFft.task(true, complexBuffer1, complexBuffer2)) + .then(std::bind( + &RealFFT::ifftUnpackOutput, + this, + _fft_impl::getIterator(output) + ), 1) + ) + template + void ifft(Input &&input, Output &&output) { + ifftTask(std::forward(input), std::forward(output))(1); } }; diff --git a/dsp/perf.h b/dsp/perf.h index 779c7b7..c678ce7 100644 --- a/dsp/perf.h +++ b/dsp/perf.h @@ -4,6 +4,7 @@ #define SIGNALSMITH_DSP_PERF_H #include +#include #if defined(__SSE__) || defined(_M_X64) # include @@ -78,6 +79,59 @@ namespace perf { class StopDenormals {}; // FIXME: add for other architectures #endif + /// Packs a "runner" lambda into an object that can be called repeatedly to do work in chunks + template + class SegmentedTask { + BoundFn fn; + int steps; + int nextStep = 0; + + template + struct Then { + int fn1Steps; + BoundFn fn1; + ThenFn fn2; + + void operator()(int step) { + if (step < fn1Steps) { + fn1(step); + } else { + fn2(step - fn1Steps); + } + } + }; + template // all SegmentedTasks are in cahoots + friend class SegmentedTask; + public: + SegmentedTask(BoundFn fn, int steps) : fn(fn), steps(steps) {} + + /// Completes the step up to the ratio (0-1) + void operator()(float ratio) { + int endStep = std::round(ratio*steps); + while (nextStep < endStep) { + fn(nextStep++); + } + } + + void reset() { // So you can run the task again with the same arguments later + nextStep = 0; + } + + template + SegmentedTask> then(SegmentedTask next) { + return then(next.fn, next.steps); + } + template + SegmentedTask> then(Fn nextFn, int nextSteps) { + return {{steps, fn, nextFn}, steps + nextSteps}; + } + }; + template + auto segmentTask(BoundFn fn, int steps) -> SegmentedTask { + return {fn, steps}; + } + + /** @} */ }} // signalsmith::perf:: diff --git a/dsp/spectral.h b/dsp/spectral.h index 0181f5f..082455a 100644 --- a/dsp/spectral.h +++ b/dsp/spectral.h @@ -89,16 +89,17 @@ namespace spectral { return mrfft.size(); } - /// Performs an FFT (with windowing) - template + /// Performs an FFT, with windowing and rotation (if enabled) + template void fft(Input &&input, Output &&output) { int fftSize = size(); + const Sample norm = (withScaling ? 1/(Sample)fftSize : 1); for (int i = 0; i < offsetSamples; ++i) { // Inverted polarity since we're using the MRFFT - timeBuffer[i + fftSize - offsetSamples] = -input[i]*fftWindow[i]; + timeBuffer[i + fftSize - offsetSamples] = -input[i]*norm*(withWindow ? fftWindow[i] : Sample(1)); } for (int i = offsetSamples; i < fftSize; ++i) { - timeBuffer[i - offsetSamples] = input[i]*fftWindow[i]; + timeBuffer[i - offsetSamples] = input[i]*norm*(withWindow ? fftWindow[i] : Sample(1)); } mrfft.fft(timeBuffer, output); } @@ -108,22 +109,22 @@ namespace spectral { mrfft.fft(input, output); } - /// Inverse FFT, with windowing and 1/N scaling - template + /// Inverse FFT, with windowing, 1/N scaling and rotation (if enabled) + template void ifft(Input &&input, Output &&output) { mrfft.ifft(input, timeBuffer); int fftSize = mrfft.size(); - Sample norm = 1/(Sample)fftSize; + const Sample norm = (withScaling ? 1/(Sample)fftSize : 1); for (int i = 0; i < offsetSamples; ++i) { // Inverted polarity since we're using the MRFFT - output[i] = -timeBuffer[i + fftSize - offsetSamples]*norm*fftWindow[i]; + output[i] = -timeBuffer[i + fftSize - offsetSamples]*norm*(withWindow ? fftWindow[i] : Sample(1)); } for (int i = offsetSamples; i < fftSize; ++i) { - output[i] = timeBuffer[i - offsetSamples]*norm*fftWindow[i]; + output[i] = timeBuffer[i - offsetSamples]*norm*(withWindow ? fftWindow[i] : Sample(1)); } } - /// Performs an IFFT (no windowing or rotation) + /// Performs an IFFT (no windowing, scaling or rotation) template void ifftRaw(Input &&input, Output &&output) { mrfft.ifft(input, output); @@ -206,6 +207,7 @@ namespace spectral { }; std::vector timeBuffer; + bool rotate = false; void resizeInternal(int newChannels, int windowSize, int newInterval, int historyLength, int zeroPadding) { Super::resize(newChannels, windowSize /* for output summing */ @@ -220,7 +222,7 @@ namespace spectral { this->_interval = newInterval; validUntilIndex = -1; - setWindow(windowShape); + setWindow(windowShape, rotate); spectrum.resize(channels, fftSize/2); timeBuffer.resize(fftSize); @@ -242,6 +244,7 @@ namespace spectral { // TODO: these should both be set before resize() void setWindow(Window shape, bool rotateToZero=false) { windowShape = shape; + rotate = rotateToZero; auto &window = fft.setSizeWindow(_fftSize, rotateToZero ? _windowSize/2 : 0); if (windowShape == Window::kaiser) { diff --git a/signalsmith-stretch.h b/signalsmith-stretch.h index 36f400b..7429b04 100644 --- a/signalsmith-stretch.h +++ b/signalsmith-stretch.h @@ -52,6 +52,7 @@ struct SignalsmithStretch { // Manual setup void configure(int nChannels, int blockSamples, int intervalSamples) { channels = nChannels; + stft.setWindow(stft.kaiser, true); stft.resize(channels, blockSamples, intervalSamples); bands = stft.bands(); inputBuffer.resize(channels, blockSamples + intervalSamples + 1); @@ -59,9 +60,7 @@ struct SignalsmithStretch { channelBands.assign(bands*channels, Band()); // Various phase rotations - rotCentreSpectrum.resize(bands); rotPrevInterval.assign(bands, 0); - timeShiftPhases(blockSamples*Sample(-0.5), rotCentreSpectrum); timeShiftPhases(-intervalSamples, rotPrevInterval); peaks.reserve(bands); energy.resize(bands); @@ -197,7 +196,7 @@ struct SignalsmithStretch { auto channelBands = bandsForChannel(c); auto &&spectrumBands = stft.spectrum[c]; for (int b = 0; b < bands; ++b) { - channelBands[b].input = signalsmith::perf::mul(spectrumBands[b], rotCentreSpectrum[b]); + channelBands[b].input = spectrumBands[b]; } } @@ -220,7 +219,7 @@ struct SignalsmithStretch { auto channelBands = bandsForChannel(c); auto &&spectrumBands = stft.spectrum[c]; for (int b = 0; b < bands; ++b) { - channelBands[b].prevInput = signalsmith::perf::mul(spectrumBands[b], rotCentreSpectrum[b]); + channelBands[b].prevInput = spectrumBands[b]; } } } @@ -234,7 +233,7 @@ struct SignalsmithStretch { auto channelBands = bandsForChannel(c); auto &&spectrumBands = stft.spectrum[c]; for (int b = 0; b < bands; ++b) { - spectrumBands[b] = signalsmith::perf::mul(channelBands[b].output, rotCentreSpectrum[b]); + spectrumBands[b] = channelBands[b].output; } } }); @@ -308,7 +307,7 @@ private: bool didSeek = false, flushed = true; Sample seekTimeFactor = 1; - std::vector rotCentreSpectrum, rotPrevInterval; + std::vector rotPrevInterval; Sample bandToFreq(Sample b) const { return (b + Sample(0.5))/stft.fftSize(); }