Кривые Безье для обработки звука

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

Оглавление

• Предыдущая статья

• Введение

• Графический эквалайзер

• Шкала частот эквалайзера

• Построение графика АЧХ

• Кривые Безье

• Общее описание алгоритма

• Создание контура

• Построение кривой Безье

• Вычисление импульсной характеристики

• Описание командной строки

• Названия выходных файлов

• Примеры запуска программы

• Заключение

• Ссылки и литература

Введение

Данный материал является продолжением предыдущей статьи. Напомню вкратце основные выводы из неё:

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

• фильтрацию звука можно осуществить с помощью процедуры линейной свёртки отсчётов входного сигнала с конечной импульсной характеристикой КИХ-фильтра

• для расчёта коэффициентов КИХ-фильтра можно использовать как готовые аналитические формулы, так и метод получения их из необходимой амплитудно-частотной характеристики (АЧХ, в общем случае – произвольной) с помощью обратного преобразования Фурье

• полученная импульсная характеристика даёт фильтр, реальная АЧХ которого отличается от желаемой АЧХ (иногда весьма значительно), но это отличие можно минимизировать с помощью различных методов (например, применением т.н. оконной функции)

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

Графический эквалайзер

В 1990-х годах, когда звуковые карты были роскошью, а запись собственного Audio CD была сродни некоторому чуду (на студенческую стипендию можно было купить всего несколько чистых болванок CD-R, а сам записывающий привод мог стоить более 1000 у.е.), я вместе с друзьями занимался записью песен на самодельный 4-дорожечный магнитофон. Тогда компьютер с процессором Pentium был далеко не на каждом столе, и о многодорожечной цифровой звукозаписи сложно было даже мечтать. Как устройства записи, так и устройства обработки звука были в основном аналоговыми. Но уже тогда можно было оцифровать звук в компьютер с помощью звуковой карты, и появлялись первые звуковые редакторы с возможностью осуществления различных видов цифровой обработки звука (т.н. эффектов). Тогда же я начал заниматься изучением этой темы и реализацией собственных программных модулей обработки.

Одной из таких программ была Cool Edit Pro (теперь она называется Adobe Audition), и в ней присутствует эффект под названием Graphic Equalizer. Он реализует как раз обработку звука эквалайзером с помощью КИХ-фильтра (фильтр с конечной импульсной характеристикой). Как уже упоминалось в предыдущей статье, такой тип фильтра выбран не просто так – он не искажает фазу обрабатываемого звука, поэтому при изменении громкости в разных полосах частот практически не возникает неприятных на слух искажений (при условии симметричности импульсной характеристики, что легко реализовать при расчёте фильтра).

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

Реализованный в программе Cool Edit Pro эквалайзер предоставляет похожий интерфейс:

Основные элементы:

• вертикальные регуляторы громкости в каждой из частотных полос

• кнопки выбора количества полос (10, 20 или 30)

• поле Gain для задания общей громкости

Отдельно отметим дополнительные элементы: поле “Accuracy“ (points) и кнопку “Graph actual response”. Здесь задаётся длина фильтра (его импульсной характеристики) в отсчётах (точках), а нажатием кнопки можно отобразить график реальной АЧХ полученного фильтра. Вспомним ещё один важный вывод, сделанный в предыдущей статье:

Чем короче КИХ-фильтр, тем хуже он выполняет свою «работу» в низких частотах. То есть, реальная АЧХ фильтра в низких частотах тем лучше соответствует желаемой АЧХ, чем длиннее фильтр.

После того, как с помощью регуляторов громкости задана нужная АЧХ, необходимо обязательно нажимать кнопку “Graph actual response” и смотреть, сильно ли отличается график реальной АЧХ (тонкая голубая линия) от желаемой (синий закрашенный фон). Попробуем проверить утверждение о длине фильтра с реализацией ВЧ-фильтра, подавляющего низкие частоты. Нажмём кнопку “Reset All to Zero” (сбросить всё в ноль), увеличим диапазон (Range) регуляторов в полосах до 100 дБ и попробуем настроить подавление звука в частотах ниже 100 Гц. После нажатия кнопки “Graph actual response” видим, что желаемая АЧХ заметно отличается от реальной – фильтр довольно плохо подавляет низкие частоты:

Попробуем увеличить длину фильтра до 2000 точек:

Уже лучше, но только если задать длину фильтра 5000 точек и более, его реальная АЧХ станет близкой к желаемой. Даже при максимально допустимой длине 16384 (больше задать нельзя) результат не идеален, особенно если в полосе непропускания хочется максимально подавить сигнал (здесь подавление увеличено до 90 дБ):

Вспомним, что обработка звука линейной свёрткой требует большого количества вычислений, поэтому увеличивать длину фильтра нужно с осторожностью. На современных процессорах с тактовой частотой ядер в несколько ГГц обработку таким фильтром вполне можно выполнять в реальном времени (во время воспроизведения), но на процессорах Pentium даже длина фильтра в 1000 точек требовала весьма долгой обработки звукового файла и была возможной только в офлайне (запустили – подождали – сохранили результат в файл). Впрочем, пока что оставим этот момент в стороне.

Шкала частот эквалайзера

Для простоты переключим эквалайзер в 10-полосный режим и обратим внимание на частоты полос, в которых можно регулировать громкость:

Как видим, в графическом интерфейсе они расположены на одинаковом расстоянии друг от друга по горизонтали, но их частоты отличаются в 2 раза относительно соседней полосы (Гц): 31, 62, 125, 250, 500, 1000, 2000 и т.д. Это связано с тем, что восприятие человеком высоты звука обладает логарифмичностью: каждое повышение тона на одинаково воспринимаемую высоту – октаву означает увеличение частоты звука в 2 раза. Ранее упоминалось, что и единицы изменения громкости – децибелы (дБ) также введены в связи с логарифмичностью восприятия громкости звука на слух. Особенности восприятия высоты звука были описаны ещё Пифагором, и желающие могут прочитать об этом в книгах (названия некоторых приведены в конце статьи).

Для нашей задачи это означает следующее:

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

Как мы помним из предыдущей статьи, можно задать график нужной АЧХ в частотной области и получить из него импульсную характеристику КИХ-фильтра с помощью обратного дискретного преобразования Фурье (БПФ). Допустим, нам нужно получить изменение громкости в соответствии с рисунком выше. Изобразим его в виде таблицы по полосам частот (верхняя строка – частота в Гц, нижняя – громкость в дБ):

Значит, нужно в частотной области «нарисовать» примерно такой график:

Но обратное преобразование Фурье производится над массивом, который содержит линейно изменяющиеся значения в частотной шкале. Поэтому нужно соответствующим образом отобразить значения из логарифмической шкалы по горизонтали в линейную. Также необходимо значения АЧХ из децибел перевести в абсолютные значения – впрочем, это делается просто:

double dB2Gain(double dB)
{
	return pow(10.0, dB / 20.0);
}

Эксперименты, проделанные выше с ВЧ-фильтром, наглядно демонстрируют утверждение, которое можно сформулировать так:

Фильтр с АЧХ, в которой нулевая ширина переходов между частотными полосами (т.н. идеальный фильтр), возможно реализовать только теоретически – для этого нужна бесконечная длина КИХ-фильтра. На практике реализовать такой фильтр невозможно.

Требуется соблюсти некоторый компромисс – например, можно сделать переходы АЧХ между полосами более плавными, без резких участков с прямыми углами. Неизвестно, как именно это реализовано разработчиками программы Cool Edit Pro, но можно предположить, что они постарались делать желаемый график АЧХ (Requested frequency response) достаточно плавным, чтобы он как можно меньше отличался от реальной АЧХ полученного фильтра.

Построение графика АЧХ

Когда я занимался задачами обработки звука более 25 лет назад и захотел реализовать аналогичный КИХ-эквалайзер, передо мной также встал вопрос:

Каким образом задать график желаемой АЧХ, чтобы в результате обратного дискретного преобразования Фурье получился фильтр, реальная АЧХ которого мало отличается от желаемой?

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

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

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

• в верхних и нижних точках полос проходить строго через точку с нужной громкостью (то есть, не должно быть никаких лишних «выбросов» за пределы нужных значений – например, такие возможны, если применять некоторые виды сплайн-интерполяции между опорными точками)

Пытаясь нарисовать различные графики на бумаге, я постепенно понял, что мне подобные кривые напоминают. В середине 1990-х годов мне приходилось выполнять некоторую работу в редакторе Corel Draw, в котором есть удобный инструмент – кривые Безье (Bézier curves – сейчас он есть уже во многих графических редакторах). С их помощью можно соединять множество опорных точек непрерывной кривой, и при помощи дополнительных направляющих точек (т.н. guide points) можно обеспечить её гладкость, изменяя форму и кривизну участков между точками:

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

Кривые Безье

Про теорию кривых Безье в сети есть множество материалов, а здесь можно «покрутить» эти кривые самостоятельно и увидеть, как выглядят кривые разного порядка и как влияют на их форму направляющие точки:

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

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

Если каждую пару опорных точек соединить кривой Безье с соответствующими направляющими точками, должна получиться примерно такая кривая:

Выглядит вполне неплохо, с этим уже можно работать. Но в каких единицах нужно задавать координаты точек, чтобы получить окончательный график АЧХ в линейной шкале частот? Над этим мне пришлось поломать голову, пока не пришла следующая идея:

1. Задается базовая частота первой опорной точки с самой низкой частотой первой полосы.

2. Ширина полос задается в октавах.

3. Частоты остальных полос рассчитываются из базовой частоты и интервала в октавах от неё.

Например, чтобы получить 10-полосный эквалайзер с частотами полос как в примерах выше, нужно взять 31.25 Гц в качестве значения базовой частоты, а ширину каждой полосы установить в одну октаву. Для 20- и 30-полосного эквалайзера ширина полос будет 1/2 и 1/3 октавы.

Соответственно, для опорных и направляющих точек кривой Безье, задающей график АЧХ, координата по горизонтали – смещение в октавах относительно базовой частоты (её координата принимается за 0). Координаты по вертикали задаются в децибелах – их нужно будет пересчитать в абсолютные значения коэффициентов АЧХ (формула была показана выше в процедуре dB2Gain). После того, как построен полный график, координаты всех точек по горизонтали нужно будет пересчитать в значения частоты относительно частоты дискретизации:

double FreqOct(double freq, double oct)
{
	return freq * exp(log(2.0) * oct);
}

const int nSR = 48000;
const auto bfreq = 31.25 / double(nSR);
double foct = 1.0;
const auto index = uint32_t(FreqOct(bfreq, foct) * double(fftsize) + 0.5);

Здесь также показано, как по значению смещения в октавах foct получить целочисленный индекс элемента в буфере, в который нужно поместить значение коэффициента АЧХ. Затем над этим буфером производится дискретное обратное преобразование Фурье длиной fftsize точек (напомню, что для вычисления быстрого преобразования Фурье его длина должна быть степенью двойки). Из результата можно получить значения импульсной характеристики КИХ-фильтра.

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

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

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

Почти 26 лет назад! Посмотрев текст кода, я решил, что выкладывать его «как есть» мне не очень хочется (думаю, многие поймут, почему ;-), и переписал его в некоторых местах. В те давние времена я не знал о существовании стандартной библиотеки, поэтому сделал свою реализацию двунаправленного списка. Наверное, можно было и это место переписать с использованием контейнеров STL, но я оставил, как есть.

Основная функция класса:

void CBezier::CalcResp(const uint32_t fftsize, float* data, const uint32_t nbands, double* dB, const double bfreq, const double step)
{
	CreateContour(nbands, dB, step);
	const double minint = step / 24.0;
	constexpr double mindb = 0.1;
	m_minlen = minint * minint + mindb * mindb;
	PerformBezier();
	CreateResp(fftsize, data, bfreq);
}

• fftsize – длина преобразования Фурье. Должна быть числом, которое является степенью двойки, и минимум в 2 раза большим, чем желаемая длина импульсной характеристики КИХ-фильтра (которая, в свою очередь, должна быть нечётным числом для симметричности характеристики)

• data – указатель на буфер, в котором осуществляется расчёт обратного дискретного преобразования Фурье. Его объем должен быть в 2 раза больше, чем fftsize, т.к. обратное преобразование осуществляется над комплексными числами. Перед началом расчёта в буфере должны быть нулевые значения

• nbands – количество полос эквалайзера (например, 10)

• dB – указатель на массив значений громкости в частотных полосах в децибелах (их количество должно равняться nbands)

• bfreq – базовая частота самой низкочастотной полосы (не в Гц, а значение относительно частоты дискретизации звука)

• step – ширина каждой полосы в октавах

По коду процедуры легко понять общую схему алгоритма:

1. Сгенерировать контур кривой Безье, состоящий из опорных и направляющих точек (функция CreateContour).

2. Получить значение минимального расстояния между точками, задающего степень детализации кривой и сохранить в переменной m_minlen.

3. Выполнить построение кривой Безье с использованием этого параметра (функция PerformBezier).

4. По полученной кривой заполнить буфер data значениями АЧХ с интерполяцией между точками кривой (функция CreateResp).

Рассмотрим части алгоритма по отдельности.

Создание контура

В функции CreateContour формируется двунаправленный список опорных точек кривой. Каждая точка хранится в экземпляре класса (здесь показаны не все функции).

class ContourNode {
public:
	ContourNode(const double x, const double y) {
		m_x = x; m_y = y;
	}
	void SetGuides(const double x1, const double y1, const double x2, const double y2) {
		m_hasGuides = true;
		m_x1 = x1;
		m_y1 = y1;
		m_x2 = x2;
		m_y2 = y2;
	}
	auto HasGuides() const {
		return m_hasGuides;
	}
	auto GetNext() const {
		return m_next;
	}
	// ...
private:
	ContourNode* m_next = nullptr;
	ContourNode* m_prev = nullptr;
	bool m_hasGuides = false;
	double m_x = -1, m_y = 0; // contour point
	double m_x1 = -1, m_y1 = 0; // first guide point
	double m_x2 = -1, m_y2 = 0; // second guide point
};

В нём хранятся координаты опорной точки и двух направляющих точек. Если в двух соседних полосах эквалайзера одинаковые значения громкости, получается горизонтальный участок АЧХ, для которого не нужно строить кривую, и в таком случае задается только опорная точка без направляющих. Функция HasGuides возвращает true, если направляющие заданы, и false в противном случае. Направляющие точки задаются с помощью функции SetGuides. Указатели m_prev и m_next указывают на предыдущую и следующую точки в двунаправленном списке.

Класс Contour хранит двунаправленный список всех опорных точек кривой (здесь показаны не все функции).

class Contour
{
public:
	~Contour() {
		Clear();
	}
	bool IsEmpty() const {
    	return (m_nodes < 1);
    }
	auto TooShort() const {
		return (m_nodes < 2);
	}
	void Clear();
	auto GetHead() const {
		return m_head;
	}
	ContourNode* AddAfter(ContourNode* curnode, const double x, const double y);
	ContourNode* AddTail(const double x, const double y);
private:
	ContourNode* m_head = nullptr;
	ContourNode* m_tail = nullptr;
	uint32_t m_nodes = 0;
};

Функции класса позволяют получить указатель на его «голову» (первую точку), а также добавить новую опорную точку в конец списка либо после точки в середине списка. Указатели m_head и m_tail содержат указатели на первую и последнюю точки, а переменная m_nodes – длину списка.

Рассмотрим код функции создания контура в классе CBezier.

void CBezier::CreateContour(const uint32_t nbands, double* dB, const double step)
{
	// create m_contour, freq in octaves from base freq
	m_contour.Clear();
	uint32_t band = 0;
	const double bw = step * 0.125;
	double f = 0, db = dB[band++], dbnext, fnext = step;
	double f1, db1, f2, db2; // guide points
	// first point at freq = 0 oct (base freq)
	auto curnode = m_contour.AddTail(f, db);
	while (band < nbands) {
		dbnext = dB[band];
		fnext = f + step;
		if (dbnext == db) // no guide points
			curnode = m_contour.AddTail(fnext, dbnext);
		else { // add a point to make bandwidth
			curnode = m_contour.AddTail(f + bw, db);
			f1 = f2 = f + step * 0.5; // middle between bands
			db1 = db; db2 = dbnext;
			// set guide points
			curnode->SetGuides(f1, db1, f2, db2);
			curnode = m_contour.AddTail(fnext - bw, dbnext);
		}
		db = dbnext; f = fnext;
		band++;
	}
	// edge of last highpass band
	m_contour.AddTail(fnext, db);
}

Параметры: nbands – количество полос, dB – указатель на массив громкостей в полосах, step – ширина полосы в октавах. Описание алгоритма:

1. Перед началом цикла текущее смещение в октавах f устанавливается нулевым.

2. В начало списка добавляется точка, соответствующая нулевому смещению в октавах, и громкостью в первой полосе, которая запоминается в переменной db. Направляющих точек у неё нет.

3. Далее запускается цикл. Пока не достигнута последняя полоса, на каждом шаге цикла получается значение следующего смещения: fnext = f + step, т.е. к текущему смещению добавляется ширина полосы в октавах. Также берется громкость следующей полосы dbnext.

4. Если громкость dbnext совпадает с громкостью в текущей полосе db, в конец списка добавляется точка со смещением fnext и громкостью dbnext, но без направляющих.

5. Если громкость dbnext отличается от db, то сначала в конец списка добавляется точка с текущим смещением f плюс 1/8 ширины полосы и громкостью db. Также добавляется вторая точка со смещением fnext минус 1/8 ширины полосы и громкостью dbnext. Для первой добавленной точки задаются две направляющих со смещением, соответствующим середине полосы: f + step * 0.5 и громкостями db и dbnext.

6. Значение dbnext присваивается переменной db, а значение смещения fnext сохраняется в переменной f.

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

Пояснение к пункту 5. Если посмотреть внимательно на график желаемой АЧХ, приведенный выше, можно заметить, что в верхних и нижних точках перегиба участки кривой получаются совсем узкими, без горизонтального участка. Мне такой вид не понравился, поэтому я решил для каждой частотной полосы добавлять не по одной, а по две опорных точки с плоским горизонтальным участком длиной 1/4 ширины полосы между ними. Не уверен, что разница в обработке была бы как-то различима на слух, но мне показалось, что так будет правильнее:

Построение кривой Безье

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

В данном случае в качестве пороговой величины используется значение:

const double minint = step / 24.0;
constexpr double mindb = 0.1;
m_minlen = minint * minint + mindb * mindb;

Оно сохраняется в переменной m_minlen и соответствует характерному «расстоянию» в 0.1 дБ по вертикали и 1/24 ширины полосы по горизонтали (подобраны опытным путём). Для проверки минимальности расстояния между двумя точками используется функция:

bool CBezier::IsSmall(const double x1, const double y1, const double x2, const double y2) const
{
	return ((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1) < m_minlen);
}

Построение кривой выполняется в функции:

void CBezier::PerformBezier()
{
	if (m_contour.IsEmpty())
		return;
	auto curnode = m_contour.GetHead();
	while (curnode != nullptr) {
		if (curnode->HasGuides())
			BezierSample(curnode);
		curnode = curnode->GetNext();
	}
}

Цикл выполняется, начиная с «головы» по всему списку точек, и для каждой точки, имеющей направляющие, вызывается рекурсивная функция BezierSample.

void CBezier::BezierSample(ContourNode* curnode)
{
	const auto next = curnode->GetNext();
	if (next == nullptr)
		return;
	// first key point
	auto x0 = curnode->GetX(), y0 = curnode->GetY();
	// second key point
	auto x3 = next->GetX(), y3 = next->GetY();
	if (IsSmall(x0, y0, x3, y3))
		return;
	// first guide point
	auto [x1, y1] = curnode->GetGuide1();
	// second guide point
	auto [x2, y2] = curnode->GetGuide2();
	// calculate first middle point
	x0 = (x0 + x1) * 0.5;
	y0 = (y0 + y1) * 0.5;
	// calculate second middle point
	x3 = (x2 + x3) * 0.5;
	y3 = (y2 + y3) * 0.5;
	// calculate guide middle point
	double x = (x1 + x2) * 0.5;
	double y = (y1 + y2) * 0.5;
	// calculate new guide points
	x1 = (x0 + x) * 0.5;
	y1 = (y0 + y) * 0.5;
	x2 = (x3 + x) * 0.5;
	y2 = (y3 + y) * 0.5;
	// calculate new middle point
	x = (x1 + x2) * 0.5;
	y = (y1 + y2) * 0.5;
	// add new node at new middle point
	auto newnode = m_contour.AddAfter(curnode, x, y);
	curnode->SetGuides(x0, y0, x1, y1);
	BezierSample(curnode);
	newnode->SetGuides(x2, y2, x3, y3);
	BezierSample(newnode);
	return;
}

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

1. Расстояние между текущей и следующей опорными точками. Если оно меньше порогового значения (IsSmall), производится выход из функции (прерывание рекурсии).

2. Координаты двух точек, находящихся посередине между первой опорной и первой направляющей (M0) и посередине между второй опорной и второй направляющей (M3).

3. Координата точки M, находящейся посередине между двумя направляющими.

4. Координаты новых направляющих G1 и G2, находящихся посередине между М0 и M и посередине между M3 и M.

5. Координаты новой опорной точки N, находящейся посередине между G1 и G2. Точка N вставляется в двунаправленный список после текущей опорной точки.

6. Для текущей опорной точки задаются новые направляющие M0 и G1 и вызывается этот же алгоритм (рекурсия 1).

7. Для новой точки N задаются направляющие точки G2 и M3 и вызывается этот же алгоритм (рекурсия 2).

В результате получается список точек, максимально близко расположенных к нужной кривой.

Вычисление импульсной характеристики

Остались последние шаги, первый из которых – заполнение буфера для дискретного обратного быстрого преобразования Фурье (БПФ). Пока что мы имеем список точек, каждая из которых содержит значение смещения в октавах от базовой частоты первой полосы эквалайзера и значение громкости в ней в децибелах. Заполнение буфера по этому списку реализовано в функции CreateResp.

Как уже было рассказано в предыдущей статье, если длина БПФ составляет fftsize точек, то нужно выделить буфер в два раза больше, т.к. преобразование выполняется над комплексными числами (с нулевой мнимой частью). Чтобы получить импульсную характеристику, нужно заполнить первую половину буфера значениями АЧХ для частот от 0 до 1/2 от частоты дискретизации, а вторую половину – зеркально отраженной копией этих же значений относительно центра. Не буду здесь приводить полный код процедуры, а отмечу только основные моменты:

1. Каждое значение АЧХ в децибелах переводится в абсолютное значение процедурой dB2Gain.

2. Для первой полосы значение АЧХ заполняется в элементы буфера с индексами от нуля до индекса, пересчитанного из базовой частоты первой полосы:

   uint32_t index = 0, index1 = uint32_t(bfreq * double(fftsize) + 0.5);

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

   index1 = uint32_t(FreqOct(bfreq, fnext) * double(fftsize) + 0.5);

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

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

6. Для последней полосы с самой высокой частотой значение АЧХ в ней заполняется до элемента с индексом, равным fftsize/2-1.

7. При заполнении каждого элемента этим же значением заполняется элемент, находящийся на таком же расстоянии от центра буфера.

После заполнения буфера вторым шагом производится обратное БПФ и вычисляется импульсная характеристика фильтра, как это было уже описано в предыдущей статье. Длина фильтра ограничивается нужным значением, и при необходимости применяется оконная функция для уменьшения пульсаций АЧХ (эффекта Гиббса). Если необходимо применить дополнительное общее изменение громкости (master gain), оно умножается на центральный отсчёт характеристики.

В конце статьи приведена ссылка на репозиторий с исходным кодом программы на языке C++. Это программа для расчёта КИХ-фильтров, уже описанная в предыдущей статье. К ней добавлен класс CBezier, в котором реализованы перечисленные выше алгоритмы расчёта АЧХ, а в классе CFIRFilter добавлена функция CalcFIREqu, реализующая весь алгоритм с окончательным расчётом импульсной характеристики фильтра и применением оконных функций при необходимости. Помимо окна Кайзера, описанного ранее, добавлены ещё несколько оконных функций (код для них вынесен в отдельный файл window.h), которые при желании можно использовать.

Описание командной строки

Программа запускается из командной строки, генерирует КИХ-фильтр по заданным параметрам и сохраняет изображение с графиком АЧХ (при необходимости – отдельный файл с импульсной характеристикой фильтра) в выходной файл. В репозитории в подпапке win прилагается собранный для Windows исполняемый файл программы filters.exe, также программу можно собрать под Linux (протестировано под Ubuntu 22.04.5).

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

filters.exe [опции]

Остальные опции такие же, как для НЧ/ВЧ-фильтров, и описаны в предыдущей статье. Запуск программы с опцией -helpequ выдает полное их описание.

Второе состояние позволяет получить изображения АЧХ и импульсных характеристик для последовательности фильтров, при этом в каждой частотной полосе производится плавное изменение громкости от первого состояния ко второму за указанное количество шагов. Это позволяет имитировать перемещение регуляторов громкости эквалайзера и оценить, как при этом изменяется реальная АЧХ.

Названия выходных файлов

Названия выходных файлов с изображениями формируются программой автоматически. Для одиночных фильтров название файла с графиком АЧХ имеет вид:

Equ10b_8191pt_fft8192_exp_dB_wnd.gif

Название файла с графиком импульсной характеристики:

Equ10b_8191pt_imp_wnd.gif

Для последовательностей фильтров с изменением длины – отдельные файлы:

Equ10b_fft4096_exp_dB_wnd_0005_3295pt.gif
Equ10b_imp_wnd_0005_3295pt.gif

Здесь 0005 – номер фильтра в последовательности, 3295pt – длина фильтра в отсчётах. В случае с анимированным GIF (один файл):

Equ10b_fft8192_exp_dB_wnd_8191-591pt_20.gif

Здесь 8191-591pt – диапазон изменения длины фильтра, 20 – количество шагов изменения.

Для последовательностей фильтров с переходом между двумя состояниями эквалайзера – отдельные файлы:

Equ10b_fft8192_exp_dB_wnd_8191pt_0005.gif
Equ10b_imp_wnd_8191pt_0005.gif

Здесь 0005 – номер фильтра в последовательности. В случае с анимированным GIF (один файл):

Equ10b_50_fft8192_exp_dB_wnd.gif
Equ10b_50_imp_wnd.gif

Здесь 50 – количество шагов изменения между двумя состояниями.

Примеры запуска программы

Некоторые из параметров запуска являются обязательными, а именно:

• длина фильтра в параметре -fltequ, четное число будет приведено к нечетному

• путь к выходной папке (-outfolder)

• громкости в частотных полосах (-equdB), не менее 2

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

filters.exe -fltequ 8191 -fltwnd -equdB 0;3;0;6;0;-3;-6;0;12;0 -outfolder D:\tmp\flt -dB -top 20 -range 40 -grid -height 400 -width 600 -exp 20 -gif

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

Saved image 600 x 400 pixels to file: D:\tmp\flt\Equ10b_8191pt_fft8192_exp_dB_wnd.gif

Это означает, что сгенерирован КИХ-фильтр 10-полосного эквалайзера длиной 8191 отсчёт, изображение АЧХ с экспоненциальной шкалой частот по горизонтали и шкалой значений по вертикали в децибелах сохранено в выходное изображение в формате GIF размерами 600х400 пикселов. При расчёте использованы БПФ длиной 8192 точек и оконная функция Кайзера. Полученное изображение:

Пример вызова для расчёта 20-полосного эквалайзера с шириной полосы 1/2 октавы:

filters.exe -fltequ 8191 -fltwnd -octstep 1/2 -equdB 0;3;0;6;0;-3;-6;0;12;0;-3;3;0;6;0;-3;-6;0;12;0 -outfolder D:\tmp\flt -fft 65536 -dB -top 20 -range 40 -grid -height 400 -width 600 -exp 20 -gif

Полученное изображение АЧХ:

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

filters.exe -fltequ 8191 -fltwnd -equdB 0;3;0;6;0;-3;-6;0;12;0 -equdB2 0;-6;0;3;0;6;-3;0;-12;0 -freqsteps 50 -delay 20 -outfolder D:\tmp\flt -dB -top 20 -range 40 -grid -height 400 -width 600 -exp 20 -gif

Полученное изображение АЧХ:

Также можно также оценить, как изменяется АЧХ при уменьшении длины фильтра:

filters.exe -fltequ 8191 -fltwnd -equdB 0;3;0;6;0;-3;-6;0;12;0 -equdB2 0;-6;0;3;0;6;-3;0;-12;0 -freqsteps 50 -delay 20 -outfolder D:\tmp\flt -dB -top 20 -range 40 -grid -height 400 -width 600 -exp 20 -gif

Полученное изображение АЧХ:

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

Заключение

Как видно из последнего примера вызова, подтверждается то, что при уменьшении длины КИХ-фильтра ухудшается соответствие желаемой и реальной АЧХ в низких частотах. Для 10-полосного эквалайзера желательно использовать длину фильтра не менее 4000 точек (лучше больше). Но не стоит забывать о том, что обработка сигнала длинным фильтром увеличивает задержку сигнала и повышает нагрузку на процессор при обработке.

Также можно сделать вывод о том, что использование кривых Безье для «рисования» плавного графика АЧХ в частотной области вполне позволяет решить задачу получения КИХ-эквалайзера с нужной АЧХ.

Желающие могут по ссылкам ниже изучить исходный код программы, запустить её самостоятельно (для Windows имеется собранный исполняемый файл), чтобы попробовать сгенерировать КИХ-эквалайзеры с нужной длиной фильтра, количеством частотных полос и громкостями в полосах. Можно посмотреть графики их АЧХ и обработать ими звуковые файлы. Также приведены ссылки на книги и материалы, использованные при написании программы. Ссылки на ранее использованные материалы и использованные библиотеки можно найти в предыдущей статье.

Ссылки и литература

1. Репозиторий с программой и исходным кодом на C++: https://github.com/Voldemar-d/filters.git

2. Предыдущая статья

3. Графические эквалайзеры

4. Варга Б., Димень Ю., Лопариц Э. «Язык, музыка, математика»

5. Арбонес Х., Милруд П. «Числа — основа гармонии. Музыка и математика» (том 12 серии «Мир математики»)

6. Статья про кривые Безье на Хабре

7. Анимированные кривые Безье