Modulating coefficients in allpass filter implementation

dsp/filters/allpass.py

@@ -11,24 +11,47 @@
 
 class AllpassFilter:
     """
-    First-order allpass filter class with an additional constraint that the
-    b0 coefficient always equals the a1 coefficient.
+    First-order allpass filter class with modulating coefficients. Also
+    maintains the constraint that the b0 coefficient always equals the a1
+    coefficient.
 
     Parameters
-    a1 : The initial value for the a1 coefficient.
+    offset      : [0.0, 1.0] : The initial value for the a1 coefficient, and the
+                    value about which the modulating operator oscillates.
+    amplitude   : [0.0, 1.0] : The amplitude of the modulating signal.
+                    Effectively the "amount" of modulation.
+    rate        : [0, 96000] : The rate (Hz) of the modulating signal.
     """
 
-    def __init__(self, a1):
+    def __init__(self, offset, amplitude, rate):
+        # The maximum and minimum values allowed in the modulating signal, so as
+        # to keep the DC delay within a reasonable range, and to avoid the pole
+        # zero cancellation at delta = 0.0.
+        self._mmax = 0.82
+        self._mmin = -0.05
+
+        self._range = self._mmax - self._mmin
+        self._offset = offset * self._range
+        self._amp = min(amplitude, self._mmax - self._offset)
+        self._rate = rate
+
+        # Initial coefficient values
         self.a0 = 1.0
-        self.a1 = a1
-        self.b0 = self.a1
         self.b1 = 1.0
+        self.b0 = self.a1 = self._mmin + self._offset
 
         # Zero pad the initial input and output buffers. The filter will be in
         # its transient state for 1 sample before reaching steady state.
         self._x = 0.0
         self._y = 0.0
 
+    def update(self, t):
+        # Note: this implementation works only for large continuous blocks; if
+        # instead it operated on discrete blocks in a continuous stream, the `t`
+        # input would yield discontinuities.
+        return self._mmin + self._offset + \
+                self._amp * np.sin(2.0 * np.pi * self._rate * t / 44100.)
+
     def process_block(self, input_buffer, output_buffer):
         _len = input_buffer.size
         _x = np.insert(input_buffer, 0, self._x)
@@ -40,6 +63,7 @@ def process_block(self, input_buffer, output_buffer):
                    - (self.a1 / self.a0) * _y[i - 1] \
 
             _y[i] = sample
+            self.b0 = self.a1 = self.update(i)
 
         np.copyto(output_buffer, _y[1:])
         self._x = input_buffer[-1]
@@ -84,7 +108,7 @@ def plot(self, ax1, ax2, color='c', alpha=1.0):
     _, (ax1, ax2) = plt.subplots(2, sharex=True)
 
     for i in np.linspace(0.0, 0.8, 4):
-        apf = AllpassFilter(i)
+        apf = AllpassFilter(0.5, 1.0, 64000)
         apf.plot(ax1, ax2, 'c', i / 0.8)
 
     plt.show()
@@ -105,7 +129,7 @@ def plot(self, ax1, ax2, color='c', alpha=1.0):
     RealTimeResamplingOscillator(WaveType.SAWTOOTH, 43.65, 3.0, 1.0).render(rs)
 
     ap = np.zeros(size, dtype='d')
-    apf = AllpassFilter(0.68)
+    apf = AllpassFilter(0.5, 1.0, 64000)
     apf.process_block(ss, ap)
 
     plt.figure()