Как создать режекторный фильтр с отдельными регуляторами ширины и центральной частоты

Регулятор центра задаёт частоту максимального подавления. Регулятор ширины задаёт расстояние между двумя границами половинной мощности, которые обычно называют точками −3 дБ. Режекторный фильтр — это узкополосная разновидность полосового заграждающего фильтра; в официальной документации SciPy добротность определяется как отношение центральной частоты к ширине этой полосы. SciPy: iirnotch

Перед настройкой регуляторов выберите единицу измерения ширины. Ширина в герцах сохраняет абсолютный частотный интервал. Ширина в октавах задаёт частотное отношение. Рецепт аудиоэквалайзера W3C явно поддерживает ширину полосы в октавах для режекторных фильтров и учитывает искажение частоты в формулах цифровых коэффициентов. Это разные варианты поведения регулятора, поэтому чётко укажите выбранную единицу. W3C: рецепт аудиоэквалайзера

В этом примере используйте герцы. Пусть центр равен 1 000 Гц, а ширина — 200 Гц: соответствующая добротность Q равна 5. Если переместить центр на 2 000 Гц, сохранив ширину 200 Гц, Q должна стать равной 10. Если оставить Q равной 5, ширина полосы составит 400 Гц. Эти наглядные расчёты объясняют, почему регулятор с фиксированной Q не решает поставленную задачу.

Используйте фильтр второго порядка с бесконечной импульсной характеристикой, также называемый биквадом. У него три коэффициента числителя и три коэффициента знаменателя. После нормализации первого коэффициента знаменателя до единицы уравнение для отсчёта имеет вид:\n\n\[\ny[n]=b_0x[n]+b_1x[n-1]+b_2x[n-2]-a_1y[n-1]-a_2y[n-2].\n\]\n\nЭта связь для биквада описана в рецепте аудиоэквалайзера. W3C: рецепт аудиоэквалайзера\n\nДля частоты дискретизации \(F_s\), центральной частоты \(f_c\) и ширины \(B\) вычислите:\n\n\[\n\omega_c=2\pi f_c/F_s,\qquad\n\beta=\tan(\pi B/F_s),\qquad\ng=\frac{1}{1+\beta}.\n\]\n\nЗатем используйте:\n\n\[\nb=g[1,-2\cos(\omega_c),1],\n\]\n\n\[\na=[1,-2g\cos(\omega_c),2g-1].\n\]\n\nЭто уравнения коэффициентов режекторного фильтра, используемые в реализации SciPy, переписанные так, чтобы ширина полосы в герцах была непосредственным входным параметром. Исходный код проектирования фильтров SciPy\n\nПри изменении каждого из регуляторов могут меняться несколько коэффициентов. Их независимость означает независимость характеристик отклика: при изменении ширины частота подавления остаётся прежней, а при изменении центра сохраняется заданная ширина полосы. Привязка каждого регулятора лишь к одному коэффициенту такого поведения не обеспечит.

Сохраните следующую полную программу в файле band_reject.py. Для неё требуется Python 3.8 или новее, дополнительные пакеты не нужны. Частота дискретизации — 48 000 Гц, а сгенерированный входной сигнал содержит три тона: 300, 1 000 и 5 000 Гц. Амплитуда каждого тона — 0,2.

Программа напрямую вычисляет частотную характеристику фильтра, методом бисекции находит обе границы половинной мощности и обрабатывает две секунды сгенерированного аудио. Для записи WAV-файлов используется документированный интерфейс Python wave, сохраняющий монофонические файлы в формате PCM с разрядностью 16 бит. Если выходные файлы уже существуют, программа выдаёт ошибку вместо перезаписи. Python: wave — чтение и запись WAV-файлов

import argparse
import cmath
import math
import struct
import wave


def coefficients(center, width, fs):
    if not all(math.isfinite(v) for v in (center, width, fs)):
        raise ValueError('Parameters must be finite')
    if not (fs > 0 and 0 < center < fs / 2 and 0 < width < fs / 2):
        raise ValueError('Require fs > 0 and center, width inside (0, fs/2)')
    beta = math.tan(math.pi * width / fs)
    g = 1 / (1 + beta)
    c = math.cos(2 * math.pi * center / fs)
    return (g, -2 * g * c, g), (1.0, -2 * g * c, 2 * g - 1)


def magnitude(f, b, a, fs):
    z = cmath.exp(-2j * math.pi * f / fs)
    return abs((b[0] + b[1]*z + b[2]*z*z) /
               (a[0] + a[1]*z + a[2]*z*z))


def edge(lo, hi, b, a, fs):
    target = 1 / math.sqrt(2)
    for _ in range(60):
        mid = (lo + hi) / 2
        if ((magnitude(lo, b, a, fs) - target) *
                (magnitude(mid, b, a, fs) - target)) > 0:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2


def process(samples, b, a):
    x1 = x2 = y1 = y2 = 0.0
    result = []
    for x in samples:
        y = b[0]*x + b[1]*x1 + b[2]*x2 - a[1]*y1 - a[2]*y2
        result.append(y)
        x2, x1, y2, y1 = x1, x, y1, y
    return result


def write_wav(path, samples, fs):
    if any(not math.isfinite(x) or abs(x) > 1 for x in samples):
        raise ValueError('Output exceeds PCM range; lower the input level')
    data = b''.join(struct.pack('<h', round(x * 32767)) for x in samples)
    with open(path, 'xb') as stream:
        with wave.open(stream, 'wb') as wav:
            wav.setnchannels(1)
            wav.setsampwidth(2)
            wav.setframerate(fs)
            wav.writeframes(data)


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('--center', type=float, default=1000)
    parser.add_argument('--width', type=float, default=200)
    parser.add_argument('--write', action='store_true')
    args = parser.parse_args()
    fs = 48000
    b, a = coefficients(args.center, args.width, fs)
    low = edge(0, args.center, b, a, fs)
    high = edge(args.center, fs / 2, b, a, fs)
    print(f'Edges: {low:.3f}, {high:.3f} Hz; width: {high-low:.3f} Hz')
    assert math.isclose(high-low, args.width, abs_tol=1e-5)
    assert magnitude(args.center, b, a, fs) < 1e-7
    tones = (300, 1000, 5000)
    for f in tones:
        db = 20 * math.log10(max(magnitude(f, b, a, fs), 1e-12))
        print(f'{f} Hz: {db:.2f} dB')
    samples = [sum(0.2 * math.sin(2*math.pi*f*n/fs) for f in tones)
               for n in range(2*fs)]
    filtered = process(samples, b, a)
    assert len(filtered) == len(samples)
    assert all(math.isfinite(y) for y in filtered)
    if args.write:
        write_wav('input-demo.wav', samples, fs)
        write_wav('filtered-demo.wav', filtered, fs)


if __name__ == '__main__':
    main()

Запустите python band_reject.py --center 1000 --width 200, чтобы проверить характеристику. Добавьте --write, чтобы создать input-demo.wav и filtered-demo.wav. Импортируйте эти файлы в аудиоредактор на отдельные дорожки, выровняйте их начала и сравните с неизменённой громкостью воспроизведения.

Я выполнил программу целиком без экспорта файлов. Проверки ширины полосы, подавления на центральной частоте, длины выходных данных и конечности значений прошли. Запись WAV-файлов при этом запуске не проверялась. Выведенное значение отклика составляло приблизительно −0,02 дБ на 300 Гц и −0,01 дБ на 5 000 Гц; значение на центральной частоте достигло нижнего предела отображения −240 дБ. Этот предел обусловлен численным представлением и не является измеренной глубиной подавления в аудиофайле.

Проверьте три набора параметров: исходный, с изменённым центром и прежней шириной, а также с изменённой шириной и прежним центром. Так можно изолированно проверить каждый регулятор. Следующие значения получены при выполнении процедур расчёта коэффициентов и поиска границ на частоте 48 кГц; они описывают математическую характеристику, а не результат измерения записи.

Настройка центра	Настройка ширины	Нижняя граница половинной мощности	Верхняя граница половинной мощности	Расчётный интервал
1 000 Гц	200 Гц	904,959 Гц	1 104,959 Гц	200,000 Гц
2 000 Гц	200 Гц	1 902,441 Гц	2 102,441 Гц	200,000 Гц
1 000 Гц	400 Гц	819,692 Гц	1 219,692 Гц	400,000 Гц

Обратите внимание: границы не совпадают точно с центр ± ширина/2. Здесь «центр» обозначает нулевой уровень подавления, а «ширина» — расстояние между пересечениями уровня половинной мощности. Расчётная характеристика показывает, что их арифметическая середина может не совпадать с частотой нулевого уровня. Не рассчитывайте границы, просто прибавляя и вычитая половину ширины.

Для дополнительной численной проверки в среде с установленным SciPy вычислите отклик по возвращённым коэффициентам с помощью scipy.signal.freqz. Передайте fs и явно задайте частоты, включая точное значение центра. В документации поясняется, что явный массив частот задаёт точки вычисления отклика; разреженная сетка по умолчанию может пропустить узкий нулевой уровень. SciPy: freqz

Возьмите короткий повторяемый фрагмент обрабатываемого звука. Сначала определите частотную область, которую нужно ослабить, затем задайте в ней центр. Начните с небольшой ширины и постепенно увеличивайте её, пока подавление не охватит нужную область. После каждого изменения сравнивайте результат с обходом фильтра, сохраняя одинаковую громкость воспроизведения, чтобы сравнение оставалось полезным.

Используйте эту таблицу, чтобы выбрать следующее изменение:

Что показывает сравнение	Следующее изменение
Нежелательный тон остаётся, но соседние составляющие меняются	Прежде чем увеличивать ширину, перепроверьте центр
Центр выбран верно, но нежелательная область шире	Постепенно увеличьте ширину
Исчезает слишком много тембра соседних составляющих	Уменьшите ширину и сравните ещё раз
Нужно ослабить несколько разнесённых тонов	Проверьте отдельные режекторные секции на этих частотах
Требуется широкая и ровная полоса подавления	Рассмотрите другой тип и порядок фильтра

Это практические рекомендации по настройке, а не измерения конкретной записи. Пример с тремя тонами намеренно прост: он позволяет убедиться в правильной работе процессора, прежде чем оценивать более сложный звук.

Один режекторный фильтр также не позволяет независимо регулировать глубину подавления. В рецепте аудиоэквалайзера приведены отдельные наборы коэффициентов для режекторных фильтров и фильтров пикового типа. Если вместо полного подавления нужно умеренное ослабление с регулировкой, стоит рассмотреть фильтр пикового типа. W3C: рецепт аудиоэквалайзера

Сохраняйте тот же контракт параметров: center_hz и width_hz одновременно передаются в функцию коэффициентов. При изменении любого параметра обновляйте весь набор коэффициентов, а при изменении частоты дискретизации обработки пересчитывайте их заново. В этом прототипе оба параметра должны быть строго больше нуля и меньше половины частоты дискретизации; при недопустимых значениях возникает исключение.

Предоставленная функция process() обрабатывает целый монофонический буфер и при каждом вызове начинает с нулевой историей. Потоковая версия должна сохранять историю задержек между блоками. В документации SciPy для lfilter описаны начальное состояние задержки и возвращаемое конечное состояние, необходимые для этого. Для каждого аудиоканала храните отдельную историю. SciPy: lfilter

Перемещение регуляторов требует отдельного этапа реализации. Этот пример проверяет фиксированные настройки и не поддерживает непрерывную автоматизацию. Один из практических способов перехода — ненадолго запустить параллельно старый и новый фильтры с фиксированными настройками, дать новому фильтру накопить историю и выполнить кроссфейд выходных сигналов. Проверяйте переход на реальном материале и с предполагаемой скоростью изменения параметров.

При последовательном соединении нескольких фильтров сохраняйте каждый из них как отдельную секцию второго порядка, не объединяя каскад в один полином высокого порядка. SciPy рекомендует секции второго порядка для большинства задач фильтрации, поскольку они создают меньше численных проблем. SciPy: lfilter

Если режекторный фильтр срабатывает не на той частоте, проверьте, совпадает ли частота дискретизации, использованная при расчёте коэффициентов, с фактической частотой дискретизации аудио. Если при перемещении центра ширина растёт, проверьте, не сохраняет ли ваша обёртка постоянную Q. Если на дисплее характеристика выглядит недостаточно глубокой, сначала вычислите отклик точно на центральной частоте, а затем меняйте конструкцию.

Если повторно обработанные буферы содержат скачки, проверьте управление состоянием: демонстрационный пример сбрасывает историю при каждом вызове process(). Если экспорт WAV завершается ошибкой диапазона, уменьшите амплитуду входного сигнала и запустите программу снова; она намеренно не допускает незаметного клиппинга значений за пределами допустимого диапазона. Перед экспортом новой пары файлов нужно переименовать или удалить уже существующие демонстрационные файлы.

Сборка завершена, когда проверка центра подтверждает подавление на заданной частоте, расстояние между двумя пересечениями границ равно заданной ширине, а обработанный звук подходит для предполагаемой правки. Сохраняйте эти проверки вместе с двумя регуляторами при переходе от автономного прототипа к интерактивному аудиоинструменту.

• Рецепт аудиоэквалайзера — W3C
• iirnotch — руководство SciPy
• _filter_design.py — исходный код SciPy версии 1.18.0
• freqz — руководство SciPy
• lfilter — руководство SciPy
• wave — чтение и запись WAV-файлов — документация Python