More FIR specializations

Author: kunitoki

examples/graphics/source/examples/FilterDemo.h

@@ -563,16 +563,8 @@ class PartitionedConvolver::Impl
                 trimmedLength = trimSilenceFromEnd (impulseResponse, length, *options.trimEndSilenceBelowDb);
 
             // Update DirectFIR in-place
-            std::vector<float> directCoefficients;
-
             const auto directCoefficientsCount = std::min (directFIRCoefficientCount, trimmedLength);
-            if (directCoefficientsCount > 0)
-            {
-                directCoefficients.reserve (directCoefficientsCount);
-                directCoefficients.assign (impulseResponse, impulseResponse + directCoefficientsCount);
-            }
-
-            newFIR.setCoefficients (std::move (directCoefficients), headroomScale);
+            newFIR.setCoefficients (impulseResponse, directCoefficientsCount, headroomScale);
 
             // Update FFT layers
             std::size_t consumed = directCoefficientsCount;

modules/yup_dsp/convolution/yup_PartitionedConvolver.cpp

@@ -630,19 +630,21 @@ int FilterDesigner<CoeffType>::designLinkwitzRiley (
 //==============================================================================
 
 template <typename CoeffType>
-std::vector<CoeffType> FilterDesigner<CoeffType>::designFIRLowpass (
+void FilterDesigner<CoeffType>::designFIRLowpass (
+    std::vector<CoeffType>& coefficients,
     int numCoefficients,
     CoeffType cutoffFreq,
     double sampleRate,
-    WindowType windowType) noexcept
+    WindowType windowType,
+    CoeffType windowParameter) noexcept
 {
     jassert (numCoefficients > 0);
     jassert (cutoffFreq > static_cast<CoeffType> (0.0));
     jassert (sampleRate > 0.0);
     jassert (cutoffFreq < static_cast<CoeffType> (sampleRate / 2.0));
 
     numCoefficients = nextOdd (numCoefficients);
-    std::vector<CoeffType> coefficients (numCoefficients);
+    coefficients.resize (numCoefficients);
 
     const auto normalizedCutoff = static_cast<CoeffType> (2.0) * cutoffFreq / static_cast<CoeffType> (sampleRate);
     const int center = (numCoefficients - 1) / 2;
@@ -664,25 +666,27 @@ std::vector<CoeffType> FilterDesigner<CoeffType>::designFIRLowpass (
     // Apply window function
     for (int i = 0; i < numCoefficients; ++i)
     {
-        const auto windowValue = WindowFunctions<CoeffType>::getValue (windowType, i, numCoefficients);
+        const auto windowValue = WindowFunctions<CoeffType>::getValue (windowType, i, numCoefficients, windowParameter);
         coefficients[i] *= windowValue;
     }
 
     // Normalization
     const auto sum = std::accumulate (coefficients.begin(), coefficients.end(), static_cast<CoeffType> (0.0));
     if (sum != static_cast<CoeffType> (0.0))
+    {
         for (auto& c : coefficients)
             c /= sum;
-
-    return coefficients;
+    }
 }
 
 template <typename CoeffType>
-std::vector<CoeffType> FilterDesigner<CoeffType>::designFIRHighpass (
+void FilterDesigner<CoeffType>::designFIRHighpass (
+    std::vector<CoeffType>& coefficients,
     int numCoefficients,
     CoeffType cutoffFreq,
     double sampleRate,
-    WindowType windowType) noexcept
+    WindowType windowType,
+    CoeffType windowParameter) noexcept
 {
     jassert (numCoefficients > 0);
     jassert (cutoffFreq > static_cast<CoeffType> (0.0));
@@ -691,7 +695,7 @@ std::vector<CoeffType> FilterDesigner<CoeffType>::designFIRHighpass (
 
     // Generate lowpass first
     numCoefficients = nextOdd (numCoefficients);
-    auto coefficients = designFIRLowpass (numCoefficients, cutoffFreq, sampleRate, windowType);
+    designFIRLowpass (coefficients, numCoefficients, cutoffFreq, sampleRate, windowType);
 
     // Convert to highpass using spectral inversion
     const int center = (numCoefficients - 1) / 2;
@@ -707,19 +711,21 @@ std::vector<CoeffType> FilterDesigner<CoeffType>::designFIRHighpass (
         hpi += coefficients[n] * ((n & 1) ? static_cast<CoeffType> (-1.0) : static_cast<CoeffType> (1.0));
 
     if (hpi != static_cast<CoeffType> (0.0))
+    {
         for (auto& c : coefficients)
             c /= hpi;
-
-    return coefficients;
+    }
 }
 
 template <typename CoeffType>
-std::vector<CoeffType> FilterDesigner<CoeffType>::designFIRBandpass (
+void FilterDesigner<CoeffType>::designFIRBandpass (
+    std::vector<CoeffType>& coefficients,
     int numCoefficients,
     CoeffType lowCutoffFreq,
     CoeffType highCutoffFreq,
     double sampleRate,
-    WindowType windowType) noexcept
+    WindowType windowType,
+    CoeffType windowParameter) noexcept
 {
     jassert (numCoefficients > 0);
     jassert (lowCutoffFreq > static_cast<CoeffType> (0.0));
@@ -728,7 +734,7 @@ std::vector<CoeffType> FilterDesigner<CoeffType>::designFIRBandpass (
     jassert (highCutoffFreq < static_cast<CoeffType> (sampleRate / 2.0));
 
     numCoefficients = nextOdd (numCoefficients);
-    std::vector<CoeffType> coefficients (numCoefficients);
+    coefficients.resize (numCoefficients);
 
     const auto normalizedLow = static_cast<CoeffType> (2.0) * lowCutoffFreq / static_cast<CoeffType> (sampleRate);
     const auto normalizedHigh = static_cast<CoeffType> (2.0) * highCutoffFreq / static_cast<CoeffType> (sampleRate);
@@ -754,20 +760,20 @@ std::vector<CoeffType> FilterDesigner<CoeffType>::designFIRBandpass (
     // Apply window function
     for (int i = 0; i < numCoefficients; ++i)
     {
-        const auto windowValue = WindowFunctions<CoeffType>::getValue (windowType, i, numCoefficients);
+        const auto windowValue = WindowFunctions<CoeffType>::getValue (windowType, i, numCoefficients, windowParameter);
         coefficients[i] *= windowValue;
     }
-
-    return coefficients;
 }
 
 template <typename CoeffType>
-std::vector<CoeffType> FilterDesigner<CoeffType>::designFIRBandstop (
+void FilterDesigner<CoeffType>::designFIRBandstop (
+    std::vector<CoeffType>& coefficients,
     int numCoefficients,
     CoeffType lowCutoffFreq,
     CoeffType highCutoffFreq,
     double sampleRate,
-    WindowType windowType) noexcept
+    WindowType windowType,
+    CoeffType windowParameter) noexcept
 {
     jassert (numCoefficients > 0);
     jassert (lowCutoffFreq > static_cast<CoeffType> (0.0));
@@ -777,7 +783,7 @@ std::vector<CoeffType> FilterDesigner<CoeffType>::designFIRBandstop (
 
     // Generate bandpass first
     numCoefficients = nextOdd (numCoefficients);
-    auto coefficients = designFIRBandpass (numCoefficients, lowCutoffFreq, highCutoffFreq, sampleRate, windowType);
+    designFIRBandpass (coefficients, numCoefficients, lowCutoffFreq, highCutoffFreq, sampleRate, windowType);
 
     // Convert to bandstop using spectral inversion
     const int center = (numCoefficients - 1) / 2;
@@ -787,8 +793,6 @@ std::vector<CoeffType> FilterDesigner<CoeffType>::designFIRBandstop (
 
     // Add unit impulse at center
     coefficients[center] += static_cast<CoeffType> (1.0);
-
-    return coefficients;
 }
 
 //==============================================================================

modules/yup_dsp/designers/yup_FilterDesigner.cpp

@@ -640,11 +640,13 @@ class FilterDesigner
 
         @returns               Vector of FIR coefficients suitable for DirectFIR
     */
-    static std::vector<CoeffType> designFIRLowpass (
+    static void designFIRLowpass (
+        std::vector<CoeffType>& coefficients,
         int numCoefficients,
         CoeffType cutoffFreq,
         double sampleRate,
-        WindowType windowType = WindowType::hann) noexcept;
+        WindowType windowType = WindowType::hann,
+        CoeffType windowParameter = CoeffType (8)) noexcept;
 
     /**
         Designs FIR highpass filter coefficients using windowed sinc method.
@@ -656,11 +658,13 @@ class FilterDesigner
 
         @returns               Vector of FIR coefficients suitable for DirectFIR
     */
-    static std::vector<CoeffType> designFIRHighpass (
+    static void designFIRHighpass (
+        std::vector<CoeffType>& coefficients,
         int numCoefficients,
         CoeffType cutoffFreq,
         double sampleRate,
-        WindowType windowType = WindowType::hann) noexcept;
+        WindowType windowType = WindowType::hann,
+        CoeffType windowParameter = CoeffType (8)) noexcept;
 
     /**
         Designs FIR bandpass filter coefficients using windowed sinc method.
@@ -673,12 +677,14 @@ class FilterDesigner
 
         @returns               Vector of FIR coefficients suitable for DirectFIR
     */
-    static std::vector<CoeffType> designFIRBandpass (
+    static void designFIRBandpass (
+        std::vector<CoeffType>& coefficients,
         int numCoefficients,
         CoeffType lowCutoffFreq,
         CoeffType highCutoffFreq,
         double sampleRate,
-        WindowType windowType = WindowType::hann) noexcept;
+        WindowType windowType = WindowType::hann,
+        CoeffType windowParameter = CoeffType (8)) noexcept;
 
     /**
         Designs FIR bandstop filter coefficients using windowed sinc method.
@@ -691,12 +697,14 @@ class FilterDesigner
 
         @returns               Vector of FIR coefficients suitable for DirectFIR
     */
-    static std::vector<CoeffType> designFIRBandstop (
+    static void designFIRBandstop (
+        std::vector<CoeffType>& coefficients,
         int numCoefficients,
         CoeffType lowCutoffFreq,
         CoeffType highCutoffFreq,
         double sampleRate,
-        WindowType windowType = WindowType::hann) noexcept;
+        WindowType windowType = WindowType::hann,
+        CoeffType windowParameter = CoeffType (8)) noexcept;
 };
 
 } // namespace yup

modules/yup_dsp/designers/yup_FilterDesigner.h

@@ -100,29 +100,17 @@ class DirectFIR : public FilterBase<SampleType, CoeffType>
 
     //==============================================================================
     /**
-        Set the FIR filter coefficients.
+        Set the FIR filter coefficients from a raw pointer.
 
-        @param coefficients  Vector containing the FIR coefficients in time order
-        @param scaling       Scaling factor to apply to all coefficients
+        @param coefficients     Span of FIR coefficients array
+        @param scaling          Scaling factor to apply to all coefficients
 
         @note This method is not real-time safe and should be called during initialization
-              or when audio processing is paused.
+              or when audio processing is paused, unless the coefficients vector has already been set with a greater or equal size.
     */
-    void setCoefficients (std::vector<CoeffType> coefficients, CoeffType scaling = CoeffType (1))
+    void setCoefficients (yup::Span<const CoeffType> coefficients, CoeffType scaling = CoeffType (1))
     {
-        currentScaling = scaling;
-        if (! approximatelyEqual (currentScaling, CoeffType (1)))
-            FloatVectorOperations::multiply (coefficients.data(), scaling, coefficients.size());
-
-        coefficientsReversed = std::move (coefficients);
-        std::reverse (coefficientsReversed.begin(), coefficientsReversed.end());
-
-        numCoefficients = coefficientsReversed.size();
-        paddedLen = (numCoefficients + 3u) & ~3u; // Round up to multiple of 4 for SIMD
-        coefficientsReversed.resize (paddedLen, 0.0f);
-
-        history.assign (2 * numCoefficients, 0.0f);
-        reset();
+        setCoefficients (coefficients.data(), coefficients.size(), scaling);
     }
 
     /**
@@ -133,7 +121,7 @@ class DirectFIR : public FilterBase<SampleType, CoeffType>
         @param scaling          Scaling factor to apply to all coefficients
 
         @note This method is not real-time safe and should be called during initialization
-              or when audio processing is paused.
+              or when audio processing is paused, unless the coefficients vector has already been set with a greater or equal size.
     */
     void setCoefficients (const CoeffType* coefficients, std::size_t numCoefficientsIn, CoeffType scaling = CoeffType (1))
     {
@@ -144,8 +132,21 @@ class DirectFIR : public FilterBase<SampleType, CoeffType>
             return;
         }
 
-        std::vector<float> coefficientsVector (coefficients, coefficients + numCoefficientsIn);
-        setCoefficients (std::move (coefficientsVector), scaling);
+        numCoefficients = numCoefficientsIn;
+        paddedLen = (numCoefficientsIn + 3u) & ~3u; // Round up to multiple of 4 for SIMD
+
+        coefficientsReversed.resize (paddedLen, 0.0f);
+
+        currentScaling = scaling;
+        if (! approximatelyEqual (currentScaling, CoeffType (1)))
+            FloatVectorOperations::copyWithMultiply (coefficientsReversed.data(), coefficients, scaling, static_cast<int> (numCoefficientsIn));
+        else
+            FloatVectorOperations::copy (coefficientsReversed.data(), coefficients, static_cast<int> (numCoefficientsIn));
+
+        std::reverse (coefficientsReversed.begin(), coefficientsReversed.begin() + numCoefficients);
+
+        history.resize (2 * numCoefficients, 0.0f);
+        writeIndex = 0;
     }
 
     /**
@@ -263,19 +264,29 @@ class DirectFIR : public FilterBase<SampleType, CoeffType>
         if (numCoefficients == 0)
             return Complex<CoeffType> (0, 0);
 
+        // ω = 2π f / Fs
         const CoeffType omega = MathConstants<CoeffType>::twoPi * frequency / static_cast<CoeffType> (this->sampleRate);
 
-        Complex<CoeffType> response (0, 0);
+        // Standard FIR frequency response: H(e^{jω}) = Σ_{n=0}^{N-1} h[n] * e^{-jωn}
+        // coefficientsReversed stores: [h[M-1], h[M-2], ..., h[1], h[0]]
+        // So coefficientsReversed[k] = h[M-1-k], and we need: Σ h[n] * e^{-jωn}
+
+        // e^{-jω}
+        const Complex<CoeffType> ez_neg { std::cos (omega), -std::sin (omega) };
+
+        // Accumulate: Σ_{n=0}^{N-1} h[n] * e^{-jωn}
+        // Since coefficientsReversed[k] = h[M-1-k], we have h[n] = coefficientsReversed[M-1-n]
+        Complex<CoeffType> sum { 0, 0 };
+        Complex<CoeffType> ez_neg_n { 1, 0 }; // e^{-jω*0} = 1
 
-        // H(e^jω) = Σ h[n] * e^(-jωn) for n = 0 to N-1
         for (std::size_t n = 0; n < numCoefficients; ++n)
         {
-            const CoeffType angle = -omega * static_cast<CoeffType> (n);
-            Complex<CoeffType> exponential (std::cos (angle), std::sin (angle));
-            response += static_cast<CoeffType> (coefficientsReversed[numCoefficients - 1 - n]) * exponential;
+            const CoeffType h_n = coefficientsReversed[numCoefficients - 1 - n];
+            sum += h_n * ez_neg_n;
+            ez_neg_n *= ez_neg;
         }
 
-        return response;
+        return sum;
     }
 
     /**