Однополюсный фильтр нижних частот

Источники голоса дают пилу, прямоугольник и шум — сигналы, в которых много высоких частот. Огибающая ADSR меняет громкость ноты во времени, но не её тембр (timbre)?тембр (timbre) — Свойство слухового ощущения, по которому различают два звука одинаковой высоты и громкости. Зависит от состава спектра и его изменения во времени.. Тембр меняет фильтр нижних частот (low-pass filter, LPF)?фильтр нижних частот (low-pass filter, LPF) — Фильтр, который пропускает низкие частоты и ослабляет высокие. Перед дискретизацией помогает ограничить спектр входного сигнала.: он пропускает низкие частоты и ослабляет высокие. Здесь разбирается его простейший цифровой вариант — однополюсный фильтр (one-pole filter)?однополюсный фильтр (one-pole filter) — Простейший рекурсивный фильтр: новый выход вычисляется из входного отсчёта и одного предыдущего выхода, y[n] = y[n−1] + a·(x[n] − y[n−1]). При коэффициенте a от 0 до 1 работает как фильтр нижних частот., которому на каждом отсчёте (sample)?отсчёт (sample) — Значение сигнала в один момент дискретизации, одно число. В звукозаписи его также называют сэмплом. нужен только собственный предыдущий выход. Третья часть одноголосного синтезатора: формула фильтра, частота среза и коэффициент, реакция на скачок, звучание и класс на C++20.

Голос и фильтр

На вход фильтра приходит сумма $x[n]$ из микшера, выход $y[n]$ идёт в VCA (voltage-controlled amplifier)?VCA (voltage-controlled amplifier) — Усилитель голоса синтезатора, коэффициент которого задаёт управляющий сигнал. В цифровом голосе — умножение каждого отсчёта на значение огибающей.. В аналоговом синтезаторе этот блок называют VCF (voltage-controlled filter): его частоту среза можно менять управляющим напряжением, например от огибающей или LFO (low-frequency oscillator, LFO)?низкочастотный генератор (low-frequency oscillator, LFO) — Генератор медленно повторяющегося управляющего сигнала, которым изменяют параметры звука, например усиление, частоту или панораму. Общей строгой верхней границы частоты у LFO нет.. Здесь частота среза постоянна в течение ноты; как менять её на ходу, сказано в разделе о коде.

Субтрактивный синтез

Субтрактивный синтез (subtractive synthesis)?субтрактивный синтез (subtractive synthesis) — Способ синтеза, при котором сначала создают сигнал с большим числом гармоник, например пилу или шум, а затем фильтром ослабляют часть его частот. сначала создаёт сигнал с большим числом составляющих, а потом фильтрами ослабляет часть из них (Ableton, Drift: Subtractive Synthesis). Название — от английского subtract, «вычитать»: звук получается тем, что из готового набора частот убирают лишнее.

Фильтр стоит после микшера, а не на каждом источнике, потому что он линеен: его отклик на сумму сигналов равен сумме откликов на каждый из них (Sundararajan, §2.3.2). Поэтому один фильтр после микшера даёт тот же выход, что четыре одинаковых фильтра на четырёх источниках, и выполняет вчетверо меньше вычислений. Частота среза у всех источников общая, и огибающей или LFO нужно управлять одним параметром. Тот же порядок блоков — микшер источников, затем фильтр и усилитель — у Roland SH-101 (статья об осцилляторах).

Фильтр как преобразование с состоянием

Фильтр получает последовательность отсчётов $x[0], x[1], x[2], \ldots$ и выдаёт последовательность $y[0], y[1], y[2], \ldots$ с той же частотой дискретизации (sample rate)?частота дискретизации (sample rate) — Число отсчётов в каждом канале за одну секунду, то же число кадров аудиоданных в секунду. Измеряется в герцах: 44 100 Гц означает 44 100 отсчётов на канал в секунду.. VCA вычисляет $y[n]$ только по $x[n]$. Фильтру нужна память: однополюсный фильтр хранит одно число — свой предыдущий выход $y[n-1]$. Это его состояние. Новый отсчёт вычисляется так:

Скобка — расстояние от прошлого выхода до нового входа. За один отсчёт выход проходит долю $a$ этого пути. Пусть $a = 0{,}25$, $y[-1] = 0$, а вход всё время равен 1. Тогда $y[0] = 0{,}25$, $y[1] = 0{,}25 + 0{,}25 \cdot 0{,}75 = 0{,}4375$, $y[2] = 0{,}578125$. Выход догоняет вход, но не перескакивает через него: каждый раз остаётся три четверти прежнего расстояния.

Если раскрыть скобки, получится эквивалентная запись:

Это рекуррентная формула (recurrence relation)?рекуррентная формула (recurrence relation) — Правило вычисления очередного значения последовательности через предыдущие значения. Для начала вычислений нужно задать начальное состояние. из математического справочника: новое значение — взвешенное среднее прошлого выхода и нового входа, веса $1 - a$ и $a$ в сумме дают 1. В справочнике вес прошлого обозначен $\alpha = 1 - a$, и той же формулой сглаживается подавление усиления в компрессоре (dynamic range compressor)?компрессор (dynamic range compressor) — Устройство или программа, которая ослабляет сигнал, когда его уровень выше порога, и тем сужает динамический диапазон.. Там на вход подаётся требуемое подавление в децибелах, здесь — звуковые отсчёты; вычисление одно и то же. Экспоненциальные участки огибающей ADSR — тоже такой фильтр, вход которого постоянен и равен цели участка.

Название «однополюсный» (one-pole) пришло из анализа фильтров через z-преобразование: у такого фильтра один полюс, и он равен $1 - a$ (J. O. Smith, One-Pole). Для расчётов в этой статье z-преобразование не нужно. В книгах встречается и название single pole recursive filter: рекурсивным (recursive) называют фильтр, который использует собственные прошлые выходы (S. W. Smith, гл. 19).

Начальное состояние и reset

Для первого отсчёта нужно $y[-1]$. Обычно берут 0: фильтр начинает с тишины. Метод reset() в коде присваивает состоянию начальное значение. Его вызывают перед новым рендером и перед нотой, если голос молчал: иначе первая нота начнётся с остатка прошлой. Сбрасывать состояние во время звучания нельзя. Если $y[n-1]$ скачком станет нулём при ненулевом входе, выход прыгнет к нулю и только потом вернётся к сигналу — это разрыв, который слышен как щелчок.

Частота среза и коэффициент

Коэффициент $a$ неудобно задавать напрямую: одно и то же значение при разной частоте дискретизации даёт разный звук. Удобнее задать частоту среза (cutoff frequency)?частота среза (cutoff frequency) — Частота, выше которой фильтр нижних частот заметно ослабляет сигнал. У однополюсного фильтра при частоте среза много ниже fs/2 синусоида этой частоты ослабляется примерно на 3 дБ. $f_c$ в герцах — границу, выше которой фильтр заметно ослабляет сигнал. Коэффициент из неё вычисляется через экспоненту (exponential function)?экспонента (exponential function) — Функция eᵘ, также записываемая exp(u), где e ≈ 2,71828. Обратна натуральному логарифму: exp(ln(a)) = a при положительном a. При отрицательном аргументе принимает значения между 0 и 1. (S. W. Smith, гл. 19, с. 324):

У S. W. Smith в той же формуле буквой $x = e^{-2\pi f_c}$ обозначен вес прошлого выхода $1 - a$, а $f_c$ записана долей частоты дискретизации, от 0 до 0,5. Связанная величина — постоянная времени (time constant)?постоянная времени (time constant) — Время τ, за которое экспоненциальный процесс проходит 1 − 1/e ≈ 63,2 % пути к установившемуся значению; у однополюсного фильтра τ = 1/(2π·fc).:

Например, при $f_c = 1000$ Гц и $f_s = 48\,000$ Гц показатель равен $-2\pi \cdot 1000/48\,000 \approx -0{,}1309$, $e^{-0{,}1309} \approx 0{,}8773$, и $a \approx 0{,}1227$: каждый отсчёт выход проходит около 12 % оставшегося пути к входу. Постоянная времени $1/(2\pi \cdot 1000) \approx 0{,}159$ мс — 7,64 отсчёта.

Если $f_c$ много меньше $f_s/2$, синусоида (sine)?синус (sine) — Вертикальная координата точки на единичной окружности под заданным углом от положительного направления горизонтальной оси. Для острого угла прямоугольного треугольника равен отношению противолежащего катета к гипотенузе. частоты $f_c$ проходит фильтр ослабленной примерно на 3 дБ: её амплитуда (amplitude)?амплитуда (amplitude) — Наибольший модуль отклонения колеблющейся величины от равновесного значения. Для синусоиды обозначается A; единицы совпадают с единицами сигнала. умножается на $1/\sqrt{2} \approx 0{,}707$. Так S. W. Smith и определяет $f_c$ в формуле (3) — частота среза по уровню −3 дБ (S. W. Smith, гл. 19, с. 324). Условие важно: в таблице 1 при 200 и 1000 Гц ослабление совпадает с 3 дБ до сотых, а при 5000 Гц оно уже 2,86 дБ, и уровень −3 дБ приходится на 5190 Гц. Чем ближе $f_c$ к $f_s/2$, тем сильнее расхождение. Синусоиды намного ниже $f_c$ проходят почти без изменений, выше $f_c$ — ослабляются тем сильнее, чем выше их частота.

Связь с компрессором видна и в числах. Его коэффициент $\alpha = e^{-\ln 9/(f_s t)}$ (формула (2) статьи о компрессоре) равен $1 - a$ при $\tau = t/\ln 9 \approx t/2{,}2$. Атака 10 мс соответствует $\tau \approx 4{,}55$ мс и $f_c \approx 35$ Гц. Сглаживание подавления — фильтр нижних частот с очень низкой частотой среза: медленные изменения уровня он пропускает, а колебания с частотой звука убирает.

Почему фильтр устойчив

Фильтр устойчив (stable), если ограниченный вход всегда даёт ограниченный выход (Sundararajan, §2.5). Для формулы (2) при $0 < a \le 1$ это видно без вычислений: $y[n]$ — взвешенное среднее $y[n-1]$ и $x[n]$ с неотрицательными весами, а такое среднее лежит между двумя усредняемыми числами. Если вход не выходит за ±1 и начальное состояние тоже, выход за ±1 не выйдет никогда. Формула (3) всегда даёт такое $a$: при любой $f_c > 0$ экспонента лежит между 0 и 1.

На краях интервала $a = 1$ означает $y[n] = x[n]$: фильтра нет. При $a$, близком к 0, выход почти не меняется. Если подставить $a > 1$, вес прошлого $1 - a$ станет отрицательным: выход начнёт перескакивать через вход с переменой знака, а при $a > 2$ колебания будут расти без предела. Общее условие устойчивости одного полюса — модуль веса прошлого выхода меньше 1 (J. O. Smith, One-Pole); фильтром нижних частот такая формула работает только при положительном весе, то есть при $0 < a \le 1$.

Реакция на скачок

Подадим на вход скачок: до начала вход и состояние равны 0, затем вход всё время равен 1. На каждом отсчёте расстояние до входа умножается на $1 - a$, поэтому после $N$ отсчётов от скачка выход равен

Вторая запись следует из формулы (4): $1 - a = e^{-1/(f_s\tau)}$. Отсчёты выхода лежат на экспоненте с постоянной времени $\tau$. При $f_c = 1000$ Гц $f_s\tau = 7{,}64$ отсчёта: уровень 63,2 % выход проходит между 7-м и 8-м отсчётом, 95 % — за $3\tau$, около 23 отсчётов или 0,48 мс, 99 % — за $4{,}6\tau$. Так же отвечает на скачок напряжения RC-цепь: S. W. Smith сравнивает переходные характеристики цифрового фильтра и RC-фильтра нижних частот (S. W. Smith, гл. 19, рис. 19-2), а Sundararajan выводит разностное уравнение первого порядка из уравнения такой цепи (Sundararajan, §2.1).

Сдвиньте ползунок к 16 кГц: постоянная времени станет меньше половины отсчёта, и выход дойдёт почти до 1 за три отсчёта. При 24,8 Гц $\tau = 6{,}4$ мс — около 308 отсчётов, и график показывает уже не отдельные точки, а кривую. Переходная характеристика не перелетает через 1: однополюсный фильтр приближается к входу монотонно.

Видимое и слышимое сглаживание

Таблица 2 показывает, что фильтр делает со скачком пилы. Для удобства счёта взято $a = 0{,}5$ — при 48 кГц это $f_c \approx 5300$ Гц — и состояние $y[-1] = 0{,}7$. Пила растёт на 0,1 за отсчёт, а между отсчётами 2 и 3 падает от 1 к −1.

Вход меняется на 2 за один отсчёт. Выход за тот же отсчёт меняется на 0,956, за следующие два — ещё на 0,43 и 0,16: скачок превратился в экспоненциальный переход. На ровном подъёме выход отстаёт от пилы почти на постоянную величину, здесь около 0,1. Чем ниже частота среза, тем длиннее переход и тем сильнее скруглены углы формы. Скачок пилы и даёт её высокие гармоники, поэтому скруглённая пила звучит темнее.

Визуализация считает фильтр тем же кодом, что приведён в разделе C++, при 48 кГц. Петля длиной 2 с проходит через фильтр дважды: первый проход приводит состояние к установившемуся, второй даёт выход, поэтому график и звук начинаются без переходного процесса. Опустите частоту среза ниже частоты тона: ослабляется уже основная частота, и пила становится похожа на синусоиду с меньшей амплитудой. Переключите «До фильтра» и «После фильтра» на белом шуме, а затем на розовом: разница в громкости и тембре у белого намного больше.

Ограничения первого фильтра

У самого простого фильтра нижних частот есть ограничения, которые заметны в синтезаторе.

Спад пологий. Далеко выше $f_c$ и далеко ниже $f_s/2$ каждое удвоение частоты — каждая октава (octave)?октава (octave) — Музыкальный интервал между нотами одного названия в соседних регистрах. Соответствует отношению частот 2:1. — ослабляет синусоиду ещё примерно на 6 дБ, то есть вдвое по амплитуде. При $f_c = 200$ Гц и 48 кГц синусоида 400 Гц ослабляется на 7,0 дБ, 800 Гц — на 12,3 дБ, 1600 Гц — на 18,1 дБ. S. W. Smith подытоживает: однополюсные фильтры хорошо работают во временной области и плохо в частотной — соседние полосы частот они разделяют слабо (S. W. Smith, гл. 19, с. 325). В синтезаторах ставят и более крутые фильтры: например, у Ableton Drift два варианта — 12 и 24 дБ на октаву (Ableton, Drift: Filter Section).

Нет резонанса. У фильтров синтезаторов есть регулятор резонанса (resonance): он поднимает частоты возле $f_c$, и при движении частоты среза эта узкая полоса слышна как отдельный тон, который скользит вверх или вниз. Однополюсному фильтру такой подъём взять неоткуда: его выход приближается к входу монотонно, как на визуализации 2. Резонанс дают фильтры второго порядка с двумя состояниями или обратная связь: в лестничном фильтре Moog четыре фильтра первого порядка охвачены общей петлёй, и она усиливает частоты возле среза (Pirkle, гл. 12).

Это не модель аналоговой схемы. С формулой (3) отсчёты переходной характеристики в точности равны значениям реакции RC-цепи с $RC = \tau$ в моменты отсчётов. Частотная характеристика совпадает с RC-цепью только при $f_c$ много ниже $f_s/2$: у RC-цепи ослабление растёт с частотой без предела, а у цифрового фильтра с $f_c = 1000$ Гц на $f_s/2 = 24$ кГц оно всего 23,7 дБ, с $f_c = 5000$ Гц — 10 дБ. Фильтры виртуального аналога (virtual analog) строят иначе, один из них разобран ниже.

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

Код на C++

Класс LowPassFilter устроен как осциллятор и огибающая серии: параметр — коэффициент a — хранится отдельно от состояния previous, то есть $y[n-1]$. prepare задаёт частоту дискретизации и сбрасывает состояние, setCutoff пересчитывает коэффициент по формуле (3), reset присваивает состоянию начальное значение, processSample вычисляет один отсчёт по формуле (1). Используется только стандартная библиотека C++20.

// Однополюсный фильтр нижних частот: коэффициент a отдельно от состояния y[n-1].
#include <cassert>
#include <cmath>
#include <numbers>

class LowPassFilter
{
public:
    void prepare (double newSampleRate) noexcept
    {
        assert (std::isfinite (newSampleRate) && newSampleRate > 0.0);
        sampleRate = newSampleRate;
        setCutoff (cutoff);
        reset();
    }

    // Параметр: a = 1 - e^(-2π·fc/fs), 0 < fc < fs/2.
    void setCutoff (double hz) noexcept
    {
        assert (hz > 0.0 && hz < sampleRate / 2.0);
        cutoff = hz;
        a = 1.0 - std::exp (-2.0 * std::numbers::pi * hz / sampleRate);
    }

    // Состояние: значение y[n-1], с которого начнётся следующий отсчёт.
    void reset (double value = 0.0) noexcept { previous = value; }

    double processSample (double x) noexcept
    {
        previous += a * (x - previous); // y[n] = y[n-1] + a·(x[n] - y[n-1])
        return previous;
    }

private:
    double sampleRate = 48000.0, cutoff = 1000.0;
    double a = 0.0;
    double previous = 0.0;
};

setCutoff вызывает std::exp один раз на смену параметра, а processSample выполняет на отсчёт одно вычитание, одно умножение и одно сложение. std::numbers::pi появилась в C++20 в заголовке <numbers>. prepare пересчитывает коэффициент: частота среза хранится в герцах, а $a$ зависит и от частоты дискретизации. Код собран clang++ -std=c++20 -Wall -Wextra -Wconversion без предупреждений. При $f_c = 1000$ Гц и скачке входа он даёт те же отсчёты, что визуализации, до 15-го знака.

Фильтр в голосе

// Голос из трёх статей серии: микшер → фильтр нижних частот → VCA.
void prepareVoice (VoiceSources& voice, LowPassFilter& filter, AdsrEnvelope& envelope,
                   double sampleRate) noexcept
{
    voice.prepare (sampleRate);
    filter.prepare (sampleRate);
    filter.setCutoff (2000.0);
    envelope.prepare (sampleRate);
}

float nextVoiceSample (VoiceSources& voice, const VoiceLevels& levels,
                       LowPassFilter& filter, AdsrEnvelope& envelope) noexcept
{
    const double mixed = voice.nextSample (levels);
    return static_cast<float> (filter.processSample (mixed) * envelope.nextSample());
}

Состояние фильтра сбрасывается там же, где фазы осцилляторов: в startNote, когда голос молчал (envelope.isActive() вернул false), вместе с voice.noteOn вызывается filter.reset(). У звучащего голоса состояние сохраняется, иначе новая нота началась бы со скачка.

Сравнение с juce::dsp::FirstOrderTPTFilter

// Готовый фильтр первого порядка JUCE: один канал, нижние частоты, 1000 Гц.
#include <juce_dsp/juce_dsp.h>

juce::dsp::FirstOrderTPTFilter<float> tpt;

void prepareFilter (double sampleRate, int maximumBlockSize)
{
    tpt.prepare ({ sampleRate, static_cast<juce::uint32> (maximumBlockSize), 1 });
    tpt.setType (juce::dsp::FirstOrderTPTFilterType::lowpass);
    tpt.setCutoffFrequency (1000.0f);
}

float filterSample (float x) { return tpt.processSample (0, x); }

Совпадает устройство. Состояние — одно число на канал, prepare принимает juce::dsp::ProcessSpec с частотой дискретизации, наибольшим размером блока и числом каналов и сбрасывает состояние, reset() обнуляет его, а reset (value) присваивает значение. setCutoffFrequency принимает частоту в герцах и проверяет, что она лежит между 0 и $f_s/2$. processSample (channel, x) обрабатывает один отсчёт, process — целый блок.

Различается расчёт. TPT (topology-preserving transform) — метод виртуального аналога Вадима Завалишина: интегратор RC-цепи заменяют цифровым, который считает площадь под сигналом по формуле трапеций (Pirkle, §12.4.1 и 12.6.1). Это билинейное преобразование (bilinear transform) аналогового фильтра (Sundararajan, §8.2). В коде JUCE коэффициенты и шаг такие:

// JUCE 9.0.3, FirstOrderTPTFilter: update() и processSample() для нижних частот, сокращённо.
auto g = std::tan (pi * cutoffFrequency / sampleRate); // частота с предыскажением
G = g / (1 + g);

auto v = G * (inputValue - s); // s — состояние канала
auto y = v + s;
s = y + v;
return y;

• Частота среза точная. Тангенс заранее учитывает искажение оси частот при билинейном преобразовании, поэтому синусоида частоты $f_c$ ослабляется ровно на 3,01 дБ при любой $f_c$ ниже $f_s/2$. У формулы (3) при 5000 Гц — 2,86 дБ.
• На $f_s/2$ выход равен нулю: у TPT-фильтра есть ноль передаточной функции на частоте Найквиста. При $f_c = 1000$ Гц синусоида 20 кГц ослабляется на 35,1 дБ, у LowPassFilter — на 23,4 дБ.
• Реакция на скачок другая уже на первом отсчёте: выход TPT-фильтра равен $G$, а не $a$, при 1000 Гц — 0,0615 против 0,1227. Полюс TPT-фильтра равен $(1 - g)/(1 + g)$, у LowPassFilter — $1 - a$. При $f_c = 1000$ Гц это 0,8770 против 0,8773: дальше обе реакции идут по почти одинаковой экспоненте, и TPT-фильтр отстаёт примерно на половину отсчёта. При 10 кГц полюсы различаются сильнее: 0,132 против 0,270.
• Та же структура даёт фильтр верхних частот (x − y) и фазовый фильтр (2y − x) — выбор через setType.
• Документация JUCE отмечает, что такой фильтр можно модулировать с высокой скоростью без слышимых артефактов, но изменения частоты среза всё равно советует сглаживать, например juce::SmoothedValue. У LowPassFilter смена a тоже не даёт скачка выхода: состояние сохраняется, меняется только доля шага. Но скачок частоты среза мгновенно меняет тембр, поэтому при модуляции её тоже меняют плавно.
• snapToZero() заменяет нулём состояния, которые после долгого затухания стали денормализованными числами (denormals) — значениями у нижней границы точности float. process вызывает его сам, а при обработке по одному отсчёту его вызывают после блока.

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

Источники

• Smith S. W. The Scientist and Engineer’s Guide to Digital Signal Processing. Гл. 19: Recursive Filters, с. 319–325. в списке литературы → — рекурсивный фильтр, однополюсный фильтр нижних частот $a_0 = 1 - x$, $b_1 = x$, связь $x$ с постоянной времени и частотой среза по уровню −3 дБ, сравнение с RC-цепью, неустойчивость при $x > 1$, слабое разделение частот.
• Smith J. O. III. Introduction to Digital Filters with Audio Applications. One-Pole. в списке литературы → — разностное уравнение и передаточная функция фильтра с одним полюсом, условие устойчивости, нижние частоты при полюсе возле 1, нормировка усиления $1 - |a_1|$.
• Smith J. O. III. Mathematics of the Discrete Fourier Transform (DFT) with Audio Applications. Exponentials. в списке литературы → — постоянная времени как время затухания в $e$ раз.
• Sundararajan D. Digital Signal Processing: An Introduction. §2.1, с. 40–43; §2.3.2, с. 53; §2.5, с. 57; §8.2, с. 299–301. в списке литературы → — разностное уравнение первого порядка из уравнения RC-цепи, начальное состояние и отклик на скачок, линейность, устойчивость по ограниченному выходу, билинейное преобразование и искажение оси частот.
• Pirkle W. C. Designing Software Synthesizer Plugins in C++: With Audio DSP. Гл. 12: Synthesizer Filters, с. 173–186. в списке литературы → — цифровой интегратор по формуле трапеций, предыскажение частоты среза, фильтр первого порядка виртуального аналога с коэффициентом $g/(1 + g)$, ноль на частоте Найквиста, резонанс лестничного фильтра Moog.
• JUCE 9.0.3. juce_FirstOrderTPTFilter.h и juce_FirstOrderTPTFilter.cpp — методы класса, расчёт G через тангенс, шаг фильтра, типы и snapToZero; ProcessSpec — поля спецификации.
• Ableton. Live 12 Manual, Drift: Subtractive Synthesis и Filter Section — определение субтрактивного синтеза, фильтры нижних частот 12 и 24 дБ на октаву, резонанс.