Signalsmith-Audio/signalsmith-stretch
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Do phase-centre rotation in the time domain

Browse Source
This commit is contained in:
Geraint 2024-03-19 07:31:37 +00:00 committed by Geraint Luff
parent 026622300e
commit 218bd0f16c
6 changed files with 245 additions and 120 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
cmd/main.cpp
View File

@ -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";
}
}

5
dsp/common.h
View File

@ -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

268
dsp/fft.h
View File

@ -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 <bool conjugateSecond, typename V>
template <bool conjugateSecond=false, typename V>
SIGNALSMITH_INLINE std::complex<V> complexMul(const std::complex<V> &a, const std::complex<V> &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<typename T>
SIGNALSMITH_AUTO_RETURN(getIterator(T &&t), GetIterator<T>::get(t))
}
/** Floating-point FFT implementation.
@ -193,7 +195,7 @@ namespace signalsmith { namespace fft {
}
}
template<bool inverse, typename RandomAccessIterator>
template<typename RandomAccessIterator>
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<inverse>(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<inverse>(working[i], twiddle);
sum += _fft_impl::complexMul(working[i], twiddle);
}
data[f*stride] = sum;
}
@ -223,7 +225,7 @@ namespace signalsmith { namespace fft {
}
}
template<bool inverse, typename RandomAccessIterator>
template<typename RandomAccessIterator>
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<inverse>(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<bool inverse, typename RandomAccessIterator>
template<typename RandomAccessIterator>
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<inverse>(data[stride], twiddles[1]);
complex C = _fft_impl::complexMul<inverse>(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<bool inverse, typename RandomAccessIterator>
template<typename RandomAccessIterator>
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<inverse>(data[stride], twiddles[2]);
complex B = _fft_impl::complexMul<inverse>(data[stride*2], twiddles[1]);
complex D = _fft_impl::complexMul<inverse>(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<!inverse>(diffAC, diffBD);
data[stride] = _fft_impl::complexAddI<true>(diffAC, diffBD);
data[stride*2] = sumAC - sumBD;
data[stride*3] = _fft_impl::complexAddI<inverse>(diffAC, diffBD);
data[stride*3] = _fft_impl::complexAddI<false>(diffAC, diffBD);
twiddles += 4;
}
@ -295,32 +297,53 @@ namespace signalsmith { namespace fft {
}
template<typename InputIterator, typename OutputIterator>
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<bool inverse, typename InputIterator, typename OutputIterator>
void run(InputIterator &&input, OutputIterator &&data) {
permute(input, data);
template<typename InputIterator, typename OutputIterator>
struct Task {
FFT &fft;
bool inverse;
InputIterator input;
OutputIterator output;
for (const Step &step : plan) {
switch (step.type) {
case StepType::generic:
fftStepGeneric<inverse>(data + step.startIndex, step);
break;
case StepType::step2:
fftStep2<inverse>(data + step.startIndex, step);
break;
case StepType::step3:
fftStep3<inverse>(data + step.startIndex, step);
break;
case StepType::step4:
fftStep4<inverse>(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<typename InputIterator, typename OutputIterator>
signalsmith::perf::SegmentedTask<Task<InputIterator, OutputIterator>> 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<typename InputIterator, typename OutputIterator>
void fft(InputIterator &&input, OutputIterator &&output) {
auto inputIter = _fft_impl::GetIterator<InputIterator>::get(input);
auto outputIter = _fft_impl::GetIterator<OutputIterator>::get(output);
return run<false>(inputIter, outputIter);
template<typename Input, typename Output>
void fft(Input &&input, Output &&output) {
return task(false, input, output)(1);
}
template<typename InputIterator, typename OutputIterator>
void ifft(InputIterator &&input, OutputIterator &&output) {
auto inputIter = _fft_impl::GetIterator<InputIterator>::get(input);
auto outputIter = _fft_impl::GetIterator<OutputIterator>::get(output);
return run<true>(inputIter, outputIter);
template<typename Input, typename Output>
void ifft(Input &&input, Output &&output) {
return task(true, input, output)(1);
}
template<typename Input, typename Output>
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<complex> twiddlesMinusI;
std::vector<complex> modifiedRotations;
FFT<V> complexFft;
template<typename Input>
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<typename Output>
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<typename Input>
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<true>(evenRotMinusI, twiddlesMinusI[i]);
complexBuffer1[i] = odd + evenI;
complexBuffer1[conjI] = conj(odd - evenI);
}
}
template<typename Output>
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<true>(v, modifiedRotations[i]);
output[2*i] = v.real();
output[2*i + 1] = v.imag();
}
}
public:
static size_t fastSizeAbove(size_t size) {
return FFT<V>::fastSizeAbove((size + 1)/2)*2;
@ -451,63 +535,43 @@ namespace signalsmith { namespace fft {
size_t size() const {
return complexFft.size()*2;
}
template<typename InputIterator, typename OutputIterator>
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<false>({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<false>(evenI, twiddlesMinusI[i]);
output[i] = odd + evenRotMinusI;
output[conjI] = conj(odd - evenRotMinusI);
}
template<typename Input, typename Output>
SIGNALSMITH_AUTO_RETURN(fftTask(Input &&input, Output &&output),
signalsmith::perf::segmentTask(std::bind(
&RealFFT::fftPackInput<decltype(_fft_impl::getIterator(input))>,
this,
_fft_impl::getIterator(input)
), 1)
.then(complexFft.task(false, complexBuffer1, complexBuffer2))
.then(std::bind(
&RealFFT::fftOutputBufferfly<decltype(_fft_impl::getIterator(output))>,
this,
_fft_impl::getIterator(output)
), 1)
)
template<typename Input, typename Output>
void fft(Input &&input, Output &&output) {
fftTask(std::forward<Input>(input), std::forward<Output>(output))(1);
}
template<typename InputIterator, typename OutputIterator>
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<true>(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<true>(v, modifiedRotations[i]);
output[2*i] = v.real();
output[2*i + 1] = v.imag();
}
template<typename Input, typename Output>
SIGNALSMITH_AUTO_RETURN(ifftTask(Input &&input, Output &&output),
signalsmith::perf::segmentTask(std::bind(
&RealFFT::ifftInputBufferfly<decltype(_fft_impl::getIterator(input))>,
this,
_fft_impl::getIterator(input)
), 1)
.then(complexFft.task(true, complexBuffer1, complexBuffer2))
.then(std::bind(
&RealFFT::ifftUnpackOutput<decltype(_fft_impl::getIterator(output))>,
this,
_fft_impl::getIterator(output)
), 1)
)
template<typename Input, typename Output>
void ifft(Input &&input, Output &&output) {
ifftTask(std::forward<Input>(input), std::forward<Output>(output))(1);
}
};

54
dsp/perf.h
View File

@ -4,6 +4,7 @@
#define SIGNALSMITH_DSP_PERF_H
#include <complex>
#include <functional>
#if defined(__SSE__) || defined(_M_X64)
# include <xmmintrin.h>
@ -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 BoundFn>
class SegmentedTask {
BoundFn fn;
int steps;
int nextStep = 0;
template<class ThenFn>
struct Then {
int fn1Steps;
BoundFn fn1;
ThenFn fn2;
void operator()(int step) {
if (step < fn1Steps) {
fn1(step);
} else {
fn2(step - fn1Steps);
}
}
};
template<typename Fn2> // 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<class Fn>
SegmentedTask<Then<Fn>> then(SegmentedTask<Fn> next) {
return then(next.fn, next.steps);
}
template<class Fn>
SegmentedTask<Then<Fn>> then(Fn nextFn, int nextSteps) {
return {{steps, fn, nextFn}, steps + nextSteps};
}
};
template<class BoundFn>
auto segmentTask(BoundFn fn, int steps) -> SegmentedTask<BoundFn> {
return {fn, steps};
}
/** @} */
}} // signalsmith::perf::

25
dsp/spectral.h
View File

@ -89,16 +89,17 @@ namespace spectral {
return mrfft.size();
}
/// Performs an FFT (with windowing)
template<class Input, class Output>
/// Performs an FFT, with windowing and rotation (if enabled)
template<bool withWindow=true, bool withScaling=false, class Input, class Output>
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<class Input, class Output>
/// Inverse FFT, with windowing, 1/N scaling and rotation (if enabled)
template<bool withWindow=true, bool withScaling=true, class Input, class Output>
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<class Input, class Output>
void ifftRaw(Input &&input, Output &&output) {
mrfft.ifft(input, output);
@ -206,6 +207,7 @@ namespace spectral {
};
std::vector<Sample> 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) {

11
signalsmith-stretch.h
View File

@ -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<true>(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<Complex> rotCentreSpectrum, rotPrevInterval;
std::vector<Complex> rotPrevInterval;
Sample bandToFreq(Sample b) const {
return (b + Sample(0.5))/stft.fftSize();
}