陷波滤波器的脉冲响应

Question

我正在尝试为赋值实现形式为 W(z) = N(z)/D(z) 的陷波滤波器，其中 N=a+bz+cz^2 和 D=1+Bz +Cz^2。为此，我将第一个滤波器应用为直接 3 项卷积，然后将第二个滤波器应用为 3 项逆卷积。

为了测试这个过滤器，我创建了一个离散化的 delta 函数并将其作为输入传递给函数。

我用于过滤器和测试的代码如下：

import numpy as np
import matplotlib.pyplot as plt

fs=12
f0=1
M=1.05
epsilon = 0.05

# Define rational filter:
def ratFilter(N, D, x):
    '''
    Apply two filters in succession to x
    :param N: 3-tuple parameters for numerator filter
    :param D: 3-tuple parameters for denominator inverse filter
    :param x:
    :return: y
    '''
    y = np.zeros(np.size(x))
    n = np.zeros(np.size(x)+np.size(D)) #middle value
    # apply first filter:
    for i in range(np.size(n)):
        for k in range(0, i+1):
            try:
                n[i]+=N[k]*x[i-k]
            except:
                pass
    for i in range(np.size(y)):
        y[i]=n[i]/D[0]
        for k in range(1, i+1):
            try:
                y[i]-=D[k]/D[0]*y[i-k]
            except:
                pass
    return y

# Impulse response:
delta = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

N = [0.952, -1.650, 0.952]
D = [1, -1.650, 0.907]

output = ratFilter(N, D, delta)

plt.plot(output)
plt.title("impulse response")
plt.show()

生成的图表具有以下形式，我认为这是不正确的，因为我认识的其他人在开始时有一个急剧上升的尖峰，然后是线性增加：

脉冲响应输出：

必须做的另一件事是对脉冲响应进行傅里叶变换以获得频率响应|W(f)|，但直觉上我不确定傅里叶变换和z之间的link转换，因此如果有任何特殊程序而不是将 np.fft.fft 应用于输出的任何结果。

Answer 1

我不确定是什么让您认为同学的结果正确而您的不正确，但看起来您的实施产生的结果与 scipy.signal.lfilter 的标准实施产生的结果相同输入。

就|W(f)|而言，z变换与np.fft.fft的关系可以从以下观察得到：

Discrete Time Fourier Transform (DTFT) 对应于在复平面的单位圆上计算的 z 变换。
在频域中，Discrete Fourier Transform (DFT) 对应于连续 DTFT 频谱的均匀采样版本。
快速傅立叶变换（由 np.fft.fft 计算）是计算离散傅立叶变换的有效算法。

因此，计算脉冲响应的 FFT 将为您提供该序列在离散频率值处的频谱。这种情况下的主要问题是，由于滤波器的递归性质，脉冲响应具有无限长度。因此，您需要确保考虑足够长的脉冲响应部分，截断误差可以忽略不计。

提示：作为替代方案，您应该考虑在评估 N 的 FFT 与 D 的 FFT 的比率时会发生什么。

陷波滤波器的脉冲响应

Impulse response of notch filter

python

filtering

fft

filter