Wave Digital Filter Bridge-T Resonator implementation, gives expected cutoff frequency but incorrect gain and roll-off

Question

Trying to implement the WDF in Fig 5(a) of this publication. The response of the ideal op-amp implementation is given by the black curve in Fig 7:

Here is the plot I get when trying to implement the circuit myself:

I used the same scattering matrix as given in the paper:

However, the paper doesn't say what port resistance was used to adapt Rc, so I used SymPy to find S using the scattering matrix equation for R-type adaptors:

and the matrix X for the circuit, given in Fig 6:

Here is my code:

from __future__ import division
from IPython.display import display
from sympy import *
#from math import sqrt

init_printing()

RA = symbols('Ra')
RB = symbols('Rb')
RC = symbols('Rc')
RD = symbols('Rd')
RE = symbols('Re')
RF = symbols('Rf')
GA = 1/RA
GB = 1/RB
GC = 1/RC
GD = 1/RD
GE = 1/RE
GF = 1/RF
I = eye(6)
zIz = zeros(6, 10).row_join(I).row_join(zeros(6,1))
R = Matrix([
        [RA, 0, 0, 0, 0, 0],
        [0, RB, 0, 0, 0, 0],
        [0, 0, RC, 0, 0, 0],
        [0, 0, 0, RD, 0, 0],
        [0, 0, 0, 0, RE, 0],
        [0, 0, 0, 0, 0, RF]
])
#datum node 1 omitted to make X invertible
X = Matrix([
    [GA, 0, 0, 0, -GA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, -1, 0, 0, 0],
    [0, 0, GB, 0, 0, -GB, 0, 0, 0, 0, 0, 0, -1, 0, -1, 0, 0],
    [0, 0, 0, GD+GE+GF, 0, 0, 0, -GD, -GE, -GF, 0, 0, 0, 0, 0, 0, 1],
    [-GA, 0, 0, 0, GA, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
    [0, 0, -GB, 0, 0, GB, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, GC, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
    [0, 0, 0, -GD, 0, 0, 0, GD, 0, 0, 0, 0, 0, 1, 0, 0, 0],
    [0, 0, 0, -GE, 0, 0, 0, 0, GE, 0, 0, 0, 0, 0, 1, 0, 0],
    [0, 0, 0, -GF, 0, 0, 0, 0, 0, GF, 0, 0, 0, 0, 0, 1, 0],
    [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, -1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, -1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
    [-1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
])
X_inv = X.inv()
S = I + 2*R*zIz*X_inv*zIz.T
S = S.applyfunc(simplify)

By setting Rc's port resistance to RbRe/(Rb + Rd + Re), and using fs=44.1kHz, S becomes the expected scattering matrix.

*the scattering matrix I found is actually almost identical to the one in the paper, but col[0],row[1:5] have signs flipped. I've tried running my circuit code with either S matrix and I end up getting the same response.

I used a WDF library from this notebook with a few modifications. Here's the code for the WDF circuit components themselves:

#One-port
class WDFOnePort(object):
    def __init__(self):
        self.a, self.b = 0, 0

    # v = (a + b)/2
    def wave_to_voltage(self):
        voltage = (self.a + self.b)/2
        return voltage

#Resistor
class Resistor(WDFOnePort):
    def __init__(self, R):
        WDFOnePort.__init__(self)
        self.Rp = R

    def get_reflected_wave(self, a):
        self.a = a
        self.b = 0 # to avoid delay-free loop, reflected wave is always zero
        return self.b

class RootResistor(WDFOnePort):
  def __init__(self, R, Rp):
    WDFOnePort.__init__(self)
    self.R = R
    self.Rp = Rp

  # instantaneous reflection permitted
  def get_reflected_wave(self, a):
    self.b = a*(self.R - self.Rp)/(self.R + self.Rp)
    self.a = a
    return self.b

# Capacitor
class Capacitor(WDFOnePort):
    def __init__(self, C, fs=44100):
        WDFOnePort.__init__(self)
        print("WDF Capacitor fs: " + str(fs))
        self.Rp = 1/(2*fs*C)

    def get_reflected_wave(self, a):
        self.b = self.a
        self.a = a
        return self.b

    def set_incident_wave(self, a):
        self.a = a

# Resistive Voltage Source
class ResistiveVoltageSource(WDFOnePort):
    def __init__(self, Rs):
        WDFOnePort.__init__(self)
        self.Rp = Rs

    def get_reflected_wave(self, a, vs=0):
        self.a = a
        self.b = vs
        return self.b

Then I used the values from the paper to instantiate these components, form the circuit and return the impulse response:

def nullorBasedBridgedTResonator_WDF(steps=2**14):
  input = np.zeros(steps)
  input[0] = 1 #unit impulse
  output = np.zeros(steps)
  fs = 44100

  #https://pureadmin.qub.ac.uk/ws/portalfiles/portal/158209014/1570255463.pdf
  S = np.array([
    [1, 0, 0, 0, 0, 0],
    [-1, -0.001, -1, 0.001, -0.001, 0],
    [-1, -0.999, 0, -0.001, 0.001, 0],
    [882, 882.998, 882, 0.002, 0.998, 0],
    [883, 881.999, 883, 1.001, -0.001, 0],
    [884, 882.998, 882, 1.002, 0.998, -1]
  ])

  Cap_value = 1e-9
  Ra_value = 1
  Rb_value = 1/(2*Cap_value*fs)
  Rc_value = 500
  Rd_value = 10000000
  Re_value = Rb_value
  Rf_value = 10000

  # root element RC = R1
  Rc_port_value = Rb_value*Re_value/(Rb_value + Rd_value + Re_value)
  print(Rc_port_value)

  Vin = ResistiveVoltageSource(Ra_value)
  C1 = Capacitor(Cap_value)
  R1 = RootResistor(Rc_value, Rc_port_value)
  R2 = Resistor(Rd_value)
  C2 = Capacitor(Cap_value)
  RL = Resistor(Rf_value)

  b = [0, 0, 0, 0, 0, 0]
  a = [0, 0, 0, 0, 0, 0]

  for i in range(steps):
    #gather leaf node incident waves
    a[5] = RL.get_reflected_wave(b[5])
    a[4] = C2.get_reflected_wave(b[4])
    a[3] = R2.get_reflected_wave(0) # don't care
    a[1] = C1.get_reflected_wave(b[1])
    a[0] = Vin.get_reflected_wave(0, input[i]) # don't care

    # wave-up
    b = np.dot(S, a)
    # root, instantaneous reflection
    a[2] = R1.get_reflected_wave(b[2])
    #wave down
    b = np.dot(S,a)

    output[i] = RL.wave_to_voltage()

    C1.set_incident_wave(b[1])
    C2.set_incident_wave(b[4])
  return output

And the code to plot the output:

x = nullorBasedBridgeTResonator_WDF()
f, h = signal.freqz(x, 1, worN=4096, fs=44100)
H = 20*np.log10(np.abs(h))
ax_mag.semilogx(f, H, label=label)

There aren't any other obvious discrepancies I can see between the methods in the paper and my methods, and the Python library I've used for these circuits has been pretty consistent with other circuits I've tried implementing.

Answer 1

I was able to run your code, reproduce your results, and it seems like the values given by that table are too highly truncated. I modified the code to calculate the scattering matrix where you have it hard coded in, and the results are much closer. I also modified it to use an impulse of 0.01, as that's what they used in the linked article. Apologies for the lack of markers, but the x axis on the frequency plot goes from 1.6kHz to 3.0kHz to match the supplied plot. You can also see the impulse response has a decently long decay. Here's the scattering matrix I got:

[[ 1.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00]
 [ 1.00000000e+00 -1.13122172e-03 -1.00000000e+00  1.13122172e-03
  -1.13122172e-03  0.00000000e+00]
 [ 1.00000000e+00 -9.98868778e-01  2.22044605e-16 -1.13122172e-03
   1.13122172e-03  0.00000000e+00]
 [-8.82000000e+02  8.82997738e+02  8.82000000e+02  2.26244344e-03
   9.97737557e-01  0.00000000e+00]
 [-8.83000000e+02  8.81998869e+02  8.83000000e+02  1.00113122e+00
  -1.13122172e-03  0.00000000e+00]
 [-8.84000000e+02  8.82997738e+02  8.82000000e+02  1.00226244e+00
   9.97737557e-01 -1.00000000e+00]]

I included the scientific format for the sake of completeness, my apologies if it's a bit tough to read, but they match the provided table.

A couple of notes:

1. The first column is inverted same as what you observed. This would invert the signal coming out of the voltage source, so it doesn't matter much but I am curious what the deal is there, since it's using their provided controlling equations.
2. That WDF library is a little wonky. The docs are great, but the way its used is a bit odd. First, the wave up/down thing is weird. You need to calculate the reflection free port first, then use that to get the other ones. The up/down thing does the trick, but it's pretty inefficient. Second, it always has a single sample delay. The incident waves are clocked in at the top of the loop, then they get recalculated, but when you get the output it still has the previous value, which I think is annoying.
3. Lastly, the magnitude values don't match, which is a bummer, but they are suspiciously off by 20, which is a factor of 10. So maybe they meant to write '100mV' instead of 10, dunno. I'd do a paper analysis on the circuit, but I'm kind of spent on it.

Hope that helps! Cheers.