SciPy-1-12-中文文档-二十九-

201 阅读35分钟

SciPy 1.12 中文文档(二十九)

原文:docs.scipy.org/doc/scipy-1.12.0/index.html

scipy.signal.firwin

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.firwin.html#scipy.signal.firwin

scipy.signal.firwin(numtaps, cutoff, *, width=None, window='hamming', pass_zero=True, scale=True, nyq=<object object>, fs=None)

使用窗口方法设计 FIR 滤波器。

此函数计算有限冲激响应滤波器的系数。滤波器将具有线性相位;如果numtaps为奇数则为 Type I,如果numtaps为偶数则为 Type II。

Type II 滤波器在奈奎斯特频率处始终具有零响应,因此如果使用numtaps为偶数且其通带右端在奈奎斯特频率处的情况下调用 firwin,则会引发 ValueError 异常。

参数:

numtaps整数

滤波器的长度(系数数量,即滤波器阶数+1)。如果通带包含奈奎斯特频率,则numtaps必须为奇数。

cutoff浮点数或 1-D 数组

滤波器的截止频率(以与fs相同的单位表示)或截止频率数组(即带边缘)。在后一种情况下,cutoff中的频率应为正且单调增加,在 0 和fs/2之间不应包括值 0 和fs/2

width浮点数或 None,可选

如果width不为 None,则假定其为过渡区域的大致宽度(以fs的相同单位表示),用于 Kaiser FIR 滤波器设计。在这种情况下,window参数将被忽略。

window字符串或字符串和参数值的元组,可选

所需使用的窗口。有关窗口和所需参数的列表,请参阅scipy.signal.get_window

pass_zero{True, False, ‘bandpass’, ‘lowpass’, ‘highpass’, ‘bandstop’},可选

如果为 True,则频率为 0 时的增益(即“直流增益”)为 1。如果为 False,则直流增益为 0。也可以是所需滤波器类型的字符串参数(相当于btype在 IIR 设计函数中的参数)。

从版本 1.3.0 开始支持字符串参数。

scale布尔值,可选

设置为 True 以使系数按比例缩放,以便频率响应在某个频率上完全为单位。该频率可以是:

  • 如果第一个通带从 0 开始(即 pass_zero 为 True),则直流(DC)为 0。

  • fs/2(奈奎斯特频率),如果第一个通带结束于fs/2(即滤波器是单通带高通滤波器);否则为第一个通带的中心

nyq浮点数,可选,已弃用

这是奈奎斯特频率。cutoff中的每个频率必须介于 0 和nyq之间。默认为 1。

自版本 1.0.0 起不推荐使用:firwin 关键字参数nyq已弃用,推荐使用fs,并将在 SciPy 1.14.0 中移除。

fs浮点数,可选

信号的采样频率。cutoff中的每个频率必须介于 0 和fs/2之间。默认为 2。

返回:

h(numtaps,)ndarray

长度为numtaps的 FIR 滤波器系数。

引发:

ValueError

如果cutoff中的任何值小于等于 0 或大于等于fs/2,如果cutoff的值不是严格单调递增,或者numtaps是偶数但通带包含奈奎斯特频率。

参见

firwin2

firls

minimum_phase

remez

示例

低通从 0 到 f:

>>> from scipy import signal
>>> numtaps = 3
>>> f = 0.1
>>> signal.firwin(numtaps, f)
array([ 0.06799017,  0.86401967,  0.06799017]) 

使用特定的窗口函数:

>>> signal.firwin(numtaps, f, window='nuttall')
array([  3.56607041e-04,   9.99286786e-01,   3.56607041e-04]) 

高通(从 0 到 f):

>>> signal.firwin(numtaps, f, pass_zero=False)
array([-0.00859313,  0.98281375, -0.00859313]) 

带通:

>>> f1, f2 = 0.1, 0.2
>>> signal.firwin(numtaps, [f1, f2], pass_zero=False)
array([ 0.06301614,  0.88770441,  0.06301614]) 

带阻:

>>> signal.firwin(numtaps, [f1, f2])
array([-0.00801395,  1.0160279 , -0.00801395]) 

多带通(通带为 [0, f1],[f2, f3] 和 [f4, 1]):

>>> f3, f4 = 0.3, 0.4
>>> signal.firwin(numtaps, [f1, f2, f3, f4])
array([-0.01376344,  1.02752689, -0.01376344]) 

多带通(通带为 [f1, f2] 和 [f3,f4]):

>>> signal.firwin(numtaps, [f1, f2, f3, f4], pass_zero=False)
array([ 0.04890915,  0.91284326,  0.04890915]) 

scipy.signal.firwin2

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.firwin2.html#scipy.signal.firwin2

scipy.signal.firwin2(numtaps, freq, gain, *, nfreqs=None, window='hamming', nyq=<object object>, antisymmetric=False, fs=None)

使用窗口方法设计 FIR 滤波器。

根据给定的频率 freq 和相应的增益 gain,此函数构造具有线性相位和(近似)给定频率响应的 FIR 滤波器。

参数:

numtapsint

FIR 滤波器中的 taps 数。numtaps 必须小于 nfreqs

freqarray_like, 1-D

频率采样点。通常为 0.0 到 1.0,其中 1.0 为奈奎斯特。奈奎斯特频率是 fs 的一半。 freq 中的值必须是非递减的。一个值可以重复一次以实现不连续性。 freq 中的第一个值必须为 0,最后一个值必须为 fs/2。值 0 和 fs/2 不得重复。

gainarray_like

频率采样点处的滤波器增益。根据滤波器类型应用某些增益值的约束条件,请参阅备注以获取详细信息。

nfreqsint, optional

用于构建滤波器的插值网格的大小。为了实现最有效的行为,这应该是一个 2 的幂加 1(例如 129, 257 等)。默认值为大于等于 numtaps 的最小 2 的幂加 1。nfreqs 必须大于 numtaps

windowstring or (string, float) or float, or None, optional

要使用的窗口函数。默认值为“hamming”。参见 scipy.signal.get_window 获取可能值的完整列表。如果为 None,则不应用窗口函数。

nyqfloat, optional, deprecated

这是奈奎斯特频率。 freq 中的每个频率必须在 0 和 nyq 之间。默认值为 1。

自 1.0.0 版本起已弃用:firwin2 关键字参数 nyq 已弃用,改为使用 fs,将在 SciPy 1.14.0 中删除。

antisymmetricbool, optional

结果脉冲响应是否对称/反对称。更多细节请参见备注。

fsfloat, optional

信号的采样频率。 cutoff 中的每个频率必须在 0 和 fs/2 之间。默认值为 2。

返回:

tapsndarray

FIR 滤波器的滤波器系数,作为长度为 numtaps 的 1-D 数组。

参见

firls

firwin

minimum_phase

remez

备注

从给定的频率和增益集合中,在频率域中构造所需的响应。将逆 FFT 应用于所需的响应以创建相关的卷积核,并返回此卷积核的前 numtaps 系数,按 window 缩放。

FIR 滤波器将具有线性相位。滤波器的类型由 numtaps 的值和 antisymmetric 标志确定。有四种可能的组合:

  • 即使 numtaps 为奇数,antisymmetric 为 False,生成类型 I 滤波器。
  • 即使 numtaps 为偶数,antisymmetric 为 False,生成类型 II 滤波器。
  • 即使 numtaps 为奇数,antisymmetric 为 True,生成类型 III 滤波器。
  • 即使 numtaps 为偶数,antisymmetric 为 True,生成类型 IV 滤波器。

除了类型 I 滤波器外,所有滤波器的幅度响应都受以下约束的影响:

  • 类型 II – 零频率处为 Nyquist 频率。
  • 类型 III – 零频率和 Nyquist 频率处为零。
  • 类型 IV – 零频率处为零。

新版本为 0.9.0。

参考文献

[1]

Oppenheim, A. V. 和 Schafer, R. W.,“Discrete-Time Signal Processing”,Prentice-Hall,Englewood Cliffs,New Jersey(1989)。(例如,参见第 7.4 节。)

[2]

Smith, Steven W.,“The Scientist and Engineer’s Guide to Digital Signal Processing”,第十七章。www.dspguide.com/ch17/1.htm

示例

一个低通 FIR 滤波器,其响应在 [0.0, 0.5] 上为 1,并且在 [0.5, 1.0] 上从 1 线性减少到 0:

>>> from scipy import signal
>>> taps = signal.firwin2(150, [0.0, 0.5, 1.0], [1.0, 1.0, 0.0])
>>> print(taps[72:78])
[-0.02286961 -0.06362756  0.57310236  0.57310236 -0.06362756 -0.02286961] 

scipy.signal.freqs

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.freqs.html#scipy.signal.freqs

scipy.signal.freqs(b, a, worN=200, plot=None)

计算模拟滤波器的频率响应。

给定模拟滤波器的 M 阶分子b和 N 阶分母a,计算其频率响应:

 b[0]*(jw)**M + b[1]*(jw)**(M-1) + ... + b[M]
H(w) = ----------------------------------------------
        a[0]*(jw)**N + a[1]*(jw)**(N-1) + ... + a[N] 

参数:

数组b

线性滤波器的分子。

数组a

线性滤波器的分母。

worN{None, int, array_like},可选

如果为 None,则在响应曲线的有趣部分周围的 200 个频率上计算(由极点零点位置决定)。如果是单个整数,则在那么多频率上计算。否则,计算在给定的角频率(例如,rad/s)处给出的响应worn

plot可调用函数,可选

接受两个参数的可调用函数。如果给定,则将返回参数wh传递给 plot。用于在freqs内部绘制频率响应。

返回:

数组w

计算h的角频率。

数组h

频率响应。

另请参见

freqz

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

注意事项

使用 Matplotlib 的“plot”函数作为plot的可调用函数会产生意外的结果,这会绘制复数传递函数的实部,而不是幅度。尝试lambda w, h: plot(w, abs(h))

示例

>>> from scipy.signal import freqs, iirfilter
>>> import numpy as np 
>>> b, a = iirfilter(4, [1, 10], 1, 60, analog=True, ftype='cheby1') 
>>> w, h = freqs(b, a, worN=np.logspace(-1, 2, 1000)) 
>>> import matplotlib.pyplot as plt
>>> plt.semilogx(w, 20 * np.log10(abs(h)))
>>> plt.xlabel('Frequency')
>>> plt.ylabel('Amplitude response [dB]')
>>> plt.grid(True)
>>> plt.show() 

../../_images/scipy-signal-freqs-1.png

scipy.signal.freqs_zpk

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.freqs_zpk.html#scipy.signal.freqs_zpk

scipy.signal.freqs_zpk(z, p, k, worN=200)

计算模拟滤波器的频率响应。

给定滤波器的零点z,极点p和增益k,计算其频率响应:

 (jw-z[0]) * (jw-z[1]) * ... * (jw-z[-1])
H(w) = k * ----------------------------------------
           (jw-p[0]) * (jw-p[1]) * ... * (jw-p[-1]) 

参数:

zarray_like

线性滤波器的零点

parray_like

线性滤波器的极点

kscalar

线性滤波器的增益

worN{None, int, array_like}, 可选

如果为 None,则在响应曲线的有趣部分周围的 200 个频率处计算(由极点位置确定)。 如果为单个整数,则计算该数量的频率。 否则,计算在worN给定的角频率(例如,rad/s)处的响应。

返回:

wndarray

在计算h时使用的角频率。

hndarray

频率响应。

另请参见

freqs

计算 TF 形式模拟滤波器的频率响应

freqz

计算 TF 形式数字滤波器的频率响应

freqz_zpk

计算 ZPK 形式数字滤波器的频率响应

注意

新版 0.19.0 中新增。

示例

>>> import numpy as np
>>> from scipy.signal import freqs_zpk, iirfilter 
>>> z, p, k = iirfilter(4, [1, 10], 1, 60, analog=True, ftype='cheby1',
...                     output='zpk') 
>>> w, h = freqs_zpk(z, p, k, worN=np.logspace(-1, 2, 1000)) 
>>> import matplotlib.pyplot as plt
>>> plt.semilogx(w, 20 * np.log10(abs(h)))
>>> plt.xlabel('Frequency')
>>> plt.ylabel('Amplitude response [dB]')
>>> plt.grid(True)
>>> plt.show() 

../../_images/scipy-signal-freqs_zpk-1.png

scipy.signal.freqz

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.freqz.html#scipy.signal.freqz

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

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

给定数字滤波器的 M 阶分子b和 N 阶分母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:{None, int, array_like},可选

如果是单个整数,则在那么多的频率上进行计算(默认值为 N=512)。这是以下便利的替代方法:

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

使用快速 FFT 计算的数字可以导致更快的计算(见备注)。

如果是 array_like,则在给定的频率上计算响应。这些频率与fs的单位相同。

whole:bool,可选

通常,频率从 0 到 Nyquist 频率 fs/2(单位圆的上半部分)计算。如果whole为 True,则从 0 到 fs 计算频率。如果 worN 是 array_like,则忽略。

plot:callable

一个接受两个参数的可调用函数。如果提供了,返回参数wh将传递给绘图函数。用于在freqz内绘制频率响应。

fs:float,可选

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

1.2.0 版的新功能。

include_nyquist:bool,可选

如果whole为 False 且worN为整数,则将include_nyquist设置为 True 将包括最后一个频率(Nyquist 频率),否则将被忽略。

1.5.0 版的新功能。

返回:

w:ndarray

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

h:ndarray

频率响应,作为复数。

另请参阅

[freqz_zpk

sosfreqz

备注

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

在满足以下条件时使用直接计算(R)FFT 来计算频率响应:

  1. worN 是一个整数值。

  2. worN 通过 FFT 计算快速(即,next_fast_len(worN) 等于 worN)。

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

  4. worN 至少与分子系数的长度相同(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() 

../../_images/scipy-signal-freqz-1_00_00.png

广播示例

假设我们有两个 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]。这两个分母的系数存储在 2-D 数组 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 是多于 1 维的。为了使其与频率广播兼容,在调用freqz时,我们通过在虚拟维度上扩展它来实现:

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

scipy.signal.freqz_zpk

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.freqz_zpk.html#scipy.signal.freqz_zpk

scipy.signal.freqz_zpk(z, p, k, worN=512, whole=False, fs=6.283185307179586)

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

给定数字滤波器的零点、极点和增益,计算其频率响应:

(H(z)=k \prod_i (z - Z[i]) / \prod_j (z - P[j]))

其中(k)为增益,(Z)为零点,(P)为极点

参数:

zarray_like

线性滤波器的零点

parray_like

线性滤波器的极点

k标量

线性滤波器的增益

worN{None, int, array_like},可选

如果是单个整数,则在该数量的频率上计算(默认值为 N=512)。

如果是 array_like,则计算给定频率处的响应。这些频率与fs具有相同的单位。

whole布尔值,可选

通常,频率从 0 到 Nyquist 频率 fs/2(单位圆的上半部分)计算。如果whole为 True,则从 0 到 fs 计算频率。如果w为 array_like,则忽略。

fs浮点数,可选

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

新版本为 1.2.0。

返回:

wndarray

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

hndarray

作为复数的频率响应。

另请参见

freqs

计算 TF 形式模拟滤波器的频率响应

freqs_zpk

计算 ZPK 形式模拟滤波器的频率响应

freqz

计算 TF 形式数字滤波器的频率响应

笔记

新版本为 0.19.0。

示例

在采样率为 1000 Hz 的系统中,设计一个 4 阶数字 Butterworth 滤波器,截止频率为 100 Hz,并绘制其频率响应:

>>> import numpy as np
>>> from scipy import signal
>>> z, p, k = signal.butter(4, 100, output='zpk', fs=1000)
>>> w, h = signal.freqz_zpk(z, p, k, fs=1000) 
>>> import matplotlib.pyplot as plt
>>> fig = plt.figure()
>>> ax1 = fig.add_subplot(1, 1, 1)
>>> 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 [Hz]')
>>> ax1.grid(True) 
>>> ax2 = ax1.twinx()
>>> angles = np.unwrap(np.angle(h))
>>> ax2.plot(w, angles, 'g')
>>> ax2.set_ylabel('Angle [radians]', color='g') 
>>> plt.axis('tight')
>>> plt.show() 

../../_images/scipy-signal-freqz_zpk-1.png

scipy.signal.sosfreqz

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.sosfreqz.html#scipy.signal.sosfreqz

scipy.signal.sosfreqz(sos, worN=512, whole=False, fs=6.283185307179586)

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

给定sos,一个形状为(n, 6)的数组,其中包含数字滤波器的二阶段。计算系统函数的频率响应:

 B0(z)   B1(z)         B{n-1}(z)
H(z) = ----- * ----- * ... * ---------
       A0(z)   A1(z)         A{n-1}(z) 

对于 z = exp(omega*1j),其中 B{k}(z)和 A{k}(z)分别是第 k 个二阶段传递函数的分子和分母。

参数:

sos 类似数组

数组的二阶滤波器系数,必须具有形状(n_sections, 6)。每行对应一个二阶段,前三列提供分子系数,后三列提供分母系数。

worN{无,整数,类似数组},可选

如果是单个整数,则计算那么多频率(默认为 N=512)。使用 FFT 计算快速的数字可以导致更快的计算(见freqz的注意事项)。

如果是类似数组,则在给定频率处计算响应(必须是 1-D)。这些单位与fs相同。

whole 布尔值,可选

通常,频率从 0 到 Nyquist 频率 fs/2 计算(单位圆的上半部分)。如果whole为 True,则从 0 到 fs 计算频率。

fs 浮点数,可选

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

新版本 1.2.0 中。

返回:

w ndarray

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

h ndarray

频率响应,作为复数。

另请参见

freqzsosfilt

注意

新版本 0.19.0 中。

示例

设计一个 15 阶带通滤波器的 SOS 格式。

>>> from scipy import signal
>>> import numpy as np
>>> sos = signal.ellip(15, 0.5, 60, (0.2, 0.4), btype='bandpass',
...                    output='sos') 

在 1500 点处从 DC 到 Nyquist 计算频率响应。

>>> w, h = signal.sosfreqz(sos, worN=1500) 

绘制响应。

>>> import matplotlib.pyplot as plt
>>> plt.subplot(2, 1, 1)
>>> db = 20*np.log10(np.maximum(np.abs(h), 1e-5))
>>> plt.plot(w/np.pi, db)
>>> plt.ylim(-75, 5)
>>> plt.grid(True)
>>> plt.yticks([0, -20, -40, -60])
>>> plt.ylabel('Gain [dB]')
>>> plt.title('Frequency Response')
>>> plt.subplot(2, 1, 2)
>>> plt.plot(w/np.pi, np.angle(h))
>>> plt.grid(True)
>>> plt.yticks([-np.pi, -0.5*np.pi, 0, 0.5*np.pi, np.pi],
...            [r'$-\pi$', r'$-\pi/2$', '0', r'$\pi/2$', r'$\pi$'])
>>> plt.ylabel('Phase [rad]')
>>> plt.xlabel('Normalized frequency (1.0 = Nyquist)')
>>> plt.show() 

../../_images/scipy-signal-sosfreqz-1_00_00.png

如果将相同的滤波器实现为单个传递函数,数值误差会损坏频率响应:

>>> b, a = signal.ellip(15, 0.5, 60, (0.2, 0.4), btype='bandpass',
...                    output='ba')
>>> w, h = signal.freqz(b, a, worN=1500)
>>> plt.subplot(2, 1, 1)
>>> db = 20*np.log10(np.maximum(np.abs(h), 1e-5))
>>> plt.plot(w/np.pi, db)
>>> plt.ylim(-75, 5)
>>> plt.grid(True)
>>> plt.yticks([0, -20, -40, -60])
>>> plt.ylabel('Gain [dB]')
>>> plt.title('Frequency Response')
>>> plt.subplot(2, 1, 2)
>>> plt.plot(w/np.pi, np.angle(h))
>>> plt.grid(True)
>>> plt.yticks([-np.pi, -0.5*np.pi, 0, 0.5*np.pi, np.pi],
...            [r'$-\pi$', r'$-\pi/2$', '0', r'$\pi/2$', r'$\pi$'])
>>> plt.ylabel('Phase [rad]')
>>> plt.xlabel('Normalized frequency (1.0 = Nyquist)')
>>> plt.show() 

../../_images/scipy-signal-sosfreqz-1_01_00.png

scipy.signal.gammatone

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.gammatone.html#scipy.signal.gammatone

scipy.signal.gammatone(freq, ftype, order=None, numtaps=None, fs=None)

Gammatone 滤波器设计。

此函数计算 FIR 或 IIR Gammatone 数字滤波器的系数 [1]

参数:

freqfloat

滤波器的中心频率(与 fs 相同的单位表示)。

ftype{‘fir’, ‘iir’}

函数生成的滤波器类型。如果是 ‘fir’,函数将生成一个 N 阶 FIR Gammatone 滤波器。如果是 ‘iir’,函数将生成一个 8 阶数字 IIR 滤波器,建模为 4 阶 Gammatone 滤波器。

orderint, optional

滤波器的阶数。仅在 ftype='fir' 时使用。默认为 4,用于模拟人类听觉系统。必须介于 0 和 24 之间。

numtapsint, optional

滤波器的长度。仅在 ftype='fir' 时使用。默认为 fs*0.015(如果 fs 大于 1000),15(如果 fs 小于或等于 1000)。

fsfloat, optional

信号的采样频率。freq 必须介于 0 和 fs/2 之间。默认为 2。

返回:

b, andarray, ndarray

滤波器的分子 (b) 和分母 (a) 多项式。

Raises:

ValueError

如果 freq 小于或等于 0 或大于或等于 fs/2,如果 ftype 不是 ‘fir’ 或 ‘iir’,如果 orderftype='fir' 时小于或等于 0 或大于 24

参见

firwin

iirfilter

参考文献

[1]

Slaney, Malcolm, “An Efficient Implementation of the Patterson-Holdsworth Auditory Filter Bank”, Apple Computer Technical Report 35, 1993, pp.3-8, 34-39.

示例

以 440 Hz 为中心的 16 采样 4 阶 FIR Gammatone 滤波器

>>> from scipy import signal
>>> signal.gammatone(440, 'fir', numtaps=16, fs=16000)
(array([ 0.00000000e+00,  2.22196719e-07,  1.64942101e-06,  4.99298227e-06,
 1.01993969e-05,  1.63125770e-05,  2.14648940e-05,  2.29947263e-05,
 1.76776931e-05,  2.04980537e-06, -2.72062858e-05, -7.28455299e-05,
 -1.36651076e-04, -2.19066855e-04, -3.18905076e-04, -4.33156712e-04]),
 [1.0]) 

以 440 Hz 为中心的 IIR Gammatone 滤波器

>>> import matplotlib.pyplot as plt
>>> import numpy as np 
>>> b, a = signal.gammatone(440, 'iir', fs=16000)
>>> w, h = signal.freqz(b, a)
>>> plt.plot(w / ((2 * np.pi) / 16000), 20 * np.log10(abs(h)))
>>> plt.xscale('log')
>>> plt.title('Gammatone filter frequency response')
>>> plt.xlabel('Frequency')
>>> plt.ylabel('Amplitude [dB]')
>>> plt.margins(0, 0.1)
>>> plt.grid(which='both', axis='both')
>>> plt.axvline(440, color='green') # cutoff frequency
>>> plt.show() 

../../_images/scipy-signal-gammatone-1.png

scipy.signal.group_delay

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.group_delay.html#scipy.signal.group_delay

scipy.signal.group_delay(system, w=512, whole=False, fs=6.283185307179586)

计算数字滤波器的组延迟。

组延迟测量信号各谱成分幅度包络被滤波器延迟多少个样本。形式上定义为连续(展开)相位的导数:

 d        jw
D(w) = - -- arg H(e)
         dw 

参数:

system数组对(b, a)

滤波器传输函数的分子和分母系数。

w{无,整数,数组形式},可选

如果是单个整数,则在那么多的频率上进行计算(默认为 N=512)。

如果是数组形式,则计算给定频率下的延迟。这些频率与fs单位相同。

whole布尔值,可选

通常,频率从 0 到奈奎斯特频率 fs/2(单位圆的上半部分)计算。如果whole为 True,则从 0 到 fs 计算频率。如果 w 是数组形式,则忽略。

fs浮点数,可选

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

新版功能于版本 1.2.0 中添加。

返回:

wndarray

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

gdndarray

组延迟。

另请参阅

[freqz

数字滤波器的频率响应

注释

MATLAB 中的类似函数称为grpdelay

如果数字系统的传输函数(H(z))在单位圆上有零点或极点,则在相应频率下的组延迟是未定义的。当出现这种情况时,会发出警告,并将组延迟设置为这些频率上的 0。

关于组延迟的数值计算的详细信息,请参考[1]

新版功能于版本 0.16.0 中添加。

参考文献

[1]

Richard G. Lyons,《理解数字信号处理,第 3 版》,第 830 页。

示例

>>> from scipy import signal
>>> b, a = signal.iirdesign(0.1, 0.3, 5, 50, ftype='cheby1')
>>> w, gd = signal.group_delay((b, a)) 
>>> import matplotlib.pyplot as plt
>>> plt.title('Digital filter group delay')
>>> plt.plot(w, gd)
>>> plt.ylabel('Group delay [samples]')
>>> plt.xlabel('Frequency [rad/sample]')
>>> plt.show() 

../../_images/scipy-signal-group_delay-1.png

scipy.signal.iirdesign

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.iirdesign.html#scipy.signal.iirdesign

scipy.signal.iirdesign(wp, ws, gpass, gstop, analog=False, ftype='ellip', output='ba', fs=None)

完整的 IIR 数字和模拟滤波器设计。

根据给定的基本类型的通带和阻带频率及增益构造模拟或数字 IIR 滤波器的最小阶数。以分子、分母(‘ba’)、极点-零点(‘zpk’)或二阶段(‘sos’)形式返回输出。

参数:

wp, wsfloat 或 array like, 形状 (2,)

通带和阻带边缘频率。可能的取值为标量(适用于低通和高通滤波器)或范围(适用于带通和带阻滤波器)。对于数字滤波器,这些频率与fs(采样频率)的单位相同。默认情况下,fs是每个样本的 2 个半周期,因此这些频率被归一化为 0 到 1,其中 1 是奈奎斯特频率。例如:

  • 低通:wp = 0.2,ws = 0.3
  • 高通:wp = 0.3,ws = 0.2
  • 带通:wp = [0.2, 0.5],ws = [0.1, 0.6]
  • 带阻:wp = [0.1, 0.6],ws = [0.2, 0.5]

对于模拟滤波器,wpws 是角频率(例如,rad/s)。注意,对于带通和带阻滤波器,通带必须严格位于阻带内,反之亦然。

gpassfloat

通带中的最大损失(dB)。

gstopfloat

在阻带中的最小衰减(dB)。

analogbool, 可选

当为 True 时,返回模拟滤波器,否则返回数字滤波器。

ftypestr, 可选

要设计的 IIR 滤波器类型:

  • Butterworth:‘butter’
  • Chebyshev I:‘cheby1’
  • Chebyshev II:‘cheby2’
  • Cauer/elliptic:‘ellip’

output{‘ba’, ‘zpk’, ‘sos’}, 可选

输出的滤波器形式:

  • 推荐的二阶段形式:‘sos’
  • 分子/分母(默认):‘ba’
  • 极点-零点:‘zpk’

一般推荐使用二阶段形式(‘sos’),因为推断分子/分母形式(‘ba’)的系数会受到数值不稳定性的影响。出于向后兼容性的考虑,默认形式是分子/分母形式(‘ba’),其中‘b’ 和 ‘a’ 分别是系数的常用名称。

注意:有时使用二阶段形式(‘sos’)会伴随额外的计算成本:因此建议对数据密集型应用进行探索,也要考虑分子/分母形式(‘ba’)。

fsfloat, 可选

数字系统的采样频率。

新版本 1.2.0 中的新增功能。

返回:

b, andarray, ndarray

IIR 滤波器的分子(b)和分母(a)多项式。仅在output='ba'时返回。

z, p, kndarray, ndarray, float

IIR 滤波器传递函数的零点、极点和系统增益。仅在output='zpk'时返回。

sosndarray

IIR 滤波器的二阶段表示。仅在output='sos'时返回。

另请参见

butter

使用阶数和临界点设计滤波器

cheby1, cheby2, ellip, bessel

buttord

根据通带和阻带规格找到阶数和关键点

cheb1ord, cheb2ord, ellipord

iirfilter

使用阶数和关键频率进行一般滤波器设计

笔记

'sos'输出参数是在 0.16.0 版本中添加的。

示例

>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> import matplotlib.ticker 
>>> wp = 0.2
>>> ws = 0.3
>>> gpass = 1
>>> gstop = 40 
>>> system = signal.iirdesign(wp, ws, gpass, gstop)
>>> w, h = signal.freqz(*system) 
>>> 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]')
>>> ax1.grid(True)
>>> ax1.set_ylim([-120, 20])
>>> 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')
>>> ax2.set_ylim([-6, 1])
>>> nticks = 8
>>> ax1.yaxis.set_major_locator(matplotlib.ticker.LinearLocator(nticks))
>>> ax2.yaxis.set_major_locator(matplotlib.ticker.LinearLocator(nticks)) 

../../_images/scipy-signal-iirdesign-1.png

scipy.signal.iirfilter

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.iirfilter.html#scipy.signal.iirfilter

scipy.signal.iirfilter(N, Wn, rp=None, rs=None, btype='band', analog=False, ftype='butter', output='ba', fs=None)

给定阶数和临界点,设计 IIR 数字和模拟滤波器。

设计一个 N 阶数字或模拟滤波器,并返回滤波器系数。

参数:

N 整数

滤波器的阶数。

Wn 类似数组

给出临界频率的标量或长度为 2 的序列。

对于数字滤波器,Wn的单位与fs相同。默认情况下,fs为每个采样 2 个半周期,因此这些值从 0 到 1 进行归一化,其中 1 为奈奎斯特频率。(Wn因此是半周期/每个样本。)

对于模拟滤波器,Wn 是一个角频率(例如,rad/s)。

当 Wn 是长度为 2 的序列时,Wn[0]必须小于Wn[1]

rp 浮点数, 可选

对于 Chebyshev 和 elliptic 滤波器,提供通带中的最大波纹。(dB)

rs 浮点数, 可选

对于 Chebyshev 和 elliptic 滤波器,提供阻带中的最小衰减。(dB)

btype {'bandpass', 'lowpass', 'highpass', 'bandstop'}, 可选

滤波器类型。默认为‘bandpass’。

analog 布尔值, 可选

当为 True 时,返回模拟滤波器,否则返回数字滤波器。

ftype 字符串, 可选

要设计的 IIR 滤波器类型:

  • Butterworth:‘butter’
  • Chebyshev I:‘cheby1’
  • Chebyshev II:‘cheby2’
  • Cauer/elliptic:‘ellip’
  • Bessel/Thomson:‘bessel’

output {'ba', 'zpk', 'sos'}, 可选

输出的过滤器形式:

  • 二阶段形式(推荐):‘sos’
  • 默认的分子/分母形式:‘ba’
  • 极点-零点形式:‘zpk’

一般推荐使用二阶段形式(‘sos’),因为推断出分子/分母形式(‘ba’)的系数会受到数值不稳定性的影响。出于向后兼容性的考虑,默认形式为分子/分母形式(‘ba’),其中‘b’和‘a’分别指代所用系数的常用名称。

注意:使用二阶段形式(‘sos’)有时会伴随额外的计算成本:因此,建议对于数据密集的用例也研究分子/分母形式(‘ba’)。

fs 浮点数, 可选

数字系统的采样频率。

1.2.0 版中新增。

返回:

b, a 数组, 数组

IIR 滤波器的分子(b)和分母(a)多项式。仅在output='ba'时返回。

z, p, k 数组, 数组, 浮点数

IIR 滤波器传递函数的零点、极点和系统增益。仅在output='zpk'时返回。

sos 数组

IIR 滤波器的二阶段形式表示。仅在output='sos'时返回。

另请参阅

butter

使用阶数和临界点进行滤波器设计。

cheby1, cheby2, ellip, bessel

buttord

从通带和阻带规范中找到阶数和临界点。

cheb1ord, cheb2ord, ellipord

iirdesign

使用通带和阻带规范进行通用滤波器设计。

注意

'sos'输出参数在 0.16.0 版本中被添加。

示例

生成一个从 50 Hz 到 200 Hz 的 17 阶 Chebyshev II 模拟带通滤波器,并绘制频率响应图:

>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt 
>>> b, a = signal.iirfilter(17, [2*np.pi*50, 2*np.pi*200], rs=60,
...                         btype='band', analog=True, ftype='cheby2')
>>> w, h = signal.freqs(b, a, 1000)
>>> fig = plt.figure()
>>> ax = fig.add_subplot(1, 1, 1)
>>> ax.semilogx(w / (2*np.pi), 20 * np.log10(np.maximum(abs(h), 1e-5)))
>>> ax.set_title('Chebyshev Type II bandpass frequency response')
>>> ax.set_xlabel('Frequency [Hz]')
>>> ax.set_ylabel('Amplitude [dB]')
>>> ax.axis((10, 1000, -100, 10))
>>> ax.grid(which='both', axis='both')
>>> plt.show() 

../../_images/scipy-signal-iirfilter-1_00_00.png

在采样率为 2000 Hz 的系统中创建具有相同特性的数字滤波器,并绘制频率响应图。(需要使用二阶段实现以确保这一阶数的滤波器稳定性):

>>> sos = signal.iirfilter(17, [50, 200], rs=60, btype='band',
...                        analog=False, ftype='cheby2', fs=2000,
...                        output='sos')
>>> w, h = signal.sosfreqz(sos, 2000, fs=2000)
>>> fig = plt.figure()
>>> ax = fig.add_subplot(1, 1, 1)
>>> ax.semilogx(w, 20 * np.log10(np.maximum(abs(h), 1e-5)))
>>> ax.set_title('Chebyshev Type II bandpass frequency response')
>>> ax.set_xlabel('Frequency [Hz]')
>>> ax.set_ylabel('Amplitude [dB]')
>>> ax.axis((10, 1000, -100, 10))
>>> ax.grid(which='both', axis='both')
>>> plt.show() 

../../_images/scipy-signal-iirfilter-1_01_00.png

scipy.signal.kaiser_atten

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.kaiser_atten.html#scipy.signal.kaiser_atten

scipy.signal.kaiser_atten(numtaps, width)

计算 Kaiser FIR 滤波器的衰减。

给定抽头数N和过渡宽度width,使用 Kaiser 公式计算衰减a,如下所示:

a = 2.285 * (N - 1) * pi * width + 7.95

Parameters:

numtapsint

FIR 滤波器的抽头数量。

widthfloat

滤波器在通带和阻带之间(或一般来说,在任何不连续处)的过渡区域的期望宽度,以奈奎斯特频率的分数形式表示。

Returns:

afloat

波纹的衰减,单位为 dB。

另请参见

kaiserord, kaiser_beta

Examples

假设我们希望使用 Kaiser 窗口方法设计一个 FIR 滤波器,该滤波器有 211 个抽头,并且在采样频率为 480 Hz 的信号中具有 9 Hz 的过渡宽度。以奈奎斯特频率的分数形式表示,宽度为 9/(0.5*480) = 0.0375. 按照以下公式计算近似衰减(以 dB 为单位):

>>> from scipy.signal import kaiser_atten
>>> kaiser_atten(211, 0.0375)
64.48099630593983 

scipy.signal.kaiser_beta

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.kaiser_beta.html#scipy.signal.kaiser_beta

scipy.signal.kaiser_beta(a)

计算 Kaiser 参数 beta,给定阻带衰减 a

参数:

afloat

在 dB 中所需的阻带衰减和通带中的最大波动。这应该是一个 正数

返回值:

betafloat

用于计算 Kaiser 窗口公式中的 beta 参数。

参考文献

Oppenheim, Schafer, “Discrete-Time Signal Processing”, p.475-476.

示例

假设我们想设计一个低通滤波器,在阻带中具有 65 dB 的衰减。通过 kaiser_beta(65) 计算的 Kaiser 窗口参数如下:

>>> from scipy.signal import kaiser_beta
>>> kaiser_beta(65)
6.20426 

scipy.signal.kaiserord

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.kaiserord.html#scipy.signal.kaiserord

scipy.signal.kaiserord(ripple, width)

确定 Kaiser 窗口方法的滤波器窗口参数。

此函数返回的参数通常用于使用窗口法创建有限冲激响应滤波器,可以使用firwinfirwin2

参数:

ripple浮点数

滤波器的频率响应的幅度与所需滤波器的频率响应的偏差的上限(不包括任何过渡区间中的频率)。也就是说,如果 w 是以 Nyquist 频率的分数表示的频率,则 A(w)是滤波器的实际频率响应,D(w)是期望的频率响应,则设计要求是:

abs(A(w) - D(w))) < 10**(-ripple/20) 

对于 0 <= w <= 1 且 w 不在过渡区间内。

width浮点数

过渡区域的宽度,标准化为对应于每个采样π弧度。也就是说,频率表示为 Nyquist 频率的分数。

返回:

numtaps整数

Kaiser 窗口的长度。

beta浮点数

Kaiser 窗口的 beta 参数。

另见

kaiser_beta, kaiser_atten

注意事项

有几种方法可以获取 Kaiser 窗口:

  • signal.windows.kaiser(numtaps, beta, sym=True)

  • signal.get_window(beta, numtaps)

  • signal.get_window(('kaiser', beta), numtaps)

Kaiser 发现的经验方程式被使用。

参考文献

Oppenheim, Schafer, “离散时间信号处理”, pp.475-476.

示例

我们将使用 Kaiser 窗口方法为 1000 Hz 采样的信号设计低通 FIR 滤波器。

我们希望在阻带中至少有 65 dB 的抑制,在通带中增益变化不超过 0.5%。

我们希望截止频率为 175 Hz,通带和阻带之间的过渡为 24 Hz。也就是说,在区间[0, 163]内,增益变化不超过 0.5%,在区间[187, 500]内,信号至少被 65 dB 衰减。

>>> import numpy as np
>>> from scipy.signal import kaiserord, firwin, freqz
>>> import matplotlib.pyplot as plt
>>> fs = 1000.0
>>> cutoff = 175
>>> width = 24 

Kaiser 方法接受一个参数来控制通带波动和阻带抑制,因此我们选择两者中较为严格的一个。在这种情况下,通带波动为 0.005,即 46.02 dB,因此我们将使用 65 dB 作为设计参数。

使用kaiserord确定滤波器的长度和 Kaiser 窗口的参数。

>>> numtaps, beta = kaiserord(65, width/(0.5*fs))
>>> numtaps
167
>>> beta
6.20426 

使用firwin创建 FIR 滤波器。

>>> taps = firwin(numtaps, cutoff, window=('kaiser', beta),
...               scale=False, fs=fs) 

计算滤波器的频率响应。w 是频率数组,h 是相应的复数频率响应数组。

>>> w, h = freqz(taps, worN=8000)
>>> w *= 0.5*fs/np.pi  # Convert w to Hz. 

计算滤波器响应的幅度与理想低通滤波器的偏差。过渡区域的数值设为nan,因此它们不会出现在绘图中。

>>> ideal = w < cutoff  # The "ideal" frequency response.
>>> deviation = np.abs(np.abs(h) - ideal)
>>> deviation[(w > cutoff - 0.5*width) & (w < cutoff + 0.5*width)] = np.nan 

绘制偏差图。仔细观察阻带左端,显示出第一个主瓣中 65 dB 的衰减要求被超过约 0.125 dB。这对于凯泽窗方法并不罕见。

>>> plt.plot(w, 20*np.log10(np.abs(deviation)))
>>> plt.xlim(0, 0.5*fs)
>>> plt.ylim(-90, -60)
>>> plt.grid(alpha=0.25)
>>> plt.axhline(-65, color='r', ls='--', alpha=0.3)
>>> plt.xlabel('Frequency (Hz)')
>>> plt.ylabel('Deviation from ideal (dB)')
>>> plt.title('Lowpass Filter Frequency Response')
>>> plt.show() 

../../_images/scipy-signal-kaiserord-1.png

scipy.signal.minimum_phase

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.minimum_phase.html#scipy.signal.minimum_phase

scipy.signal.minimum_phase(h, method='homomorphic', n_fft=None)

将线性相位 FIR 滤波器转换为最小相位

参数:

h数组

线性相位 FIR 滤波器系数。

method{‘hilbert’, ‘homomorphic’}

使用的方法:

‘homomorphic’(默认)

此方法[4] [5] 最适用于具有奇数抽头数的滤波器,并且生成的最小相位滤波器的幅度响应近似于原始滤波器幅度响应的平方根。

‘hilbert’

此方法[1]设计用于等波纹滤波器(例如来自remez的滤波器),具有单位或零增益区域。

n_fft整数

FFT 使用的点数。应至少比信号长度大几倍(见注释)。

返回:

h_minimum数组

滤波器的最小相位版本,长度为(length(h) + 1) // 2

另请参阅

firwin

firwin2

remez

注释

希尔伯特[1]或同态[4] [5]方法都需要选择 FFT 长度以估算滤波器的复合倒谱。

在希尔伯特方法中,偏离理想频谱的epsilon与阻带零点数n_stop和 FFT 长度n_fft有关:

epsilon = 2. * n_stop / n_fft 

例如,有 100 个阻带零点和 FFT 长度为 2048 时,epsilon = 0.0976。如果我们保守地假设阻带零点数比滤波器长度少一个,我们可以将 FFT 长度取为满足epsilon = 0.01的下一个 2 的幂:

n_fft = 2 ** int(np.ceil(np.log2(2 * (len(h) - 1) / 0.01))) 

n_fft=None时使用的值,此方法对希尔伯特和同态方法都给出了合理的结果。

还存在其他方法创建最小相位滤波器,包括零反转[2]和频谱因子分解[3] [4]。更多信息请参见:

dspguru.com/dsp/howtos/how-to-design-minimum-phase-fir-filters

参考文献

[1] (1,2)

N. Damera-Venkata 和 B. L. Evans,“实数和复数最小相位数字 FIR 滤波器的最优设计”,声学、语音和信号处理,1999 年国际会议记录,凤凰城,AZ,1999,第 1145-1148 卷 3。DOI:10.1109/ICASSP.1999.756179

[2]

X. Chen 和 T. W. Parks,《通过直接因式分解设计最优最小相位 FIR 滤波器》,《信号处理》,第 10 卷,第 4 期,pp. 369-383,1986 年 6 月。

[3]

T. Saramaki,《有限冲激响应滤波器设计》,《数字信号处理手册》,第四章,纽约:Wiley-Interscience,1993 年。

[4] (1,2,3)

J. S. Lim,《信号处理的高级主题》。新泽西州恩格尔伍德克利夫斯:普林斯顿大厅,1988 年。

[5] (1,2)

A. V. Oppenheim, R. W. Schafer, 和 J. R. Buck,《离散时间信号处理》,第二版。新泽西州,上班顶部:普林斯顿大厅,1999 年。

例子

创建一个最优的线性相位滤波器,然后将其转换为最小相位:

>>> import numpy as np
>>> from scipy.signal import remez, minimum_phase, freqz, group_delay
>>> import matplotlib.pyplot as plt
>>> freq = [0, 0.2, 0.3, 1.0]
>>> desired = [1, 0]
>>> h_linear = remez(151, freq, desired, fs=2.) 

将其转换为最小相位:

>>> h_min_hom = minimum_phase(h_linear, method='homomorphic')
>>> h_min_hil = minimum_phase(h_linear, method='hilbert') 

比较这三个滤波器:

>>> fig, axs = plt.subplots(4, figsize=(4, 8))
>>> for h, style, color in zip((h_linear, h_min_hom, h_min_hil),
...                            ('-', '-', '--'), ('k', 'r', 'c')):
...     w, H = freqz(h)
...     w, gd = group_delay((h, 1))
...     w /= np.pi
...     axs[0].plot(h, color=color, linestyle=style)
...     axs[1].plot(w, np.abs(H), color=color, linestyle=style)
...     axs[2].plot(w, 20 * np.log10(np.abs(H)), color=color, linestyle=style)
...     axs[3].plot(w, gd, color=color, linestyle=style)
>>> for ax in axs:
...     ax.grid(True, color='0.5')
...     ax.fill_between(freq[1:3], *ax.get_ylim(), color='#ffeeaa', zorder=1)
>>> axs[0].set(xlim=[0, len(h_linear) - 1], ylabel='Amplitude', xlabel='Samples')
>>> axs[1].legend(['Linear', 'Min-Hom', 'Min-Hil'], title='Phase')
>>> for ax, ylim in zip(axs[1:], ([0, 1.1], [-150, 10], [-60, 60])):
...     ax.set(xlim=[0, 1], ylim=ylim, xlabel='Frequency')
>>> axs[1].set(ylabel='Magnitude')
>>> axs[2].set(ylabel='Magnitude (dB)')
>>> axs[3].set(ylabel='Group delay')
>>> plt.tight_layout() 

../../_images/scipy-signal-minimum_phase-1.png

scipy.signal.savgol_coeffs

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.savgol_coeffs.html#scipy.signal.savgol_coeffs

scipy.signal.savgol_coeffs(window_length, polyorder, deriv=0, delta=1.0, pos=None, use='conv')

计算 1-D Savitzky-Golay FIR 滤波器的系数。

参数:

window_length:整数

滤波器窗口的长度(即系数的数量)。

polyorder:整数

用于拟合样本的多项式的顺序。polyorder必须小于window_length

deriv:整数,可选

要计算的导数阶数。这必须是非负整数。默认值为 0,表示在不进行微分的情况下过滤数据。

delta:浮点数,可选

要应用滤波器的样本的间距。仅当 deriv > 0 时使用。

pos:整数或 None,可选

如果 pos 不为 None,则指定窗口内的评估位置。默认值为窗口的中间。

use:字符串,可选

‘conv’或‘dot’。此参数选择系数的顺序。默认值为‘conv’,表示系数按卷积使用的顺序排列。使用‘dot’时,顺序反转,因此通过将系数与数据集点乘来应用滤波器。

返回:

coeffs:1-D ndarray

滤波器系数。

参见

savgol_filter

注意事项

新版本 0.14.0 中引入。

参考文献

A. Savitzky, M. J. E. Golay, 简化最小二乘法的平滑和微分数据处理。分析化学,1964 年,36(8),第 1627-1639 页。罗建文,应奎,白静。2005 年。用于偶数数据的 Savitzky-Golay 平滑和微分滤波器。信号处理。85,7(2005 年 7 月),第 1429-1434 页。

示例

>>> import numpy as np
>>> from scipy.signal import savgol_coeffs
>>> savgol_coeffs(5, 2)
array([-0.08571429,  0.34285714,  0.48571429,  0.34285714, -0.08571429])
>>> savgol_coeffs(5, 2, deriv=1)
array([ 2.00000000e-01,  1.00000000e-01,  2.07548111e-16, -1.00000000e-01,
 -2.00000000e-01]) 

注意,use='dot'仅简单地反转系数。

>>> savgol_coeffs(5, 2, pos=3)
array([ 0.25714286,  0.37142857,  0.34285714,  0.17142857, -0.14285714])
>>> savgol_coeffs(5, 2, pos=3, use='dot')
array([-0.14285714,  0.17142857,  0.34285714,  0.37142857,  0.25714286])
>>> savgol_coeffs(4, 2, pos=3, deriv=1, use='dot')
array([0.45,  -0.85,  -0.65,  1.05]) 

x包含从抛物线 x = t**2 采样的数据,采样点为 t = -1, 0, 1, 2, 3。c保存了在最后一个位置计算导数的系数。当与x点乘时,结果应为 6。

>>> x = np.array([1, 0, 1, 4, 9])
>>> c = savgol_coeffs(5, 2, pos=4, deriv=1, use='dot')
>>> c.dot(x)
6.0 

scipy.signal.remez

原文:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.remez.html#scipy.signal.remez

scipy.signal.remez(numtaps, bands, desired, *, weight=None, Hz=<object object>, type='bandpass', maxiter=25, grid_density=16, fs=None)

使用 Remez 交换算法计算极小值最优滤波器。

使用 Remez 交换算法在指定频段中最小化所需增益与实际增益之间的最大误差,计算有限脉冲响应(FIR)滤波器的滤波器系数。

参数:

numtaps整数

滤波器中所需的阶数。阶数是滤波器中的项数,或者是滤波器阶数加一。

bands数组型

包含带边缘的单调序列。所有元素必须为非负且小于由 fs 给出的采样频率的一半。

desired数组型

包含每个指定频段中所需增益的带的一半大小的序列。

weight数组型,可选

给每个带区域赋予的相对权重。weight 的长度必须是 bands 长度的一半。

Hz标量,可选,已弃用

信号的采样频率(单位:赫兹)。默认为 1。

自 1.0.0 版起弃用:remez关键字参数Hz,将被fs取代,并将在 SciPy 1.14.0 中删除。

type{‘bandpass’, ‘differentiator’, ‘hilbert’},可选

滤波器类型:

  • ‘bandpass’:带通响应。这是默认设置。
  • ‘differentiator’:频率比例响应带。
  • ‘hilbert’滤波器,具有奇对称性,即类型 III
  • (对于偶数阶)或类型 IV(对于奇数阶)线性相位滤波器。

maxiter整数,可选

算法的最大迭代次数。默认为 25。

grid_density整数,可选

网格密度。在remez中使用的密集网格大小为(numtaps + 1) * grid_density。默认为 16。

fs浮点数,可选

信号的采样频率。默认为 1。

返回:

outndarray

一个秩为 1 的数组,包含最优(在最小最大意义上)滤波器的系数。

另见

firls

firwin

firwin2

minimum_phase

参考文献

[1]

J. H. McClellan 和 T. W. Parks,“一种统一的最优 FIR 线性相位数字滤波器设计方法”,IEEE Trans. Circuit Theory, vol. CT-20, pp. 697-701, 1973 年。

[2]

J. H. McClellan, T. W. Parks 和 L. R. Rabiner,“用于设计最优 FIR 线性相位数字滤波器的计算机程序”,IEEE Trans. Audio Electroacoust., vol. AU-21, pp. 506-525, 1973 年。

示例

在这些示例中,remez 用于设计低通、高通、带通和带阻滤波器。定义每个滤波器的参数包括滤波器阶数、频带边界、边界过渡宽度、每个频带中的期望增益以及采样频率。

在所有示例中,我们将使用 22050 Hz 的采样频率。在每个示例中,每个频带中的期望增益为 0(阻带)或 1(通带)。

freqz 用于计算每个滤波器的频率响应,下面定义的实用函数 plot_response 用于绘制响应。

>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt 
>>> fs = 22050   # Sample rate, Hz 
>>> def plot_response(w, h, title):
...     "Utility function to plot response functions"
...     fig = plt.figure()
...     ax = fig.add_subplot(111)
...     ax.plot(w, 20*np.log10(np.abs(h)))
...     ax.set_ylim(-40, 5)
...     ax.grid(True)
...     ax.set_xlabel('Frequency (Hz)')
...     ax.set_ylabel('Gain (dB)')
...     ax.set_title(title) 

第一个示例是一个低通滤波器,截止频率为 8 kHz。滤波器长度为 325,从通带到阻带的过渡宽度为 100 Hz。

>>> cutoff = 8000.0    # Desired cutoff frequency, Hz
>>> trans_width = 100  # Width of transition from pass to stop, Hz
>>> numtaps = 325      # Size of the FIR filter.
>>> taps = signal.remez(numtaps, [0, cutoff, cutoff + trans_width, 0.5*fs],
...                     [1, 0], fs=fs)
>>> w, h = signal.freqz(taps, [1], worN=2000, fs=fs)
>>> plot_response(w, h, "Low-pass Filter")
>>> plt.show() 

../../_images/scipy-signal-remez-1_00_00.png

此示例展示了一个高通滤波器:

>>> cutoff = 2000.0    # Desired cutoff frequency, Hz
>>> trans_width = 250  # Width of transition from pass to stop, Hz
>>> numtaps = 125      # Size of the FIR filter.
>>> taps = signal.remez(numtaps, [0, cutoff - trans_width, cutoff, 0.5*fs],
...                     [0, 1], fs=fs)
>>> w, h = signal.freqz(taps, [1], worN=2000, fs=fs)
>>> plot_response(w, h, "High-pass Filter")
>>> plt.show() 

../../_images/scipy-signal-remez-1_01_00.png

此示例展示了一个带通滤波器,通带从 2 kHz 到 5 kHz。过渡宽度为 260 Hz,滤波器长度为 63,比其他示例中的要小:

>>> band = [2000, 5000]  # Desired pass band, Hz
>>> trans_width = 260    # Width of transition from pass to stop, Hz
>>> numtaps = 63         # Size of the FIR filter.
>>> edges = [0, band[0] - trans_width, band[0], band[1],
...          band[1] + trans_width, 0.5*fs]
>>> taps = signal.remez(numtaps, edges, [0, 1, 0], fs=fs)
>>> w, h = signal.freqz(taps, [1], worN=2000, fs=fs)
>>> plot_response(w, h, "Band-pass Filter")
>>> plt.show() 

../../_images/scipy-signal-remez-1_02_00.png

低阶数导致更高的波动和更缓的过渡。

接下来的示例展示了一个带阻滤波器。

>>> band = [6000, 8000]  # Desired stop band, Hz
>>> trans_width = 200    # Width of transition from pass to stop, Hz
>>> numtaps = 175        # Size of the FIR filter.
>>> edges = [0, band[0] - trans_width, band[0], band[1],
...          band[1] + trans_width, 0.5*fs]
>>> taps = signal.remez(numtaps, edges, [1, 0, 1], fs=fs)
>>> w, h = signal.freqz(taps, [1], worN=2000, fs=fs)
>>> plot_response(w, h, "Band-stop Filter")
>>> plt.show() 

../../_images/scipy-signal-remez-1_03_00.png

scipy.signal.unique_roots

原文链接:docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.signal.unique_roots.html#scipy.signal.unique_roots

scipy.signal.unique_roots(p, tol=0.001, rtype='min')

从根列表中确定唯一根及其重数。

参数:

parray_like

根的列表。

tolfloat,可选

两个根被认为相等的公差。默认值为 1e-3。有关根分组细节,请参阅备注。

rtype{‘max’, ‘maximum’, ‘min’, ‘minimum’, ‘avg’, ‘mean’},可选

如果多个根在tol范围内,则如何确定返回的根。

  • ‘max’、‘maximum’:选择这些根中的最大值
  • ‘min’、‘minimum’:选择这些根中的最小值
  • ‘avg’、‘mean’:取这些根的平均值

在找到复根的最小或最大值时,首先比较实部,然后再比较虚部。

返回:

uniquendarray

唯一根的列表。

multiplicityndarray

每个根的重数。

备注

如果我们有根abc,使得a接近b,而b接近c(距离小于tol),则并不一定意味着a接近c。这意味着根分组不是唯一的。在此函数中,我们使用“贪婪”分组,按照输入p中给定的顺序遍历根。

此实用函数不专门用于根,而是可用于需要确定唯一性和重数的任何值序列。有关更通用的程序,请参阅numpy.unique

示例

>>> from scipy import signal
>>> vals = [0, 1.3, 1.31, 2.8, 1.25, 2.2, 10.3]
>>> uniq, mult = signal.unique_roots(vals, tol=2e-2, rtype='avg') 

检查具有大于 1 的重数的根:

>>> uniq[mult > 1]
array([ 1.305])