Adding a low-pass filter before downsample reduces aliasing quite dramatically when using non-linear processing at higher rate (96k or 192k)

dsp/RecursiveLinearFilter.cpp

@@ -162,3 +162,20 @@ void recursive_linear_filter::HighShelf::SetParams(const recursive_linear_filter
 
   this->_AssignCoefficients(a0, a1, a2, b0, b1, b2);
 }
+
+void recursive_linear_filter::LowPassBiquad::SetParams(const recursive_linear_filter::BiquadParams& params)
+{
+    const double omega0 = params.GetOmega0();
+    const double cosw0 = std::cos(omega0);
+    const double alpha = params.GetAlpha(omega0);
+
+    const double b0 = (1.0 - cosw0) / 2.0;
+    const double b1 = 1.0 - cosw0;
+    const double b2 = (1.0 - cosw0) / 2.0;
+    const double a0 = 1.0 + alpha;
+    const double a1 = -2.0 * cosw0;
+    const double a2 = 1.0 - alpha;
+
+    // Normalize the coefficients by a0 and assign them
+    this->_AssignCoefficients(a0, a1, a2, b0, b1, b2);
+}

dsp/RecursiveLinearFilter.h

@@ -202,4 +202,11 @@ class LowPass : public Base
   }
 };
 
+class LowPassBiquad : public Biquad
+{
+public:
+    LowPassBiquad() : Biquad() {}
+    void SetParams(const BiquadParams& params) override;
+};
+
 }; // namespace recursive_linear_filter

dsp/ResamplingContainer/Dependencies/LanczosResampler.h

@@ -208,7 +208,7 @@ class LanczosResampler
     return static_cast<size_t>(std::max(res + 1.0, 0.0));
   }
 
-  inline void PushBlock(T** inputs, size_t nFrames)
+  virtual inline void PushBlock(T** inputs, size_t nFrames)
   {
     for (auto s = 0; s < nFrames; s++)
     {

dsp/ResamplingContainer/Dependencies/LanczosResamplerWithLPF.h

@@ -0,0 +1,66 @@
+// File: ResamplingContainer.h
+// Created Date: Monday February 5th 2024
+// Author: Mikko Honkala (mikko.honkala@gmail.com)
+
+// A container for real-time resampling using a Low pass filtered Lanczos anti-aliasing filter
+
+#pragma once
+
+#include <iostream>
+#include <functional>
+#include <cmath>
+
+#include "LanczosResampler.h"
+#include "AudioDSPTools/dsp/RecursiveLinearFilter.h"
+
+
+namespace dsp
+{
+/**
+ * Extends the LanczosResampler class by integrating a LowPassBiquad filter to reduce aliasing in
+ * the case of downsampilng.
+ * This subclass applies a low-pass filter to the input audio signal before performing
+ * the Lanczos downsampling process. In case of upsamping, the filter is not used.
+ *
+ * Template parameters:
+ * - T: The data type of the audio samples (e.g., float or double).
+ * - NCHANS: The number of audio channels to process (e.g., 1 for mono, 2 for stereo).
+ * - A: The filter size parameter of the Lanczos resampler, affecting quality and latency.
+ */
+template<typename T = double, int NCHANS = 2, size_t A = 12>
+class LanczosResamplerWithLPF : public LanczosResampler<T, NCHANS, A> {
+public:
+    recursive_linear_filter::LowPassBiquad lowPassFilter;
+    bool applyLPF = false; // Conditionally apply LPF
+
+    // Constructor
+    LanczosResamplerWithLPF(float inputRate, float outputRate, float cutoffRatioOfNyquist = 0.9)
+    : LanczosResampler<T, NCHANS, A>(inputRate, outputRate) {
+        // Compute cutoff frequency based on the output rate, set to just below Nyquist
+        float cutoffFrequency = outputRate / 2.0 * cutoffRatioOfNyquist; // e.g., 90% of Nyquist frequency
+
+        // Only initialize and apply LPF if downsampling
+        if (outputRate < inputRate) {
+            double qualityFactor = 0.707; // Common choice for a Butterworth filter
+            double gainDB = 0.0; // No gain change for a low-pass filter
+            recursive_linear_filter::BiquadParams params(inputRate, cutoffFrequency, qualityFactor, gainDB);
+            lowPassFilter.SetParams(params);
+            applyLPF = true;
+        }
+    }
+
+    // Override the PushBlock method to conditionally apply LPF
+    void PushBlock(T** inputs, size_t nFrames) override {
+        if (applyLPF) {
+            // Apply the low-pass filter to the inputs before resampling
+            DSP_SAMPLE** dspInputs = reinterpret_cast<DSP_SAMPLE**>(inputs);
+            DSP_SAMPLE** filteredOutputs = lowPassFilter.Process(dspInputs, NCHANS, nFrames);
+            LanczosResampler<T, NCHANS, A>::PushBlock(reinterpret_cast<T**>(filteredOutputs), nFrames);
+        } else {
+            // Directly call base class PushBlock without LPF
+            LanczosResampler<T, NCHANS, A>::PushBlock(inputs, nFrames);
+        }
+    }
+};
+
+}; // namespace dsp

dsp/ResamplingContainer/ResamplingContainer.h

@@ -54,10 +54,12 @@ iPlug 2 includes the following 3rd party libraries (see each license info):
 #include "Dependencies/WDL/ptrlist.h"
 
 #include "Dependencies/LanczosResampler.h"
+#include "AudioDSPTools/dsp/RecursiveLinearFilter.h"
+
+#include "Dependencies/LanczosResamplerWithLPF.h"
 
 namespace dsp
 {
-
 /** A multi-channel real-time resampling container that can be used to resample
  * audio processing to a specified sample rate for the situation where you have
  * some arbitary DSP code that requires a specific sample rate, then back to
@@ -86,7 +88,7 @@ class ResamplingContainer
 {
 public:
   using BlockProcessFunc = std::function<void(T**, T**, int)>;
-  using LanczosResampler = LanczosResampler<T, NCHANS, A>;
+  using LanczosResamplerWithLPF = LanczosResamplerWithLPF<T, NCHANS, A>;
 
   // :param renderingSampleRate: The sample rate required by the code to be encapsulated.
   ResamplingContainer(double renderingSampleRate)
@@ -130,9 +132,9 @@ class ResamplingContainer
     }
 
     {
-      mResampler1 = std::make_unique<LanczosResampler>(mInputSampleRate, mRenderingSampleRate);
-      mResampler2 = std::make_unique<LanczosResampler>(mRenderingSampleRate, mInputSampleRate);
-
+      mResampler1 = std::make_unique<LanczosResamplerWithLPF>(mInputSampleRate, mRenderingSampleRate);
+      mResampler2 = std::make_unique<LanczosResamplerWithLPF>(mRenderingSampleRate, mInputSampleRate);
+      
       // Zeroes the scratch pointers so that we warm up with silence.
       ClearBuffers();
 
@@ -182,6 +184,7 @@ class ResamplingContainer
       {
         throw std::runtime_error("Got more encapsulated samples than the encapsulated DSP is prepared to handle!");
       }
+
       func(mEncapsulatedInputPointers.GetList(), mEncapsulatedOutputPointers.GetList(), (int)populated1);
       // And push the results into the second resampler so that it has what the external context requires.
       mResampler2->PushBlock(mEncapsulatedOutputPointers.GetList(), populated1);
@@ -210,6 +213,7 @@ class ResamplingContainer
   int GetLatency() const { return mLatency; }
 
 private:
+
   static inline int LinearInterpolate(T** inputs, T** outputs, int inputLen, double ratio, int maxOutputLen)
   {
     // FIXME check through this!
@@ -329,7 +333,10 @@ class ResamplingContainer
   // The sample rate required by the DSP that this object encapsulates
   const double mRenderingSampleRate;
   // Pair of resamplers for (1) external -> encapsulated, (2) encapsulated -> external
-  std::unique_ptr<LanczosResampler> mResampler1, mResampler2;
+  std::unique_ptr<LanczosResamplerWithLPF> mResampler1, mResampler2;
+
 };
 
+
+
 }; // namespace dsp