Python SciPy signal.freqz用法及代码示例

本文简要介绍 python 语言中 scipy.signal.freqz 的用法。

用法:

scipy.signal.freqz(b, a=1, worN=512, whole=False, plot=None, fs=6.283185307179586, include_nyquist=False)#

计算数字滤波器的频率响应。

给定数字滤波器的M-order 分子 b 和 N-order 分母 a，计算其频率响应：

jw                 -jw              -jwM
   jw    B(e  )    b[0] + b[1]e    + ... + b[M]e
H(e  ) = ------ = -----------------------------------
            jw                 -jw              -jwN
         A(e  )    a[0] + a[1]e    + ... + a[N]e

参数 ：：

b： array_like

线性滤波器的分子。如果b维度大于 1，假设系数存储在第一个维度中，并且b.shape[1:],a.shape[1:]，并且频率阵列的形状必须与广播兼容。

a： array_like

线性滤波器的分母。如果b维度大于 1，假设系数存储在第一个维度中，并且b.shape[1:],a.shape[1:]，并且频率阵列的形状必须与广播兼容。

worN： {无，int，数组}，可选

如果是单个整数，则以那么多频率进行计算(默认为 N=512)。这是一个方便的替代方案：

np.linspace(0, fs if whole else fs/2, N, endpoint=include_nyquist)

对 FFT 计算使用快速的数字可以提高计算速度(请参阅注释)。

如果是 数组，计算给定频率的响应。这些与 fs 的单位相同。

whole： 布尔型，可选

通常，频率的计算范围是从 0 到奈奎斯特频率 fs/2(upper-half of unit-circle)。如果整体为真，则计算从 0 到 fs 的频率。如果 worN 是 数组，则忽略。

plot： 可调用的

带有两个参数的可调用对象。如果给定，返回参数w和h被传递给情节。用于绘制内部的频率响应freqz.

fs： 浮点数，可选

数字系统的采样频率。默认为 2*pi 弧度/样本(所以 w 是从 0 到 pi)。

include_nyquist： 布尔型，可选

如果整体为 False 且 worN 为整数，则将 include_nyquist 设置为 True 将包括最后一个频率(奈奎斯特频率)，否则将被忽略。

返回 ：：

w： ndarray

计算 h 的频率，单位与 fs 相同。默认情况下，w 被归一化为范围 [0, pi)(弧度/样本)。

h： ndarray

频率响应，作为复数。

注意：

使用 Matplotlib 的matplotlib.pyplot.plot作为可调用函数阴谋产生意想不到的结果，因为这绘制了复数传递函数的实部，而不是幅度。尝试lambda w, h: plot(w, np.abs(h)).

当满足以下条件时，通过 (R)FFT 的直接计算用于计算频率响应：

1. 给 worN 一个整数值。

2. 用过通过 FFT 计算速度很快(即，scipy.fft.next_fast_len等于用过)。

3. 分母系数是单个值(a.shape[0] == 1)。

4. 用过至少与分子系数一样长 (worN >= b.shape[0])。

5. 如果 b.ndim > 1 ，那么 b.shape[-1] == 1 。

对于长 FIR 滤波器，FFT 方法可以具有更低的误差并且比等效的直接多项式计算要快得多。

例子：

>>> from scipy import signal
>>> import numpy as np
>>> b = signal.firwin(80, 0.5, window=('kaiser', 8))
>>> w, h = signal.freqz(b)

>>> import matplotlib.pyplot as plt
>>> fig, ax1 = plt.subplots()
>>> ax1.set_title('Digital filter frequency response')

>>> ax1.plot(w, 20 * np.log10(abs(h)), 'b')
>>> ax1.set_ylabel('Amplitude [dB]', color='b')
>>> ax1.set_xlabel('Frequency [rad/sample]')

>>> ax2 = ax1.twinx()
>>> angles = np.unwrap(np.angle(h))
>>> ax2.plot(w, angles, 'g')
>>> ax2.set_ylabel('Angle (radians)', color='g')
>>> ax2.grid(True)
>>> ax2.axis('tight')
>>> plt.show()

广播示例

假设我们有两个 FIR 滤波器，它们的系数存储在形状为 (2, 25) 的数组的行中。对于这个演示，我们将使用随机数据：

>>> rng = np.random.default_rng()
>>> b = rng.random((2, 25))

要通过一次调用 freqz 来计算这两个滤波器的频率响应，我们必须传入 b.T ，因为 freqz 期望第一个轴保存系数。然后，我们必须将形状扩展为长度为 1 的简单维度，以允许使用频率阵列进行广播。也就是说，我们传入 b.T[..., np.newaxis] ，其形状为 (25, 2, 1)：

>>> w, h = signal.freqz(b.T[..., np.newaxis], worN=1024)
>>> w.shape
(1024,)
>>> h.shape
(2, 1024)

现在，假设我们有两个传递函数，具有相同的分子系数b = [0.5, 0.5].两个分母的系数存储在二维数组的第一维中a：

a = [   1      1  ]
    [ -0.25, -0.5 ]

>>> b = np.array([0.5, 0.5])
>>> a = np.array([[1, 1], [-0.25, -0.5]])

仅有的a大于一维。为了使其与频率广播兼容，我们在调用中将其扩展为一个微不足道的维度freqz：

>>> w, h = signal.freqz(b, a[..., np.newaxis], worN=1024)
>>> w.shape
(1024,)
>>> h.shape
(2, 1024)