Python中的IIR响应

6 投票

3 回答

7340 浏览

提问于 2025-04-18 06:31

我有一个新问题，跟这个帖子直接相关。在Python中，我创建了一个二阶的IIR带通滤波器，具有特定的特性[以下代码是故意使用习惯用法的]：

fs = 40e6           # 40 MHz f sample frequency
fc = 1e6/fs         # 1 MHz center cutoff
BW = 20e3/fs        # 20 kHz bandwidth 
fl = (fc - BW/2)/fs # 0.99 MHz f lower cutoff
fh = (fc + BW/2)/fs # 1.01 MHz f higher cutoff

这段代码给出了滤波器的系数：

R  = 1 - (3*BW)
K  = (1 - 2*R*np.cos(2*np.pi*fc) + (R*R)) / (2 - 2*np.cos(2*np.pi*fc))

a0 = 1 - K                       # a0 =  0.00140
a1 = 2*(K-R)*np.cos(2*np.pi*fc)  # a1 =  0.00018 
a2 = (R*R) - K                   # a2 = -0.00158

b1 = 2*R*np.cos(2*np.pi*fc)      # b1 =  1.97241
b2 = -(R*R)                      # b2 = -0.99700

正如ukrutt在之前的帖子中提到的，我使用了scipy.signal.freqz，但遗憾的是没有得到我想要的响应。不过，我相信这个滤波器是按预期工作的（代码在下面）。这是freqz的结果：

在这里输入图片描述

我的问题是：我该如何生成一个更接近预期响应的图表？

代码：

a = [0.0014086232031758072, 0.00018050359364826498, -0.001589126796824103]
b = [1.9724136161684902, -0.9970022500000001]

w,h  = signal.freqz(a, b)
h_dB = 20 * np.log10(np.abs(h))
plt.plot(w/np.max(w),h_dB)
plt.grid()

3 个回答

-1

问题在于，signal.freqz 返回的是半圆上的点……所以你不能直接扩展到更大的 x 范围，除非通过 signal.freqz 来实现。我试着研究了一下，发现可以给 signal.freqz 传递 whole=True，这样你会得到上面那样的结果，但是关于负 x 的镜像。所以这不是解决办法。不过，还有另一个参数可以让你传入想要 signal.freqz 计算的 x 点……我试着用 np.arange(-5., 5., 0.1)，结果看起来和右边的图完全不一样——看起来像是原始图的多个反射。这让我在想……也许你右边的图和左边的图有不同的坐标轴？具体来说，一个是角频率，另一个是普通的频率？

进一步研究后，发现 signal.freqz 返回的是 w,h，其中 w 是以弧度/样本为单位的归一化角频率。所以，你在代码中不需要通过 np.max(w) 来进行归一化来生成图。不过，这仍然没有解决问题。你右边的图似乎是以 fc 为单位，而 fc 是以 MHz 为单位（例如 1/样本）。

所以，要让左边的图和右边的图匹配，我想这意味着你需要“去归一化”你的 x 轴，然后需要把角频率的单位转换成 MHz。

或者更可能的是，使用一个不同于 signal.freqz 的函数。

回答于 2025-04-18 由 Python大师

分享举报

我觉得问题不在于你绘制响应的方式，而是在于你选择的滤波器。你现在使用的是一个低阶的IIR滤波器，这样会导致你得到的响应范围很窄。我认为你要么需要使用一个高阶的滤波器，要么就放宽一些限制。

举个例子，下面这个使用了巴特沃斯滤波器，它是一个IIR滤波器，能够提供一个更接近你想要的响应形状。当然，要达到你期望的滤波特性，还需要做更多的工作。

    b, a = signal.butter(4, [1.0/4-1.0/2e2,1.0/4+1.0/2e2], 'bandpass', analog=False)
    w, h = signal.freqs(b, a)

    import matplotlib.pyplot as plt
    fig = plt.figure()
    plt.title('Digital filter frequency response')
    ax1 = fig.add_subplot(111)
    plt.semilogy(w, np.abs(h), 'b')

    plt.ylabel('Amplitude (dB)', color='b')
    plt.xlabel('Frequency (rad/sample)')
    ax2 = ax1.twinx()

    angles = np.unwrap(np.angle(h))
    plt.plot(w, angles, 'g')
    plt.ylabel('Angle (radians)', color='g')
    plt.grid()
    plt.axis('tight')
    plt.show()

这样就得到了：

滤波器响应图

回答于 2025-04-18 由 Python大师

分享举报

用线性的x轴比例，你不会看到什么好看的东西。我不太懂numpy，但我对matlab比较熟悉，matlab里有一些函数可以用来做对数图。你可以试试用x对数比例来绘图，代码如下：

import matplotlib.pyplot  as pyplot
fig = pyplot.figure()
ax = fig.add_subplot(2,1,1)    
line, = ax.plot(w/np.max(w), h_dB, color='blue', lw=2)
ax.set_xscale('log')

show()

顺便说一下，我还没测试过，因为我没有安装python :(

编辑：

我试着在matlab里模拟一个巴特沃斯滤波器，做了一个4阶的IIR滤波器和一个20阶的IIR滤波器。

%!/usr/local/bin/matlab

%% Inputs
fs = 40e6;
fc = 1e6;
BW = 20e3;
fl = (fc - BW/2);
fh = (fc + BW/2);

%% Build bandpass filter IIR Butterworth order 4
N   = 4;        % Filter Order
h  = fdesign.bandpass('N,F3dB1,F3dB2', N, fl, fh, fs);
Hd1 = design(h, 'butter');

%% Build bandpass filter IIR Butterworth order 50
N   = 20;        % Filter Order
h  = fdesign.bandpass('N,F3dB1,F3dB2', N, fl, fh, fs);
Hd2 = design(h, 'butter');

%% Compare
fvtool(Hd1,Hd2);

两个滤波器

稍微放大一下

这是第一个滤波器的系数A和B：

FilterStructure: 'Direct-Form II Transposed'                                               
A: [2.46193004641106e-06 0 -4.92386009282212e-06 0 2.46193004641106e-06]     
B: [1 -3.94637005453608 5.88902106889851 -3.93761314372475 0.995566972065978]

如果有时间的话，我会尝试用numpy做同样的事情！

回答于 2025-04-18 由 Python大师

分享举报

Python中的IIR响应

3 个回答

撰写回答