Програмування&Музика: Частотний фільтр Баттервота. Частина 3

Всім привіт!
Ви читаєте третю частину статті про створення VST-синтезатора на З#. У попередніх частинах було розглянуто SDK і бібліотеки для створення VST плагінів, розглянуто програмування осцилятора і ADSR-огинаючої для керування амплітудою сигналу.
У цій частині я розповім, як розрахувати і закодить фільтр частот, без якого не обходиться жоден синтезатор. А без еквалайзера немислима обробка звуку)
Буде розглянуто вихідний код і застосування еквалайзера з бібліотеки NAudio (бібліотека для роботи зі звуком .NET).
Увага — буде багато матана — будемо розраховувати формули для коефіцієнтів фільтра.
Вихідний код написаного мною синтезатора доступний на GitHub'е.

Скріншот VST плагін-еквалайзера Fab Filter Pro Q

Цикл статей
  1. Розуміємо і пишемо VSTi синтезатор на C# WPF
  2. ADSR-обвідна сигналу
  3. Частотний фільтр Баттервота
  4. далі буде...
Зміст
  1. Еквалайзер
  2. Фільтрація через перетворення Фур'є
  3. Цифрові фільтри
  4. Чому фільтр Баттервота?
  5. Виведення формули фільтра НЧ
  6. Виведення формули фільтра ВЧ і смугового фільтра
  7. Програмування фільтра Баттервота за отриманими формулами
  8. Смуговий еквалайзер в бібліотеці NAudio
  9. Програми для розрахунку фільтрів
  10. Список літератури

Еквалайзер
Часто при обробці звуку ми хочемо змінити його характер/забарвлення/тембр. Зробити звук більш басовим, прибрати верхні частоти, або навпаки, зробити звук "прозорим", залишивши лише середину і верху. Впевнений, багато людей, які не працювали з обробкою звуку, знають що таке еквалайзер — вони є акустичних колонках, музичних центрах, магнітофонах, плеєрах, і т. д. Еквалайзер — це набір фільтрів, кожен з яких змінює амплітуду сигналу в обраній смузі частот. На побутових колонках це, зазвичай, 2-3 крутилки — низькі частоти, середні і верху, з фіксованими смугами частот.
В Winamp ' вському еквалайзері вже є 10 заздалегідь визначених смуг.

Скріншот еквалайзера в плеєрі Winamp
У світі обробки звуку існує безліч плагінів-еквалайзерів, на будь-який смак і колір. Плагін Fab Filter Pro Q (скріншот на початку статті) є параметричним еквалайзером, який дозволяє створити велику кількість смуг, редагувати їх параметри.
Кожна смуга в еквалайзері — це, по суті, фільтр частот. Фільтри частот змінюють тембральні/частотні характеристики сигналу. В електроніці існує багато типів і класифікацій фільтрів, з відповідними характеристиками і параметрами — дивимося вікіпедію.
Ми розглянемо і запрограммируем найпростіші фільтри: НЧ, ВЧ і смуговий фільтри.

Фільтрація через перетворення Фур'є
За ідеєю, вам ніхто не заважає робити з сигналом дискретне перетворення Фур'є, обробити частоти і потім зробити зворотне перетворення.
Якщо не думати над реалізацією ДПФ, то такий підхід я б назвав досить інтуїтивним і простим у програмуванні (знову ж таки, якщо взяти ДПФ з якої-небудь либы і не кодити самому).
Мінуси підходу — по-перше, ДПФ приймає на вхід масив з семплів, розмір якого є ступенем двійки. Це означає, що вихідний сигнал вже буде з затримкою. По-друге, кожен 512-й семпл ми будемо виробляти даний алгоритм: ДПФ, обробка частот сигналу, зворотне ДФП. Це не малі обчислення. По-третє, є ще мінуси і тонкощі, які знають адепти цифрової обробки сигналів.
Ми не будемо розглядати застосування ДПФ, а заглянемо в теорію цифрових фільтрів; напишемо фільтр, який обробляє значення семплів і має лінійну обчислювальну складність в залежності від довжини вхідного масиву семплів.

Цифрові фільтри
Більшу частину інформації і вивід формул я взяв з книги Digital Signal Processing: A Practical Approach — дуже рекомендую, вона є в російській версії — Цифрова обробка сигналів. Практичний підхід, зацікавлені знайдуть PDF в мережі.
Хочу зробити важливе зауваження. Тема побудови і розрахунку фільтра дійсно дуже складна, містить масу тонкощів і нюансів, вимагає знання та розуміння теорії. У цій статті я покажу, як розрахувати формули фільтра Баттервота, щоб у читача виникло розуміння, звідки виводяться ці формули. Але чому саме такі вихідні формули, чому саме такі заміни — можна зрозуміти лише занурившись в глибоку теорію цифрової обробки сигналів.
Коли я починав гуглити код фільтрів, я відразу знаходив безліч незрозумілого математичного коду, і хотілося хоч трохи зрозуміти, звідки беруться такі розрахункові формули. Осцилятори, обвідна, ділей — розуміння і програмування роботи цих складових особисто мені здається інтуїтивним, але тільки не фільтрів. Цією статтею я хочу пробудити інтерес до цифрової обробки сигналів) Буду радий, якщо виникне бажання розібратися в цій темі більш грунтовно.
Вам треба знати (хоча б трошки) такі терміни як згортка, імпульсна характеристика фільтра, передатна функція фільтра.

АЧХ фільтра (картинка з радянського підручника, не знайшов оригінал)
Фільтр змінює сигнал, "прибираючи" у ньому вибрані частоти. Існуючі фільтри не ідеальні. Смуга пропущення — це смуга частот, яку фільтр "не зачіпає" (на графіку є деяка зміни — особливість неідеального представленого фільтра). Смуга придушення — смуга небажаних частот. У смузі переходу відбувається спад частот. Природно, фільтр ближче до "ідеальному" тим, наскільки менше він спотворює смугу пропускання, наскільки сильно він пригнічує частоти в смузі придушення і наскільки вузька смуга переходу. Є різні "наближення" фільтрів — фільтр Чебишева, Баттервота, і так далі — їх ви знайдете в книжках та на просторах мережі.

Чому фільтр Баттервота?
Все дуже просто, АЧХ фільтра Баттерворта максимально гладка на частотах смуги пропускання — імхо, важливіше за все не зіпсувати сигнал в смузі пропускання.

Логарифмічна АЧХ для фільтрів нижніх частот Баттерворта різних порядків (скріншот взято з Вікіпедії)

Виведення формули фільтра НЧ
Передатна функція для фільтра Баттервота на s-площині записується наступними формулами:
при парних n
при непарних n
Тут n — порядок фільтра: амплітуда на частоті зрізу w дорівнює -3n dB, амплітудно-частотна характеристика затухає на −6n дб на октаву.
Щоб отримати згорткові коефіцієнти, потрібно отримати передавальну функцію на z-площині у вигляді

Знайдемо передавальну функцію фільтра другого порядку (загасання на -6 дб на октаву), підставляємо в формулу для H(s) n = 2:

Тоді згортка для фільтра буде виглядати так:

Хай нам задані частота зрізу w (на якій амплітуда сигналу буде -3 dB) і Fc — частота дискретизації (кількість семплів в секунду), в герцах.
У формулі потрібно використовувати денормированные частоти, тобто провести заміну (в смуговому фільтрі будуть дві частоти w1 і w2, що визначають смугу пропускання):

Якщо ми хочемо розраховувати НЧ-фільтр, то треба зробити перетворення — зробити заміну параметра s в передатної функції:

Для розрахунку інших типів фільтрів (ВЧ, смугового, режекторного) потрібно робити інші заміни. Вони розглянуті в книзі Цифрова обробка сигналів. Практичний підхід в частині 8.8.2 і далі в статті.
Далі, для переходу в z-площину робимо заміну:

Для аналітичних розрахунків я використовував пакет Mathematica.

Треба отримати чисельник і знаменник у вигляді поліномів від z. Наведемо складові знаменника H(z) до спільного знаменника. Для цього знайдемо найбільший загальний дільник (НОД, GCD) доданків знаменника і поділимо на нього чисельник і знаменник вихідної функції H(z).

Знайдемо коефіцієнти при степенях у отриманих поліномів, використовуючи функцію CoefficientList:

Якщо все робити чесно, то, за умовою, a0 повинен дорівнювати 1 — поділимо на всі коефіцієнти a0 (для кодинга будемо використовувати попередні формули без поділу):


Виведення формули фільтра ВЧ і смугового фільтра
Вивід формул для ВЧ-фільтра аналогічний НЧ-фільтр з іншим перетворенням:


Для виводу формул смугового фільтра застосовується перетворення:

Якщо робити заміну, то ступінь поліномів у чисельнику та знаменнику H(z) подвоїться (заміни є s^2), значить порядок фільтра збільшиться вдвічі. Тому спочатку використовуємо функцію H(s) для n = 1:


Програмування фільтра Баттервота за отриманими формулами
Фільтр буде мати 2 параметри: тип фільтра (НЧ, ВЧ, смуговий) і частота зрізу w. Для смугового фільтра будемо розглядати частоту зрізу частоту посередині смуги пропускання. Смугу пропускання ж визначимо як відрізок частот [w — w/4, w + w/4] (тут можна придумати більш складний і логічний тут логарифмічний закон, на ваш розсуд).
Нехай ми визначили коефіцієнти b0, b1, b2, a1 і a2 (a0 за умовою дорівнює 1) за розрахованими формулами. Алгоритм роботи фільтра зводиться до згортку, яка робиться послідовно для кожного семпла:

y(n) — це нове значення семпла, яке потрібно розрахувати. x(n) — поточне значення семпла, відповідно y(n-1) і y(n-2) — попередні 2 розрахованих семпла, а x(n-1) і x(n-2) — попередні вхідні значення семплів.
Потрібно організувати запам'ятовування попередніх семплом.
Не будемо мудрувати з циклічними буферами, зробимо просто і зрозуміло: два масиву з трьох елементів. Кожен раз будемо "проштовхувати" нові значення в цей масив, послідовно копіюючи більш старі значення семплів.
Отримуємо простий клас:
class BiquadConvolutionTable
{
public double B0, B1, B2, A1, A2;

private readonly double[] _x = new double[3];
private readonly double[] _y = new double[3];

public double Process(double s)
{
// "зсовуємо" попередні семпли
_x[2] = _x[1];
_x[1] = _x[0];
_x[0] = s;

_y[2] = _y[1];
_y[1] = _y[0];

// згортка
_y[0] = B0 * _x[0] + B1 * _x[1] + B2 * _x[2] - A1 * _y[1] - A2 * _y[2];

return _y[0];
}
}

Напишемо каркасний клас для фільтра (дивись архітектуру синтезатора першої статті).
Клас BiquadConvolutionTable працює з одним сигналом, тобто з одним каналом — моно. Тому нам потрібні дві BiquadConvolutionTable — для лівого і правого каналів.
Щоб коректно застосувати фільтр, потрібно послідовно, для всіх семплів вхідної послідовності застосувати функцію BiquadConvolutionTable.Process і заповнити результуючий масив семплів.
Розрахунком коефіцієнтів для BiquadConvolutionTable буде займатися функція CalculateCoefficients.
public enum EFilterPass
{
None,

LowPass,
HiPass,
BandPass
}

public class ButterworthFilter : SyntageAudioProcessorComponentWithparameters<AudioProcessor>, IProcessor
{
private readonly BiquadConvolutionTable _tablel;
private readonly BiquadConvolutionTable _tabler;

public EnumParameter<EFilterPass> FilterType { get; private set; }
public FrequencyParameter CutoffFrequency { get; private set; }

public ButterworthFilter(AudioProcessor audioProcessor) : base(audioProcessor)
{
_tablel = new BiquadConvolutionTable();
_tabler = new BiquadConvolutionTable();
}

public override IEnumerable<Parameter> CreateParameters(string parameterPrefix)
{
FilterType = new EnumParameter<EFilterPass>(parameterPrefix + "Pass", "Filter Type", "Filter", 'false');
CutoffFrequency = new FrequencyParameter(parameterPrefix + "Cutoff", "Filter Cutoff Frequency", "Cutoff");

return new List<Parameter> {FilterType, CutoffFrequency};
}

public void Process(IAudioStream stream)
{
if (FilterType.Value == EFilterPass.None)
return;

var count = Processor.CurrentStreamLenght;
var lc = stream.Channels[0];
var rc = stream.Channels[1];
for (int i = 0; i < count; ++i)
{
var cutoff = CutoffFrequency.Value;
CalculateCoefficients(cutoff);

var ls = _tablel.Process(lc.Samples[i]);
lc.Samples[i] = ls;

var rs = _tabler.Process(rc.Samples[i]);
rc.Samples[i] = rs;
}
}

private void CalculateCoefficients(double cutoff)
{
...
}
}

Функція CalculateCoefficients викликається щоразу в циклі — навіщо? У наступній статті я розповім про модуляцію (зміна в часі) параметрів, і тому, частота зрізу може змінюватися, а значить, потрібно перерозраховувати коефіцієнти. Звичайно, по-трушному, потрібно підписатися на зміну частоти зрізу і вже в оброблювачі розраховувати коефіцієнти. Але в цих статтях я оптимизами займатися не буду, мета — закодить фільтр.
Залишилося закодить функцію CalculateCoefficients по розрахованим формул для коефіцієнтів.
Згадаймо, що потрібно використовувати денормированные частоти, тобто зробити заміну:

Списуємо все формули для коефіцієнтів b0, b1, b2, a0, a1, a2. Після розрахунків потрібно поділити всі коефіцієнти на a0, щоб a0 стало дорівнює 1.
private double TransformFrequency(double w)
{
return Math.Tan(Math.PI * w / Processor.SampleRate);
}

private void CalculateCoefficients(double cutoff)
{
double b0, b1, b2, a0, a1, a2;

switch (FilterType.Value)
{
case EFilterPass.LowPass:
{
var w = TransformFrequency(cutoff);

a0 = 1 + Math.Sqrt(2) * w + w * w;
a1 = -2 + 2 * w * w;
a2 = 1 - Math.Sqrt(2) * w + w * w;

b0 = w * w;
b1 = 2 * w * w;
b2 = w * w;
}

break;

case EFilterPass.HiPass:
{
var w = TransformFrequency(cutoff);

a0 = 1 + Math.Sqrt(2) * w + w * w;
a1 = -2 + 2 * w * w;
a2 = 1 - Math.Sqrt(2) * w + w * w;

b0 = 1;
b1 = -2;
b2 = 1;
}

break;

case EFilterPass.BandPass:
{
var w = cutoff;
var d = w / 4;

// визначимо смугу фільтра [w * 3 / 4, w * 5 / 4]
var w1 = Math.Max(w - d, CutoffFrequency.Min);
var w2 = Math.Min(w + d, CutoffFrequency.Max);
w1 = TransformFrequency(w1);
w2 = TransformFrequency(w2);

var w0Sqr = w2 * w1; // w0^2
var wd = w2 - w1; // W

a0 = -1 - wd - w0Sqr;
a1 = 2 - 2 * w0Sqr;
a2 = -1 + wd - w0Sqr;
b0 = -wd;
b1 = 0;
b2 = wd;
}

break;

default:
throw new ArgumentOutOfRangeException();
}

_tablel.B0 = _tabler.B0 = b0 / a0;
_tablel.B1 = _tabler.B1 = b1 / a0;
_tablel.B2 = _tabler.B2 = b2 / a0;
_tablel.A1 = _tabler.A1 = a1 / a0;
_tablel.A2 = _tabler.A2 = a2 / a0;
}

Повний код класу ButterworthFilter
public enum EFilterPass
{
None,

LowPass,
HiPass,
BandPass
}

public class ButterworthFilter : SyntageAudioProcessorComponentWithparameters<AudioProcessor>, IProcessor
{
private class BiquadConvolutionTable
{
public double B0, B1, B2, A1, A2;

private readonly double[] _x = new double[3];
private readonly double[] _y = new double[3];

public double Process(double s)
{
// "зсовуємо" попередні семпли
_x[2] = _x[1];
_x[1] = _x[0];
_x[0] = s;

_y[2] = _y[1];
_y[1] = _y[0];

// згортка
_y[0] = B0 * _x[0] + B1 * _x[1] + B2 * _x[2] - A1 * _y[1] - A2 * _y[2];

return _y[0];
}
}

private readonly BiquadConvolutionTable _tablel;
private readonly BiquadConvolutionTable _tabler;

public EnumParameter<EFilterPass> FilterType { get; private set; }
public FrequencyParameter CutoffFrequency { get; private set; }

public ButterworthFilter(AudioProcessor audioProcessor) : base(audioProcessor)
{
_tablel = new BiquadConvolutionTable();
_tabler = new BiquadConvolutionTable();
}

public override IEnumerable<Parameter> CreateParameters(string parameterPrefix)
{
FilterType = new EnumParameter<EFilterPass>(parameterPrefix + "Pass", "Filter Type", "Filter", 'false');
CutoffFrequency = new FrequencyParameter(parameterPrefix + "Cutoff", "Filter Cutoff Frequency", "Cutoff");

return new List<Parameter> {FilterType, CutoffFrequency};
}

public void Process(IAudioStream stream)
{
if (FilterType.Value == EFilterPass.None)
return;

var count = Processor.CurrentStreamLenght;
var lc = stream.Channels[0];
var rc = stream.Channels[1];
for (int i = 0; i < count; ++i)
{
var cutoff = CutoffFrequency.Value;
CalculateCoefficients(cutoff);

var ls = _tablel.Process(lc.Samples[i]);
lc.Samples[i] = ls;

var rs = _tabler.Process(rc.Samples[i]);
rc.Samples[i] = rs;
}
}

private double TransformFrequency(double w)
{
return Math.Tan(Math.PI * w / Processor.SampleRate);
}

private void CalculateCoefficients(double cutoff)
{
double b0, b1, b2, a0, a1, a2;

switch (FilterType.Value)
{
case EFilterPass.LowPass:
{
var w = TransformFrequency(cutoff);

a0 = 1 + Math.Sqrt(2) * w + w * w;
a1 = -2 + 2 * w * w;
a2 = 1 - Math.Sqrt(2) * w + w * w;

b0 = w * w;
b1 = 2 * w * w;
b2 = w * w;
}

break;

case EFilterPass.HiPass:
{
var w = TransformFrequency(cutoff);

a0 = 1 + Math.Sqrt(2) * w + w * w;
a1 = -2 + 2 * w * w;
a2 = 1 - Math.Sqrt(2) * w + w * w;

b0 = 1;
b1 = -2;
b2 = 1;
}

break;

case EFilterPass.BandPass:
{
var w = cutoff;
var d = w / 4;

// визначимо смугу фільтра [w * 3 / 4, w * 5 / 4]
var w1 = Math.Max(w - d, CutoffFrequency.Min);
var w2 = Math.Min(w + d, CutoffFrequency.Max);
w1 = TransformFrequency(w1);
w2 = TransformFrequency(w2);

var w0Sqr = w2 * w1; // w0^2
var wd = w2 - w1; // W

a0 = -1 - wd - w0Sqr;
a1 = 2 - 2 * w0Sqr;
a2 = -1 + wd - w0Sqr;
b0 = -wd;
b1 = 0;
b2 = wd;
}

break;

default:
throw new ArgumentOutOfRangeException();
}

_tablel.B0 = _tabler.B0 = b0 / a0;
_tablel.B1 = _tabler.B1 = b1 / a0;
_tablel.B2 = _tabler.B2 = b2 / a0;
_tablel.A1 = _tabler.A1 = a1 / a0;
_tablel.A2 = _tabler.A2 = a2 / a0;
}
}


Смуговий еквалайзер в бібліотеці NAudio
Для роботи з звуком, звуковими файлами різних форматів .NET є хороша бібліотека NAudio.
Вона містить простір імен NAudio.Dsp з функціоналом для фільтрації, згортки, гейта, огинаючої, ШПФ та інших цікавих штук.
Розглянемо клас Equalizer (з прикладів, простір імен NAudioWpfDemo.EqualizationDemo), що дозволяє виробляти эквализацию сигналу. Клас реалізує ISampleProvider, який у функції Read(float[] buffer, int offset, int count) обробляє (змінює) масив семплів buffer.
Конструктор приймає масив структури EqualizerBand, які описують "смуги" еквалайзера:
class EqualizerBand
{
public float Frequency { get; set; }
public float Gain { get; set; }
public float Bandwidth { get; set; }
}

Тут Frequency — центральна частота смуги з параметром Q (Bandwidth, добротність фільтру), c посиленням Gain dB.
Якщо заглянути в реалізацію, то кожній смузі EqualizerBand відповідає клас BiQuadFilter який використовується як смуговий (Peaking) фільтр. Всі фільтри змінюють сигнал використовуються послідовно.
Клас EqualizerBand є реалізацією фільтра Баттервота, з великим виборів типів фільтрів і параметрами. Якщо подивитися реалізацію, можна побачити схожі формули і коефіцієнти.
Приклад використання класу Equalizer ви знайдете в проекті NAudioWpfDemo в класі EqualizationDemoViewModel.

Програми для розрахунку фільтрів
Прабатьком цифрових фільтрів були аналогові фільтри. Теорія для аналогових схем і аналогової обробки сигналів в подальшому переросла в теорію цифрової обробки сигналів.
Для розглянутого фільтра Баттервота і розрахованих формул для коефіцієнтів можна зібрати аналогову схему.
Є багато програм для розрахунку і побудови схем, параметрів елементів для них, згорткових коефіцієнтів різних фільтрів.
Можете погуглити "filter calculation software".

Iowa Hills Software RF Filter Designer
У наступній статті я розповім про ефект delay, distortion і модуляцію параметрів.
Всім добра!
Удачі в програмуванні!


Список літератури
Не забувайте дивитися списки статей і книг в попередніх статтях.
  1. Digital Signal Processing: A Practical Approach. Російська версія — Цифрова обробка сигналів. Практичний підхід.
  2. Хабр-стаття "Простими словами про перетворення Фур'є"
  3. Фільтр Баттерворта на вікі
  4. Github-репозитарій з кодом для розрахунку фільтрів Баттервота
  5. DISCRETE SIGNALS AND SYSTEMS, Chapter 7. FIR and IIR Filters
  6. How does a low-pass filter programmatically work? (dsp.stackexchange.com/)
  7. Designing Digital Баттерворта and Чебишева Filters
  8. Audio EQ Cookbook
  9. Iowa Hills Software — Digital and Analog Filters
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.