Poly octave working with new filter bank

CMakeLists.txt

@@ -4,9 +4,9 @@ if(WIN32)
     set(CMAKE_SYSTEM_VERSION 7.1 CACHE STRING INTERNAL FORCE) # Windows SDK for Windows 7 and up
 endif()
 if(APPLE AND (CMAKE_SYSTEM_NAME STREQUAL "iOS"))
-    project(BYOD VERSION 2.1.0)
+    project(BYOD VERSION 2.1.1)
 else()
-    project(BYOD VERSION 1.3.0)
+    project(BYOD VERSION 1.3.1)
 endif()
 set(CMAKE_CXX_STANDARD 20)
 

src/processors/other/poly_octave/PolyOctave.cpp

@@ -13,11 +13,11 @@ const String v1Tag = "v1_mode";
 
 PolyOctave::PolyOctave (UndoManager* um)
     : BaseProcessor (
-        "Poly Octave",
-        createParameterLayout(),
-        BasicInputPort {},
-        OutputPort {},
-        um)
+          "Poly Octave",
+          createParameterLayout(),
+          BasicInputPort {},
+          OutputPort {},
+          um)
 {
     using namespace ParameterHelpers;
     const auto setupGainParam = [this] (const juce::String& paramID,
@@ -70,10 +70,18 @@ void PolyOctave::prepare (double sampleRate, int samplesPerBlock)
     upOctaveBuffer_double.setMaxSize (2, 2 * samplesPerBlock); // allocate extra space for SIMD
     downOctaveBuffer_double.setMaxSize (2, samplesPerBlock);
 
-    poly_octave_v1::designFilterBank (octaveUpFilterBank, 2.0, 5.0, 6.0, sampleRate);
-    poly_octave_v1::designFilterBank (octaveUp2FilterBank, 3.0, 7.0, 4.0, sampleRate);
-
+    poly_octave_v2::design_filter_bank<poly_octave_v2::N1> (octaveUpFilterBank, 2.0, 5.0, 5.0, sampleRate);
+    poly_octave_v2::design_filter_bank<poly_octave_v2::N1> (octaveUp2FilterBank, 3.0, 7.0, 6.0, sampleRate);
     for (auto& shifter : downOctavePitchShifters)
+    {
+        shifter.prepare (sampleRate);
+        shifter.set_pitch_factor (0.5f);
+    }
+
+    poly_octave_v1::designFilterBank (octaveUpFilterBank_v1, 2.0, 5.0, 6.0, sampleRate);
+    poly_octave_v1::designFilterBank (octaveUp2FilterBank_v1, 3.0, 7.0, 4.0, sampleRate);
+
+    for (auto& shifter : downOctavePitchShifters_v1)
     {
         shifter.prepare (sampleRate);
         shifter.set_pitch_factor (0.5);
@@ -86,9 +94,9 @@ void PolyOctave::prepare (double sampleRate, int samplesPerBlock)
     }
 
     mixOutBuffer.setSize (2, samplesPerBlock);
-    up1OutBuffer.setSize (2, samplesPerBlock);
-    up2OutBuffer.setSize (2, samplesPerBlock);
-    down1OutBuffer.setSize (2, samplesPerBlock);
+    up1OutBuffer.setSize (2, 4 * samplesPerBlock + 8); // padding for SIMD
+    up2OutBuffer.setSize (2, 4 * samplesPerBlock + 8); // padding for SIMD
+    down1OutBuffer.setSize (2, 4 * samplesPerBlock + 8); // padding for SIMD
 }
 
 void PolyOctave::processAudio (AudioBuffer<float>& buffer)
@@ -99,7 +107,62 @@ void PolyOctave::processAudio (AudioBuffer<float>& buffer)
         return;
     }
 
+    const auto numChannels = buffer.getNumChannels();
+    const auto numSamples = buffer.getNumSamples();
+
+    mixOutBuffer.setSize (numChannels, numSamples, false, false, true);
+    up1OutBuffer.setSize (numChannels, numSamples, false, false, true);
+    up2OutBuffer.setSize (numChannels, numSamples, false, false, true);
+    down1OutBuffer.setSize (numChannels, numSamples, false, false, true);
+
+    // "down" processing
+    for (auto [ch, data_in, data_out] : chowdsp::buffer_iters::zip_channels (std::as_const (buffer), down1OutBuffer))
+    {
+        for (auto [x, y] : chowdsp::zip (data_in, data_out))
+            y = downOctavePitchShifters[(size_t) ch].process_sample (x);
+    }
+    downOctaveGain.process (numSamples);
+    chowdsp::BufferMath::applyGainSmoothedBuffer (down1OutBuffer, downOctaveGain);
+
+    // "up1" processing
+    for (auto [ch, data_in, data_out] : chowdsp::buffer_iters::zip_channels (std::as_const (buffer), up1OutBuffer))
+    {
+        poly_octave_v2::process<1> (octaveUpFilterBank[ch],
+                                    data_in.data(),
+                                    data_out.data(),
+                                    numSamples);
+    }
+    upOctaveGain.process (numSamples);
+    chowdsp::BufferMath::applyGainSmoothedBuffer (up1OutBuffer, upOctaveGain);
+
+    // "up2" processing
+    for (auto [ch, data_in, data_out] : chowdsp::buffer_iters::zip_channels (std::as_const (buffer), up2OutBuffer))
+    {
+        poly_octave_v2::process<2> (octaveUp2FilterBank[ch],
+                                    data_in.data(),
+                                    data_out.data(),
+                                    numSamples);
+    }
+    up2OctaveGain.process (numSamples);
+    chowdsp::BufferMath::applyGainSmoothedBuffer (up2OutBuffer, up2OctaveGain);
 
+    chowdsp::BufferMath::copyBufferData (buffer, mixOutBuffer);
+    dryGain.process (numSamples);
+    chowdsp::BufferMath::applyGainSmoothedBuffer (mixOutBuffer, dryGain);
+
+    chowdsp::BufferMath::addBufferData (up1OutBuffer, mixOutBuffer);
+    chowdsp::BufferMath::addBufferData (up2OutBuffer, mixOutBuffer);
+    chowdsp::BufferMath::addBufferData (down1OutBuffer, mixOutBuffer);
+
+    dcBlocker[MixOutput].processBlock (mixOutBuffer);
+    dcBlocker[Up1Output].processBlock (up1OutBuffer);
+    dcBlocker[Up2Output].processBlock (up2OutBuffer);
+    dcBlocker[Down1Output].processBlock (down1OutBuffer);
+
+    outputBuffers.getReference (MixOutput) = &mixOutBuffer;
+    outputBuffers.getReference (Up1Output) = &up1OutBuffer;
+    outputBuffers.getReference (Up2Output) = &up2OutBuffer;
+    outputBuffers.getReference (Down1Output) = &down1OutBuffer;
 }
 
 void PolyOctave::processAudioV1 (AudioBuffer<float>& buffer)
@@ -118,7 +181,7 @@ void PolyOctave::processAudioV1 (AudioBuffer<float>& buffer)
     for (auto [ch, data_in, data_out] : chowdsp::buffer_iters::zip_channels (std::as_const (doubleBuffer), downOctaveBuffer_double))
     {
         for (auto [x, y] : chowdsp::zip (data_in, data_out))
-            y = downOctavePitchShifters[(size_t) ch].process_sample (x);
+            y = downOctavePitchShifters_v1[(size_t) ch].process_sample (x);
     }
 
     // "up" processing
@@ -137,8 +200,8 @@ void PolyOctave::processAudioV1 (AudioBuffer<float>& buffer)
         jassert (juce::snapPointerToAlignment (up2DataSIMD, xsimd::default_arch::alignment()) == up2DataSIMD);
         std::fill (up2DataSIMD, up2DataSIMD + numSamples, float_2 {});
 
-        auto& upFilterBank = octaveUpFilterBank[static_cast<size_t> (ch)];
-        auto& up2FilterBank = octaveUp2FilterBank[static_cast<size_t> (ch)];
+        auto& upFilterBank = octaveUpFilterBank_v1[static_cast<size_t> (ch)];
+        auto& up2FilterBank = octaveUp2FilterBank_v1[static_cast<size_t> (ch)];
         static constexpr auto eps = std::numeric_limits<double>::epsilon();
 
         for (size_t k = 0; k < poly_octave_v1::ComplexERBFilterBank::numFilterBands; k += float_2::size)
@@ -232,7 +295,6 @@ void PolyOctave::processAudioV1 (AudioBuffer<float>& buffer)
     outputBuffers.getReference (Down1Output) = &down1OutBuffer;
 }
 
-
 void PolyOctave::processAudioBypassed (AudioBuffer<float>& buffer)
 {
     const auto numSamples = buffer.getNumSamples();

src/processors/other/poly_octave/PolyOctave.h

@@ -40,14 +40,19 @@ class PolyOctave : public BaseProcessor
     chowdsp::SmoothedBufferValue<float> up2OctaveGain {};
     chowdsp::SmoothedBufferValue<float> downOctaveGain {};
 
+    // V2 processing stuff...
+    std::array<poly_octave_v2::ComplexERBFilterBank<poly_octave_v2::N1>, 2> octaveUpFilterBank;
+    std::array<poly_octave_v2::ComplexERBFilterBank<poly_octave_v2::N1>, 2> octaveUp2FilterBank;
+    std::array<pitch_shift::Processor<float>, 2> downOctavePitchShifters;
+
     // V1 processing stuff...
     chowdsp::Buffer<double> doubleBuffer;
     chowdsp::Buffer<double> upOctaveBuffer_double;
     chowdsp::Buffer<double> up2OctaveBuffer_double;
     chowdsp::Buffer<double> downOctaveBuffer_double;
-    std::array<poly_octave_v1::ComplexERBFilterBank, 2> octaveUpFilterBank;
-    std::array<poly_octave_v1::ComplexERBFilterBank, 2> octaveUp2FilterBank;
-    std::array<pitch_shift::Processor<double>, 2> downOctavePitchShifters;
+    std::array<poly_octave_v1::ComplexERBFilterBank, 2> octaveUpFilterBank_v1;
+    std::array<poly_octave_v1::ComplexERBFilterBank, 2> octaveUp2FilterBank_v1;
+    std::array<pitch_shift::Processor<double>, 2> downOctavePitchShifters_v1;
 
     std::array<chowdsp::FirstOrderHPF<float>, (size_t) numOutputs> dcBlocker;
 

src/processors/other/poly_octave/PolyOctaveFilterBandHelpers.h

@@ -159,9 +159,178 @@ static void designFilterBank (std::array<ComplexERBFilterBank, 2>& filterBank,
         }
     }
 }
+} // namespace poly_octave_v1
 
-void process_filter_bank()
+namespace poly_octave_v2
 {
+// Reference for filter-bank design and octave shifting:
+// https://aaltodoc.aalto.fi/server/api/core/bitstreams/ff9e52cf-fd79-45eb-b695-93038244ec0e/content
+
+inline std::pair<T, T> process_sample (const T& x,
+                                       const std::array<T, 3>& b_shared_coeffs,
+                                       const std::array<T, 3>& b_real_coeffs,
+                                       const std::array<T, 3>& b_imag_coeffs,
+                                       const std::array<T, 3>& a_coeffs,
+                                       std::array<T, 3>& z_shared,
+                                       std::array<T, 3>& z_real,
+                                       std::array<T, 3>& z_imag)
+{
+    const auto y_shared = z_shared[1] + x * b_shared_coeffs[0];
+    z_shared[1] = z_shared[2] + x * b_shared_coeffs[1] - y_shared * a_coeffs[1];
+    z_shared[2] = x * b_shared_coeffs[2] - y_shared * a_coeffs[2];
+
+    const auto y_real = z_real[1] + y_shared; // for the real filter, we know that b[0] == 1
+    z_real[1] = z_real[2] + y_shared * b_real_coeffs[1] - y_real * a_coeffs[1];
+    z_real[2] = y_shared * b_real_coeffs[2] - y_real * a_coeffs[2];
+
+    const auto y_imag = z_imag[1]; // for the imaginary filter, we know that b[0] == 0
+    z_imag[1] = z_imag[2] + y_shared * b_imag_coeffs[1] - y_imag * a_coeffs[1];
+    z_imag[2] = y_shared * b_imag_coeffs[2] - y_imag * a_coeffs[2];
+
+    return { y_real, y_imag };
+}
+
+static constexpr auto q_c = 4.0;
+static auto design_erb_filter (size_t erb_index,
+                               double gamma,
+                               double erb_start,
+                               double q_ERB,
+                               double sample_rate,
+                               double (&b_coeffs_cplx_shared)[3],
+                               double (&b_coeffs_cplx_real)[3],
+                               double (&b_coeffs_cplx_imag)[3],
+                               double (&a_coeffs_cplx)[3])
+{
+    const auto q_PS = gamma;
+
+    const auto z = erb_start + static_cast<double> (erb_index) * (q_c / q_ERB);
+    const auto center_target_freq = 228.7 * (std::pow (10.0, z / 21.3) - 1.0);
+    const auto filter_q = (1.0 / (q_PS * q_ERB)) * (24.7 + 0.108 * center_target_freq);
+
+    double a_coeffs_proto[3];
+    chowdsp::CoefficientCalculators::calcSecondOrderBPF (b_coeffs_cplx_shared,
+                                                         a_coeffs_proto,
+                                                         center_target_freq / gamma,
+                                                         filter_q * 0.5,
+                                                         sample_rate);
+
+    auto pole = (std::sqrt (std::pow (std::complex { a_coeffs_proto[1] }, 2.0) - 4.0 * a_coeffs_proto[2]) - a_coeffs_proto[1]) / 2.0;
+    if (std::imag (pole) < 0.0)
+        pole = std::conj (pole);
+    const auto pr = std::real (pole);
+    const auto pi = std::imag (pole);
+
+    // a[] = 1 - 2 pr z + (pi^2 + pr^2) z^2
+    a_coeffs_cplx[0] = 1.0;
+    a_coeffs_cplx[1] = -2.0 * pr;
+    a_coeffs_cplx[2] = pi * pi + pr * pr;
+
+    // b_real[] = 1 - 2 pr z + (-pi^2 + pr^2) z^2
+    b_coeffs_cplx_real[0] = 1.0;
+    b_coeffs_cplx_real[1] = -2.0 * pr;
+    b_coeffs_cplx_real[2] = -pi * pi + pr * pr;
+
+    // b_imag[] = 2 pi z - 2 pi pr z^2
+    b_coeffs_cplx_imag[0] = 0.0;
+    b_coeffs_cplx_imag[1] = 2.0 * pi;
+    b_coeffs_cplx_imag[2] = -2.0 * pi * pr;
+
+    return center_target_freq;
+}
+
+template <size_t N>
+static void design_filter_bank (std::array<ComplexERBFilterBank<N>, 2>& filter_bank,
+                                double gamma,
+                                double erb_start,
+                                double q_ERB,
+                                double sample_rate)
+{
+    for (size_t kiter = 0; kiter < ComplexERBFilterBank<N>::num_filter_bands; ++kiter)
+    {
+        double b_coeffs_cplx_shared_double[3] {};
+        double b_coeffs_cplx_real_double[3] {};
+        double b_coeffs_cplx_imag_double[3] {};
+        double a_coeffs_cplx_double[3] {};
+        design_erb_filter (kiter,
+                           gamma,
+                           erb_start,
+                           q_ERB,
+                           sample_rate,
+                           b_coeffs_cplx_shared_double,
+                           b_coeffs_cplx_real_double,
+                           b_coeffs_cplx_imag_double,
+                           a_coeffs_cplx_double);
+
+        for (auto& bank : filter_bank)
+        {
+            const auto k_div = kiter / T::size;
+            const auto k_off = kiter - (k_div * T::size);
+
+            bank.erb_filter_complex[k_div].z_shared = {};
+            bank.erb_filter_complex[k_div].z_real = {};
+            bank.erb_filter_complex[k_div].z_imag = {};
+
+            for (size_t i = 0; i < 3; ++i)
+            {
+                reinterpret_cast<S*> (&bank.erb_filter_complex[k_div].b_shared_coeffs[i])[k_off] = static_cast<S> (b_coeffs_cplx_shared_double[i]);
+                reinterpret_cast<S*> (&bank.erb_filter_complex[k_div].b_real_coeffs[i])[k_off] = static_cast<S> (b_coeffs_cplx_real_double[i]);
+                reinterpret_cast<S*> (&bank.erb_filter_complex[k_div].b_imag_coeffs[i])[k_off] = static_cast<S> (b_coeffs_cplx_imag_double[i]);
+                reinterpret_cast<S*> (&bank.erb_filter_complex[k_div].a_coeffs[i])[k_off] = static_cast<S> (a_coeffs_cplx_double[i]);
+            }
+        }
+    }
+}
+
+template <size_t num_octaves_up, size_t N>
+static void process (ComplexERBFilterBank<N>& filter_bank,
+                     const float* buffer_in,
+                     float* buffer_out,
+                     int num_samples) noexcept
+{
+    // buffer_out size is padded by 4x
+    static constexpr auto eps = std::numeric_limits<T>::epsilon();
+    auto* buffer_out_simd = juce::snapPointerToAlignment (reinterpret_cast<T*> (buffer_out), xsimd::default_arch::alignment());
+    std::fill (buffer_out_simd, buffer_out_simd + num_samples, T {});
+    for (size_t k = 0; k < N; k += T::size)
+    {
+        const auto filter_idx = k / T::size;
+        auto& cplx_filter = filter_bank.erb_filter_complex[filter_idx];
+        chowdsp::ScopedValue z_shared { cplx_filter.z_shared };
+        chowdsp::ScopedValue z_re { cplx_filter.z_real };
+        chowdsp::ScopedValue z_im { cplx_filter.z_imag };
+        for (int n = 0; n < num_samples; ++n)
+        {
+            const auto x_in = static_cast<T> (buffer_in[n]);
+            const auto [x_re, x_im] = process_sample (x_in,
+                                                      cplx_filter.b_shared_coeffs,
+                                                      cplx_filter.b_real_coeffs,
+                                                      cplx_filter.b_imag_coeffs,
+                                                      cplx_filter.a_coeffs,
+                                                      z_shared.get(),
+                                                      z_re.get(),
+                                                      z_im.get());
+
+            auto x_re_sq = x_re * x_re;
+            auto x_im_sq = x_im * x_im;
+            auto x_abs_sq = x_re_sq + x_im_sq;
+
+            if constexpr (num_octaves_up == 1)
+            {
+                auto x_abs_r = xsimd::select (x_abs_sq > eps, xsimd::rsqrt (x_abs_sq), {});
+                buffer_out_simd[n] += (x_re_sq - x_im_sq) * x_abs_r;
+            }
+            else if constexpr (num_octaves_up == 2)
+            {
+                auto x_abs_sq_r = xsimd::select (x_abs_sq > eps, xsimd::reciprocal (x_abs_sq), {});
+                buffer_out_simd[n] += x_re * (x_re_sq - (S) 3 * x_im_sq) * x_abs_sq_r;
+            }
+        }
+    }
+
+    for (int n = 0; n < num_samples; ++n)
+        buffer_out[n] = xsimd::reduce_add (buffer_out_simd[n]);
 
+    static constexpr auto norm_gain = 2.0f / static_cast<float> (N);
+    juce::FloatVectorOperations::multiply (buffer_out, norm_gain, num_samples);
 }
-} // namespace FilterBankHelpers
+} // namespace poly_octave_v2

src/processors/other/poly_octave/PolyOctaveFilterBankTypes.h

@@ -11,3 +11,28 @@ struct ComplexERBFilterBank
     std::array<chowdsp::IIRFilter<4, float_2>, numFilterBands / float_2::size> erbFilterReal, erbFilterImag;
 };
 }
+
+namespace poly_octave_v2
+{
+using T = xsimd::batch<float>;
+using S = float;
+static constexpr auto N1 = 32;
+template <size_t N>
+struct ComplexERBFilterBank
+{
+    static constexpr size_t num_filter_bands = N;
+
+    struct ComplexFilter
+    {
+        std::array<T, 3> b_shared_coeffs {};
+        std::array<T, 3> b_real_coeffs {};
+        std::array<T, 3> b_imag_coeffs {};
+        std::array<T, 3> a_coeffs {};
+        std::array<T, 3> z_shared {};
+        std::array<T, 3> z_real {};
+        std::array<T, 3> z_imag {};
+    };
+
+    std::array<ComplexFilter, num_filter_bands / T::size> erb_filter_complex;
+};
+}