From f33f1c4b02b2f974ea1531189250208eada22ffd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ju=CC=88rgen=20Hock?= Date: Sun, 12 Nov 2023 10:36:58 +0100 Subject: [PATCH] Remove SDFT prototype #5 --- cpp/StftPitchShift/LibStftPitchShift.cmake | 1 - cpp/StftPitchShift/SDFT.h | 226 --------------------- cpp/StftPitchShift/StftPitchShift.cpp | 45 +--- cpp/StftPitchShift/main.cpp | 2 +- cpp/StftPitchShift/wasm.cpp | 2 +- python/stftpitchshift/main.py | 10 - python/stftpitchshift/sdft.py | 96 --------- 7 files changed, 12 insertions(+), 370 deletions(-) delete mode 100644 cpp/StftPitchShift/SDFT.h delete mode 100644 python/stftpitchshift/sdft.py diff --git a/cpp/StftPitchShift/LibStftPitchShift.cmake b/cpp/StftPitchShift/LibStftPitchShift.cmake index e6b97d1..24dcf37 100644 --- a/cpp/StftPitchShift/LibStftPitchShift.cmake +++ b/cpp/StftPitchShift/LibStftPitchShift.cmake @@ -37,7 +37,6 @@ set(HEADERS "${CMAKE_CURRENT_LIST_DIR}/Pitcher.h" "${CMAKE_CURRENT_LIST_DIR}/Resampler.h" "${CMAKE_CURRENT_LIST_DIR}/RFFT.h" - "${CMAKE_CURRENT_LIST_DIR}/SDFT.h" "${CMAKE_CURRENT_LIST_DIR}/STFT.h" "${CMAKE_CURRENT_LIST_DIR}/StftPitchShift.h" "${CMAKE_CURRENT_LIST_DIR}/StftPitchShiftCore.h" diff --git a/cpp/StftPitchShift/SDFT.h b/cpp/StftPitchShift/SDFT.h deleted file mode 100644 index bbdd8f3..0000000 --- a/cpp/StftPitchShift/SDFT.h +++ /dev/null @@ -1,226 +0,0 @@ -#pragma once - -#include -#include -#include -#include -#include -#include - -#include - -namespace stftpitchshift -{ - template - class SDFT - { - - public: - - SDFT(const size_t framesize, const double latency, const bool chronometry = false) : - framesize(framesize), - dftsize(framesize / 2), - latency(1), - chronometry(chronometry) - { - analysis.roi = { 0, dftsize }; - synthesis.roi = { 0, dftsize }; - - analysis.twiddles.resize(dftsize); - synthesis.twiddles.resize(dftsize); - - analysis.cursor = 0; - analysis.input.resize(dftsize * 2); - analysis.accoutput.resize(dftsize); - analysis.auxoutput.resize(dftsize + 2); - analysis.fiddles.resize(dftsize, 1); - - const double pi = -2.0 * std::acos(-1.0) / (dftsize * 2); - const double weight = 2.0 / (1.0 - std::cos(pi * dftsize * latency)); - - for (size_t i = 0; i < dftsize; ++i) - { - analysis.twiddles[i] = std::polar(1.0, pi * i); - synthesis.twiddles[i] = std::polar(weight, pi * i * dftsize * latency); - } - } - - void operator()(const std::vector& input, std::vector& output, const std::function>& dft)> callback) - { - (*this)(input.size(), input.data(), output.data(), callback); - } - - void operator()(const size_t size, const T* input, T* const output, const std::function>& dft)> callback) - { - std::vector> dft(framesize / 2 + 1); - - if (chronometry) - { - struct - { - Timer analysis; - Timer synthesis; - Timer callback; - Timer loop; - } - timers; - - timers.loop.tic(); - for (size_t i = 0; i < size; ++i) - { - timers.analysis.tic(); - sdft(input[i], dft); - timers.analysis.toc(); - - timers.callback.tic(); - callback(dft); - timers.callback.toc(); - - timers.synthesis.tic(); - output[i] = isdft(dft); - timers.synthesis.toc(); - } - timers.loop.toc(); - - std::cout << "analysis " << timers.analysis.str() << std::endl; - std::cout << "synthesis " << timers.synthesis.str() << std::endl; - std::cout << "callback " << timers.callback.str() << std::endl; - std::cout << "loop " << timers.loop.str() << std::endl; - } - else - { - for (size_t i = 0; i < size; ++i) - { - sdft(input[i], dft); - callback(dft); - output[i] = isdft(dft); - } - } - } - - private: - - const size_t framesize; - const size_t dftsize; - const double latency; - const bool chronometry; - - struct - { - std::pair roi; - std::vector> twiddles; - - size_t cursor; - std::vector input; - std::vector> accoutput; - std::vector> auxoutput; - std::vector> fiddles; - } - analysis; - - struct - { - std::pair roi; - std::vector> twiddles; - } - synthesis; - - inline static std::complex window(const std::complex& left, - const std::complex& middle, - const std::complex& right, - const double weight) - { - // the factor 1/4 is already included in the weight - - const std::complex value = - /* 0.25 */ ((middle + middle) - (left + right)) * weight; - - return std::complex( - static_cast(value.real()), - static_cast(value.imag())); - } - - inline static T exchange(T& oldvalue, const T newvalue) - { - const T oldcopy = oldvalue; - - oldvalue = newvalue; - - return oldcopy; - } - - void sdft(const T sample, std::vector>& dft) - { - // actually the weight denominator needs to be dftsize*2 to get proper magnitude scaling, - // but then requires a correction by factor 2 in synthesis and is therefore omitted - - const double weight = 0.25 / dftsize; // incl. factor 1/4 for windowing - - const double delta = sample - exchange(analysis.input[analysis.cursor], sample); - - for (size_t i = analysis.roi.first, j = i + 1; i < analysis.roi.second; ++i, ++j) - { - const std::complex oldfiddle = analysis.fiddles[i]; - const std::complex newfiddle = oldfiddle * analysis.twiddles[i]; - - analysis.fiddles[i] = newfiddle; - - analysis.accoutput[i] += delta * oldfiddle; - analysis.auxoutput[j] = analysis.accoutput[i] * std::conj(newfiddle); - } - - // theoretically the DFT periodicity needs to be preserved for proper windowing, - // but the both outer bins seem to be noisy for an unclear reason - // and will be suppressed anyway after windowing - - // analysis.auxoutput[0] = analysis.auxoutput[dftsize]; - // analysis.auxoutput[dftsize + 1] = analysis.auxoutput[1]; - - for (size_t i = analysis.roi.first, j = i + 1; i < analysis.roi.second; ++i, ++j) - { - dft[i] = window(analysis.auxoutput[j - 1], - analysis.auxoutput[j], - analysis.auxoutput[j + 1], - weight); - } - - // finally suppress outer DFT bins as announced in the comment above - - dft[0] = dft[dftsize - 1] = 0; - - if (++analysis.cursor >= analysis.input.size()) - { - analysis.cursor = 0; - - // TODO std::fill(analysis.fiddles.begin(), analysis.fiddles.end(), 1); - } - } - - T isdft(const std::vector>& dft) - { - double sample = 0; - - if (latency == 1) - { - for (size_t i = synthesis.roi.first; i < synthesis.roi.second; ++i) - { - const double value = dft[i].real(); - - sample += value * (i % 2 ? -1 : +1); - } - } - else - { - for (size_t i = synthesis.roi.first; i < synthesis.roi.second; ++i) - { - const std::complex value = { dft[i].real(), dft[i].imag() }; - - sample += (value * synthesis.twiddles[i]).real(); - } - } - - return static_cast(sample); - } - - }; -} diff --git a/cpp/StftPitchShift/StftPitchShift.cpp b/cpp/StftPitchShift/StftPitchShift.cpp index b55e769..de7bfa3 100644 --- a/cpp/StftPitchShift/StftPitchShift.cpp +++ b/cpp/StftPitchShift/StftPitchShift.cpp @@ -3,7 +3,6 @@ #include #include -#include #include #include @@ -162,24 +161,12 @@ void StftPitchShift::shiftpitch( core.distortion(distortion); core.normalization(normalization); - if (hopsize == 1) - { - SDFT sdft(framesize, /* latency */ 1, chronometry); - - sdft(size, input, output, [&](std::vector>& dft) - { - core.shiftpitch(dft); - }); - } - else - { - STFT stft(fft, framesize, hopsize, chronometry); + STFT stft(fft, framesize, hopsize, chronometry); - stft(size, input, output, [&](std::vector>& dft) - { - core.shiftpitch(dft); - }); - } + stft(size, input, output, [&](std::vector>& dft) + { + core.shiftpitch(dft); + }); } void StftPitchShift::shiftpitch( @@ -197,22 +184,10 @@ void StftPitchShift::shiftpitch( core.distortion(distortion); core.normalization(normalization); - if (hopsize == 1) - { - SDFT sdft(framesize, /* latency */ 1, chronometry); - - sdft(size, input, output, [&](std::vector>& dft) - { - core.shiftpitch(dft); - }); - } - else - { - STFT stft(fft, framesize, hopsize, chronometry); + STFT stft(fft, framesize, hopsize, chronometry); - stft(size, input, output, [&](std::vector>& dft) - { - core.shiftpitch(dft); - }); - } + stft(size, input, output, [&](std::vector>& dft) + { + core.shiftpitch(dft); + }); } diff --git a/cpp/StftPitchShift/main.cpp b/cpp/StftPitchShift/main.cpp index 780c932..6ce65b9 100644 --- a/cpp/StftPitchShift/main.cpp +++ b/cpp/StftPitchShift/main.cpp @@ -66,7 +66,7 @@ int main(int argc, char** argv) StftPitchShift stft( cli.framesize, - cli.hoprate ? cli.framesize / cli.hoprate : 1, + cli.framesize / cli.hoprate, samplerate, cli.normalization, cli.chronometry); diff --git a/cpp/StftPitchShift/wasm.cpp b/cpp/StftPitchShift/wasm.cpp index 06bdd1c..d96fc0c 100644 --- a/cpp/StftPitchShift/wasm.cpp +++ b/cpp/StftPitchShift/wasm.cpp @@ -119,7 +119,7 @@ extern "C" { StftPitchShift stft( cli.framesize, - cli.hoprate ? cli.framesize / cli.hoprate : 1, + cli.framesize / cli.hoprate, samplerate, cli.normalization, cli.chronometry); diff --git a/python/stftpitchshift/main.py b/python/stftpitchshift/main.py index 20aeb6d..d9a3742 100644 --- a/python/stftpitchshift/main.py +++ b/python/stftpitchshift/main.py @@ -2,7 +2,6 @@ from stftpitchshift import __version__ as version from stftpitchshift.io import read, write from stftpitchshift.plot import spectrogram -from stftpitchshift.sdft import sdft from stftpitchshift.stft import stft import click @@ -55,18 +54,9 @@ def cent(value): return pow(2, float(re.match('([+,-]?\\d+){1}([+,-]\\d+){0,1}', for channel in range(channels): - # STFT framesX = stft(x[:, channel], framesize, hopsize) framesY = stft(y[:, channel], framesize, hopsize) - # SDFT - # framesX = sdft(x[:, channel], framesize // 2) - # framesY = sdft(y[:, channel], framesize // 2) - - # TEST - # framesX = stft(x[:, channel], framesize, hopsize) - # framesY = sdft(x[:, channel], framesize // 2) - figure = plot.figure(f'Channel {channel+1}/{channels}') spectrogramX = figure.add_subplot(2, 1, 1, title='Input Spectrogram') diff --git a/python/stftpitchshift/sdft.py b/python/stftpitchshift/sdft.py deleted file mode 100644 index 60945a5..0000000 --- a/python/stftpitchshift/sdft.py +++ /dev/null @@ -1,96 +0,0 @@ -import numpy as np - - -def window(x): - - M, N = x.shape - - left = x[:, :-2] - right = x[:, +2:] - middle = x[:, +1:-1] - - y = ((middle + middle) - (left + right)) / (N * 4) - - y = np.pad(y, ((0, 0), (1, 1))) - - return y - - -def sdft(samples, dftsize): - - M = samples.size - N = dftsize - - m = np.arange(M + 1)[:, None] - n = np.arange(N) - - twiddles = np.exp(-2j * np.pi * m * n / (N * 2)) - - delayline = np.resize(np.pad(samples, (N * 2, 0)), samples.shape) - - deltas = samples - delayline - - data = deltas[:, None] * twiddles[:-1] - - # TODO loop vs. accumulate - # for i in range(1, M): data[i] += data[i - 1] - np.add.accumulate(data, axis=0, out=data) - - data *= np.conj(twiddles[1:]) - - dfts = window(data) - - return dfts - - -def isdft(dfts, latency=1): - - M, N = dfts.shape - - twiddles = np.array([-1 if n % 2 else +1 for n in np.arange(N)]) \ - if latency == 1 else \ - np.exp(-1j * np.pi * latency * np.arange(N)) - - weight = 2 / (1 - np.cos(np.pi * latency)) - - samples = np.sum(np.real(dfts * twiddles * weight), axis=-1) - - return samples - - -if __name__ == '__main__': - - import matplotlib.pyplot as plot - - sr = 44100 - - f = 1000 - t = np.arange(1 * sr) / sr - x = np.sin(2 * np.pi * f * t) - - dftsize = 1024 - latency = 1 - - dfts = sdft(x, dftsize) - samples = isdft(dfts, latency) - - assert x.shape[0] == dfts.shape[0] - assert x.shape == samples.shape - - print(dfts.shape, samples.shape) - - with np.errstate(divide='ignore', invalid='ignore'): - db = 20 * np.log10(np.abs(dfts)) - - roi = (np.min(t), np.max(t), 0, sr / (dftsize * 2)) - - plot.figure() - plot.imshow(db.T, extent=roi, cmap='inferno', aspect='auto', interpolation='nearest', origin='lower') - plot.clim(-120, 0) - - plot.figure() - plot.plot(t, x, alpha=0.5) - plot.plot(t, samples, alpha=0.5) - plot.xlim(0, 0.1) - - plot.show()