Dry/wet（即包括部分原始信号）IIR 陷波滤波器在与其他滤波器卷积之前的参数

Question

我有一组 scipy.signal.iirnotch 过滤器相互卷积，我希望每个过滤器都为每个过滤器实现独立的 dry/wet 值。我知道我可以通过将它们彼此相加而不是对它们执行卷积来实现一个伪 dry/wet 参数，但这是简陋且难以控制的。

我才刚刚开始学习 DSP，但我还没有找到任何 documentation/information 类似的东西，所以除了查看系数之外，我不知道从哪里开始（无济于事）。到目前为止，这是孤立的相关代码：

from scipy import signal
notch_freqs = (220, 440, 880) # could be any frequencies, these are just placeholders
notch_coeffs = [signal.iirnotch(freq, BANDWIDTH, samplerate) for freq in notch_freqs]

# here would be the step in applying the dry/wet to each coefficient in notch_coeffs

# and then the convolution would happen...
def convolve_notches():  # I know this can probably be optimized but it works for now
    filter_conv_a = signal.fftconvolve(*[filt[0] for filt in notch_coeffs[:2]])
    filter_conv_b = signal.fftconvolve(*[filt[1] for filt in notch_coeffs[:2]])

    if len(notch_coeffs) > 2:
         for filt in notch_coeffs[:2]:
            filter_conv_a = signal.fftconvolve(filter_conv_a, filt[0])
            filter_conv_b = signal.fftconvolve(filter_conv_b, filt[1])
    return filter_conv_a, filter_conv_b

# ... then apply using signal.filtfilt

这不是实时音频应用程序，所以我不太在意计时。重读这篇文章，我想澄清一下dry/wet；这是一个 ASCII 表示：

Wet: 25%       50%         100%
---\    /----\     /----\       /----
    \__/      \   /      \     /
               \_/        \   /
                           | |
                           | |

Answer 1

这是您可以执行此操作的一种方法。首先，一些观察结果：

要包含“湿”因素，我们可以使用一点代数。 IIR 滤波器可以表示为有理函数
```
        P_B(z)
H(z) =  ------
        P_A(z)
```
其中 P_A(z) 和 P_B(z) 是 z⁻¹ 中的多项式。陷波滤波器在中心频率处的增益为 0，但您想要混合“湿”和“干”（即过滤和未过滤）信号。令 w 为“湿”分数，0 <= w <= 1，因此 w = 0 表示没有滤波，w=0.25 表示在陷波频率处增益为 0.75，等这样的过滤器可以表示为
```
                        P_B(z)   (1 - w)*P_A(z) + w*P_B(z)
H(z) =  (1 - w) + (w) * ------ = -------------------------
                        P_A(z)            P_A(z)
```
就调用这样的函数返回的系数数组而言作为b, a = iirnotch(w0, Q, fs)，修改后的系数过滤器是 b_mod = (1-w)*a + w*b 和 a_mod = a.
iirnotch(w0, Q, fs) returns 两个长度为 3 的一维数组。这些是“双二阶”陷波滤波器的系数。去创造一个新的滤波器，它是几个陷波滤波器的级联，你可以简单地将 iirnotch 返回的数组堆叠成一个形状为 (n, 6) 的数组；这是 SOS（二阶部分）过滤器的格式。这种格式实际上是推荐的过滤器的格式超过几个订单，因为它是数字比 (b, a) 格式更健壮（需要高阶多项式）。

这是一个演示如何使用这些点的脚本实现你想要的级联陷波滤波器每个档位的“湿度”：

import numpy as np
from scipy.signal import iirnotch, sosfreqz
import matplotlib.pyplot as plt


samplerate = 4000
notch_freqs = (220, 440, 880)
Q = 12.5

filters = [iirnotch(freq, Q, samplerate) for freq in notch_freqs]

# Stack the filter coefficients into an array with shape
# (len(notch_freqs, 6)).  This array of "second order sections"
# can be used with sosfilt, sosfreqz, etc.

# This version, `sos`, is the combined filter without taking
# into account the desired "wet" factor.  This is created just
# for comparison to the filter with "wet" factors.
sos = np.block([[b, a] for b, a in filters])

# sos2 includes the desired "wet" factor of the notch filters.
wet = (0.25, 0.5, 1.0)
sos2 = np.block([[(1-w)*a + w*b, a] for w, (b, a) in zip(wet, filters)])

# Compute the frequency responses of the two filters.
w, h = sosfreqz(sos, worN=8000, fs=samplerate)
w2, h2 = sosfreqz(sos2, worN=8000, fs=samplerate)

plt.subplot(2, 1, 1)
plt.plot(w, np.abs(h), '--', alpha=0.75, label='Simple cascade')
plt.plot(w2, np.abs(h2), alpha=0.8, label='With "wet" factor')
plt.title('Frequency response of cascaded notch filters')
for yw in wet:
    plt.axhline(1 - yw, color='k', alpha=0.4, linestyle=':')
plt.grid(alpha=0.25)
plt.ylabel('|H(f)|')
plt.legend(framealpha=1, shadow=True, loc=(0.6, 0.25))

plt.subplot(2, 1, 2)
plt.plot(w, np.angle(h), '--', alpha=0.75, label='Simple cascade')
plt.plot(w2, np.angle(h2), alpha=0.8, label='With "wet" factor')
plt.xlabel('frequency')
plt.grid(alpha=0.25)
plt.ylabel('∠H(f) (radians)')
plt.legend(framealpha=1, shadow=True, loc=(0.6, 0.25))
plt.show()

这是脚本生成的情节：

Dry/wet（即包括部分原始信号）IIR 陷波滤波器在与其他滤波器卷积之前的参数

Dry/wet (i.e. include partial raw signal) parameter on IIR notch filter before convolution with other filters

python

signal-processing

scipy