Обработка цифрового звука фильтрами

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

Оглавление

• Введение

• Начнем с простого

• Частотные фильтры и АЧХ

• Реализация цифровых фильтров сигнала

• Оконная функция

• Расчет фильтров с произвольной АЧХ

• Описание программы

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

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

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

• Некоторые наблюдения и выводы

• Заключение

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

Введение

Цифровой звук в компьютере представлен в виде последовательности отсчётов сигнала, сгенерированных с помощью синтезаторов либо записанных с каких-либо источников звука (микрофонов, электронных музыкальных инструментов и др.). Преобразование звука в цифровой вид производится с помощью аналогово-цифровых преобразователей (АЦП), которые сегодня есть в практически любой звуковой карте. При оцифровке звуковой сигнал записывается с периодичностью, определяемой частотой дискретизации (sampling frequency, далее Fs). Каждый отсчёт сигнала хранится в виде числа с какой-либо точностью. Самыми популярными форматами таких чисел являются знаковое целочисленное 16-битное значение (например, так хранится звук на Audio CD), а также числа с плавающей точкой типа float (точность представления 24 бит) – в таком формате звук удобно хранить и обрабатывать в памяти.

Частота дискретизации задает диапазон частот, которые возможно записать в оцифрованном сигнале. Согласно теореме Котельникова, при значении частоты дискретизации Fs из цифрового сигнала возможно восстановить аналоговый, в котором присутствуют частоты не более Fs/2. Например, в звуке, записанном на Audio CD с частотой дискретизации 44100 Гц, могут присутствовать звуки частотой не выше 22050 Гц.

Начнем с простого

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

#define DB_INF -999

double dB2Gain(double dB) {
	if (dB < DB_INF) return 0;
	return pow(10.0, dB / 20.0);
}

double Gain2dB(double g) {
	if (g <= 0) return DB_INF;
	double dB = double(20.0 * log10(g));
	if (dB < DB_INF) dB = DB_INF;
	return dB;
}

Изменение амплитуды в 2 раза соответствует примерно 6 дБ. То есть, чтобы сделать звук громче в 2 раза, нужно изменить громкость на 6 дБ, а чтобы уменьшить в 2 раза – изменить на -6 дБ. Константа DB_INF используется для обозначения «громкости» полной тишины – как известно, значение логарифма стремится к минус бесконечности при стремлении аргумента к нулю.

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

std::vector<int16_t> samplesInt;
std::vector<float> samplesFlt;

Процедура изменения громкости звука в контейнере:

template<typename T>
void ApplyGain(std::vector<T>& samples, const double dB, const size_t count) {
	const float fg = float(dB2Gain(dB));
	const size_t nsmp = std::min(count, samples.size());
	for (size_t i = 0; i < nsmp; i++) {
		if constexpr (std::is_floating_point<T>)
			samples[i] *= fg;
		else if constexpr (std::is_same_v<T, int16_t>) {
			const float s = float(samples[i]) * fg;
			if (s < SHRT_MIN)
				samples[i] = SHRT_MIN;
			else if (s > SHRT_MAX)
				samples[i] = SHRT_MAX;
			else
				samples[i] = int16_t(s);
		}
		else
			static_assert(false, "T must be float or int16_t!");
	}
}

// значение громкости может быть взято из файла настроек, введено пользователем и т.д.
const double dB = 6.0;
// количество отсчётов, которые должны быть обработаны
size_t count = 1024;
ApplyGain(samplesInt, dB, count);
ApplyGain(samplesFlt, dB, count);

Переменная count задает количество отсчётов звука, которые нужно обработать – это может пригодиться в случае, если контейнер заполнен не весь. Например, звук может считываться из файла порциями разной длины в заранее выделенный буфер заведомо достаточного размера. При желании можно расширить данную функцию для обработки других форматов несжатого звука (например, 24/32-битный целочисленный).

Частотные фильтры и АЧХ

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

Как можно реализовать подобную обработку в своей программе?

Решение подобных задач изучается в области цифровой обработки сигналов, или DSP (digital signal processing). Для изменения громкости разных частотных компонент применяются цифровые фильтры сигнала. Чтобы оценить, каким образом фильтр влияет на громкость разных частотных компонент, нужно посмотреть на график его амплитудно-частотной характеристики (АЧХ). Пример подобного графика приведен в заглавном изображении к данной статье – это фильтр низких частот, или НЧ-фильтр (правильнее его называть «фильтр, пропускающий низкие частоты», или «low-pass filter»).

Чтобы построить график АЧХ, нужно сначала обработать фильтром цифровой сигнал, в первом отсчёте которого находится значение 1, а во всех остальных – 0, то есть, полная тишина (в математике это соответствует функции Дирака, или дельта-функции). В результате такой обработки получится импульсная характеристика фильтра. Затем нужно произвести над ней преобразование Фурье (ПФ), которое отображает сигнал из временной области в частотную, и из полученного набора чисел можно получить АЧХ и изобразить его в виде графика.

Здесь есть два важных момента: ПФ выдает набор комплексных чисел с действительной и мнимой частями, и отображение производится в диапазон частот от 0 Гц до частоты дискретизации Fs. Выше мы уже упоминали, что в сигнале при этом могут присутствовать частотные компоненты только до значения Fs/2. Например, если мы обработаем импульсную характеристику длиной 1024 отсчёта, то в результате ПФ такой же длины получим только 512 отсчётов полезной информации для частот от 0 до Fs/2. Оставшаяся половина с частотами выше Fs/2 будет зеркальным отображением, её нужно просто отбросить.

Для вычисления преобразования Фурье в данном материале используется одна из многих реализаций быстрого преобразования Фурье (БПФ), или Fast Fourier Transform (FFT). Имеется класс TFFTF, который на входе берет массив из N отсчётов сигнала, каждый из которых должен быть дополнен нулевой мнимой частью, а N для быстроты вычислений должно быть степенью двойки. В классе имеется функция SetSize, задающая длину преобразования, а также функция CDFT, собственно вычисляющая комплексное БПФ над массивом данных (результат размещается в нем же):

class TFFTF {
public:
	void SetSize(int N);
	void CDFT(float* data);
	void CDFTI(float* data);
};

Допустим, также имеется функция Process, осуществляющая обработку сигнала цифровым фильтром по какому-то алгоритму. Тогда алгоритм вычисления АЧХ можно представить в таком виде:

1. Выделить массив длиной N отсчётов сигнала в формате float (числа с плавающей точкой). В самый первый элемент массива поместить значение 1, остальные заполнить нулями.

2. Обработать массив вызовом функции Process, в результате в массиве будет импульсный отклик фильтра.

3. Найти минимальную длину БПФ, необходимую для обработки – минимальное значение степени двойки fsz, которое больше либо равно N. Чтобы частотная выборка АЧХ была достаточно подробной, это значение лучше делать не меньше какого-то значения – например, 1024.

4. Выделить второй массив, длина которого должна составлять 2 * fsz. В каждый нечетный элемент скопировать значения из массива с импульсным откликом, а все нечетные заполнить нулями, чтобы получить комплексные числа. Все незадействованные в копировании элементы также заполнить нулями.

5. Выполнить комплексное БПФ над вторым массивом. В нем же будут содержаться значения комплексных чисел с результатом вычисления.

6. Из каждого комплексного числа нужно взять модуль, т.е. квадратный корень из суммы квадратов действительной и мнимой частей. Полученные значения и будут значениями АЧХ фильтра в интервалах, соответствующих разбиению частотного диапазона от 0 до Fs/2 на количество частей, равное половине длины БПФ, т.е. fsz/2.

7. Далее полученный массив можно изобразить в виде графика, в котором минимальное значение 0, а максимальное 1 (иногда несколько больше, об этом будет сказано ниже). Также можно каждую величину массива перевести в децибелы, чтобы получить более информативный график АЧХ.

Изобразим этот алгоритм в виде кода:

template<bool dB>
void GetFreqResp(const int nSamples, float* imp, std::vector<float>& resp) {
	// найти минимальную длину БПФ
	int fsz = 2;
	while (fsz < nSamples)
		fsz <<= 1;
	fsz = std::min(fsz, 1024);
	// преобразовать импульсный отклик в комплексные числа для БПФ
	std::vector<float> cimp;
	cimp.resize(size_t(fsz) * 2, 0.0f);
	for (int i = 0; i < nSamples; i++)
		cimp[i * 2] = imp[i];
	auto fdata = cimp.data();
	// выполнить БПФ нужной длины
	TFFTF fft;
	fft.SetSize(fsz);
	fft.CDFT(fdata);
	// посчитать значения АЧХ до частоты Fs/2 из комплексных чисел
	const int nf = int(fsz / 2);
	double re, im, val;
	resp.resize(nf);
	for (int i = 0; i < nf; i++) {
		re = fdata[i * 2];
		im = fdata[i * 2 + 1];
		// модуль комплексного числа
		val = sqrt(re * re + im * im);
		if constexpr(dB)
			resp[i] = (float)Gain2dB(val);
		else
			resp[i] = (float)val;
	}
}

auto Process = [](const size_t nSamples, std::vector<float>& signal) {
	for (size_t i = 0; i < nSamples; i++) {
		//signal[i] = ...
	}
};

// сигнал с единичным первым отсчётом
const int nSamples = 1000;
std::vector<float> imp;
imp.resize(1000, 0.0f);
imp[0] = 1.0f;
// обработать сигнал фильтром
Process(imp.size(), imp);
// получить АЧХ из импульсной характеристики
std::vector<float> resp;
GetFreqResp<false>(nSamples, imp.data(), resp);
// получить АЧХ в децибелах
GetFreqResp<true>(nSamples, imp.data(), resp);

В данном примере функция Process ничего не делает. Запустив этот код под отладчиком, легко увидеть, что в результате его работы весь массив resp будет заполнен значениями 1.0 после первого вызова и значениями 0 дБ после второго. Поскольку сигнал никак не обрабатывается, АЧХ совершенно ровная, и показывает, что сигнал во всех частотах никак не изменяется.

Как реализовать функцию Process, чтобы она изменяла громкость частотных составляющих сигнала нужным образом?

Реализация цифровых фильтров сигнала

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

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

КИХ-фильтр при обработке использует только набор коэффициентов h[k] и отсчёты входного сигнала:

Здесь y[n] — отсчёт выходного сигнала в момент времени n; x[n-k] — отсчёт входного сигнала на k позиций раньше текущего; h[k] — коэффициент фильтра с индексом k; N — длина фильтра (количество коэффициентов).

Набор коэффициентов КИХ-фильтра h[k] является его импульсной характеристикой. Если обработать таким фильтром импульсный сигнал, рассмотренный выше, то на выходе будет получен тот же набор коэффициентов фильтра. Обработку сигнала КИХ-фильтром также называют сверткой сигнала с импульсной характеристикой фильтра. Здесь N – длина импульсной характеристики, для простоты её также называют просто длиной КИХ-фильтра. Наглядная схема такой обработки приведена в статье (знаками X показаны умножители, знаками + - сумматоры):

Из формулы видно, что для получения каждого отсчёта выходного сигнала y нужно использовать текущий отсчёт входного сигнала x и N-1 предыдущих отсчётов входного сигнала. То есть, на каждом шаге обработки нужно использовать и обновлять массив чисел с «историей», в которой нужно хранить предыдущие отсчёты сигнала. Существуют разные реализации таких алгоритмов. Самая простая и «наивная» - на каждом шаге обработки сдвигать в массиве N-1 чисел назад, а текущий отсчёт сохранять в последнем элементе.

Более правильный и быстрый способ – использовать циклический массив, не сдвигая данные, а обновляя только текущую позицию в нем, с проверкой, что эта позиция не вышла за границы массива. Такой алгоритм реализован в классе CFIRFilter в процедуре Apply (ссылка на репозиторий с кодом в конце статьи).

В нем же реализованы процедуры для получения импульсных характеристик фильтров: НЧ-фильтр (пропускающий низкие частоты), ВЧ-фильтр (пропускающий высокие частоты). Рассмотрим для начала расчет НЧ-фильтра. График АЧХ идеального НЧ-фильтра, начиная с нулевой частоты, выглядит как горизонтальная линия со значением 1.0 (или 0 дБ - фильтр пропускает низкие частоты без изменений) – этот диапазон частот называют полосой пропускания. Далее, начиная с некоторой частоты значения АЧХ падают нуля – фильтр подавляет высокие частоты. Основным параметром такого фильтра является значение частоты среза – это граничная частота, «в районе» которой начинается диапазон (или полоса) подавления.

Для длины фильтра N, частоты среза F Гц и частоты дискретизации Fs Гц алгоритм расчета коэффициентов НЧ-фильтра следующий:

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

int CalcLowPassFIR(const int N, const int F, const int Fs, std::vector<float>& hImp) {
	int n = N;
	if (0 == n % 2) n--;
	const int Np = n / 2;
	hImp.resize(n);
	auto h = hImp.data();
	const double freq = double(F) / double(Fs);
	const double ratio = 0.5 / freq;
	const double pi_r = M_PI / ratio;
	double ck, dk;
	int i, j, k;
	for (k = 0, i = j = Np; k <= Np; k++, i++, j--) {
		dk = double(k);
		if (0 == k)
			ck = 1.0;
		else
			ck = ratio * sin(dk * pi_r) / (dk * M_PI);
		h[i] = h[j] = float(ck / ratio);
	}
	return n;
}

Он реализует формулу расчета коэффициентов идеального НЧ-фильтра:

Из алгоритма видно, что длина фильтра всегда нечетная, и импульсная характеристика фильтра симметрична относительно середины. Симметричный фильтр не искажает фазу обрабатываемого сигнала, а также это позволяет ускорить процедуру обработки за счет сокращения количества умножений почти в 2 раза. В классе CFIRFilter это реализовано в процедуре ApplySym.

Оконная функция

Попробуем произвести расчет фильтра по приведенному алгоритму и посмотреть, как выглядит его АЧХ. Например, для НЧ-фильтра с частотой среза 5000 Гц длиной 255 отсчётов:

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

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

Рассмотрим изменение АЧХ такого фильтра для длины от 63 до 1023 отсчётов:

Здесь видно, что с увеличением длины фильтра ширина выбросов пульсации уменьшается, но их максимальный уровень остаётся примерно тем же, и их амплитуда не уменьшается. Для борьбы с этим эффектом импульсную характеристику дополнительно умножают на некоторую функцию, значение которой в середине равно 1, а к краям плавно уменьшается до 0. Это так называемая оконная функция, самый простой вид которой – треугольная (окно Бартлетта):

Существуют и другие виды оконных функций: Хэмминга, Ханна (Хеннинга), Блэкмана-Харриса, Гаусса, но одной из самых оптимальных оконных функций является окно Кайзера:

Формула расчета данной функции:

Здесь n — индекс отсчёта (0 ≤ n ≤ N−1, N — длина окна); β — параметр, который позволяет управлять компромиссом между шириной главного лепестка спектра и уровнем боковых лепестков; I₀ — модифицированная функция Бесселя первого рода нулевого порядка.

Для расчета значения модифицированной функции Бесселя I₀ с достаточной нам точностью можно использовать следующий алгоритм, использующий разложение функции в ряд Тейлора:

double izero(double y) {
	double s = 1.0, ds = 1.0, d = 0.0;
	do {
		d += 2.0;
		ds = ds * (y * y) / (d * d);
		s = s + ds;
	} while (ds > 1E-7 * s);
	return s;
}

Рассмотрим АЧХ тех же фильтров, рассчитанных с использованием окна Кайзера (от 63 до 1023 отсчётов):

Как видно, проявления эффекта Гиббса практически сошли на нет. Может показаться, что их не стало совсем, но правильнее смотреть АЧХ не в виде абсолютных значений, а в децибелах – пример для 255 отсчётов:

В районе полосы пропускания колебания АЧХ остались, но совсем небольшие – разницу в громкости в доли децибела человеческое ухо различить не способно. В полосе подавления остались выбросы, пики которых ниже уровня -90 дБ. Это очень хороший результат – полный динамический диапазон 16-битного звука, как на Audio CD, составляет около 96 дБ. Это означает, что если данным фильтром обработать 16-битный цифровой звук, то выше примерно 5 с небольшим килогерц в сигнале будет практически полный ноль, т.е. тишина.

Расчет фильтров с произвольной АЧХ

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

Откуда она взялась? Вообще говоря, она получена именно так: задан график АЧХ идеального НЧ-фильтра, который до определенной частоты должен не изменять сигнал, а начиная с частоты среза должен полностью его подавлять. То есть, в нижних частотах задано значение АЧХ, равное 1, а начиная с частоты среза и до Fs/2 оно равно 0. Затем над этой АЧХ произведено обратное дискретное преобразование Фурье.

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

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

1. Для получения КИХ-фильтра длиной N отсчётов нужно сначала определить минимальное значение длины БПФ: взять двойную длину фильтра и найти значение степени двойки, которое больше либо равно 2*N. Для этого в классе есть функция getFFTsize, результат сохраним в переменной fftSz.

2. Поскольку нужно произвести комплексное БПФ, необходимо выделить массив m_buf длиной еще в 2 раза больше, то есть 2*fftSz, и заполнить его нулями.

3. Для частоты среза nFltFreq и частоты дискретизации nSampleRate определить номер отсчёта в АЧХ, соответствующий частоте среза:

   int nF = int(double(m_buf.size()) * double(nFreq) / double(nSampleRate));

4. Комплексные числа в буфере m_buf, начиная с нулевого и вплоть до числа с индексом nF, заполнить действительным значением 1 (остальные в п.2 были заполнены нулями).

5. Как мы помним, в результате прямого БПФ значения АЧХ зеркально отображаются из диапазона от 0 до половины частоты дискретизации Fs/2 в диапазон от Fs/2 до Fs, поэтому при построении АЧХ мы также должны зеркально «размножить» АЧХ в верхнюю часть буфера – первое комплексное число скопировать в последнее, второе – в предпоследнее и т.д.

6. Выполнить обратное комплексное БПФ над буфером m_buf. В классе TFFTF для этого предназначена функция CDFTI.

7. Выделить буфер m_flt для импульсной характеристики фильтра длиной N отсчётов. Также рассчитать нормирующий коэффициент, обратный длине БПФ:

   double R = 1.0 / double(nFFTsz);

8. Из полученных в буфере m_buf комплексных чисел взять первое, его действительную часть умножить на R и поместить в центральный отсчёт фильтра с индексом N/2 (как мы помним, длина фильтра – нечетное число).

9. Далее, начиная со второго комплексного числа в буфере m_buf, в цикле скопировать N/2-1 действительных частей, умноженных на R, в отсчёты фильтра симметрично от центра к его краям.

10. При необходимости можно дополнительно применить окно Кайзера (или другое) к рассчитанной импульсной характеристике фильтра.

В упомянутом выше классе CFIRFilter основная часть алгоритма (пп.5-10) реализована в функции CalcFilter. Она позволяет рассчитать фильтр с произвольной АЧХ, если предварительно заполнить буфер m_buf нужными значениями коэффициентов усиления в действительной части чисел. Полный расчет, включающий весь алгоритм, а также применение окна Кайзера при необходимости, реализован в функции InvLowHighPass – в шаблонном параметре bLowPass нужно указать true для расчета НЧ-фильтра и false для ВЧ-фильтра.

Описание программы

Все изображения АЧХ и импульсных характеристик в данной статье получены в программе, ссылка на исходный код которой на языке C++ приведена в конце статьи. Программа запускается из командной строки и позволяет:

• рассчитывать НЧ/ВЧ КИХ-фильтры заданной длины

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

• применять либо не применять окно Кайзера к полученному фильтру

• задать путь к выходной папке для сохранения файлов (изображений и звуковых)

• сохранять АЧХ полученного фильтра в изображения формата BMP или GIF

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

• формировать график АЧХ с линейной либо экспоненциальной (логарифмической) шкалой частот

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

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

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

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

Упомянутые выше классы CFIRFilter и TFFTF являются частями программы. Предназначение всех классов программы следующее:

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

В репозитории в подпапке win прилагается собранный для Windows исполняемый файл программы filters.exe, также программу можно собрать под Linux (протестировано под Ubuntu 22.04.5). Запуск программы из командной строки:

filters.exe [опции]

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

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

LowPass_511pt_1000hz_fft4096_exp_abs_wnd.gif

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

LowPass_255pt_5000hz_imp_wnd.gif

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

HighPass_5000hz_fft4096_exp_abs_wnd_0005_255pt.gif

HighPass_5000hz_imp_wnd_0005_255pt.gif

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

HighPass_5000hz_fft4096_exp_dB_wnd_127-415pt_10.gif

HighPass_5000hz_imp_wnd_127-415pt_10.gif

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

Пример изображения импульсной характеристики ВЧ-фильтров длиной от 127 до 415 отсчётов:

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

HighPass_fft4096_exp_dB_wnd_127pt_0010_5900hz.gif

HighPass_imp_wnd_127pt_0010_5900hz.gif

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

HighPass_5000-5900hz_10_fft4096_exp_dB_wnd_127pt.gif

HighPass_5000-5900hz_10_imp_wnd_127pt.gif

Здесь 5000-5900hz – диапазон изменения частоты среза, 10 – количество шагов изменения.

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

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

• тип и длина фильтра (-flthp/-fltlp), четное число будет приведено к нечетному

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

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

filters.exe -outfolder d:\tmp\flt -fltlp 127

Программа по завершению работы выдаст следующую информацию:

Saved image 2048 x 480 pixels to file: d:\tmp\flt\LowPass_127pt_1000hz_fft4096_lin_abs_nownd.bmp

Это означает, что сгенерирован НЧ-фильтр длиной 127 отсчётов и частотой среза 1000 Гц, изображение АЧХ с линейной шкалой частот по горизонтали и абсолютной шкалой значений по вертикали (не в децибелах) сохранено в выходное изображение в формате BMP размерами 2048х480 пикселов. Ширина не указана, поэтому сохранена полная АЧХ длиной 2048 точек – взята ½ длины БПФ (4096 точек). При расчете использована формула идеального НЧ-фильтра без оконной функции, поэтому наблюдаются проявления эффекта Гиббса – пульсации АЧХ как в полосе пропускания, так и в полосе подавления. Приведем фрагмент этого изображения в районе частоты среза:

Чтобы получить изображение, приведенное в заголовке статьи, добавим опции, чтобы включить сетку на изображении, вертикальную шкалу в децибелах, экспоненциальную шкалу частот, а также применение оконной функции Кайзера; также увеличим длину фильтра до 511 точек, чтобы уменьшить ширину переходной полосы в районе частоты среза, и сохраним результат в формате GIF:

filters.exe -outfolder d:\tmp\flt -fltlp 511 -fltwnd -dB -top 10 -range 200 -exp 100 -width 600 -height 400 -gif

Теперь проверим обработку звукового файла с помощью полученного фильтра. Для примера в репозитории имеется файл sweep.wav, в котором сгенерирован специальный звук: синусоидальный сигнал, частота которого в течение 10 секунд меняется от 20 Гц до 20 кГц. Спектральная картинка такого сигнала выглядит так:

Поместим этот файл в папку с программой и обработаем полученным выше фильтром:

filters.exe -outfolder d:\tmp\flt -fltlp 511 -fltwnd -dB -top 10 -range 200 -exp 100 -gif -width 600 -height 400 -wav ./sweep.wav

Программа по завершению работы выдаст информацию также о формате и обработке WAV-файла:

Saved image 600 x 400 pixels to file: d:\tmp\flt\LowPass_511pt_1000hz_fft4096_exp_dB_wnd.gif
Loaded input WAV file: ./sweep.wav
|======================================|
Num Channels: 1
Num Samples Per Channel: 480000
Sample Rate: 48000
Bit Depth: 16
Length in Seconds: 10
|======================================|
Processed 480000 samples in 1 channels
Saved output WAV file: d:\tmp\flt\sweep.wav

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

Видно, что от частоты среза фильтра 1 кГц и выше звук постепенно переходит в тишину, но реально звука практически нет только на частоте примерно 1.5 кГц. Всё в соответствии с графиком АЧХ: по нему видно, что подавление звука на -100 дБ и более происходит как раз примерно в частотах от 1.5 кГц и выше.

Естественно, обрабатывать можно не только конкретно этот, но и любые другие WAV-файлы с несжатым звуком и любым количеством каналов – каждый канал будет обработан по отдельности.

Некоторые наблюдения и выводы

Снова вернемся к вопросу: возможно ли получить любой КИХ-фильтр по заданному графику АЧХ? Выше мы уже отмечали, что можно, но с ограничениями:

• возникают неравномерности АЧХ в виде «выбросов» (эффект Гиббса), которые можно частично (но не до конца) подавить применением оконных функций

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

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

Насколько существенны эти ограничения в реальных задачах обработки звука? Попробуем рассмотреть изменение АЧХ НЧ-фильтра длиной 255 отсчётов при уменьшении частоты среза от 1 кГц до 50 Гц. Для этого запустим программу следующим образом:

filters.exe -outfolder d:\tmp\flt -fltlp 255 -fltwnd -dB -top 10 -range 200 -exp 100 -gif -width 600 -height 400 -stepsize -50 -freqsteps 20 -delay 50

На первый взгляд, для НЧ-фильтра всё выглядит вполне приемлемо. А теперь рассмотрим такое же изменение ВЧ-фильтра, подавляющего низкие частоты:

Видно, что чем ниже частота среза, тем хуже ВЧ-фильтр работает - на частотах ниже 500 Гц очень плохо подавляет звук в полосе непропускания. Конечно, всё относительно, но подавление на 10…40 дБ – это не идёт в сравнение с подавлением на 100 дБ, когда сигнал превращается в почти полную тишину.

Попробуем увеличить длину фильтра до 4095 отсчётов и снова посмотреть изменение АЧХ:

Видно, что здесь ситуация гораздо лучше: ширина переходной полосы уже небольшая (десятки Гц), и даже на частоте среза ниже 200 Гц подавление сигнала около 100 дБ. Можно и дальше продолжить подобные эксперименты, но теперь уже можно сформулировать еще одну и, пожалуй, главную особенность КИХ-фильтров:

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

Насколько это плохо? Попробуем прикинуть, какая длина фильтра нужна, чтобы подавить ВЧ-фильтром частоты ниже 100 Гц на 100 дБ. Чтобы не утомлять читателя, приведу результат следующего эксперимента. Я запустил генерирование нескольких сотен ВЧ-фильтров с частотой среза 100 Гц, начиная с длины 2047 отсчётов, с увеличением длины на 100 отчетов на каждом шаге. Затем просмотрел изображения полученных АЧХ (чтобы можно было их получше рассмотреть, увеличил размеры картинок до 2000х1000 пикселей). В результате могу сделать примерно такой вывод: чтобы в полосе непропускания подавление таким фильтром было гарантированно не хуже 100 дБ, нужна длина фильтра 7000-8000 отсчётов и более.

Много это или мало? У КИХ-фильтров есть ещё два недостатка, которые до сих пор не упоминались. Во-первых, такой фильтр выдает обработанный сигнал с задержкой, равной ½ длины фильтра. Для фильтра длиной 8000 отсчётов задержка составит 4000 отсчётов, при частоте дискретизации звука 48000 Гц это составляет 1/12 секунды. Если обрабатывать таким фильтром звук в реальном времени (например, слушать результат обработки прямо во время записи звука), такая задержка уже ощутимо заметна на слух.

Во-вторых, как мы помним, обработка звука КИХ-фильтром означает, что при обработке каждого отсчёта входного сигнала мы должны выполнить линейную свертку с «историей» их отсчётов входного сигнала в почти полную длину фильтра. Даже при симметричной импульсной характеристике фильтра длиной 8000 отсчётов нужно выполнить 8000 сложений и 4000 умножений для получения каждого отсчёта обработанного сигнала. Итак, мы можем сформулировать ещё два недостатка КИХ-фильтров:

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

• требуется высокая вычислительная мощность при такой обработке, если хочется получить фильтр с АЧХ, максимально близкой к той, что необходима для обработки звука в низких частотах

Может сложиться впечатление, что недостатков у КИХ-фильтров как-то многовато получается – а в чем же тогда их преимущества? Можно сформулировать следующие:

1. Возможность получить фильтр с практически произвольной АЧХ (пусть и с ограничениями). Вполне можно создать программу, в которой можно «нарисовать» нужный график АЧХ и получить фильтр, который довольно близко её реализует.

2. Симметричная импульсная характеристика фильтра означает, что при обработке не искажается фаза обрабатываемого сигнала. На первый взгляд, не очень понятно, чем это хорошо. Для обработки звука это означает следующее: если с помощью такого фильтра изменять громкость разных частотных компонент, это будет звучать именно как просто изменение громкости (нужные частоты как бы «выходят на первый план») – при этом не возникает эффектов «звука как из бочки», «гулкости» и прочих подобных «некрасивых» искажений.

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

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

Заключение

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

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

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

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

2. Дельта-функция Дирака

3. Цифровая фильтрация на ПЛИС – Часть 2

4. Расчет КИХ-фильтров и эффект Гиббса

5. Оконные функции

6. Paul M. Embree "C algorithms for real-time DSP", Prentice Hall

7. Emmanuel C. Ifeachor, Barrie W. Jervis "Digital Signal Processing: A Practical Approach", Addison-Wesley; "Цифровая обработка сигналов: практический подход" Эммануэль С. Айфичер, Барри У. Джервис, "Вильямс"

8. Библиотека fmt

9. Вычисление БПФ

10. Библиотека ezdib

11. Библиотека gif-h

12. Библиотека AudioFile