From 096027468b156618a511ce535abe7dd63ed1e4e3 Mon Sep 17 00:00:00 2001 From: Paul Date: Fri, 24 Nov 2023 17:26:54 -0500 Subject: [PATCH] Add a Hilbert Transform operator (#72) Add a Hilbert Transform operator, both in float -> complex and stereo pair -> sse vector form. Add tests and update documentation. Thanks to Sean Costello for providing the weights and idea to use serial biquads. --- README.md | 3 +- .../sst/basic-blocks/dsp/HilbertTransform.h | 292 ++++++++++++++++++ .../sst/basic-blocks/dsp/LanczosResampler.h | 2 +- .../basic-blocks/modulators/FXModControl.h | 8 +- tests/dsp_tests.cpp | 164 +++++++++- tests/modulator_tests.cpp | 8 +- 6 files changed, 469 insertions(+), 8 deletions(-) create mode 100644 include/sst/basic-blocks/dsp/HilbertTransform.h diff --git a/README.md b/README.md index c06063e..3374b99 100644 --- a/README.md +++ b/README.md @@ -12,4 +12,5 @@ where they are (1) standalone (or can be trivially made standalone by removing a individual header files are also available to use in an MIT license context. Those header files are explicitly marked in the text of the header file and are listed here also -- include/sst/basic-blocks/dsp/LanczosResampler.h \ No newline at end of file +- include/sst/basic-blocks/dsp/LanczosResampler.h +- include/sst/basic-blocks/dsp/HilbertTransform.h \ No newline at end of file diff --git a/include/sst/basic-blocks/dsp/HilbertTransform.h b/include/sst/basic-blocks/dsp/HilbertTransform.h new file mode 100644 index 0000000..5ad5ea2 --- /dev/null +++ b/include/sst/basic-blocks/dsp/HilbertTransform.h @@ -0,0 +1,292 @@ +/* + * sst-basic-blocks - an open source library of core audio utilities + * built by Surge Synth Team. + * + * Provides a collection of tools useful on the audio thread for blocks, + * modulation, etc... or useful for adapting code to multiple environments. + * + * Copyright 2023, various authors, as described in the GitHub + * transaction log. Parts of this code are derived from similar + * functions original in Surge or ShortCircuit. + * + * sst-basic-blocks is released under the GNU General Public Licence v3 + * or later (GPL-3.0-or-later). The license is found in the "LICENSE" + * file in the root of this repository, or at + * https://www.gnu.org/licenses/gpl-3.0.en.html. + * + * A very small number of explicitly chosen header files can also be + * used in an MIT/BSD context. Please see the README.md file in this + * repo or the comments in the individual files. Only headers with an + * explicit mention that they are dual licensed may be copied and reused + * outside the GPL3 terms. + * + * All source in sst-basic-blocks available at + * https://github.com/surge-synthesizer/sst-basic-blocks + */ + +#ifndef INCLUDE_SST_BASIC_BLOCKS_DSP_HILBERTTRANSFORM_H +#define INCLUDE_SST_BASIC_BLOCKS_DSP_HILBERTTRANSFORM_H + +/* + * This file is available for use under the MIT license or + * the GPL3 license, as described in README. + * + * Thanks to Sean Costello for the conversation which led + * to BaconPaul understanding the hilbert transform in the + * context of pitch shifters, and from sharing this + * serial-biquad implementation with allpass coefficients. + * The coefficients first appear in one of the + * Electronotes series papers by Bernie Hutchins. + * + * The hilbert transform takes a real valued signal + * and returns a complex valued signal which allows you + * to create the "analytic signal" - namely the signal with + * only positive frequencies. That is, take a real signal + * u(t) = cos(wt) and imagine it as a complex signal + * U(t) = 0.5 (e^iwt + e^-iwt). The modified signal + * v(t) = H(u(t)) contains only the positive + * frequency. Basically at this point go read wikipedia + * and stuff, but that gives you loads of cool properties + * for modifying phase and pitch in the complex plane. + * + * The actual calculation of the hilbert transform numerically + * is tricky, which is why we are particularly grateful to + * Sean for sharing this implementation (which he shared + * in float form, and which is presented here both in + * mono-float form and in stereo-SSE form). + */ + +#include +#include + +namespace sst::basic_blocks::dsp +{ +struct HilbertTransformMonoFloat +{ + struct BQ + { + float a1{1}, a2{0}, b0{1}, b1{0}, b2{0}, reg0{0}, reg1{0}; + inline void reset() + { + reg0 = 0; + reg1 = 0; + } + inline void setCoefs(float _a1, float _a2, float _b0, float _b1, float _b2) + { + a1 = _a1; + a2 = _a2; + b0 = _b0; + b1 = _b1; + b2 = _b2; + } + + inline float step(float input) + { + double op; + + op = input * b0 + reg0; + reg0 = input * b1 - a1 * op + reg1; + reg1 = input * b2 - a2 * op; + + return (float)op; + } + } allpass[2][3]; + + float sampleRate{0}; + void setSampleRate(float sr) + { + sampleRate = sr; + setHilbertCoefs(); + } + + float hilbertCoefs[12]; + void setHilbertCoefs() + { + assert(sampleRate); + float a1, a2, b0, b1, b2; + + // phase difference network normalized pole frequencies + // first six numbers are for real, the rest for imaginary + float poles[12] = {.3609f, 2.7412f, 11.1573f, 44.7581f, 179.6242f, 798.4578f, + 1.2524f, 5.5671f, 22.3423f, 89.6271f, 364.7914f, 2770.1114f}; + + float minFreq = 25.0f; // minimum frequency of phase differencing network, in Hertz + + // bilinear z-transform for calculating coefficients for 1st order allpasses + for (int j = 0; j < 12; j++) + hilbertCoefs[j] = (1.0f - (minFreq * M_PI * poles[j]) / sampleRate) / + (1.0f + (minFreq * M_PI * poles[j]) / sampleRate); + + // calculate biquad coeffcients for 2nd order allpasses, using the first order allpasses + // in the following order: highest (in a branch) by lowest, next highest by next lowest, + // next highest by next lowest. + + for (int j = 0; j < 3; j++) + { + allpass[0][j].reset(); // clears states + a1 = -(hilbertCoefs[j] + hilbertCoefs[5 - j]); + a2 = hilbertCoefs[j] * hilbertCoefs[5 - j]; + b0 = hilbertCoefs[j] * hilbertCoefs[5 - j]; + b1 = -(hilbertCoefs[j] + hilbertCoefs[5 - j]); + b2 = 1.0f; + allpass[0][j].setCoefs(a1, a2, b0, b1, b2); + } + + for (int j = 0; j < 3; j++) + { + allpass[1][j].reset(); // clears states + a1 = -(hilbertCoefs[j + 6] + hilbertCoefs[11 - j]); + a2 = hilbertCoefs[j + 6] * hilbertCoefs[11 - j]; + b0 = hilbertCoefs[j + 6] * hilbertCoefs[11 - j]; + b1 = -(hilbertCoefs[j + 6] + hilbertCoefs[11 - j]); + b2 = 1.0f; + allpass[1][j].setCoefs(a1, a2, b0, b1, b2); + } + } + + std::pair stepPair(float in) + { + float im{in}, re{in}; + + for (int i = 0; i < 3; ++i) + { + re = allpass[0][i].step(re); + im = allpass[1][i].step(im); + } + return {re, im}; + } + + std::complex stepComplex(float in) + { + auto [r, i] = stepPair(in); + return {r, i}; + } +}; + +struct HilbertTransformStereoSSE +{ + /* + * Same as above but where we use SSE so the allpass is a 4-wide single + * step with layout reL, imL, reR, imR + */ + struct BQ + { + __m128 a1{1}, a2{0}, b0{1}, b1{0}, b2{0}, reg0{0}, reg1{0}; + inline void reset() + { + reg0 = _mm_setzero_ps(); + reg1 = _mm_setzero_ps(); + } + inline void setOne(int idx, __m128 &on, float f) + { + float r alignas(16)[4]; + _mm_store_ps(r, on); + r[idx] = f; + on = _mm_load_ps(r); + } + inline void setCoefs(int idx, float _a1, float _a2, float _b0, float _b1, float _b2) + { + setOne(idx, a1, _a1); + setOne(idx, a2, _a2); + setOne(idx, b0, _b0); + setOne(idx, b1, _b1); + setOne(idx, b2, _b2); + } + + inline __m128 step(__m128 input) + { + auto op = _mm_add_ps(_mm_mul_ps(input, b0), reg0); + reg0 = _mm_add_ps(_mm_sub_ps(_mm_mul_ps(input, b1), _mm_mul_ps(a1, op)), reg1); + reg1 = _mm_sub_ps(_mm_mul_ps(input, b2), _mm_mul_ps(a2, op)); + + return op; + } + } allpassSSE[3]; + + float sampleRate{0}; + void setSampleRate(float sr) + { + sampleRate = sr; + setHilbertCoefs(); + } + + float hilbertCoefs[12]; + void setHilbertCoefs() + { + assert(sampleRate); + float a1, a2, b0, b1, b2; + + // phase difference network normalized pole frequencies + // first six numbers are for real, the rest for imaginary + float poles[12] = {.3609f, 2.7412f, 11.1573f, 44.7581f, 179.6242f, 798.4578f, + 1.2524f, 5.5671f, 22.3423f, 89.6271f, 364.7914f, 2770.1114f}; + + float minFreq = 25.0f; // minimum frequency of phase differencing network, in Hertz + + // bilinear z-transform for calculating coefficients for 1st order allpasses + for (int j = 0; j < 12; j++) + hilbertCoefs[j] = (1.0f - (minFreq * M_PI * poles[j]) / sampleRate) / + (1.0f + (minFreq * M_PI * poles[j]) / sampleRate); + + // calculate biquad coeffcients for 2nd order allpasses, using the first order allpasses + // in the following order: highest (in a branch) by lowest, next highest by next lowest, + // next highest by next lowest. + + for (int j = 0; j < 3; ++j) + allpassSSE[j].reset(); + + for (int j = 0; j < 3; j++) + { + a1 = -(hilbertCoefs[j] + hilbertCoefs[5 - j]); + a2 = hilbertCoefs[j] * hilbertCoefs[5 - j]; + b0 = hilbertCoefs[j] * hilbertCoefs[5 - j]; + b1 = -(hilbertCoefs[j] + hilbertCoefs[5 - j]); + b2 = 1.0f; + // Real is 0, 2 + allpassSSE[j].setCoefs(0, a1, a2, b0, b1, b2); + allpassSSE[j].setCoefs(2, a1, a2, b0, b1, b2); + } + + for (int j = 0; j < 3; j++) + { + a1 = -(hilbertCoefs[j + 6] + hilbertCoefs[11 - j]); + a2 = hilbertCoefs[j + 6] * hilbertCoefs[11 - j]; + b0 = hilbertCoefs[j + 6] * hilbertCoefs[11 - j]; + b1 = -(hilbertCoefs[j + 6] + hilbertCoefs[11 - j]); + b2 = 1.0f; + allpassSSE[j].setCoefs(1, a1, a2, b0, b1, b2); + allpassSSE[j].setCoefs(3, a1, a2, b0, b1, b2); + } + } + + // Returns reL, imL, reR, imR + __m128 stepStereo(float L, float R) + { + float r alignas(16)[4]{L, L, R, R}; + auto in = _mm_load_ps(r); + for (int i = 0; i < 3; ++i) + { + in = allpassSSE[i].step(in); + } + return in; + } + + std::pair, std::complex> stepToComplex(float L, float R) + { + auto v = stepStereo(L, R); + float r alignas(16)[4]; + _mm_store_ps(r, v); + return {{r[0], r[1]}, {r[2], r[3]}}; + } + + std::pair, std::pair> stepToPair(float L, float R) + { + auto v = stepStereo(L, R); + float r alignas(16)[4]; + _mm_store_ps(r, v); + return {{r[0], r[1]}, {r[2], r[3]}}; + } +}; +} // namespace sst::basic_blocks::dsp + +#endif diff --git a/include/sst/basic-blocks/dsp/LanczosResampler.h b/include/sst/basic-blocks/dsp/LanczosResampler.h index 2ea43ae..ca1a302 100644 --- a/include/sst/basic-blocks/dsp/LanczosResampler.h +++ b/include/sst/basic-blocks/dsp/LanczosResampler.h @@ -236,7 +236,7 @@ template bool LanczosResampler::tablesInitialized{false}; template inline size_t LanczosResampler::populateNext(float *fL, float *fR, size_t max) { - int populated = 0; + size_t populated = 0; while (populated < max && (phaseI - phaseO) > A + 1) { read((phaseI - phaseO), fL[populated], fR[populated]); diff --git a/include/sst/basic-blocks/modulators/FXModControl.h b/include/sst/basic-blocks/modulators/FXModControl.h index 15cb51a..d413f6c 100644 --- a/include/sst/basic-blocks/modulators/FXModControl.h +++ b/include/sst/basic-blocks/modulators/FXModControl.h @@ -12,7 +12,13 @@ * sst-basic-blocks is released under the GNU General Public Licence v3 * or later (GPL-3.0-or-later). The license is found in the "LICENSE" * file in the root of this repository, or at - * https://www.gnu.org/licenses/gpl-3.0.en.html + * https://www.gnu.org/licenses/gpl-3.0.en.html. + * + * A very small number of explicitly chosen header files can also be + * used in an MIT/BSD context. Please see the README.md file in this + * repo or the comments in the individual files. Only headers with an + * explicit mention that they are dual licensed may be copied and reused + * outside the GPL3 terms. * * All source in sst-basic-blocks available at * https://github.com/surge-synthesizer/sst-basic-blocks diff --git a/tests/dsp_tests.cpp b/tests/dsp_tests.cpp index c6463bd..c9e7452 100644 --- a/tests/dsp_tests.cpp +++ b/tests/dsp_tests.cpp @@ -30,6 +30,7 @@ #include "sst/basic-blocks/dsp/BlockInterpolators.h" #include "sst/basic-blocks/dsp/QuadratureOscillators.h" #include "sst/basic-blocks/dsp/LanczosResampler.h" +#include "sst/basic-blocks/dsp/HilbertTransform.h" #include "sst/basic-blocks/dsp/FastMath.h" #include "sst/basic-blocks/dsp/Clippers.h" #include "sst/basic-blocks/mechanics/block-ops.h" @@ -252,7 +253,7 @@ TEST_CASE("LanczosResampler", "[dsp]") } float outBlock alignas(16)[64], outBlockR alignas(16)[64]; - int q, gen; + int gen; dp /= 88100.0 / 48000.0; phase = 0; @@ -306,7 +307,6 @@ TEST_CASE("LanczosResampler", "[dsp]") if (spool > 10) { auto expL = std::sin((phase - 2 * dp) * 2.0 * M_PI); - auto expR = phase * 2 - 1; diff += fabs(expL - L); maxDiff = std::max((float)fabs(expL - L), maxDiff); @@ -725,7 +725,7 @@ TEST_CASE("lipol_ps class", "[dsp]") assert(mypol.blockSize == nfloat); mypol.store_block(storeTarget); - for (auto i = 0; i < nfloat; ++i) + for (auto i = 0U; i < nfloat; ++i) REQUIRE(storeTarget[i] == prevtarget); // should be constant in the first instance for (int i = 0; i < 10; ++i) @@ -738,7 +738,7 @@ TEST_CASE("lipol_ps class", "[dsp]") REQUIRE(storeTarget[nfloat - 1] == Approx(target)); float dy = storeTarget[1] - storeTarget[0]; - for (auto j = 1; j < nfloat; ++j) + for (auto j = 1U; j < nfloat; ++j) { REQUIRE(storeTarget[j] - storeTarget[j - 1] == Approx(dy).epsilon(1e-3)); } @@ -748,3 +748,159 @@ TEST_CASE("lipol_ps class", "[dsp]") prevtarget = target; } } + +TEST_CASE("Hilbert Float", "[dsp]") +{ + SECTION("Hilbert has Positive Frequencies") + { + /* + * This test does the hibert transform of real sin + * and shows that the angle of the resulting complex + * tracks only positive frequencies. Which is basically + * the definition. + */ + auto sr = 48000; + sst::basic_blocks::dsp::QuadratureOscillator src, qoP, qoN; + sst::basic_blocks::dsp::HilbertTransformMonoFloat hilbert; + hilbert.setSampleRate(sr); + + src.setRate(2.0 * M_PI * 440 / sr); + qoP.setRate(2.0 * M_PI * 440 / sr); + qoN.setRate(-2.0 * M_PI * 440 / sr); + + auto accDP{0.f}, accDN{0.f}, dP{0.f}, dN{0.f}; + constexpr int nSteps = 2500, warmUp = 100; + for (auto i = 0; i < nSteps; ++i) + { + auto sC = hilbert.stepComplex(src.u); + auto qP = std::complex(qoP.u, -qoP.v); + auto qN = std::complex(qoN.u, -qoN.v); + + auto aC = std::arg(sC); + auto aP = std::arg(qP); + auto aN = std::arg(qN); + + if (i == warmUp) + { + dP = aC - aP; + dN = aC - aN; + } + if (i > warmUp) + { + if (aC < aP) + { + aP -= 2.0 * M_PI; + } + if (aC < aN) + { + aN -= 2.0 * M_PI; + } + accDN += std::abs(aC - aN - dN); + accDP += std::abs(aC - aP - dP); + } + src.step(); + qoP.step(); + qoN.step(); + } + + accDN /= (nSteps - warmUp); + accDP /= (nSteps - warmUp); + REQUIRE(accDP < 0.1); + REQUIRE(accDN > 1); + } +} + +TEST_CASE("Hilbert SSE", "[dsp]") +{ + SECTION("Hilbert SSE has Positive Frequencies") + { + /* + * This test does the hibert transform of real sin + * and shows that the angle of the resulting complex + * tracks only positive frequencies. Which is basically + * the definition. + */ + auto sr = 48000; + sst::basic_blocks::dsp::QuadratureOscillator src, qoP, qoN; + sst::basic_blocks::dsp::QuadratureOscillator srcR, qoPR, qoNR; + sst::basic_blocks::dsp::HilbertTransformStereoSSE hilbert; + hilbert.setSampleRate(sr); + + src.setRate(2.0 * M_PI * 440 / sr); + qoP.setRate(2.0 * M_PI * 440 / sr); + qoN.setRate(-2.0 * M_PI * 440 / sr); + + srcR.setRate(2.0 * M_PI * 762 / sr); + qoPR.setRate(2.0 * M_PI * 762 / sr); + qoNR.setRate(-2.0 * M_PI * 762 / sr); + + auto accDP{0.f}, accDN{0.f}, dP{0.f}, dN{0.f}; + auto accDPR{0.f}, accDNR{0.f}, dPR{0.f}, dNR{0.f}; + constexpr int nSteps = 2500, warmUp = 100; + for (auto i = 0; i < nSteps; ++i) + { + auto [sC, sR] = hilbert.stepToComplex(src.u, srcR.u); + auto qP = std::complex(qoP.u, -qoP.v); + auto qN = std::complex(qoN.u, -qoN.v); + + auto qPR = std::complex(qoPR.u, -qoPR.v); + auto qNR = std::complex(qoNR.u, -qoNR.v); + + auto aC = std::arg(sC); + auto aP = std::arg(qP); + auto aN = std::arg(qN); + + auto aCR = std::arg(sR); + auto aPR = std::arg(qPR); + auto aNR = std::arg(qNR); + + if (aC < aP) + { + aP -= 2.0 * M_PI; + } + if (aC < aN) + { + aN -= 2.0 * M_PI; + } + if (aCR < aPR) + { + aPR -= 2.0 * M_PI; + } + if (aCR < aNR) + { + aNR -= 2.0 * M_PI; + } + + if (i == warmUp) + { + dP = aC - aP; + dN = aC - aN; + dPR = aCR - aPR; + dNR = aCR - aNR; + } + if (i > warmUp) + { + accDN += std::abs(aC - aN - dN); + accDP += std::abs(aC - aP - dP); + + accDNR += std::abs(aCR - aNR - dNR); + accDPR += std::abs(aCR - aPR - dPR); + } + src.step(); + srcR.step(); + qoP.step(); + qoN.step(); + qoPR.step(); + qoNR.step(); + } + + accDN /= (nSteps - warmUp); + accDP /= (nSteps - warmUp); + accDNR /= (nSteps - warmUp); + accDPR /= (nSteps - warmUp); + REQUIRE(accDP < 0.1); + REQUIRE(accDN > 1); + REQUIRE(accDPR < 0.1); + REQUIRE(accDNR > 1); + } +} \ No newline at end of file diff --git a/tests/modulator_tests.cpp b/tests/modulator_tests.cpp index 44c029d..53672fb 100644 --- a/tests/modulator_tests.cpp +++ b/tests/modulator_tests.cpp @@ -12,7 +12,13 @@ * sst-basic-blocks is released under the GNU General Public Licence v3 * or later (GPL-3.0-or-later). The license is found in the "LICENSE" * file in the root of this repository, or at - * https://www.gnu.org/licenses/gpl-3.0.en.html + * https://www.gnu.org/licenses/gpl-3.0.en.html. + * + * A very small number of explicitly chosen header files can also be + * used in an MIT/BSD context. Please see the README.md file in this + * repo or the comments in the individual files. Only headers with an + * explicit mention that they are dual licensed may be copied and reused + * outside the GPL3 terms. * * All source in sst-basic-blocks available at * https://github.com/surge-synthesizer/sst-basic-blocks