Спектр сигнала: ДПФ и БПФ

Осциллограмма показывает, как сигнал меняется во времени. О звуке часто спрашивают другое: какие частоты в нём есть и насколько сильна каждая. Ответ даёт спектр (spectrum)?спектр (spectrum) — Представление сигнала через его частотные составляющие: их частоты, амплитуды и фазы. Для отсчётов спектр вычисляют дискретным преобразованием Фурье.. По отсчётам (sample)?отсчёт (sample) — Значение сигнала в один момент дискретизации, одно число. В звукозаписи его также называют сэмплом. его вычисляет дискретное преобразование Фурье, а на практике — быстрый алгоритм того же преобразования. Статья выводит формулу из сравнения сигнала с пробной синусоидой, разбирает разрешение по частоте, утечку и оконные функции, устройство анализатора спектра и заканчивается кодом на C++23 и JUCE. Понадобятся комплексные числа и формула Эйлера.

Что показывает спектр

Спектр отвечает на вопрос «сколько каждой частоты» в сигнале. Для каждой частоты он даёт амплитуду (amplitude)?амплитуда (amplitude) — Наибольший модуль отклонения колеблющейся величины от равновесного значения. Для синусоиды обозначается A; единицы совпадают с единицами сигнала. и фазу (phase)?фаза (phase) — Положение внутри периода колебания, выраженное углом: полный период соответствует 360°, или 2π радиан. Начальная фаза — её значение в момент t = 0. синусоиды этой частоты, а сумма всех таких синусоид снова даёт сигнал. На графике спектра по горизонтали откладывают частоту, по вертикали — амплитуду или уровень в децибелах. Фазу обычно не рисуют: по ней трудно что-то увидеть, хотя для обратного перехода к сигналу она нужна.

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

Сравнение с пробной синусоидой

Возьмём кадр из $N$ отсчётов $x[0], x[1], \ldots, x[N-1]$ и проверим, есть ли в нём синусоида заданной частоты. Для этого каждый отсчёт умножают на отсчёт пробной синусоиды этой частоты и всё складывают (S. W. Smith, гл. 8, с. 157–159). Если в сигнале есть такая синусоида, произведения почти всё время одного знака, и сумма растёт. Если нет, положительные и отрицательные произведения уравновешивают друг друга.

Пусть $N = 64$, а сигнал — косинусоида с амплитудой 1, которая делает в кадре 3 периода: $x[n] = \cos(2\pi \cdot 3n/64)$. Пробный косинус тех же трёх периодов даёт произведения $\cos^2$, они неотрицательны, а их среднее равно 1/2. Сумма 64 произведений равна 32, то есть $N/2$. Пробный косинус с пятью периодами даёт сумму ровно 0. Так получается для любых двух разных целых чисел периодов от 1 до $N/2 - 1$: сумма произведения таких косинусоид по кадру равна нулю. Это свойство называют ортогональностью (orthogonality); на нём основано ДПФ (Sundararajan, §3.1).

Если амплитуда сигнала равна $A$, сумма равна $AN/2$. Значит, амплитуду можно найти, умножив сумму на $2/N$.

Одной пробной синусоиды мало. Сдвинем фазу сигнала на четверть периода: $x[n] = \cos(2\pi \cdot 3n/64 + \pi/2) = -\sin(2\pi \cdot 3n/64)$. Частота та же, а сумма с пробным косинусом теперь равна нулю: произведения синуса на косинус одной частоты тоже уравновешиваются. Поэтому берут две пробы, косинус и синус, и две суммы:

Для $x[n] = A\cos(2\pi kn/N + \varphi)$ при целом $k$ от 1 до $N/2 - 1$ получается $C = (AN/2)\cos\varphi$ и $S = -(AN/2)\sin\varphi$. Фаза делит длину $AN/2$ между двумя суммами, как угол делит радиус между проекциями на единичной окружности. Отсюда амплитуда $A = (2/N)\sqrt{C^2 + S^2}$, а фаза — угол точки $(C, -S)$.

Пару $(C, -S)$ удобно записать одним комплексным числом $C - jS$. По формуле Эйлера $\cos\theta - j\sin\theta = e^{-j\theta}$, поэтому обе суммы — это одна сумма с комплексной пробой:

Модуль результата пропорционален амплитуде, аргумент равен фазе. Это комплексная амплитуда синусоиды, умноженная на $N/2$.

В начальном состоянии визуализации 2 сигнал и пробный косинус делают по 3 периода, и накопленная сумма растёт до 32. Сдвиньте фазу сигнала на 90°: сумма с косинусом станет нулём, сумма с синусом — −32, а модуль и амплитуда не изменятся. Поставьте пробную частоту 5: обе суммы обнулятся. Наконец, поставьте частоту сигнала 3,5 периода: ни на одной пробной частоте сумма не обнуляется, при $k = 3$ и $k = 4$ оценки амплитуды 0,59 и 0,68 вместо 1. К этому случаю вернёмся в разделе об утечке.

Дискретное преобразование Фурье

Если повторить комплексную сумму для каждого целого числа периодов $k$ от 0 до $N - 1$, получится дискретное преобразование Фурье (discrete Fourier transform, DFT)?дискретное преобразование Фурье (discrete Fourier transform, DFT) — Преобразование N отсчётов в N комплексных чисел X[k] = Σ x[n]·e^(−j2πkn/N). Модуль и аргумент X[k] дают амплитуду и фазу составляющей с частотой k·fₛ/N., сокращённо ДПФ (J. O. Smith, DFT Definition; Sundararajan, §3.1, формула 3.3):

ДПФ превращает $N$ отсчётов в $N$ комплексных чисел. Для каждого $k$ выполняется $N$ комплексных умножений, всего $N^2$. Формула (1) не требует, чтобы сигнал был синусоидой: она применима к любому кадру, а результат описывает этот кадр как сумму синусоид с целым числом периодов.

Бины и частоты

Каждое значение $k$ называют бином (frequency bin)?бин (frequency bin) — Одна частотная точка ДПФ с номером k: частота k·fₛ/N, шаг до соседней fₛ/N. У вещественного сигнала независимы бины от 0 до N/2.. Проба бина $k$ делает $k$ периодов за кадр, а кадр длится $N/f_s$ секунд, где $f_s$ — частота дискретизации (sample rate)?частота дискретизации (sample rate) — Число отсчётов в каждом канале за одну секунду, то же число кадров аудиоданных в секунду. Измеряется в герцах: 44 100 Гц означает 44 100 отсчётов на канал в секунду.. Отсюда частота бина:

При $N = 1024$ и $f_s = 48$ кГц бин 1 соответствует 46,875 Гц, бин 21 — 984,4 Гц, бин 512 — 24 кГц. Синусоида 1000 Гц делает в таком кадре $1000/46{,}875 \approx 21{,}33$ периода и попадает между бинами 21 и 22. Шаг между бинами одинаков на всех частотах (Sundararajan, §3.1, с. 75).

Амплитуда и фаза

Пусть $x[n] = A\cos(2\pi k_0 n/N + \varphi)$, где $k_0$ — целое от 1 до $N/2 - 1$. По разделу о пробе $X[k_0] = (AN/2)\,e^{j\varphi}$, а остальные бины от 1 до $N/2$ равны нулю. Амплитуду и фазу дают модуль и аргумент:

Множитель $2/N$ не действует для двух бинов. Бин 0 — сумма отсчётов, и $X[0]/N$ — их среднее, постоянная составляющая. Бин $N/2$ — сумма отсчётов с чередующимся знаком, и $X[N/2]/N$ — амплитуда последовательности $+1, -1, +1, \ldots$. Для этих бинов амплитуда равна $|X[k]|/N$. Нормировку разные книги и библиотеки выбирают по-разному: Bellanger относит множитель $1/N$ к прямому преобразованию (Bellanger, §2.2), Sundararajan, J. O. Smith и JUCE — к обратному. Перед расчётом уровней нужно проверить соглашение своей библиотеки.

Таблица 1 — пример из учебника Sundararajan: четыре отсчёта сигнала $x[n] = 1 + \cos(2\pi n/4 + \pi/3) - \cos(\pi n)$ и их ДПФ (Sundararajan, §3.1, с. 70–72).

Для $k = 1$ сумма с пробой $e^{-j2\pi n/4}$, то есть с числами $1, -j, -1, j$, равна $0{,}5 - j1{,}134 + 0{,}5 + j2{,}866 = 1 + j1{,}732$. Её модуль $\sqrt{1 + 3} = 2$ и аргумент $\operatorname{atan2}(1{,}732;\ 1) = 60°$ восстанавливают амплитуду 1 и фазу $\pi/3$ второго слагаемого.

Симметрия и частота Найквиста

Новых сведений бины выше $N/2$ не несут, и анализаторы показывают только бины от 0 до $N/2$ — это $N/2 + 1$ чисел (S. W. Smith, гл. 8). Верхний бин $N/2$ соответствует $f_s/2$, частоте Найквиста (Nyquist frequency)?частота Найквиста (Nyquist frequency) — Половина частоты дискретизации, fₛ/2. В принятой низкочастотной модели точное восстановление по отсчётам требует, чтобы все частоты исходного сигнала были строго ниже этой границы.. Частот выше неё в отсчётах нет: если до дискретизации они были, то при ней уже отразились ниже $f_s/2$, и ДПФ покажет их на отражённой частоте.

Разрешение по частоте

Шаг между бинами задаёт разрешение по частоте:

Шаг по частоте обратен длительности кадра. Bellanger называет это соотношение принципом неопределённости для спектрального анализа: точнее различить частоты можно, только дольше наблюдая сигнал (Bellanger, §2.4.2). Таблица 2 даёт числа для 48 кГц.

Ноты E1 и F1 нижнего регистра бас-гитары звучат на частотах 41,20 и 43,65 Гц, между ними 2,45 Гц. Чтобы они попали хотя бы в разные бины, нужен шаг меньше 2,45 Гц, то есть $N \ge 48\,000/2{,}45 \approx 19\,600$: ближайшая степень двойки — 32 768 отсчётов, 0,68 с. За это время нота может закончиться, а удар барабана в таком кадре не локализовать точнее, чем в пределах всего кадра. Короткий кадр из 256 отсчётов покажет начало удара с точностью 5 мс, но вся полоса от 0 до 187 Гц уложится в один шаг между бинами. Подробная картина по частоте требует длинного кадра, а подробная во времени — короткого.

Можно дополнить кадр нулями до большего $N$. Бинов станет больше, но новые значения лишь интерполируют спектр между прежними: два близких пика, которые сливались, так и останутся слитыми (S. W. Smith, гл. 9, с. 174). Разрешение определяет длительность настоящего сигнала в кадре, а не длина преобразования.

Утечка и окна

Синусоида 1000 Гц делает в кадре из 1024 отсчётов 21,33 периода: кадр обрывается посреди периода. ДПФ считает кадр одним периодом бесконечно повторяющегося сигнала (Sundararajan, §3.1). Конец кадра стыкуется с началом со скачком, и чтобы описать скачок, нужны все частоты. Уровень появляется во всех бинах, а не только у частоты синусоиды. Это утечка спектра (spectral leakage)?утечка спектра (spectral leakage) — Появление уровня в бинах вдали от частоты синусоиды, когда в кадре не укладывается целое число её периодов. Причина — обрыв сигнала на краях кадра. (J. O. Smith, Frequencies in the «Cracks»).

Утечка видна и по числам. Синусоида с амплитудой 0,5 дала бы на бине −6,02 дБ. При 1000 Гц бин 21 показывает −7,64 дБ, бин 22 — −13,74 дБ, бин на 10 выше пика — −37,4 дБ, на 50 выше — −51,7 дБ, и до нуля уровень не опускается нигде. Слабый тон 1600 Гц на 40 дБ тише, то есть на уровне −46 дБ, в 13 бинах от сильного тонет в этом основании. Если частота попадает точно посередине между бинами, пик проседает сильнее всего: в коде этой статьи — на 3,82 дБ. Увеличение $N$ сужает основание в герцах, но уровни боковых выбросов в децибелах остаются прежними (J. O. Smith).

Окно Ханна

Скачок на стыке можно убрать, если перед ДПФ умножить кадр на функцию, которая плавно сводит его края к нулю. Такой множитель называют оконной функцией (window function)?оконная функция (window function) — Множитель, на который поотсчётно умножают кадр перед ДПФ. Сводит края кадра к нулю и уменьшает утечку спектра ценой более широкого пика; пример — окно Ханна., или окном. Без него кадр тоже умножен на окно — прямоугольное, из одних единиц. Одно из самых распространённых плавных окон — окно Ханна (Hann window) (Sundararajan, §5.2.2):

Так окно Ханна заполняет juce::dsp::WindowingFunction в JUCE (juce_Windowing.cpp). Окно меняет форму пика. Основание пика синусоиды — главный лепесток (main lobe) — у окна Ханна вдвое шире: ±2 бина вместо ±1. Зато боковые лепестки (side lobes) падают: первый — до −31,5 дБ от пика вместо −13 дБ, а дальше они спадают на 18 дБ за октаву (J. O. Smith, Hann Window). Таблица 3 сравнивает оба окна на синусоиде 1000 Гц.

S. W. Smith описывает это как обмен: окно уменьшает утечку, то есть высоту «хвостов» пика, ценой разрешения, то есть ширины пика (S. W. Smith, гл. 9, с. 174–175). Прямоугольное окно лучше различает два тона одинаковой силы на соседних частотах. Окно Ханна лучше показывает слабый тон рядом с сильным; его же использует анализатор из урока JUCE.

Окно уменьшает отсчёты: среднее значение окна Ханна около 1/2, поэтому синусоида на бине без поправки покажет уровень на 6 дБ ниже. JUCE по умолчанию (normalise = true) умножает окно на $N/\sum w[n]$, то есть на $2N/(N-1)$. Тогда масштаб $2/N$ из формулы (3) верен для любого окна.

Визуализация 3 считает спектр начала сигнала тем же способом, что анализатор из раздела о коде: окно с нормировкой, БПФ и масштаб $2/N$. Звучит тот же сигнал. В начальном состоянии синусоида 1000 Гц попадает между бинами, и пик −6,65 дБ стоит на бине 21. Переключите окно на прямоугольное: основание расширится, через 10 бинов уровень будет −37 дБ вместо −76. Поставьте 1500 Гц — ровно 32 бина при $N = 1024$: прямоугольное окно даст одну линию −6,02 дБ, окно Ханна — пик −6,02 дБ и два соседних бина по −12,04 дБ. На сигнале «два тона» слабый тон 1600 Гц с прямоугольным окном виден только при $N = 4096$, а с окном Ханна — уже при 512. Выберите белый шум и меняйте $N$: уровень бинов при каждом удвоении падает на 3 дБ, а уровень синусоиды не меняется.

Быстрое преобразование Фурье

По формуле (1) кадр из 1024 отсчётов требует $1024^2 = 1\,048\,576$ комплексных умножений. Анализатор с перекрытием кадров на половину считает при 48 кГц около 94 кадров в секунду — около 98 миллионов комплексных умножений в секунду на один канал.

Быстрое преобразование Фурье (fast Fourier transform, FFT)?быстрое преобразование Фурье (fast Fourier transform, FFT) — Алгоритм, который вычисляет то же ДПФ за число операций порядка N·log₂N вместо N². Распространённые реализации требуют, чтобы N было степенью двойки., БПФ, вычисляет те же $X[k]$ с точностью до округления, но за число операций порядка $N\log_2 N$, а не $N^2$ (запись O). Простейший вариант делит кадр на чётные и нечётные отсчёты, вычисляет два ДПФ половинной длины и собирает из них результат, добавляя $N/2$ умножений. Половины делятся дальше, пока не останутся пары отсчётов. Ступеней $\log_2 N$, на каждой $N/2$ умножений, всего $(N/2)\log_2 N$. Bellanger уточняет, что преобразования из двух отсчётов, до которых доходит деление, умножений не требуют, и получает $(N/2)\log_2(N/2)$ (Bellanger, §2.2.1). Элементарное действие алгоритма, которое из двух чисел делает два новых, называют бабочкой (butterfly) по виду его схемы (S. W. Smith, гл. 12).

У S. W. Smith БПФ на 1024 точки на компьютере 1990-х годов работает в 300 с лишним раз быстрее, чем ДПФ по формуле (S. W. Smith, гл. 12, с. 237). Выигрыш растёт с $N$. Деление пополам требует, чтобы $N$ было степенью двойки: juce::dsp::FFT принимает показатель степени order, и размер кадра равен $2^{order}$ (juce_FFT.h). Для других длин существуют свои алгоритмы, но обычно кадр дополняют нулями до степени двойки (Sundararajan, гл. 10, с. 357).

Анализатор спектра

Уровни и ось частот

Анализатор переводит амплитуду бина в децибелы: $L_k = 20\lg\bigl(2|X[k]|/N\bigr)$. При таком масштабе синусоида с амплитудой 1 на бине показывает 0 дБ, а с амплитудой 0,5 — −6,02 дБ, как на рисунке 1. Опорный уровень выбирают по-разному. В уроке JUCE об анализаторе из уровня вычитают gainToDecibels (fftSize), то есть делят модуль на $N$, а не на $N/2$, и синусоида полной шкалы там показывает −6 дБ (JUCE, Spectrum analyser tutorial). Анализатор полезно проверить тестовой синусоидой известной амплитуды.

Перекрытие и усреднение

Звук идёт непрерывно, и анализатор берёт кадр за кадром. Окно Ханна почти обнуляет края кадра, поэтому отсчёты у границ кадров почти не учитываются. Чтобы не терять их, кадры перекрывают: новый кадр начинается через $N/2$ отсчётов после предыдущего. Перекрытие на половину и больше — обычный выбор (J. O. Smith, Classic Spectrograms). При $N = 1024$ и 48 кГц новый спектр появляется каждые 512 отсчётов, 93,75 раза в секунду. Экран обновляется 60 раз в секунду, и между кадрами экрана анализатор показывает последний готовый спектр.

Спектрограмма

Если спектры кадров выстроить по времени, получится спектрограмма (spectrogram): время по горизонтали, частота по вертикали, уровень — яркостью (J. O. Smith, Classic Spectrograms). На рисунке 4 один и тот же сигнал разобран кадрами двух размеров.

При $N = 256$ кадр длится 5,3 мс: смена нот и щелчок видны точно по времени, но тоны 440 и 494 Гц, между которыми 54 Гц, сливаются в одну широкую полосу — шаг бинов 187,5 Гц. При $N = 4096$ кадр длится 85 мс, шаг 11,7 Гц, и два тона — отдельные линии. Зато щелчок длиной 20 мс растянут на 0,1 с, а смена нот размыта. Это компромисс из раздела о разрешении.

Обратное ДПФ

Из $N$ чисел $X[k]$ отсчёты восстанавливаются точно (Sundararajan, §3.1, формула 3.2):

От формулы (1) она отличается знаком в показателе и множителем $1/N$. Проверим на таблице 1 при $n = 0$: все экспоненты равны 1, и $x[0] = (4 + 1 + j1{,}732 - 4 + 1 - j1{,}732)/4 = 0{,}5$. Обратное преобразование складывает все синусоиды бинов с их амплитудами и фазами.

Пара прямого и обратного ДПФ позволяет обрабатывать звук в частотной области: вычислить спектр кадра, изменить $X[k]$ и вернуться к отсчётам. Простое обнуление бинов работает плохо: ДПФ считает кадр периодическим, и изменения переходят через края кадра. Правильная схема опирается на свёртку: умножение спектров соответствует свёртке сигналов. Через БПФ свёртка длинных сигналов выполняется за $O(N\log N)$ операций вместо $O(N^2)$ (Sundararajan, §3.3.9, с. 90), так работает свёрточный ревербератор. Свёртку и фильтры разбирает статья о цифровых фильтрах. В JUCE обратное преобразование выполняет performRealOnlyInverseTransform; множитель $1/N$ она применяет сама.

Код на C++

Наивное ДПФ

// ДПФ по определению: N² комплексных умножений; эталон для проверки быстрых реализаций.
#include <complex>
#include <numbers>
#include <span>
#include <vector>

std::vector<std::complex<double>> naiveDft (std::span<const float> x)
{
    const auto size = x.size();
    std::vector<std::complex<double>> spectrum (size);

    for (std::size_t k = 0; k < size; ++k)
    {
        for (std::size_t n = 0; n < size; ++n)
        {
            // Угол по остатку k·n mod N: при больших N произведение не теряет точность.
            const auto turn = static_cast<double> ((k * n) % size) / static_cast<double> (size);
            const auto angle = -2.0 * std::numbers::pi * turn;
            spectrum[k] += static_cast<double> (x[n]) * std::polar (1.0, angle);
        }
    }

    return spectrum;
}

std::polar (1.0, angle) строит $e^{j\cdot angle}$ — точку единичной окружности (математический справочник). Угол берётся по остатку $kn \bmod N$: проба периодична с периодом $N$, а аргумент синуса и косинуса остаётся меньше $2\pi$. Сумма накапливается в double. Код собран clang++ -std=c++20 -Wall -Wextra -Wconversion без предупреждений. На косинусоиде с тремя периодами в 64 отсчётах он даёт $|X[3]| = 32$, то есть амплитуду 1, а на отсчётах таблицы 1 — те же $X[k]$.

Анализатор на JUCE

Анализатор делит работу между двумя потоками, как в статье об аудиопотоке. Аудиопоток в processBlock только копирует отсчёты в очередь без блокировок на juce::AbstractFifo. Поток сообщений забирает их, собирает кадр, применяет окно, вычисляет БПФ и уровни.

// Анализатор: аудиопоток кладёт отсчёты в очередь, поток сообщений считает спектр кадра.
#include <juce_dsp/juce_dsp.h>

class SpectrumAnalyzer
{
public:
    static constexpr int fftOrder = 10;
    static constexpr int fftSize = 1 << fftOrder; // N = 1024
    static constexpr int hopSize = fftSize / 2;   // кадры перекрываются на 50 %

    // Аудиопоток: без блокировок и выделений; если очередь полна, отсчёты теряются.
    void pushSamples (const float* samples, int count) noexcept
    {
        const auto scope = fifo.write (count);
        std::copy_n (samples, scope.blockSize1, queue.begin() + scope.startIndex1);
        std::copy_n (samples + scope.blockSize1, scope.blockSize2,
                     queue.begin() + scope.startIndex2);
    }

    // Поток сообщений: true, если пришло hopSize новых отсчётов и спектр пересчитан.
    bool update()
    {
        if (fifo.getNumReady() < hopSize)
            return false;

        while (fifo.getNumReady() >= hopSize)
        {
            std::copy (history.begin() + hopSize, history.end(), history.begin());
            const auto scope = fifo.read (hopSize);
            auto* tail = history.data() + (fftSize - hopSize);
            std::copy_n (queue.begin() + scope.startIndex1, scope.blockSize1, tail);
            std::copy_n (queue.begin() + scope.startIndex2, scope.blockSize2,
                         tail + scope.blockSize1);
        }

        std::copy (history.begin(), history.end(), fftData.begin());
        window.multiplyWithWindowingTable (fftData.data(), fftSize);
        fft.performFrequencyOnlyForwardTransform (fftData.data(), true); // |X[k]|, k = 0 … N/2

        for (int k = 0; k <= fftSize / 2; ++k)
        {
            const auto bin = static_cast<std::size_t> (k);
            const float scale = (k == 0 || k == fftSize / 2) ? 1.0f / fftSize : 2.0f / fftSize;
            levels[bin] = juce::Decibels::gainToDecibels (fftData[bin] * scale, -140.0f);
        }
        return true;
    }

    const auto& getLevels() const noexcept { return levels; } // дБ, бины 0 … N/2

private:
    juce::dsp::FFT fft { fftOrder };
    juce::dsp::WindowingFunction<float> window { fftSize,
                                                 juce::dsp::WindowingFunction<float>::hann };
    juce::AbstractFifo fifo { 8 * fftSize };
    std::array<float, 8 * fftSize> queue {};
    std::array<float, fftSize> history {};
    std::array<float, 2 * fftSize> fftData {};
    std::array<float, fftSize / 2 + 1> levels {};
};

performFrequencyOnlyForwardTransform требует массив из $2N$ чисел: первая половина — отсчёты, вторая — рабочее место. С аргументом true она вычисляет только бины от 0 до $N/2$ и записывает их модули $|X[k]|$ в начало массива (juce_FFT.h). Прямое преобразование в JUCE не делит на $N$: это сумма формулы (1). WindowingFunction с окном hann по умолчанию нормирует окно, поэтому масштаб $2/N$ даёт амплитуду синусоиды. Бины 0 и $N/2$ получают масштаб $1/N$. gainToDecibels с порогом −140 дБ возвращает −140 для нулевого модуля, а не $-\infty$.

Код собран с JUCE 9.0.3 в консольном приложении с флагами предупреждений JUCE и -Wconversion. На синусоиде 1500 Гц с амплитудой 0,5 при 48 кГц он даёт −6,02 дБ на бине 32 и −12,03 дБ на соседних, как визуализация 3 с окном Ханна.

В уроке JUCE об анализаторе спектра кадр копируется в общий массив в аудиопотоке, а готовность отмечает поле bool nextFFTBlockReady, которое таймер читает в потоке сообщений (JUCE, Spectrum analyser tutorial). Обычный bool и общий массив, которые один поток пишет, а другой читает, по правилам C++ — гонка данных, а значит, неопределённое поведение. Очередь AbstractFifo упорядочивает запись и чтение атомарными операциями, а БПФ остаётся в потоке сообщений.

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

Источники

• Sundararajan D. Digital Signal Processing: An Introduction. §3.1, с. 67–75; §3.3.9, с. 90; §5.2.1–5.2.2, с. 165–171; гл. 10, с. 357–358. в списке литературы → — ДПФ и обратное ДПФ, ортогональность проб, пример из четырёх отсчётов, сопряжённая симметрия, шаг $f_s/N$, свёртка через ДПФ, периодограмма и метод Уэлча, окна и таблица уровней боковых лепестков, сложность $O(N\log_2 N)$ и длины-степени двойки.
• Bellanger M. Digital Signal Processing: Theory and Practice. §2.2, с. 38–40; §2.4, с. 46–49. в списке литературы → — матричная запись ДПФ с множителем $1/N$, БПФ с прореживанием по времени и число умножений $(N/2)\log_2(N/2)$, ДПФ как банк фильтров, разрешение $1/(NT)$, утечка между бинами и окна.
• Smith S. W. The Scientist and Engineer’s Guide to Digital Signal Processing. Гл. 8: The Discrete Fourier Transform, с. 157–160, гл. 9: Applications of the DFT, с. 169–175, гл. 12: The Fast Fourier Transform, с. 225–237. в списке литературы → — ДПФ через корреляцию с пробными синусоидами, $N/2 + 1$ частот, усреднение спектров кадров, разрешение и дополнение нулями, утечка и окна, бабочка и сравнение скорости ДПФ и БПФ.
• Smith J. O. III. Mathematics of the Discrete Fourier Transform (DFT) with Audio Applications. DFT Definition и Frequencies in the «Cracks». в списке литературы → — определение ДПФ и частоты бинов, утечка для частот между бинами, независимость уровня утечки от $N$.
• Smith J. O. III. Spectral Audio Signal Processing. Hann or Hanning or Raised Cosine и Classic Spectrograms. в списке литературы → — главный и боковые лепестки окна Ханна, спектрограмма как последовательность спектров перекрывающихся кадров.
• JUCE 9.0.3. juce_FFT.h, juce_FFT.cpp, juce_Windowing.cpp — размер $2^{order}$, performFrequencyOnlyForwardTransform и массив $2N$, прямое преобразование без деления на $N$, окно Ханна с делителем $N-1$ и нормировка normalise.
• JUCE. Tutorial: Visualise the frequencies of a signal in real time — анализатор из урока JUCE: размер кадра, окно Ханна, логарифмическая ось и пересчёт в децибелы.