信息处理 - 波形数字滤波器 Bridge-T 谐振器实现，给出了预期的截止频率，但增益和滚降不正确 - 吾爱随笔录

波形数字滤波器 Bridge-T 谐振器实现，给出了预期的截止频率，但增益和滚降不正确

信息处理过滤器 Python 数字的矩阵海浪

2022-02-23 18:07:09

尝试实现本出版物图 5(a) 中的 WDF 。理想运算放大器实现的响应由图 7 中的黑色曲线给出：

这是我尝试自己实现电路时得到的图：

我使用了与论文中给出的相同的散射矩阵：

但是，论文并没有说明使用什么端口电阻来适配 Rc，所以我使用 SymPy 使用 R 型适配器的散射矩阵方程来找到 S：

以及电路的矩阵X，如图 6 所示：

这是我的代码：

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)

通过将 Rc 的端口电阻设置为 RbRe/(Rb + Rd + Re)，并使用 fs=44.1kHz，S 成为预期的散射矩阵。

*我发现的散射矩阵实际上与论文中的几乎相同，但是 col[0],row[1:5] 有符号翻转。我尝试使用任一 S 矩阵运行我的电路代码，最终得到相同的响应。

我使用了此笔记本中的 WDF 库，并进行了一些修改。这是 WDF 电路组件本身的代码：

#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

然后我使用论文中的值来实例化这些组件，形成电路并返回脉冲响应：

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

以及绘制输出的代码：

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)

在论文中的方法和我的方法之间，我没有看到任何其他明显的差异，而且我用于这些电路的 Python 库与我尝试实现的其他电路非常一致。

1个回答

我能够运行您的代码，重现您的结果，并且该表给出的值似乎被截断得太高了。我修改了代码来计算你硬编码的散射矩阵，结果更接近。我还修改了它以使用 0.01 的脉冲，因为这就是他们在链接文章中使用的。抱歉缺少标记，但频率图上的 x 轴从 1.6kHz 到 3.0kHz 以匹配提供的图。您还可以看到脉冲响应具有相当长的衰减。这是我得到的散射矩阵：

[[ 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]]

为了完整起见，我包括了科学格式，如果阅读起来有点困难，我深表歉意，但它们与提供的表格相匹配。

几点注意事项：

第一列与您观察到的相同。这会反转来自电压源的信号，所以没关系，但我很好奇那里有什么交易，因为它使用了他们提供的控制方程。
那个 WDF 库有点奇怪。文档很棒，但它的使用方式有点奇怪。首先，波浪向上/向下的事情很奇怪。您需要先计算无反射端口，然后使用它来获取其他端口。向上/向下的东西可以解决问题，但效率很低。其次，它总是有一个单一的样本延迟。入射波被记录在循环的顶部，然后它们被重新计算，但是当你得到输出时，它仍然具有以前的值，我认为这很烦人。
最后，幅度值不匹配，这很糟糕，但它们可疑地偏离了 20，这是 10 的因数。所以也许他们打算写“100mV”而不是 10，不知道。我会对电路进行论文分析，但我有点花在它上面。

希望有帮助！干杯。

其它你可能感兴趣的问题

上一篇验证 FPGA 信号处理链中的中间输出下一篇以图形方式理解 LTI 系统