陷波滤波器的脉冲响应

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变换之间的联系，因此如果有是任何特殊程序，而不是将 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 的比率时会发生什么。