3
$\begingroup$

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:

enter image description here

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

enter image description here

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

enter image description here

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:

enter image description here

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

enter image description here

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)

enter image description here

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.

$\endgroup$
2
  • $\begingroup$ The Q is 72? That should have a really long impulse response. Have you tried plotting the impulse response to make sure you’re getting enough of it? Similarly, is worN big enough to get the whole thing? $\endgroup$
    – Dan Szabo
    Apr 20 at 13:44
  • $\begingroup$ So I just plotted the impulse response and not only did I capture the complete response, but it settles to zero very quickly (after ~100 steps). I don't know much about the relationship between Q and impulse response length but that does seem odd to me $\endgroup$
    – Nick Nagy
    Apr 22 at 3:45
1
$\begingroup$

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. WDF Test 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.

$\endgroup$
1
  • 1
    $\begingroup$ Wow. Sure enough, when I substituted in your (un-rounded) scattering matrix, the expected results - or at least, your results - came through. I did not expect rounding to have that big of an effect. The WDF library definitely feels a little off to me, but like you said, the documentation is nice, and at least when testing impulse/frequency responses (ie not in real time) I don't mind the redundancies. Thanks for the help! $\endgroup$
    – Nick Nagy
    Apr 22 at 22:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.