Fix Comb filter

Author: kunitoki

examples/graphics/source/examples/FilterDemo.h

@@ -531,7 +531,10 @@ class FrequencyResponsePlot : public yup::Component
             return;
         }
 
-        const int numPoints = 512;
+        const int numPoints = isCombFilter() ? 4096 : 512;
+
+        minDb = isCombFilter() ? -80.0 : -60.0;
+        maxDb = isCombFilter() ? 40.0 : 20.0;
 
         responseData.clear();
         responseData.resize (numPoints);
@@ -629,7 +632,9 @@ class FrequencyResponsePlot : public yup::Component
         }
 
         // Horizontal dB lines
-        for (double db = -60.0; db <= 20.0; db += 20.0)
+        const auto gridStepDb = isCombFilter() ? 40.0 : 20.0;
+
+        for (double db = minDb; db <= maxDb; db += gridStepDb)
         {
             float y = dbToY (db, bounds);
             g.strokeLine ({ bounds.getX(), y }, { bounds.getRight(), y });
@@ -700,7 +705,9 @@ class FrequencyResponsePlot : public yup::Component
         }
 
         // dB labels
-        for (double db = -60.0; db <= 20.0; db += 20.0)
+        const auto gridStepDb = isCombFilter() ? 40.0 : 20.0;
+
+        for (double db = minDb; db <= maxDb; db += gridStepDb)
         {
             float y = dbToY (db, bounds);
             yup::String label = yup::String (db, 0) + " dB";
@@ -720,6 +727,11 @@ class FrequencyResponsePlot : public yup::Component
         return static_cast<float> (bounds.getBottom() - ratio * bounds.getHeight());
     }
 
+    bool isCombFilter() const
+    {
+        return dynamic_cast<yup::CombFilter<float>*> (filter.get()) != nullptr;
+    }
+
     std::shared_ptr<yup::FilterBase<float, double>> filter;
     std::vector<yup::Complex<float>> responseData;
     std::vector<yup::Complex<float>> phaseData;

modules/yup_dsp/filters/yup_CombFilter.h

@@ -209,13 +209,13 @@ class CombFilter : public FilterBase<SampleType, CoeffType>
     */
     Complex<CoeffType> getComplexResponse (CoeffType responseFrequency) const override
     {
-        const auto omega = frequencyToAngular (responseFrequency, static_cast<CoeffType> (this->sampleRate));
-        const auto delay = getDelayForResponse();
-        const auto delayed = polar (static_cast<CoeffType> (1), -omega * delay);
+        const auto delayed = getFractionalDelayResponse (responseFrequency);
         const auto denominator = Complex<CoeffType> (static_cast<CoeffType> (1)) - delayLineFeedback * delayed;
 
         if (std::abs (denominator) <= std::numeric_limits<CoeffType>::epsilon())
-            return Complex<CoeffType> (static_cast<CoeffType> (1));
+            return Complex<CoeffType> (static_cast<CoeffType> (1))
+                 + (static_cast<CoeffType> (0.5) * delayed)
+                       / (denominator + Complex<CoeffType> (std::numeric_limits<CoeffType>::epsilon()));
 
         return Complex<CoeffType> (static_cast<CoeffType> (1))
              + (static_cast<CoeffType> (0.5) * delayed) / denominator;
@@ -359,6 +359,26 @@ class CombFilter : public FilterBase<SampleType, CoeffType>
                        static_cast<CoeffType> (this->sampleRate) / frequency);
     }
 
+    Complex<CoeffType> getFractionalDelayResponse (CoeffType responseFrequency) const noexcept
+    {
+        const auto delay = getDelayForResponse();
+        const auto readDelay = static_cast<CoeffType> (std::ceil (delay));
+        const auto fraction = readDelay - delay;
+        const auto fraction2 = fraction * fraction;
+        const auto fraction3 = fraction2 * fraction;
+
+        const auto weight0 = static_cast<CoeffType> (-0.5) * fraction + fraction2 - static_cast<CoeffType> (0.5) * fraction3;
+        const auto weight1 = static_cast<CoeffType> (1) - static_cast<CoeffType> (2.5) * fraction2 + static_cast<CoeffType> (1.5) * fraction3;
+        const auto weight2 = static_cast<CoeffType> (0.5) * fraction + static_cast<CoeffType> (2) * fraction2 - static_cast<CoeffType> (1.5) * fraction3;
+        const auto weight3 = static_cast<CoeffType> (-0.5) * fraction2 + static_cast<CoeffType> (0.5) * fraction3;
+        const auto omega = frequencyToAngular (responseFrequency, static_cast<CoeffType> (this->sampleRate));
+
+        return weight0 * polar (static_cast<CoeffType> (1), -omega * (readDelay + static_cast<CoeffType> (1)))
+             + weight1 * polar (static_cast<CoeffType> (1), -omega * readDelay)
+             + weight2 * polar (static_cast<CoeffType> (1), -omega * (readDelay - static_cast<CoeffType> (1)))
+             + weight3 * polar (static_cast<CoeffType> (1), -omega * (readDelay - static_cast<CoeffType> (2)));
+    }
+
     //==============================================================================
     std::size_t delayLineSizeInSamples = defaultDelayLineSize;
     std::size_t delayLineMask = defaultDelayLineSize - 1;

modules/yup_dsp/yup_dsp.h

@@ -139,10 +139,6 @@
 #include "frequency/yup_FFTProcessor.h"
 #include "frequency/yup_SpectrumAnalyzerState.h"
 
-// Wavetable and antialiased oscillators (before FilterDesigner for AA-IIR2 coefficient types)
-#include "oscillators/yup_Wavetable.h"
-#include "oscillators/yup_AaIir2Oscillator.h"
-
 // Base filter interfaces and common structures
 #include "base/yup_FilterMode.h"
 #include "base/yup_FilterBase.h"

tests/yup_dsp/yup_CombFilter.cpp

@@ -124,3 +124,22 @@ TEST (CombFilterTests, ComplexResponseIsFinite)
     EXPECT_TRUE (std::isfinite (response.real()));
     EXPECT_TRUE (std::isfinite (response.imag()));
 }
+
+TEST (CombFilterTests, ComplexResponseMatchesIntegerDelayFeedForwardComb)
+{
+    CombFilter<double> filter;
+    filter.prepare (sampleRate, blockSize);
+    filter.setParameters (sampleRate / 8.0, 0.0, 0.0, sampleRate);
+
+    EXPECT_NEAR (std::abs (filter.getComplexResponse (0.0)), 1.5, 1e-9);
+    EXPECT_NEAR (std::abs (filter.getComplexResponse (sampleRate / 16.0)), 0.5, 1e-9);
+}
+
+TEST (CombFilterTests, ComplexResponseIncludesFeedbackResonance)
+{
+    CombFilter<double> filter;
+    filter.prepare (sampleRate, blockSize);
+    filter.setParameters (sampleRate / 8.0, 0.81, 0.0, sampleRate);
+
+    EXPECT_NEAR (std::abs (filter.getComplexResponse (sampleRate / 8.0)), 6.0, 1e-9);
+}