Семплы. Анализ спектра оцифрованного звука

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

   ПРОСТОЕ РЕШЕНИЕ

   Поскольку программа «Пианола_М» использует в своём синтаксисе нотные знаки, то для анализа спектра я не стал заморачиваться с преобразованиями Фурье, и поступил гораздо проще, я связал распознавание нот с физической моделью «груза на пружинке» - каждая нота снабжается своим собственным распознавателем, груз на пружинке имеет массу равную единице, а жёсткость пружины такова, что резонанс колебаний этой системы равен частоте ноты.

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

1. Задаём коэффициент жёсткости пружины по формуле K=(2пи*F)^2 и
  устанавливаем груз в начальное положение, то есть, полагаем
  его координату и скорость равными нулю - X=0: V=0
  и, выбрав на семпле место для начала анализа, запускаем процесс.

2. Читаем амплитуду семпла и используем её как дробинку, которая
  стреляет в груз, передавая ему свой импульс А. Одновременно и
  пружина, если она сжата или растянута, следуя второму закону
  Ньютона, изменяет импульс груза на величину X*K^2
  В результате груз получает приращение скорости -
    dV = ( A - X*K^2 )*dT

3. Реализуем полученное приращение, изменяя скорость, и, вслед за
  скоростью, изменяем координату X -
    V = V + dV: X = X + V*dT

4. Анализируем полученное Х - если оно по абсолютной величине
  превосходит предыдущие, то запоминаем его.
  Переходим к следующему семплу и выясняем, не последний ли это
 семпл в нашем анализируемом интервале времени. Если последний, то
 заканчиваем процесс, а если нет, то тогда возвращаемся на пункт 2.

   Описанный процесс можно проделать поочерёдно для каждого датчика из набора, или же включить его внутренним циклом и перебирать 96 датчиков каждый раз при продвижении на 1 семпл. Результат будет один и тот же - мы получим 96 значений максимального отклонения в колебаниях наших датчиков. Это и есть искомый спектр звука.

   ОДНАКО ИМЕЮТСЯ НЮАНСЫ

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

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

   Об этом свидетельствует первая диаграмма, приведённая на иллюстрации. Спектр частот на ней относится к нотам октав О4-О5 (первая и вторая в обычном понимании). На диаграмме анализируется составленный из синусоид звук аккорда D:F:F#:G:B:c в котором ноты F и G имеют громкость 1/2, а нота  B имеет громкость 1/10 относительно громкости опорной ноты D.
   По рисунку видно, что эта пропорция выдержана правильно, так же, как можно заметить одинаковость амплитуд для синусоид разных частот одинаковой единичной громкости.

   На следующий шести диаграммах приведены спектры аккорда, составленного из нот хроматической гаммы. Для такого аккорда в идеале, при большом отрезке анализируемого звука, должен получаться частокол из одинаковых тонких линий.
   Однако условия у нас не идеальные.
   Диаграммы отличаются длиной интервала - в первой анализируется звук на отрезке 1/32 сек, на второй - 1/16, далее 1/8, 1/4, 1 сек и 4 сек. И да, при 4-х секундах тонкие ровные линии мы получаем (такие же линии получаются и в более низких октавах), но в 5-й октаве становится заметным плавное падение чувствительности наших датчиков по мере увеличения частоты.

   Исследования показали, что наблюдаемое падение чувствительности на самом деле кажущееся - дело в том, что если отрезок для анализа взят большой, то спектральное разрешение датчиков увеличивается, и они перестают замечать те колебания, частота которых находится между частотами той гребёнки, на которые датчики настроены.
   Однако, возможно ли это? Можно ли промахнуться мимо точно заданных частот, если сами датчики на эти частоты настраиваются?
   Оказывается, что можно.
   Причина такого явления кроется в дискретности оцифровки, и хотя время раздачи составляет 32000 семпла в секунду, а частота ноты Фа первой октавы в сто раз меньше, но уже начиная с этой ноты наблюдается расхождение между той частотой, на которую мы датчик настроили, и той резонансной частотой, которая получилась реально.
   А реальная частота оказывается больше, при небольших частотах это расхождение небольшое - для ноты до# второй октавы оно составляет 0.06%, однако октавой выше это уже 0.8% и далее по нарастающей.

   Избежать этого неприятного явления можно, если настраивать датчики не на теоретическую, а несколько меньшую частоту, используя для этого корректирующий мержитель:
    F = F_теоретич. * Mf    -здесь множитель Mf зависит от частоты и он меньше единицы.

   Возможность коррекции спектрограммы на высоких частотах показана на следующих двух рисунках - слева не корректированный спектр, а справа - спектр, полученный от датчиков с корректированными частотами. На этих диаграммах представлены октавы О2-О5 и пара нот из октавы О6, охват частот более широкий, поэтому и линии идут чаще.
   Видно, что коррекция даёт отличные результаты - частокол линий хроматического аккорда получается очень ровным.

   Практическое применение спектрального анализа продемонстрировано в нижней части иллюстрации.
   После входа в редактор выбранный семпл отображается на нижней части клавиатуры. В нашем случае это семпл певческого голоса сопрано длительностью около 0.8 секунд. После клика по тому месту, которое указывается поинтером-стрелочкой, под осциллограммой появляется красная полоска, обозначающая отрезок, на котором проводится спектральный анализ.
 
   Длина этой полоски может быть разной, от 1/32 секунды до 4-х секунд. Используя полоску короткой длины, удобно исследовать особенности спектра в различных частях семпла и находить изменения спектрального состава звука. Клик по спектральной диаграмме выведет сообщение о ноте, находящейся на указанной кликом линии.
   В семпле певческого голоса, приведённом на иллюстрации, в начале звука доминируют ноты f-e октавы О5 и имеется низкий призвук ноты E в октаве О4. В середине звука нота f становится более явной, а низкий призвук исчезает.

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

_________
3.04.2025    
 
ПРИЛОЖЕНИЕ
Два варианта использования датчиков. Во втором варианте индексированные переменные не используются во внутреннем цикле, и поэтому он работает быстрее. Зато первый вариант можно использовать в режиме анализа звука онлайн, подключив анализатор непосредственно к потоку семплов.

 
1 For I = 1 To 96: aaX(I) = 0: aaV(I) = 0: aaXM(I) = 0: Next I 'исходные условия
  For I = Io To Im: A = ArS(I, 15): A = A * 100:
  For J = 1 To 96:
    dV = (A - aaKm(J) * aaX(J)) * DT:
    aaV(J) = aaV(J) + dV: aaX(J) = aaX(J) + aaV(J) * DT:
    If aaX(J) > aaXM(J) Then aaXM(J) = aaX(J):
  Next J: Next I:
 
    'более быстрый алгоритм
2 For J = 1 To 96: aKm = aaKm(J): aX = 0: AV = 0: aXM = 0:
  For I = Io To Im: A = ArS(I, 15): A = A * 100:
    dV = (A - aKm * aX) * DT:
    AV = AV + dV: aX = aX + AV * DT: 'B = aaXM(J):
    If aX > aXM Then aXM = aX:
  Next I: aaXM(J) = aXM: Next J: