将滤波器方程转换为频率响应

Question

我被告知这对于滤波器设计者来说“应该很容易”。虽然我不太确定，但我还是会问:)

我需要根据这些方程构造一个函数，该函数采用一个频率范围并产生滤波器在该范围内的积分响应：

float filter_response(float f_low, float f_high) {
  return {magic};
}

滤波器的频率响应已知：

并且从表格中可以得知所需函数的预期输出以进行验证：

𝑓c	𝑓低	𝑓高	输出
4	3.5	4.5	0.375
5	4.4	5.6	0.545
6.3	5.6	7.0	0.727
8	7.1	8.9	0.873
10	8.8	11.2	0.951
12.5	11.1	13.9	0.958
16	14.2	17.8	0.896
20	17.7	22.3	0.782
25	22.1	27.9	0.647
31.5	27.9	35.1	0.519
40	35.4	44.6	0.411
50	44.2	55.8	0.324
63	55.8	70.2	0.256
80	70.8	89.2	0.202
100	88.5	111.5	0.16
125	110.6	139.4	0.127
160	141.6	178.4	0.101
200	177.0	223.0	0.0799
250	221.2	278.8	0.0634
315	278.8	351.2	0.0503
400	354.0	446.0	0.0398
500	442.5	557.5	0.0315
630	557.5	702.5	0.0245
800	708.0	892.0	0.0186
1000	885.0	1115.0	0.0135
1250	1106.2	1393.8	0.00894
1600	1416.0	1784.0	0.00538
2000	1770.0	2230.0	0.00295

其中𝑓c为中心频率，𝑓low-𝑓high频率范围，输出期望的输出。只有 𝑓c 最初是从表中抄录的：𝑓low-𝑓high 由我计算，代表给定“中”频率的三分之一倍频程的频率范围。为了清楚起见，我提供了这些内容，因为它们与所需的功能直接相关。

这是我拥有的数据，虽然我确信这是可以做到的，但我也明白这可能超出了在线获取帮助的预期范围（除非它真的像向我提供的那样简单）。

无论如何，任何正确方向的指示都将不胜感激。

谢谢你<3

由于我既不是滤波器设计者，也不精通数学，所以我主要研究从图表中提取数据或通过在表中插入给定值。虽然这目前有效，但它并没有给我带来我所寻求的准确性或可信度。我阅读几乎任何语言的代码，并且对 FIR 和 IIR 滤波器非常了解。

我尝试根据给定的滤波器常量使用 Python、SciPy、Numpy 设计一个滤波器，看看是否可以复制给定的响应曲线，但没有成功。

我还尝试过人工智能资源来帮助我简化这些公式，为我提供一些适合我在该领域有限知识的术语来描述过滤器，但它们似乎完全无法处理这个问题。

最后，我希望能够通过几种不同的方式来解决这个问题。虽然用数值计算函数会很棒，但从理论上讲，我还可以分析 FFT/FIR 滤波器以生成具有所需范围和精度的新表。

更新：

import math
import numpy as np
import matplotlib.pyplot as plt

# constants
pi, f_1, f_2, Q_1, f_3, f_4, Q_2, K = math.pi, 6.31, 1258.9, 0.71, 15.915, 15.915, 0.64, 1    

# band-limiting transfer function
def band_limit(s):
    return ((s**2)*4*(pi**2)*(f_2**2))/(((s**2)+2*pi*f_1*s/Q_1+4*(pi**2)*(f_1**2))*((s**2)+2*pi*f_2*s/Q_1+4*(pi**2)*(f_2**2)))

# frequency-weighting transfer function
def frequency_weight(s):
    return ((s+2*pi*f_3)*2*pi*K*(f_4**2))/(((s**2)+2*pi*f_4*s/Q_2+4*(pi**2)*(f_4**2))*f_3)

px, py, py1, py2 = [], [], [], []

for x in np.arange(1,1000,0.1):
    w = x*2*pi # angular frequency
    px.append(x)
    py1.append(band_limit(w))
    py2.append(frequency_weight(w))
    py.append(band_limit(w)*frequency_weight(w))

plt.plot(px, py1, label = "band-limit filter")
plt.plot(px, py2, label = "frequency-weighting filter")
plt.plot(px, py, label = "total")
plt.yscale("log");
plt.xscale("log");
plt.ylim(0.01, 1);
plt.xlim(1, 1000);
plt.legend()
plt.show()

这个 python 程序生成以下图：

我已经覆盖了原始图进行比较。

情节不匹配，但它们看起来并不完全不同，这告诉我我走在正确的道路上。是否可以使用纯实值函数来计算每个频率的幅度响应，或者我得到的错误是由于不执行复数计算而导致的？

Answer 1

看起来不错，虽然我不知道你的容忍度是多少：

from io import StringIO

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt


def Hb(
    s: np.ndarray,
    f1: float = 6.310,
    f2: float = 1258.9,
    Q1: float = 0.71,
) -> np.ndarray:
    from numpy import pi

    return (
        s * 2*pi * f2
    )**2 / (
        (
            s**2 + 2*pi*f1*s / Q1
            + (2*pi*f1)**2
        ) * (
            s**2 + 2*pi*f2*s / Q1
            + (2*pi*f2)**2
        )
    )


def Hw(
    s: np.ndarray,
    f3: float = 15.915,
    f4: float = 15.915,
    Q2: float = 0.64,
    K: int = 1,
) -> np.ndarray:
    from numpy import pi

    return (s + 2*pi*f3) * (
        2*pi * K * f4**2
    ) / (
        (
            s**2 + 2*pi*f4*s / Q2
            + (2*pi*f4)**2
        ) * f3
    )


def verify_band_limiting() -> None:
    content = '''
fc  flow    fhigh   output
4   3.5     4.5     0.375
5   4.4     5.6     0.545
6.3     5.6     7.0     0.727
8   7.1     8.9     0.873
10  8.8     11.2    0.951
12.5    11.1    13.9    0.958
16  14.2    17.8    0.896
20  17.7    22.3    0.782
25  22.1    27.9    0.647
31.5    27.9    35.1    0.519
40  35.4    44.6    0.411
50  44.2    55.8    0.324
63  55.8    70.2    0.256
80  70.8    89.2    0.202
100     88.5    111.5   0.16
125     110.6   139.4   0.127
160     141.6   178.4   0.101
200     177.0   223.0   0.0799
250     221.2   278.8   0.0634
315     278.8   351.2   0.0503
400     354.0   446.0   0.0398
500     442.5   557.5   0.0315
630     557.5   702.5   0.0245
800     708.0   892.0   0.0186
1000    885.0   1115.0  0.0135
1250    1106.2  1393.8  0.00894
1600    1416.0  1784.0  0.00538
2000    1770.0  2230.0  0.00295
'''
    with StringIO(content) as f:
        expected = pd.read_csv(f, delim_whitespace=True)

    s = 2j*np.pi * expected['fc']
    out = Hb(s) * Hw(s)
    out_mag = np.abs(out)

    assert np.allclose(
        a=out_mag,
        b=expected['output'],
        rtol=0, atol=1e-2,
    )

    fig, (ax_transfer, ax_err) = plt.subplots(ncols=2)
    ax_transfer.loglog(expected['fc'], out_mag, label='actual')
    ax_transfer.loglog(expected['fc'], expected['output'], label='expected')
    ax_transfer.set_xlabel('Frequency (Hz)')
    ax_transfer.set_ylabel('Magnitude')
    ax_transfer.legend()
    ax_err.semilogx(expected['fc'], out_mag - expected['output'])
    ax_err.set_xlabel('Frequency (Hz)')
    ax_err.set_ylabel('Magnitude error')
    
    plt.show()


if __name__ == '__main__':
    verify_band_limiting()