Цифровые фильтры: импульсная характеристика, КИХ и БИХ

Однополюсный фильтр (one-pole filter)?однополюсный фильтр (one-pole filter) — Простейший рекурсивный фильтр: новый выход вычисляется из входного отсчёта и одного предыдущего выхода, y[n] = y[n−1] + a·(x[n] − y[n−1]). При коэффициенте a от 0 до 1 работает как фильтр нижних частот. описан в своей статье во времени: новый выход вычисляется из входного отсчёта (sample)?отсчёт (sample) — Значение сигнала в один момент дискретизации, одно число. В звукозаписи его также называют сэмплом. и прошлого выхода. По такой формуле трудно сразу ответить, насколько фильтр ослабит синусоиду (sine)?синус (sine) — Вертикальная координата точки на единичной окружности под заданным углом от положительного направления горизонтальной оси. Для острого угла прямоугольного треугольника равен отношению противолежащего катета к гипотенузе. 5 кГц и на сколько её задержит. Ответ дают средства, общие для всех линейных фильтров: импульсная характеристика и свёртка, частотная характеристика, передаточная функция с нулями и полюсами, групповая задержка. С их помощью разобраны фильтры с конечной и бесконечной импульсной характеристикой и биквад — фильтр второго порядка, из которого собирают эквалайзеры. В конце — класс Biquad на C++23 и фильтры JUCE. Понадобятся комплексные числа и формула Эйлера.

Линейность и инвариантность во времени

Фильтр превращает последовательность входных отсчётов $x[n]$ в выходные $y[n]$. Всё, что сказано в статье дальше, верно для фильтров, которые линейны и инвариантны во времени.

Линейность (linearity): если вход $x_1$ даёт выход $y_1$, а вход $x_2$ — выход $y_2$, то вход $c_1x_1 + c_2x_2$ даёт $c_1y_1 + c_2y_2$ при любых числах $c_1$ и $c_2$ (Sundararajan, §2.3.2). Вдвое более громкий вход даёт вдвое более громкий выход, а сумма двух сигналов проходит фильтр так, будто каждый прошёл отдельно и выходы сложили. Инвариантность во времени (time invariance): если вход $x[n]$ даёт $y[n]$, то задержанный вход $x[n - D]$ даёт тот же выход, задержанный на $D$ отсчётов (там же, §2.3.3). Коэффициенты фильтра не меняются, и ему всё равно, когда пришёл сигнал. Системы с обоими свойствами называют линейными инвариантными во времени (linear time-invariant, LTI).

Линейность проверяется на числах. Однополюсный фильтр с $a = 0{,}5$ из нулевого состояния превращает вход 0,2; 0,2; 0,2 в 0,1; 0,15; 0,175, а вход 0,4; 0,4; 0,4 — в 0,2; 0,3; 0,35, ровно вдвое больше. Возведение в квадрат $y[n] = x^2[n]$ так себя не ведёт: 0,2 даёт 0,04, а 0,4 — 0,16, вчетверо больше. Квадрат — нелинейная операция.

Импульсная характеристика и свёртка

Единичный импульс $\delta[n]$ — последовательность, в которой $\delta[0] = 1$, а все остальные отсчёты равны нулю: цифровой щелчок. Импульсная характеристика (impulse response, IR)?импульсная характеристика (impulse response, IR) — Отклик линейной инвариантной во времени системы на единичный импульс из нулевого состояния. Полностью описывает такую систему: выход на любой вход — свёртка входа с импульсной характеристикой. $h[n]$ — выход фильтра на такой вход из нулевого состояния. У однополюсного фильтра $y[n] = a\,x[n] + (1 - a)\,y[n-1]$ её легко посчитать по шагам: $h[0] = a$, $h[1] = a(1 - a)$, $h[2] = a(1 - a)^2$, и вообще $h[n] = a(1 - a)^n$. При частоте среза 1000 Гц и 48 кГц коэффициент $a \approx 0{,}1227$, и первые отсчёты равны 0,1227; 0,1076; 0,0944; 0,0828 — каждый в $1 - a \approx 0{,}877$ раза меньше предыдущего.

Свёртка

Любой вход можно записать как сумму задержанных импульсов с весами — значениями отсчётов. Вход 2, 1, 3, 4 — это $2\delta[n] + \delta[n-1] + 3\delta[n-2] + 4\delta[n-3]$. Инвариантность во времени говорит, что импульс $\delta[n - k]$ даёт отклик $h[n - k]$, а линейность — что отклик на сумму равен сумме откликов с теми же весами. Отсюда выход для любого входа (Sundararajan, §2.2, с. 47):

Такая операция называется свёрткой (convolution)?свёртка (convolution) — Операция y[n] = Σ h[k]·x[n−k]: каждый выходной отсчёт — сумма текущего и прошлых входных отсчётов с весами из импульсной характеристики. Так линейный инвариантный во времени фильтр вычисляет выход. и записывается $y = x * h$. Обе формы суммы дают одно и то же: от перестановки последовательностей свёртка не меняется. Вторая форма прямо задаёт вычисление: выход — это текущий и прошлые входы, взвешенные отсчётами импульсной характеристики.

Пример из учебника Sundararajan: $h$ = 3, 1, 2, 4 и $x$ = 2, 1, 3, 4. Отсчёт $y[3]$ собирает четыре произведения: $x[0]h[3] + x[1]h[2] + x[2]h[1] + x[3]h[0] = 2 \cdot 4 + 1 \cdot 2 + 3 \cdot 1 + 4 \cdot 3 = 25$. Весь выход — 6, 5, 14, 25, 14, 20, 16: свёртка последовательностей длиной $N$ и $M$ содержит $N + M - 1$ отсчётов (Sundararajan, §2.2, с. 47–48). Проверка: сумма выхода, 100, равна произведению сумм входа и характеристики, $10 \cdot 10$.

В визуализации 1 характеристика развёрнута во времени: значение $h[m]$ стоит под входом $x[n - m]$. С ростом $n$ она сдвигается вправо, и под ней оказываются всё новые входные отсчёты. В начальном состоянии $n = 3$, перекрываются все четыре отсчёта, и сумма равна 25. Выберите вход «щелчок»: выход повторит $h$, потому что импульсная характеристика и есть отклик на щелчок. С характеристикой «эхо» каждый входной отсчёт возвращается через три отсчёта вдвое тише — так работает линия задержки (delay line)?линия задержки (delay line) — Блок обработки, который хранит последние отсчёты сигнала и выдаёт вход на заданное число отсчётов позже: y[n] = x[n − D]. На нём построены эхо, хорус, фленжер, реверберация и компенсация задержки. с прямой связью. Вход «скачок» с характеристикой «разность» даёт 1 в начале, −1 в конце и нули посередине: разность реагирует только на изменения.

КИХ-фильтры

Если импульсная характеристика конечна, то есть после $M$ отсчётов идут только нули, формула (1) содержит $M$ слагаемых, и фильтр вычисляется прямо по ней: $y[n] = b_0x[n] + b_1x[n-1] + \ldots + b_{M-1}x[n - M + 1]$. Такой фильтр называют КИХ-фильтром (finite impulse response filter, FIR)?КИХ-фильтр (finite impulse response filter, FIR) — Фильтр с конечной импульсной характеристикой: выход — взвешенная сумма текущего и нескольких прошлых входов, без обратной связи. Всегда устойчив; при симметричных коэффициентах задерживает все частоты одинаково., фильтром с конечной импульсной характеристикой. Его коэффициенты и есть отсчёты $h$: если подать щелчок, на выход по очереди выйдут $b_0$, $b_1$ и так далее (Sundararajan, §7.1, с. 230–231).

Простейший КИХ-фильтр — скользящее среднее $M$ отсчётов: все коэффициенты равны $1/M$. Среднее двух соседних отсчётов, $y[n] = (x[n] + x[n-1])/2$, сглаживает сигнал: медленные изменения проходят почти без потерь, а чередование знаков +1, −1, +1 даёт на выходе ноль. Насколько ослабляются промежуточные частоты, покажет частотная характеристика.

КИХ-фильтр всегда устойчив: выход — сумма $M$ ограниченных слагаемых, и расти без предела ему нечем (Sundararajan, §7.1, с. 232). Зато ему нужно $M$ умножений на каждый отсчёт. Документация JUCE советует считать FIR::Filter по формуле, пока коэффициентов меньше 128, а более длинные характеристики — через БПФ в классе Convolution (juce_FIRFilter.h).

Если коэффициенты симметричны, $b_k = b_{M-1-k}$, как у скользящего среднего, фильтр задерживает все частоты на одно и то же время, $(M - 1)/2$ отсчётов. Это свойство называют линейной фазой; оно разобрано в разделе о задержке.

БИХ-фильтры

Однополюсный фильтр устроен иначе: его выход зависит от прошлого выхода. Общая запись фильтра с обратной связью (feedback)?обратная связь (feedback) — Возврат части выходного сигнала на вход обработки. В delay повторно направляет копию через линию задержки и создаёт цепочку повторений. — разностное уравнение (Sundararajan, §2.3.1 и §2.7):

Если хоть один $a_k$ отличен от нуля, импульсная характеристика бесконечна: каждый выход возвращается во все следующие. Такой фильтр называют БИХ-фильтром (infinite impulse response filter, IIR)?БИХ-фильтр (infinite impulse response filter, IIR) — Фильтр с бесконечной импульсной характеристикой: в разностное уравнение входят прошлые выходы. Даёт крутой спад малым числом коэффициентов, устойчив, только если все полюсы внутри единичной окружности., фильтром с бесконечной импульсной характеристикой (Sundararajan, §2.7, с. 60). Знак минус перед второй суммой — соглашение, с которым проще выглядит передаточная функция; так записывают коэффициенты Sundararajan, RBJ Audio EQ Cookbook и JUCE. Однополюсный фильтр — БИХ-фильтр первого порядка с $b_0 = a$ и $a_1 = -(1 - a)$. Здесь $a$ без индекса — коэффициент из статьи об однополюсном фильтре, а $a_1$ — коэффициент обратной связи; это разные величины.

Обратная связь экономит вычисления. Импульсная характеристика однополюсного фильтра с частотой среза 1000 Гц опускается на 60 дБ ниже первого отсчёта только через 53 отсчёта. КИХ-фильтру с той же характеристикой, обрезанной на этом уровне, понадобилось бы 53 умножения на отсчёт, а однополюсному хватает двух. При 200 Гц КИХ-фильтру понадобилось бы 264 умножения. Sundararajan подводит итог: фильтрам с обратной связью для тех же требований нужно меньше коэффициентов, но у них появляется вопрос устойчивости (Sundararajan, §2.7, с. 59).

Устойчивость

Фильтр устойчив, если ограниченный вход всегда даёт ограниченный выход (определение из статьи об однополюсном фильтре). Через импульсную характеристику условие записывается так: сумма модулей $\sum_n |h[n]|$ конечна (Sundararajan, §2.5, с. 57). У КИХ-фильтра она конечна всегда. У однополюсного фильтра это геометрическая прогрессия $a\,(1 + |1 - a| + |1 - a|^2 + \ldots)$, и она конечна, только если $|1 - a| < 1$. При $1 - a = 1{,}1$ каждый отсчёт характеристики на 10 % больше предыдущего, и через 1000 отсчётов он больше $10^{41}$. Общее правило для БИХ-фильтров любого порядка дают полюсы.

Частотная характеристика

Подадим на LTI-фильтр комплексную синусоиду $x[n] = e^{j\omega n}$, где $\omega = 2\pi f/f_s$ — частота в радианах на отсчёт, а $f_s$ — частота дискретизации (sample rate)?частота дискретизации (sample rate) — Число отсчётов в каждом канале за одну секунду, то же число кадров аудиоданных в секунду. Измеряется в герцах: 44 100 Гц означает 44 100 отсчётов на канал в секунду. (вращающийся вектор). По формуле (1) выход равен $\sum_k h[k]\,e^{j\omega(n-k)} = e^{j\omega n}\sum_k h[k]\,e^{-j\omega k}$. Вторая сумма не зависит от $n$: это комплексное число (complex number)?комплексное число (complex number) — Число вида a + jb, где a и b вещественные, а мнимая единица j удовлетворяет j² = −1. Изображается точкой или вектором на плоскости; модуль и аргумент хранят длину и угол, например амплитуду и фазу синусоиды., одно для каждой частоты (Sundararajan, §2.4, с. 56–57).

Синусоида выходит из фильтра той же частоты, умноженная на $H(e^{j\omega})$. Для вещественной синусоиды это значит: $\cos(\omega n)$ превращается в $|H|\cos(\omega n + \arg H)$ — амплитуда умножается на модуль, к фазе прибавляется аргумент (Sundararajan, §2.4, с. 57). Модуль $|H(e^{j\omega})|$ как функцию частоты называют амплитудно-частотной характеристикой (magnitude response)?амплитудно-частотная характеристика (magnitude response) — Зависимость модуля частотной характеристики |H(e^jω)| от частоты: во сколько раз фильтр меняет амплитуду синусоиды каждой частоты. Обычно рисуется в дБ., АЧХ (magnitude response), аргумент — фазочастотной, ФЧХ (phase response). АЧХ обычно рисуют в децибелах (decibel)?децибел (decibel) — Логарифмическая единица отношения двух величин: для амплитуд L = 20 lg(A/A₀), для мощностей — 10 lg(P/P₀). Смысл числа задаёт опорная величина.: $20\lg|H|$.

Отсюда следует, что LTI-фильтр не создаёт новых частот. Каждая синусоида входа выходит той же частоты, меняются только её амплитуда (amplitude)?амплитуда (amplitude) — Наибольший модуль отклонения колеблющейся величины от равновесного значения. Для синусоиды обозначается A; единицы совпадают с единицами сигнала. и фаза (phase)?фаза (phase) — Положение внутри периода колебания, выраженное углом: полный период соответствует 360°, или 2π радиан. Начальная фаза — её значение в момент t = 0.. Новые частоты появляются только в нелинейной обработке: срез добавляет гармоники и комбинационные частоты.

У среднего двух отсчётов $h$ = 0,5; 0,5, и по формуле (3) $H = (1 + e^{-j\omega})/2$. Если вынести $e^{-j\omega/2}$, в скобках останется полусумма $e^{j\omega/2} + e^{-j\omega/2}$, то есть $\cos(\omega/2)$ по формуле Эйлера: $H = e^{-j\omega/2}\cos(\omega/2)$. АЧХ — $\cos(\omega/2)$: 1 на нулевой частоте, $1/\sqrt{2} \approx 0{,}707$ (−3,01 дБ) на $f_s/4$ — 12 кГц при 48 кГц — и 0 на $f_s/2$. ФЧХ — прямая $-\omega/2$: все частоты сдвинуты на половину отсчёта.

Проверка однополюсного фильтра

У однополюсного фильтра $h[k] = a(1 - a)^k$, и формула (3) — сумма геометрической прогрессии со знаменателем $(1 - a)e^{-j\omega}$. Её модуль меньше 1, и сумма равна $a/\bigl(1 - (1 - a)e^{-j\omega}\bigr)$. Квадрат модуля:

В статье об однополюсном фильтре сказано: если частота среза много ниже $f_s/2$, синусоида этой частоты ослабляется примерно на 3 дБ. Теперь это можно проверить. Обозначим $p = 1 - a = e^{-\omega_c}$, где $\omega_c = 2\pi f_c/f_s$. Знаменатель формулы (4) равен $(1 - p)^2 + 2p\,(1 - \cos\omega)$. При малой $\omega_c$ выполняется $1 - p = a \approx \omega_c$ и $1 - \cos\omega_c \approx \omega_c^2/2$, а $p \approx 1$. На частоте среза знаменатель близок к $a^2 + \omega_c^2 \approx 2a^2$, и $|H|^2 \approx 1/2$: ослабление 3 дБ. Приближение требует $\omega_c \ll 1$, то есть $f_c \ll f_s/(2\pi)$. Таблица 1 даёт точные значения по формуле (4).

На $f_s/2$, где $\cos\omega = -1$, формула (4) даёт $|H| = a/(2 - a)$: при частоте среза 1000 Гц ослабление всего 23,7 дБ, как и сказано в разделе об ограничениях. Рисунок 3 сравнивает оба фильтра.

Передаточная функция и z-преобразование

Частотная характеристика (3) — сумма степеней $e^{-j\omega}$. Если обозначить $z = e^{j\omega}$, задержка на $k$ отсчётов превращается в множитель $z^{-k}$. На этом построено z-преобразование (z-transform)?z-преобразование (z-transform) — Запись последовательности x[n] выражением X(z) = Σ x[n]·z⁻ⁿ, в которой задержка на один отсчёт становится умножением на z⁻¹. Превращает разностное уравнение фильтра в передаточную функцию.: последовательности $x[n]$ сопоставляют выражение $X(z) = \sum_n x[n]\,z^{-n}$, где $z$ — комплексное число (Sundararajan, §6.1, с. 196–197). Для фильтров важно одно его свойство: задержке на один отсчёт соответствует умножение на $z^{-1}$. Поэтому элемент задержки на схемах фильтров подписывают $z^{-1}$. В учебниках к z-преобразованию добавляют условие, при каких $z$ сумма сходится; для чтения передаточных функций оно не понадобится.

Применим это к разностному уравнению (2): каждое слагаемое $x[n-k]$ превращается в $z^{-k}X(z)$, каждое $y[n-k]$ — в $z^{-k}Y(z)$. Перенесём выходы влево и поделим на $X(z)$. Отношение называют передаточной функцией (transfer function)?передаточная функция (transfer function) — Отношение H(z) = Y(z)/X(z) многочленов от z⁻¹ с коэффициентами разностного уравнения фильтра. Корни числителя — нули, корни знаменателя — полюсы; на единичной окружности z = e^jω она равна частотной характеристике. (Sundararajan, §6.3.2, с. 217):

Числитель и знаменатель — многочлены, коэффициенты которых совпадают с коэффициентами разностного уравнения. На единичной окружности, при $z = e^{j\omega}$, передаточная функция совпадает с частотной характеристикой (3) (Sundararajan, §6.1, с. 196).

Нули и полюсы

Корни числителя называют нулями (zeros): в них $H(z) = 0$. Корни знаменателя называют полюсами (poles): знаменатель там обращается в ноль, и $H(z)$ бесконечна. Их рисуют на комплексной плоскости $z$: нуль — кружком, полюс — крестиком. Точки единичной окружности соответствуют частотам: угол 0 — постоянной составляющей, $\pi/2$ — частоте $f_s/4$, $\pi$ — частоте $f_s/2$, частоте Найквиста (Nyquist frequency)?частота Найквиста (Nyquist frequency) — Половина частоты дискретизации, fₛ/2. В принятой низкочастотной модели точное восстановление по отсчётам требует, чтобы все частоты исходного сигнала были строго ниже этой границы.. Нижняя половина окружности — отрицательные частоты. У фильтра с вещественными коэффициентами комплексные нули и полюсы стоят парами, симметрично относительно вещественной оси (комплексные корни квадратного уравнения).

У среднего двух отсчётов $H(z) = (1 + z^{-1})/2 = (z + 1)/(2z)$: нуль в точке $z = -1$, то есть на $f_s/2$, где АЧХ обнулилась, и полюс в начале координат. До начала координат от любой точки окружности расстояние 1, и на АЧХ он не влияет. У однополюсного фильтра $H(z) = a/\bigl(1 - (1 - a)z^{-1}\bigr) = az/\bigl(z - (1 - a)\bigr)$: полюс на вещественной оси в точке $z = 1 - a$, при частоте среза 1000 Гц — в точке 0,877. Он стоит у угла 0, рядом с низкими частотами, поэтому фильтр их пропускает лучше, чем высокие.

Отсюда общее условие устойчивости: фильтр устойчив, если все его полюсы лежат внутри единичной окружности, $|p| < 1$ (Sundararajan, §6.3.2, с. 218). У однополюсного фильтра это $|1 - a| < 1$ — то же условие, что получилось из прогрессии. У КИХ-фильтра знаменатель равен 1, все его полюсы в начале координат, и он устойчив всегда.

Фазовая и групповая задержка

ФЧХ даёт фазовый сдвиг в радианах, а знать удобнее задержку во времени. Для этого есть две величины (J. O. Smith, Phase and Group Delay):

Фазовая задержка показывает, на сколько отсчётов сдвинута синусоида частоты $\omega$. Групповая задержка (group delay)?групповая задержка (group delay) — Производная фазочастотной характеристики по частоте с обратным знаком, −dφ/dω. Показывает, на сколько задерживается огибающая сигнала в узкой полосе около данной частоты. — производная фазы по частоте с обратным знаком. Она показывает задержку огибающей (envelope)?огибающая (envelope) — Плавная кривая, которая следует за размахом волны и показывает, как уровень сигнала меняется во времени.: короткой пачки колебаний около частоты $\omega$, «группы» близких синусоид (J. O. Smith, Group Delay). Удар барабана и начало ноты — именно такие огибающие.

Если ФЧХ — прямая $\varphi(\omega) = -D\omega$, обе задержки равны $D$ отсчётам на всех частотах. Тогда составляющие сигнала в полосе пропускания задержаны одинаково, и форма сигнала сохраняется: он просто приходит позже (Sundararajan, §7.2, с. 237). У среднего двух отсчётов $\varphi = -\omega/2$, задержка — половина отсчёта, при 48 кГц это 10,4 мкс.

Линейная фаза

У КИХ-фильтра с симметричными коэффициентами длины $M$ ФЧХ — прямая с наклоном $-(M - 1)/2$, а групповая задержка постоянна и равна $(M - 1)/2$ отсчётов (Sundararajan, §7.2, с. 235–237). На частотах, где АЧХ проходит через ноль, фаза скачком меняется на $\pi$, но прямой наклон и задержка от этого не меняются. У скользящего среднего восьми отсчётов задержка 3,5 отсчёта. Такие фильтры называют линейно-фазовыми.

У БИХ-фильтров ФЧХ — не прямая, и групповая задержка зависит от частоты. У однополюсного фильтра с частотой среза 1000 Гц она равна 7,15 отсчёта на нулевой частоте, 3,33 на 1 кГц, 1,04 на 2 кГц, а на 10 кГц становится даже немного отрицательной, −0,41. Групповая задержка описывает огибающую узкой полосы частот, а не момент, когда на выходе появляется сигнал, поэтому бывает и отрицательной; такие значения Bellanger отмечает и у секции второго порядка (Bellanger, §6.2, с. 131). Низкие частоты такой фильтр задерживает сильнее высоких.

Биквад

Фильтр второго порядка помнит два прошлых входа и два прошлых выхода:

Полюсы и добротность

Пусть полюсы — пара сопряжённых чисел $re^{\pm j\theta}$. По формуле корней квадратного уравнения знаменатель $z^2 + a_1z + a_2$ задают коэффициенты (Bellanger, §6.2, с. 128–129):

Угол $\theta$ задаёт частоту резонанса, а радиус $r$ — его остроту. Подъём АЧХ на резонансе обратно пропорционален расстоянию от полюса до окружности, $1 - r$ (Bellanger, §6.2, с. 129–130). Ширина резонанса по уровню −3 дБ при $r$, близком к 1, равна $B \approx -\ln r \cdot f_s/\pi \approx (1 - r)f_s/\pi$ (J. O. Smith, Resonator Bandwidth). Добротность (quality factor, Q)?добротность (quality factor, Q) — Отношение центральной частоты полосы к её ширине между частотами среза. Чем больше добротность, тем уже полоса. — отношение частоты к ширине полосы, $Q = f_0/B$. Подставив $f_0 = \theta f_s/(2\pi)$, получим приближение в формуле (8): чем ближе полюсы к окружности, тем выше добротность.

Коэффициенты по частоте среза $f_0$ и добротности даёт кулинарная книга Роберта Бристоу-Джонсона, RBJ Audio EQ Cookbook (W3C). Для фильтра нижних частот (low-pass filter, LPF)?фильтр нижних частот (low-pass filter, LPF) — Фильтр, который пропускает низкие частоты и ослабляет высокие. Перед дискретизацией помогает ограничить спектр входного сигнала.:

Для фильтра верхних частот (high-pass filter, HPF)?фильтр верхних частот (high-pass filter, HPF) — Фильтр, который пропускает частоты выше частоты среза и ослабляет частоты ниже неё. В сведении убирает гул и низкочастотные помехи. меняется только числитель: $b_0 = b_2 = (1 + \cos\omega_0)/2$, $b_1 = -(1 + \cos\omega_0)$. Числитель ФНЧ равен $b_0(1 + z^{-1})^2$: два нуля в точке $z = -1$, поэтому на $f_s/2$ АЧХ обнуляется. У ФВЧ оба нуля в точке $z = 1$, на нулевой частоте. Формулы получены билинейным преобразованием аналогового фильтра с поправкой на искажение оси частот, поэтому частота $f_0$ стоит точно там, где задана (W3C, разд. 2). Формулы колокола и полок той же книги нужны эквалайзеру, их разберёт отдельная статья.

На частоте среза АЧХ ФНЧ RBJ равна $Q$: $20\lg 0{,}707 = -3{,}01$ дБ, $20\lg 5 = 13{,}98$ дБ. При $Q = 5$ и 10 ширина по радиусу полюса воспроизводит добротность: $1000/199 \approx 5{,}0$ и $1000/99{,}7 \approx 10{,}0$. При $Q = 0{,}707$ полюсы далеко от окружности, пика нет: АЧХ ровная почти до частоты среза, а за ней спадает примерно на 12 дБ на октаву (octave)?октава (octave) — Музыкальный интервал между нотами одного названия в соседних регистрах. Соответствует отношению частот 2:1. — от 2 до 4 кГц на 12,1 дБ. Однополюсный фильтр спадает на 6 дБ на октаву.

Прямая форма I и транспонированная прямая форма II

Одно и то же разностное уравнение (7) можно вычислять по-разному. Прямая форма I (direct form I, DF-I) считает его как написано и хранит четыре состояния: $x[n-1]$, $x[n-2]$, $y[n-1]$, $y[n-2]$. Задержек вдвое больше необходимого, зато в целочисленной арифметике с дополнительным кодом внутренних переполнений не бывает, пока выход в допустимых пределах (J. O. Smith, Direct-Form I). Кулинарная книга RBJ называет её лучшим и простейшим способом для процессоров с аккумулятором двойной длины (W3C, разд. 2).

Транспонированная прямая форма II (transposed direct form II, TDF-II) обходится двумя состояниями (Sundararajan, §8.4, с. 317):

Транспонирование схемы не меняет передаточной функции (J. O. Smith, Transposed Direct-Forms). Подставив $s_1$ и $s_2$ в первую строку, можно убедиться, что $y[n]$ совпадает с формулой (7). Для точности у TDF-II есть преимущество: нули в ней фактически стоят раньше полюсов, и большое усиление, которое дают одни полюсы, нули частично компенсируют ещё до них (J. O. Smith, Numerical Robustness of TDF-II). По этой схеме работает juce::dsp::IIR::Filter (juce_IIRFilter.h).

Смена коэффициентов на ходу

При медленной модуляции три структуры дают одинаковый выход: в каждый момент фильтр почти неподвижен, и плавная автоматизация не вредит ни одной из них. При модуляции 500 Гц выход прямых форм в 4–5 раз выше, чем у фильтра переменных состояний, а при 2000 Гц и $Q = 20$ он уходит в бесконечность, хотя каждый «замороженный» фильтр устойчив. Для быстрой модуляции нужна другая структура.

Фильтр переменных состояний

Без модуляции ФНЧ переменных состояний с теми же частотой среза и $Q$ имеет ту же передаточную функцию, что ФНЧ RBJ: на белом шуме их выходы в коде из раздела о JUCE расходятся не больше чем на $9 \cdot 10^{-8}$. Разница появляется только при смене параметров. У полюсов сохраняется смысл формулы (8), поэтому параметр resonance в JUCE — это и есть $Q$: на частоте среза АЧХ равна $1/(2R)$, а JUCE вычисляет $2R$ как 1 / resonance (Zavalishin, §4.1, с. 97; juce_StateVariableTPTFilter.h).

Код на C++

Класс Biquad

Класс считает формулу (10) и хранит пять коэффициентов отдельно от двух состояний. Функция makeRbjFilter вычисляет коэффициенты ФНЧ и ФВЧ по формуле (9), а метод response — частотную характеристику по формуле (7) при $z = e^{j\omega}$. Нужна только стандартная библиотека C++23.

// Биквад в транспонированной прямой форме II и коэффициенты ФНЧ и ФВЧ по RBJ Audio EQ Cookbook.
#include <cmath>
#include <complex>
#include <numbers>

struct BiquadCoefficients
{
    double b0 = 1.0; // числитель
    double b1 = 0.0;
    double b2 = 0.0;
    double a1 = 0.0; // знаменатель, a0 = 1
    double a2 = 0.0;
};

// ФНЧ (lowPass = true) или ФВЧ второго порядка; все коэффициенты делятся на a0.
inline BiquadCoefficients makeRbjFilter (double sampleRate, double cutoffHz, double q, bool lowPass)
{
    const auto w0 = 2.0 * std::numbers::pi * cutoffHz / sampleRate;
    const auto cosW0 = std::cos (w0);
    const auto alpha = std::sin (w0) / (2.0 * q);
    const auto a0 = 1.0 + alpha;
    const auto edge = (lowPass ? 1.0 - cosW0 : 1.0 + cosW0) / 2.0; // b0 = b2

    return { .b0 = edge / a0,
             .b1 = (lowPass ? 2.0 : -2.0) * edge / a0,
             .b2 = edge / a0,
             .a1 = -2.0 * cosW0 / a0,
             .a2 = (1.0 - alpha) / a0 };
}

class Biquad
{
public:
    void setCoefficients (const BiquadCoefficients& newCoefficients) noexcept { c = newCoefficients; }
    void reset() noexcept { s1 = s2 = 0.0; }

    // y[n] = b0·x[n] + s1; s1 и s2 хранят ещё не сложенные части будущих выходов.
    float processSample (float input) noexcept
    {
        const auto x = static_cast<double> (input);
        const auto y = c.b0 * x + s1;
        s1 = c.b1 * x - c.a1 * y + s2;
        s2 = c.b2 * x - c.a2 * y;
        return static_cast<float> (y);
    }

    // H(e^jω) = (b0 + b1·e^−jω + b2·e^−2jω) / (1 + a1·e^−jω + a2·e^−2jω)
    std::complex<double> response (double frequencyHz, double sampleRate) const
    {
        const auto z1 = std::polar (1.0, -2.0 * std::numbers::pi * frequencyHz / sampleRate);
        const auto z2 = z1 * z1;
        return (c.b0 + c.b1 * z1 + c.b2 * z2) / (1.0 + c.a1 * z1 + c.a2 * z2);
    }

private:
    BiquadCoefficients c;   // параметры фильтра
    double s1 = 0.0;        // состояние, меняется каждый отсчёт
    double s2 = 0.0;
};

Код собран clang++ -std=c++23 -Wall -Wextra -pedantic -Wconversion без предупреждений. Для ФНЧ 1 кГц с $Q = 0{,}7071$ при 48 кГц response даёт 0,00 дБ на 100 Гц, −3,01 дБ на 1 кГц, −12,37 дБ на 2 кГц и −24,48 дБ на 4 кГц, а первые отсчёты импульсной характеристики — 0,003916; 0,014941; 0,027785. ФВЧ с теми же параметрами на 100 Гц даёт −40,03 дБ, а при $Q = 5$ ФНЧ поднимает частоту среза на 13,98 дБ, как в таблице 2.

Фильтры JUCE

В модуле juce_dsp есть КИХ- и БИХ-фильтры и фильтр переменных состояний. Класс ниже готовит биквад для двух каналов, КИХ-фильтр скользящего среднего и фильтр переменных состояний; процессор вызывает prepare из prepareToPlay и process из processBlock.

// Три фильтра JUCE в одном классе: биквад для двух каналов, КИХ и фильтр переменных состояний.
#include <juce_dsp/juce_dsp.h>

class FilterSection
{
public:
    void prepare (double sampleRate, int maxBlockSize)
    {
        const juce::dsp::ProcessSpec stereo { sampleRate, static_cast<juce::uint32> (maxBlockSize), 2 };
        setLowPass (sampleRate, 1000.0f, 0.707f);
        lowPass.prepare (stereo);

        // Скользящее среднее четырёх отсчётов; FIR::Filter обрабатывает только один канал.
        const std::array<float, 4> taps { 0.25f, 0.25f, 0.25f, 0.25f };
        average.coefficients = new juce::dsp::FIR::Coefficients<float> (taps.data(), taps.size());
        average.prepare ({ sampleRate, static_cast<juce::uint32> (maxBlockSize), 1 });

        sweep.setType (juce::dsp::StateVariableTPTFilterType::lowpass);
        sweep.setCutoffFrequency (1000.0f);
        sweep.setResonance (0.707f); // в JUCE resonance — это Q
        sweep.prepare (stereo);
    }

    // Новые коэффициенты без выделения памяти: массив копируется в готовый объект Coefficients.
    void setLowPass (double sampleRate, float cutoffHz, float q)
    {
        *lowPass.state = juce::dsp::IIR::ArrayCoefficients<float>::makeLowPass (sampleRate, cutoffHz, q);
    }

    void process (juce::AudioBuffer<float>& buffer)
    {
        juce::dsp::AudioBlock<float> block { buffer };
        lowPass.process (juce::dsp::ProcessContextReplacing<float> { block });
    }

    juce::dsp::ProcessorDuplicator<juce::dsp::IIR::Filter<float>, juce::dsp::IIR::Coefficients<float>> lowPass;
    juce::dsp::FIR::Filter<float> average;
    juce::dsp::StateVariableTPTFilter<float> sweep;
};

Код собран с JUCE 9.0.3 в консольном приложении с флагами предупреждений JUCE и -Wconversion без предупреждений. На 1 с белого шума выход класса Biquad расходится с IIR::Filter не больше чем на $1{,}6 \cdot 10^{-6}$ — JUCE хранит состояния в float, а Biquad в double, — а с ФНЧ StateVariableTPTFilter — не больше чем на $9 \cdot 10^{-8}$. АЧХ коэффициентов lowPass на 1 кГц, по getMagnitudeForFrequency, равна −3,01 дБ.

Самопроверка

Источники

• Sundararajan D. Digital Signal Processing: An Introduction. §2.2–2.5, с. 46–57; §2.7, с. 59–60; §6.1, с. 196–197; §6.3.2, с. 217–218; §7.1–7.2, с. 230–237; §8.4, с. 315–318. в списке литературы → — свёртка как сумма откликов на импульсы и пример 4 ∗ 4, длина N + M − 1, линейность и инвариантность во времени, отклик на комплексную экспоненту, устойчивость через сумму модулей, КИХ и БИХ, z-преобразование, передаточная функция и устойчивость по полюсам, КИХ-фильтры, линейная фаза и групповая задержка, каноническая и транспонированная формы второго порядка, разложение на секции.
• Bellanger M. Digital Signal Processing: Theory and Practice. §6.2, с. 127–131. в списке литературы → — чисто рекурсивная секция второго порядка: коэффициенты b₁ = −2r·cos θ и b₂ = r² (в обозначениях книги), АЧХ через расстояния до полюсов, подъём на резонансе обратно пропорционален расстоянию до окружности, ширина (1 − r)/π, групповая задержка и её отрицательные значения.
• Smith J. O. III. Introduction to Digital Filters with Audio Applications. Phase and Group Delay, Group Delay, BiQuad Section, Resonator Bandwidth in Terms of Pole Radius, Direct-Form I, Transposed Direct-Forms, Numerical Robustness of TDF-II. в списке литературы → — фазовая и групповая задержка, задержка огибающей узкополосной группы, биквад и полюсы через радиус и угол, ширина резонанса −ln r/(πT), прямая форма I без внутренних переполнений, транспонирование без изменения передаточной функции, нули перед полюсами в TDF-II.
• Bristow-Johnson R. Audio EQ Cookbook / ed. by R. Toy. W3C Working Group Note, 2021 — передаточная функция биквада, прямая форма I, коэффициенты ФНЧ и ФВЧ, учёт искажения частот билинейного преобразования.
• Zavalishin V. The Art of VA Filter Design. §3.11, с. 79–81; §4.1, с. 95–97. в списке литературы → — прямые формы и их поведение при модуляции, метод TPT, точность прямых форм высокого порядка, фильтр переменных состояний с выходами ФВЧ, полосового и ФНЧ, АЧХ 1/(2R) на частоте среза.
• JUCE 9.0.3. juce_IIRFilter.h, juce_IIRFilter.cpp, juce_IIRFilter_Impl.h — TDF-II, ArrayCoefficients и Coefficients::makeLowPass, нормировка на a₀, резерв на восемь коэффициентов, совет для быстрой модуляции.
• JUCE 9.0.3. juce_FIRFilter.h — КИХ-фильтр для одного канала, граница 128 коэффициентов и класс Convolution.
• JUCE 9.0.3. juce_StateVariableTPTFilter.h и juce_StateVariableTPTFilter.cpp — шаг фильтра переменных состояний, resonance и сглаживание частоты среза.
• FabFilter. Pro-Q 4 Help: Processing mode — задержка и предзвон линейно-фазового режима.