В задачах обработки сигналов часто возникает необходимость фильтрации сигналов, когда сигнал разбивается на узкополосные диапазоны. В бытовом плане мы с этим сталкиваемся при воспроизведении музыки через акустические системы, в которых каждый громкоговоритель (динамик) воспроизводит свою полосу частот, которых обычно три — низкие (НЧ), средние (СЧ) и высокие (ВЧ); для воспроизведения сверхнизких частот иногда выделяют отдельную акустическую систему под названием «сабвуфер». Конкретные границы частот зависят от реализации и ориентировочно находятся на границах 100 Гц, 1 кГц и 5 кГц. Для того, чтобы не было резких скачков громкости между динамиками, используют частичное перекрытие — когда амплитуда воспроизводимой полосы частот плавно спадает на одном, одновременно нарастая на другом.

Наиболее популярными фильтрами для такого разбиения являются фильтры Линквитца-Рейли 4-го порядка, представляющих из себя два последовательно соединённых фильтра Баттерворта, изображение АЧХ которых многим хорошо знакомо:

Популярность их обусловлена простотой схемотехнической реализации и гладкой АЧХ с отсутствием пульсаций. Порядок фильтра, обуславливающий скорость затухания, обычно выбирают 4-ым.

У всех аналоговых фильтров есть недостатки — нестабильность параметров и сложность в настройке. Другой немаловажным недостатком является неустранимый сдвиг фаз. Реализовать фазолинейный фильтр в аналоге невозможно, потому что это нарушает причинно-следственную связь выходного сигнала от входного.

Появление цифровой техники в целом и цифровых источников сигналов в частности позволило избавиться от недостатков аналоговых фильтров путём реализации их непосредственно программным образом, в цифре. Такой подход требует «мультиампинг» — когда используется отдельный усилитель мощности для каждой полосы частот — в отличие от классического подхода, когда сначала усиливается широкополосный сигнал, который разбивается на частотные полосы (обычно) пассивными фильтрами. Применение нескольких усилителей вместо одного очевидным образом удорожает технику, поэтому такое решение выбирают для особо качественных систем.

▍ Краткое введение в цифровые фильтры

Цифровой фильтр — это функция, которую можно рассматривать как во временной, так и в частотной области. Во временной области она определяет импульсную характеристику, получаемой при реакции некоторого устройства или его мат. модели на единичный импульс (на слух воспринимаемый как щелчок). В частотной области она определяет затухание или усиление амплитуды и сдвиг фазы на отдельно взятой частоте. Переход между этими областями (в англоязычной литературе используется слово domain) осуществляется через преобразование Фурье, которое между функциями во временной и в частотной области задаёт однозначное соответствие.

Разделяют два типа цифровых фильтров:

IIR (Infinite Impulse Response) — фильтры с бесконечной импульсной характеристикой. По сути, представляют собой всё те же аналоговые фильтры, но «моделируемые» в цифре — математический аппарат в обоих случаях идентичен (в основе которого лежит преобразование Лапласа). Отсюда следует бесконечность импульсной характеристики, которую можно представить в виде суммы экспоненциально затухающих синусоид.

Вот так выглядит (обрезанная справа) импульсная характеристика рекурсивного фильтра низких частот с высокой добротностью



Их преимуществом является низкая вычислительная сложность — каждое новое значение вычисляется рекурсивно в зависимости от предыдущих. Они также чувствительны к погрешностям вычисления и ошибкам округления — поэтому вопрос устойчивости (гарантированного затухания при отсутствии сигнала на входе) таких фильтров в теории рассматривается отдельно.

Наиболее простым является фильтр низких частот первого порядка:

FIR (Finite Impulse Response) — фильтры с конечной импульсной характеристикой. Самый известный фильтр такого рода — это скользящее среднее, а сама фильтрация осуществляется посредством линейной свёртки отсчётов фильтра с отсчётами входного сигнала. Здесь сложность уже квадратичная — однако её можно уменьшить до , если использовать алгоритм быстрой свёртки с использованием быстрого преобразования Фурье (FFT). Конечность количества отсчётов накладывает свои ограничения на форму АЧХ, приводя к пульсациям.


А так выглядит импульсная характеристика фазолинейного FIR-фильтра низких частот



Довольно популярной практикой при проектировании фильтров является линеаризация — когда берётся АЧХ известного IIR фильтра и формируется из неё FIR с линейной фазой.

▍ Постановка задачи

Итак, нашей задачей является сформировать пару фильтров, делящих полосу частот на две с перекрытием. При этом они должны обладать следующими характеристиками:

1. Гладкая АЧХ без пульсаций;
2. Фазолинейность (сдвиг фаз на всех частотах равен нулю);
3. Симметричность АЧХ в логарифмическом масштабе частот.

Cимметричность (зеркальность) АЧХ решает две задачи:
1) ВЧ-составляющую сигнала можно получить вычитанием из исходного НЧ-составляющей — как в частотной, так и во временной области;
2) избавляет от мучительного выбора, для какого фильтра — НЧ или ВЧ — предпочтительнее более пологий или крутой спад АЧХ.

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

▍ Общий алгоритм построения

Существует несколько подходов к проектированию КИХ-фильтров, из которых наиболее простым и интуитивно-понятным является следующий:

1. В частотном домене задаётся желаемая АЧХ фильтра;
2. Производится обратное преобразование Фурье:
3. Накладывается оконная функция.

Его можно выполнить как чисто аналитически, так и численно, используя дискретное преобразование Фурье (далее FFT). Аналитическое решение сложнее и имеет очевидное ограничение на использование только тех функций, для которых известны Фурье-образы.

▍ Аналитическое решение

Озвученным выше требованиям гладкости и симметрии (но не фазолинейности) соответствует фильтр Линквитца-Рейли. Зная формулу АЧХ , где — порядок фильтра, формулу импульсной характеристики можно получить аналитически через интеграл Фурье, а конкретные отсчёты FIR фильтра считать непосредственным её вычислением.

График фильтра в логарифмическом масштабе:



Лучший способ для вычисления интеграла Фурье — это использовать какую-нибудь систему компьютерной алгебры. Например, в Wolfram Mathematica это будет выглядеть так:

InverseFourierTransform[1/(1 + w^2), w, x] // FullSimplify
↓


График получившейся функции, она же импульсная характеристика:



Как видно из формулы и графика, импульсная характеристика бесконечна — стремится к нулю, но не достигает его; а также симметрична относительно нуля (за счёт фазолинейности). Ограничение её во времени осуществляется путём умножения на оконную функцию, за счёт чего значения функции за пределами окна приобретают нулевые значения, не требующих участия в расчётах. Побочным эффектом этого становится искажения исходного спектра фильтра — поскольку при умножении функций их спектры сворачиваются. Итоговый спектр можно увидеть также через преобразование Фурье. Например, для прямоугольного окна получим


в результате чего АЧХ фильтра изменится на



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

Помножив её на импульсную характеристику фильтра, получим:


а АЧХ примет вид



Как видно, величина пульсаций значительно уменьшилась. По графику видно, что «полка» фильтра слегка опустилась — в идеале это тоже нужно учитывать и компенсировать. В данном случае она составила (на нулевой частоте):

FourierTransform[E^-Abs[x] Sqrt[Pi/2] NuttallWindow[x/20], x, w] /.w->0.0
↓


Мнимая часть здесь получилась вследствие погрешности при численных вычислениях и её можно смело отбрасывать (об этом говорит значение , поскольку точность вычислений в формате double и представляет примерно 16 десятичных цифр).

Аналогичным образом можно посчитать импульсные характеристики и для более высоких (но только целых) порядков, например:

InverseFourierTransform[1/(1+w^4), w, x] // FullSimplify
↓


Как видно, автоматического упрощения уже недостаточно и здесь требуется ручная работа по приведению формулы к удобочитаемому виду. В итоге получается следующая формула (в комплексных числах) для импульсной характеристики фильтра Линквитца-Рейли произвольного порядка:

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

▍ Переход в дискретную область

Синтез фильтра в дискретном виде осуществляется сходным образом. Принципиальное отличие состоит в том, что:

1. Частотный диапазон фильтра ограничен частотой дискретизации по определению;
2. На выходе мы получаем периодический сигнал — в отличие от аналитического решения, в которой функция во времени затухает к бесконечности. Этот момент особенно важен, когда проектируются фазосдвигающие фильтры. Из этого также следует, что даже если мы не будем явно накладывать оконную функцию, по факту у нас в любом случае будет наложено прямоугольное окно.

Попробуем реализовать тот же фильтр НЧ. Первым делом необходимо определиться с частотой дискретизации и размером массива. Размер массива ограничивает нижнюю частоту фильтра (но не постоянную составляющую). Если не ставить задачей минимизацию размера, то 2048 отсчётов обычно достаточно. Размер выбирается кратным степени двойки исходя из ограничений стандартной реализации БПФ.


Также имеются различия (но не результат) при работе с классическим (комплексным) FFT, Real-FFT, который работает только с половиной спектра, предполагая его симметрию, и FHT (быстрое преобразование Хартли), в котором изначально и данные, и спектр определены в поле действительных чисел. Здесь будет использоваться классический FFT.

Итак, определим формулу для фильтра — для примера 4-го порядка:

FilterAmps[w_] := 1/(1+w^8);

Выберем частоту дискретизации (в герцах), стандартную для современных ЦАП:

samplerate = 48000;

Выберем частоту среза, на которой подавления будет составлять 0.5 (-6 дБ) и соответствовать нормированной частоте w=1:

fc = 1000;

Выберем размер массива (маленький, для наглядности):

size = 128;

Далее необходимо заполнить массив. Каждому элементу в нём будет соответствовать частота от 0 (постоянная составляющая) до (частота Найквиста), равномерно распределённых в линейном масштабе. Определим функцию, которая в зависимости от индекса будет считать нормированную частоту:

w[index_] := (index*samplerate)/(size*fc)

Заполняем первую половину (положительную часть спектра) массива:

halfspectrumdata =
Table[FilterAmps[w[index]]*(-1)^index, {index, 0, size/2}];

Здесь фазу на каждой нечётной частоте мы повернули на 180° для того чтобы после БПФ максимум импульса был расположен по центру массива.

Вторую половину (отрицательную часть спектра) получаем реверсом первой, за исключением крайних частот (нулевой — постоянной составляющей и максимальной — частоты Найквиста, потому что для их описания достаточно одного числа).

spectrumdata =
Join[halfspectrumdata, Reverse[halfspectrumdata[[2;;size/2]]]];

Получили следующее:


Теперь делаем обратное преобразование Фурье:

wavedata = InverseFourier[spectrumdata];


И накладываем оконную функцию — Нутталла, как и в прошлый раз:


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

Fourier[wavedata][[1]]
↓


и затем просто поделить на него значение каждого отсчёта в импульсе

wavedata /= Fourier[wavedata][[1]];

Сделав FFT ещё раз, можно убедиться, что постоянная составляющая стала равной единице:

Fourier[wavedata][[1]]
↓


Если бы мы проектировали фильтр высоких частот, но нормализацию нужно было бы делать по амплитуде не на нулевой частоте, а на частоте Найквиста.

▍ Кусочно-непрерывные фильтры

Исходя из изначально поставленной задачи, описывать АЧХ удобнее кусочно-непрерывной функцией, в которой полосы пропускания и подавления константны и равны 1 и 0 соответственно, заданной в логарифмическом масштабе. Выводу таких функций была посвящена отдельная статья, здесь для примера мы возьмём стандартную smooth-step функцию, построенной из кубического полинома 3-го порядка:

Из её графика в линейном масштабе симметрия, обеспечивающая суммирование в единицу, хорошо видна:


А вот в логарифмическом масштабе — уже нет, зато хорошо видно главное отличие от фильтра Линквитца-Рейли — спад не пологий, а имеет выраженную границу справа, за которой громкость (в децибелах) равна минус бесконечности. Дополнение до единицы даёт ВЧ-фильтр (жёлтым цветом), который и обладает желаемой симметрией:


На большем диапазоне громкости это свойство фильтров выглядит ещё более наглядным:


В отличие от классических фильтров, где параметризация крутизны спада АЧХ через порядок обусловлена их схемотехнической реализацией, здесь мы можем оперировать непосредственно шириной полосы сопряжения, задавая её в октавах от частоты раздела в -6 дБ. Для этого потребуется дополнительная функция для перевода логарифмического масштаба в линейный, значение которой будет передаваться в smooth-step функцию для определения амплитуды на частоте . Её можно определить несколькими путями:

1) через граничные частоты

2) через центральную частоту и ширину (в октавах)

3) через центральную и верхнюю граничную частоту

4) через центральную и нижнюю граничную частоту

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

SymmetricFilterImpulseResponse[f0_, width_, size_, samplerate_,
ClipFunction_, WindowFunction_] :=
Table[ClipFunction[(2 Log[(index samplerate)/(f0 size)])/
(width Log[2]) ] (-1)^index, {index, 0, size/2}] //
Join[#, Reverse[#[[2 ;; size/2]]]] & // InverseFourier[#]
Table[WindowFunction[(index-1)/size-1/2], {index, 1, size}] & // #/
Fourier[#, FourierParameters -> {1, -1}][[1]] & // Re


Передав в эту функцию наш кубический фильтр и окно Нуттала получим

SymmetricFilterImpulseResponse[1000, 0.5, 512, 48000, FilterCubicClip, NuttallWindow]
↓

▍ Несколько интересных фильтр-функций

Помимо уже рассмотренной кубической, можно придумать множество и других функций для фильтров, основанных на какой-нибудь математической идее. Также имеет значение, насколько легко может быть вычислена (и может ли быть вычислена вообще) обратная функция.

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

Кубическая. Обратная к ней функция уже немного контринтуитивна (а если вам интересно, откуда здесь взялись тригонометрические функции — тогда сюда):

Используя полином 5-ой степени, с двумя нулевыми производными на границах сопряжения. Обратная функция здесь уже в элементарных функциях не выражается (потому и не приведена):

Используя полином 13-ой степени — с дополнительным «выпрямлением»:

Используя рациональный полином — даёт более гладкую характеристику и более простую обратную функцию, чем у кубической:

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

Линейно спадающая в логарифмическом масштабе:

Более плавный вариант предыдущего с возможностью регулировки «жёсткости». Здесь уже обратная функция для произвольного n в элементарных функциях не выражается:

С бесконечным количеством нулевых производных в точках стыковки, что обеспечивает идеальное сопряжение. А здесь обратная функция легко находится:

С максимально быстро затухающей импульсной характеристикой — что при достаточно больших () позволяет обойтись без оконной функции вообще. Сама формула получена модуляцией аргумента кубической функции, т.е. :

Здесь — это специальная функция (функция ошибок), определяемая как интеграл от гауссианы. Вместо неё можно использовать функцию гиперболического тангенса — импульсная характеристика также будет достаточно быстро затухать, хоть и не так быстро, как в предыдущем случае. Причём при больших n визуально она будет походить на Линквитца-Рейли и это совсем не случайно —

▍ Заключение

Несмотря на очевидность постановки задачи, мне не встречалось ничего подобного ни в специализированных решениях, ни в мат. пакетах общего назначения. Возможно, что это связано с инерцией мышления в стиле «всё, что можно придумать, уже придумано». Однако это не так хотя бы потому, что сейчас доступны технологии и вычислительные мощности, о которых всего лишь 50 лет назад инженеры и математики не могли и мечтать.

В качестве proof of concept фильтры по такой методике были реализованы в плагине для популярного в аудиофильских кругах плеера foobar, извлекающего канал сабвуфера из стереосигнала (официальный репозиторий); и ещё в одном, просто для демонстрации звучания (неофициальный репозиторий). Для реализации таких фильтров в «мультиампинг»-системах можно использовать готовые решения, позволяющие загружать заранее посчитанные импульсные файлы. При отсутствии мат.пакетов их можно посчитать даже в excel-е (используя преобразование Фурье из пакета анализа).

Насколько такие фильтры «звучат» лучше классических решений — вопрос уже аудиофильский, но, как минимум, разницу в звучании обеспечивают вполне объективную и обеспечивают возможность для более точной настройки.