Алиасинг в нелинейной обработке: интермодуляция и передискретизация

По материалам статьи Боба Каца (Bob Katz), инженера мастеринга студии Digital Domain, «Analog Saved My Life: The Good, The Bad and The Ugly» на сайте digido.com. Условия и результаты измерений взяты оттуда и сверены с учебниками по цифровой обработке; оценки звучания — позиция Каца как практика, а не стандарт.

Аналоговый компрессор (dynamic range compressor)?компрессор (dynamic range compressor) — Устройство или программа, которая ослабляет сигнал, когда его уровень выше порога, и тем сужает динамический диапазон. с ламповым каскадом и его цифровая модель могут давать одинаковые гармоники, но в записи цифровой модели появляются ещё и составляющие, не связанные с нотами. Их источник — граница $f_s/2$, половина частоты дискретизации (sample rate)?частота дискретизации (sample rate) — Число отсчётов в каждом канале за одну секунду, то же число кадров аудиоданных в секунду. Измеряется в герцах: 44 100 Гц означает 44 100 отсчётов на канал в секунду.: гармоники выше неё не исчезают, а отражаются в слышимую полосу. Ниже разобрано, куда попадают такие отражения, как отличить их от интермодуляционных искажений (intermodulation distortion)?интермодуляционные искажения (intermodulation distortion) — Составляющие на суммах и разностях частот входных тонов, |m·f₁ + k·f₂| при ненулевых m и k, которые нелинейная обработка добавляет к сигналу из нескольких тонов. В отличие от гармоник, они обычно не кратны частоте ни одного из тонов., которые даёт и аналоговая схема, что показал тест двумя тонами 19 и 20 кГц и как отражения ослабляют передискретизацией (oversampling)?передискретизация (oversampling) — Работа с частотой дискретизации выше необходимой для полезной полосы. В АЦП уменьшает долю шума квантования в полезной полосе и смягчает требования к входному фильтру, в ЦАПе отодвигает спектральные копии, вокруг нелинейной обработки оставляет место для новых составляющих, чтобы они не отражались в полезную полосу. Новых сведений об исходном звуке не добавляет. и методом антипроизводной. В конце — код с juce::dsp::Oversampling.

Новые частоты на выходе кривой

Сумма $f_1 + f_2$ и разность $f_2 - f_1$ — составляющие второго порядка, $2f_1 - f_2$ и $2f_2 - f_1$ — третьего. Симметричная кривая, например гиперболический тангенс (hyperbolic tangent, tanh)?гиперболический тангенс (hyperbolic tangent, tanh) — Функция tanh u = (eᵘ − e⁻ᵘ)/(eᵘ + e⁻ᵘ). Около нуля почти совпадает с u, при росте |u| плавно приближается к ±1, не достигая их. $\tanh$, даёт только нечётные порядки, асимметричная — все. Чем выше порядок, тем слабее составляющая, но у кривой, которая не записывается многочленом, порядки не кончаются.

В аналоговой схеме высокие порядки ничем не мешают: их ослабляют ограниченная полоса самих элементов и фильтры, а то, что остаётся выше 20 кГц, не слышно (Reiss, McPherson, с. 197). В цифровой обработке кривая получает только отсчёты, и порядок «сначала дискретизация, потом кривая» даёт те же числа, что «сначала кривая, потом дискретизация». Значит, цифровой выход — это отсчёты непрерывного сигнала $f\bigl(G\,x(t)\bigr)$, полоса которого не ограничена $f_s/2$ (Zölzer, §10.4, с. 350). Все его составляющие выше частоты Найквиста (Nyquist frequency)?частота Найквиста (Nyquist frequency) — Половина частоты дискретизации, fₛ/2. В принятой низкочастотной модели точное восстановление по отсчётам требует, чтобы все частоты исходного сигнала были строго ниже этой границы. переносятся наложением спектров (aliasing)?наложение спектров (aliasing) — Наложение частотных составляющих при дискретизации: разные непрерывные сигналы дают одинаковые отсчёты. Частота выше границы Найквиста может проявиться как более низкая. в полосу от 0 до $f_s/2$ по формуле статьи о дискретизации: частота $f$ превращается в $|f - k f_s|$ с ближайшим к $f/f_s$ целым $k$.

Отражённая составляющая не кратна частоте ноты и не совпадает с сочетаниями входных частот. Отфильтровать её после кривой нельзя: она уже лежит в слышимой полосе среди полезных составляющих (Reiss, McPherson, с. 197–198). На слух алиасинг описывают как размытость, лёгкую пелену и резкость звучания.

Тест двумя тонами

Для проверки на вход подают два тона у верхнего края полосы, например 19 и 20 кГц. Их гармоники и суммы лежат выше 22,05 кГц. Все сочетания этих тонов кратны 1 кГц, и ниже 17 кГц сильнее всего разностный тон $f_2 - f_1 = 1$ кГц и его кратные: это интермодуляция чётных порядков, её уровень говорит о том, насколько асимметрична нелинейность. Остальные сочетания в этой области — высоких порядков: 16 кГц — уже седьмого, $4f_1 - 3f_2$. Таблица 1 перечисляет сочетания до третьего порядка.

При 44,1 кГц отражения смещены на 100 Гц относительно сетки 1 кГц, на которой лежат все настоящие сочетания: 4,1, 5,1, 12,9 кГц. Их видно по одной частоте. При 48 кГц отражения попадают на ту же сетку: $3f_1 = 57$ кГц переходит в 9 кГц, где могла бы стоять и настоящая составляющая. Отличить их можно только сравнением с выходом, который считался при более высокой частоте. Высокие порядки при 44,1 кГц ложатся на любые частоты, кратные 100 Гц: например, $4f_1 + 5f_2 = 176$ кГц девятого порядка оказывается на 400 Гц.

Измерения Каца

Кац сравнивал компрессоры, которые работают как модели компрессора LA-2A или похожи на него. Сигнал 19 + 20 кГц уходил из Pro Tools в прибор по AES/EBU в 24 битах и возвращался в анализатор SpectraPlus. Для рисунка 1 оригинала уровень был около −20 dBFS (decibels relative to full scale)?dBFS (decibels relative to full scale) — Уровень цифрового сигнала в децибелах относительно полной шкалы. Для пиковых значений 0 dBFS соответствует её границе; формат с плавающей точкой допускает значения выше неё. Соглашение об опоре RMS указывают отдельно. RMS при ослаблении компрессором около 1 дБ; на рисунках 1–3 частота дискретизации 44,1 кГц.

Опыт: два тона и отражения

Визуализация 1 повторяет тест на модели. Тоны 19 и 20 кГц по 0,25 (пик суммы −6 dBFS) проходят кривую $\tanh$, кривую $\tanh$ со смещением рабочей точки или жёсткое ограничение. График показывает спектр выхода при выбранной $f_s$: линии входных тонов, настоящих сочетаний и отражений различаются цветом, толщиной, пунктиром и метками, а крестики отмечают уровни, которые дала бы та же кривая в аналоговой цепи перед идеальным АЦП.

Начальное состояние — $\tanh$ с усилением +6 дБ при 44,1 кГц без передискретизации. Настоящие сочетания 18 и 21 кГц лежат на −33,3 dBFS, отражения 13,9 и 14,9 кГц — на том же уровне: это сочетания $2f_1 + f_2$ и $f_1 + 2f_2$ того же третьего порядка. Суммарная мощность отражений −29,9 дБ относительно синусоиды 0 dBFS. Если выбрать кривую со смещением, появится разностный тон 1 кГц на −30 dBFS, и он останется при любой кратности передискретизации, а отражения второго порядка 4,1, 5,1 и 6,1 кГц исчезнут уже при 2×. При 48 кГц отражения встают на сетку 1 кГц, и отличить их от сочетаний можно только по крестикам.

Передискретизация

Составляющая частоты $f$ при $L f_s$ отражается в $L f_s - f$, если $f$ лежит между $L f_s/2$ и $L f_s$. В исходную полосу она попадёт, только если $L f_s - f < f_s/2$. Всё, что ниже границы из формулы (2), фильтр после кривой уберёт.

При 44,1 кГц граница равна 66,15 кГц для 2×, 154,35 кГц для 4× и 330,75 кГц для 8×. У тонов 19 и 20 кГц сочетание порядка $p$ не выше $20p$ кГц, поэтому 2× убирает все отражения до третьего порядка, 4× — до седьмого, 8× — до шестнадцатого. В опыте с $\tanh$ при +6 дБ суммарный уровень отражений падает с −29,9 дБ до −52,6 при 2× и −96,8 при 4×, а при 8× опускается ниже −180 дБ. При 2× самое сильное отражение — 8,8 кГц от $3f_1 + 2f_2 = 97$ кГц пятого порядка, при 4× — 400 Гц от сочетания девятого порядка.

Насколько помогает каждое удвоение, зависит от того, как быстро убывают составляющие с ростом частоты (Zölzer, §10.4, с. 351). У гладкой $\tanh$ они убывают быстро. У жёсткого ограничения излом формы волны даёт медленный спад: в опыте при +12 дБ и 44,1 кГц отражения остаются на −49 дБ даже при 8×. Поэтому мягкую кривую и передискретизацию применяют вместе: у мягкой кривой высокие порядки слабее, а повышенная частота уводит их за границу.

Частота проекта 96 кГц сама работает как двукратная передискретизация относительно 48 кГц: в опыте отражения ниже 20 кГц при 96 кГц лежат на −53 дБ против −30 при 48 кГц. Тот же сигнал на 96 кГц обходится всей цепочке плагинов вдвое большим числом отсчётов, а передискретизация внутри одного эффекта повышает частоту только там, где стоит кривая.

Фильтры и задержка

Выигрыш из опыта получен с идеальными фильтрами. Реальный фильтр интерполяции или децимации имеет переходную полосу и конечное подавление, и составляющие из переходной полосы проходят децимацию и отражаются. Pirkle рекомендует крутые фильтры с линейной фазой длиной от 256 до 1024 коэффициентов и напоминает, что любая обработка внутри передискретизации, например фильтр, рассчитывается для частоты $L f_s$, а не для $f_s$ (Pirkle, гл. 22, с. 624–626).

Таблица 3 — измерения того же теста с juce::dsp::Oversampling в режиме максимального качества, с целой задержкой. В режиме 2× результат совпадает с идеальной моделью. При 4× и 8× отражения упираются в фильтры: остаток около −92…−98 дБ лежит на 21,1 кГц. Это сочетание $4f_2 - 3f_1 = 23$ кГц седьмого порядка: оно не отражается при повышенной частоте, но попадает в переходную полосу фильтра децимации и после понижения частоты переходит в $44{,}1 - 23 = 21{,}1$ кГц. Для этого теста 8× поэтому почти ничего не добавляет к 4×.

Без округления задержка КИХ-варианта 49, 59,5 и 64,25 отсчёта, БИХ-варианта — 3,14, 4,43 и 4,95. Хосту плагин сообщает целое число отсчётов, поэтому с useIntegerLatency = true класс добавляет в конце дробную задержку и сообщает целую. При 44,1 кГц 61 отсчёт — 1,4 мс. Кривая вычисляется в $L$ раз чаще, и к ней добавляются фильтры интерполяции и децимации каждой ступени в каждом канале.

Другие приёмы

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

Метод антипроизводной

Идеальное решение — восстановить непрерывный вход, применить к нему кривую, отфильтровать результат до $f_s/2$ и снова взять отсчёты. Метод антипроизводной (antiderivative antialiasing, ADAA) заменяет каждый шаг грубым приближением. Между соседними отсчётами вход считается отрезком прямой, а вместо фильтра выход усредняется за один период дискретизации (Zölzer, §10.4, с. 351–352). Среднее функции на отрезке выражается через её первообразную $F$, и получается формула (3).

Числитель и знаменатель дроби при почти равных соседних входах стремятся к нулю, и дробь теряет точность. В пределе она равна $f(u[n])$; Zölzer берёт среднее значений кривой в двух точках, потому что для прямой $f(u) = u$ дробь даёт ровно его (Zölzer, §10.4, с. 352–353). Порог $\varepsilon$ выбирают по точности чисел: во float вычитание $F$ теряет значащие цифры раньше, чем в double. Прямое вычисление std::log (std::cosh (u)) во float переполняется уже при $|u|$ около 89: $\cosh u$ превышает наибольшее число float. Поэтому $\ln\cosh u$ вычисляют как $|u| + \ln\bigl(1 + e^{-2|u|}\bigr) - \ln 2$.

Пример: соседние входы $u[n-1] = 0{,}2$ и $u[n] = 1$. Обычная кривая выдаёт $\tanh 1 \approx 0{,}762$. Метод антипроизводной выдаёт $(\ln\cosh 1 - \ln\cosh 0{,}2)/0{,}8 \approx (0{,}4338 - 0{,}0199)/0{,}8 \approx 0{,}517$ — среднее $\tanh$ на всём отрезке от 0,2 до 1.

Усреднение за период — фильтр нижних частот с задержкой на половину отсчёта. Для прямой $f(u) = u$ формула (3) даёт $(u[n] + u[n-1])/2$, и амплитуда синусоиды частоты $f$ умножается на $\cos(\pi f/f_s)$ (Zölzer, §10.4, с. 353; вывод множителя — в разделе о комплексных амплитудах). На низких частотах множитель близок к 1, а у верхнего края полосы фильтр срезает и полезный сигнал: при 44,1 кГц для 19 кГц множитель равен 0,216, то есть −13,3 дБ. Составляющие, которые отражаются в низкие частоты, метод подавляет сильно, а те, что лежат сразу выше $f_s/2$ и отражаются чуть ниже неё, — всего примерно на 3 дБ (Zölzer, §10.4, с. 353).

На тесте 19 + 20 кГц это видно по модели опыта. При 44,1 кГц без передискретизации суммарные отражения падают с −29,9 до −44,2 дБ, но и тоны теряют 12,5 и 15,9 дБ, так что относительно сигнала выигрыша почти нет; отражение на 400 Гц ослабевает с −101,4 до −120,3 dBFS. С двукратной передискретизацией метод даёт −65,5 дБ против −52,6 у одной передискретизации, а тоны теряют около 2 дБ. При 4× отражения −118,1 дБ против −96,8, потеря тонов 0,4 дБ. Поэтому метод применяют вместе с невысокой передискретизацией: она отодвигает $f_s/2$ от слышимой полосы, а метод подавляет то, что всё же отразилось (Zölzer, §10.4, с. 354).

Вычислений метод добавляет немного: если $F$ вычисляется так же быстро, как $f$, и $F(u[n])$ сохраняется для следующего шага, нужны два вычитания и одно деление (Zölzer, §10.4, с. 352). Ограничение — только кривые без памяти; для схем с фильтрами и обратной связью внутри нелинейности он в простом виде не подходит (Zölzer, §10.4, с. 354).

Код на C++

Класс ниже выполняет насыщение $\tanh$ при $4 f_s$. Конструктор juce::dsp::Oversampling принимает число каналов, показатель степени двойки для кратности (2 — две ступени по 2×, то есть 4×), вид фильтров, режим максимального качества и признак целой задержки. Число каналов известно только в prepare, поэтому объект создаётся там. processSamplesUp возвращает блок при повышенной частоте, кривая меняет его на месте, processSamplesDown фильтрует результат и записывает его обратно в исходный блок (JUCE 9.0.3, juce_Oversampling.h).

// Насыщение tanh с передискретизацией: кривая работает при 4·fs, фильтры JUCE
// убирают продукты выше исходной fs/2 до возврата к fs.
class OversampledSaturator
{
public:
    void prepare (const juce::dsp::ProcessSpec& spec)
    {
        oversampling = std::make_unique<juce::dsp::Oversampling<float>> (
            spec.numChannels, factorLog2,
            juce::dsp::Oversampling<float>::filterHalfBandPolyphaseIIR, true, true);
        oversampling->initProcessing (spec.maximumBlockSize);
    }

    void reset() noexcept { oversampling->reset(); }

    // Задержка в отсчётах исходной частоты: её сообщают хосту и учитывают в dry/wet.
    int getLatencySamples() const noexcept
    {
        return static_cast<int> (oversampling->getLatencyInSamples());
    }

    void process (juce::dsp::AudioBlock<float> block) noexcept
    {
        auto up = oversampling->processSamplesUp (block);

        for (size_t channel = 0; channel < up.getNumChannels(); ++channel)
            for (auto& x : std::span { up.getChannelPointer (channel), up.getNumSamples() })
                x = std::tanh (drive * x);

        oversampling->processSamplesDown (block);
    }

    float drive = 2.0f; // входное усиление, +6 дБ

private:
    static constexpr size_t factorLog2 = 2; // 2² = 4×
    std::unique_ptr<juce::dsp::Oversampling<float>> oversampling;
};

В процессоре плагина prepareToPlay вызывает prepare и reset, затем setLatencySamples (saturator.getLatencySamples()); processBlock передаёт в process блок буфера. Если у эффекта есть регулятор dry/wet, juce::dsp::DryWetMixer получает ту же задержку через setWetLatency, иначе смесь даст гребенчатую фильтрацию (comb filtering)?гребенчатая фильтрация (comb filtering) — Ряд чередующихся подъёмов и провалов частотной характеристики, возникающий, например, при сложении сигнала с его задержанной копией. (статья о dry/wet). Внутри передискретизации всё, что зависит от частоты дискретизации, — фильтры, сглаживание параметров — настраивают на $4 f_s$.

Код собран с JUCE 9.0.3 в консольном приложении с -Wall -Wextra -Wconversion -Wshadow без предупреждений. С БИХ-фильтрами и 4× класс сообщает задержку 6 отсчётов, как в таблице 3. Строка for (auto& x : std::span { … }) использует std::span (C++20) с выводом типа по аргументам конструктора (C++17).

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

Источники

• Katz B. Analog Saved My Life: The Good, The Bad and The Ugly. в списке литературы → — тест двумя тонами 19 и 20 кГц, условия измерения, рисунки 1–4: цифровая модель LA-2A, McDSP APB-16, копия LA-2A от ADL, Pendulum OCL-2, Tone Projects Unisum при 96 и 44,1 кГц; описание звучания алиасинга, многополосная и параллельная обработка, оговорка о связи измерений со звучанием.
• Zölzer U. Digital Audio Signal Processing. §10.1, с. 342–343; §10.4: Aliasing and its Mitigation, с. 350–354. в списке литературы → — сочетания частот нескольких тонов, перестановка дискретизации и кривой без памяти, передискретизация с фильтрами интерполяции и децимации, зависимость от спада гармоник, метод антипроизводной: вывод, формулы (10.25)–(10.27), случай близких отсчётов, задержка на половину отсчёта, подавление вблизи fs/2, сочетание с передискретизацией, ограничение кривыми без памяти.
• Reiss J. D., McPherson A. P. Audio Effects: Theory, Implementation and Application. Гл. 8: Aliasing and Oversampling, Filtering, с. 197–198. в списке литературы → — бесконечный ряд гармоник, отсутствие проблемы в аналоговой схеме, отражения, которые нельзя отфильтровать, передискретизация вставкой нулей и фильтрами, фильтр нижних частот перед кривой против интермодуляции.
• Pirkle W. C. Designing Audio Effect Plugins in C++: For AAX, AU, and VST3 With DSP Theory. Гл. 22: Sample Rate Conversion, с. 624–626. в списке литературы → — интерполятор и дециматор вокруг нелинейной обработки, крутые фильтры с линейной фазой длиной 256–1024, расчёт обработки внутри передискретизации для частоты L·fs.
• JUCE 9.0.3. juce_Oversampling.h и juce_Oversampling.cpp — параметры конструктора, КИХ- и БИХ-фильтры полуполосных ступеней и их задержка, initProcessing, processSamplesUp, processSamplesDown, getLatencyInSamples и целая задержка.
• JUCE 9.0.3. juce_DryWetMixer.h — setWetLatency: задержка обработанного сигнала, на которую смеситель задерживает исходный.