如何手动绘制Butterworth过滤器？

Question

我想问一些关于Python中Butterworth过滤器的问题。我知道这是如何绘制低通Butterworth过滤器的图形的：

  from scipy import signal
  import matplotlib.pyplot as plt

  b,a=signal.butter(N,fCut,'low',analog=True)
  w,h=signal.freqs(b,a)
  plt.plot(w,(abs(h)))

其中N是阶数，fCut是截止频率。 但是我在尝试手动操作时遇到了问题，没有使用'signal.freqs'。 这意味着我必须计算模拟传递函数： H（s）=1/Σ(s-sk），k=1..N，sk=e^（（j*pi）（（2*k+N-1）/2*N）） 在python中有这样做的方法吗

Answer 1

下面是与公式对应的Python代码：

import numpy as np
def butterTF(s, N, fCut):
    sk = np.exp(0.5j*np.pi*(2*np.arange(N) + N + 1)/N)
    return 1/np.multiply.reduce(np.subtract.outer(1.j*s/fCut, sk), axis=1)

为了测试它：

from scipy import signal
import matplotlib.pyplot as plt

b,a=signal.butter(N,fCut,'low',analog=True)
w,h=signal.freqs(b,a)
plt.figure()
plt.plot(w, np.abs(h))
plt.plot(w, np.abs(butterTF(w, N, fCut)), '.')
plt.grid()
plt.show()

一些解释：

计算sk很简单，只是我们需要稍微修改一下公式，因为np.arange(N)给出的是0, 1, ... , N-1，而不是1, 2, ... , N

现在棘手的部分：np.subtract.outer形成一个矩阵len(s) x N，其元素i, k是1.j*s[i]/fCut - sk[k]的结果。可能不是一种非常有效的内存计算方法，特别是当N很大时，它是一种相对快速的计算方法

np.multiply.reduce(..., axis=1)计算矩阵每行的乘积