Анализ

Одни говорят, что это камень,
другие - что это птица.
В действительности это яйцо.

Ланца дель Васто

Cпектральный анализ, в частности, основанный на дискретном преобразовании Фурье (ДПФ), является исключительно эффективным инструментом. Однако следует отдавать себе отчет, что это не ковер-самолет с Хоттабычем в придачу. Применение спектрального анализа требует аккуратной настройки множества параметров в зависимости от конкретной задачи и специфики анализируемых сигналов, а для правильной интерпретации результатов нужен некоторый практический опыт.

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

Пожалуй, наиболее окупаемые области применения спектрального анализа следующие:

- беспристрастное сравнение звучания акустических систем, музыкальных инструментов, воспроизведения звука магнитофонами, CD- и DVD-устройствами и т. п.;

- ремастеринг старых записей;

- прогнозирование неисправностей по появлению специфических шумов, звуков, вибраций (например, в винчестерах).

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

В первой области - свои нюансы, но их гораздо меньше. Причем каждый меломан, владеющий азами программирования, будет счастлив получить свои законные спектры с желаемым разрешением по частоте и сравнить качество звука. Можно запастись терпением и, изнасиловав компьютер, вычислить спектр, не используя алгоритм быстрого преобразования Фурье (БПФ), в лоб, по формуле ДПФ [1]:

width=192

,

где k=0, ..., N/2 - индекс спектральных компонент (другими словами, частотных составляющих); xi - блок (массив), вырезанный из анализируемого сигнала, объемом N дискрет; wi - функция временного окна (в его отсутствие wi=1), width=57

. Ясно, что результатом вычислений будет комплексное число. Возведя в квадрат и сложив его мнимую и реальную части, получим мощность соответствующей спектральной составляющей, а извлекая корень квадратный от последней - условную амплитуду, иногда называемую магнитудой. Пробежав в цикле по всем k, получим спектр. При k=0 имеем дело с условно нулевой частотой, то есть со средним значением амплитуды анализируемого сигнала, которая на практике чаще всего отлична от нуля из-за DC-эффекта, часто называемого "смещением нуля" - неизбежного спутника аналоговой аппаратуры и особенно - фильтров. Величина разрешения (resolution) по частоте определяется тривиально:

width=56

,

где f - частота дискретизации. Чем больше объем блока в дискретах, тем лучше разрешение по частоте и тем больше деталей в частотной области может быть выяснено. Например, для 44100 Гц и блока объемом 2048 получим, что каждый частотный компонент "отвечает" за диапазон частот 21,5 Гц. В этом случае, если говорят о спектральной составляющей с частотой 21,5 Гц, подразумевают, что данная составляющая занимает частотный диапазон от 10,75 до 32,25 Гц включительно [1].

БПФ даст точно такой же результат, что и ДПФ, но при существенно меньшем количестве вычислений (в N/log2N раз). При этом необходимо, чтобы N составляло 2m, например 2048. Существует несколько различных алгоритмов БПФ, но не будем изобретать велосипед (при желании исходные тексты БПФ на распространенных языках программирования можно найти в Интернете).

Спектр сигнала, а точнее, спектральная плотность мощности (power spectral density), вычисленная напрямую по формуле ДПФ с нахождением модуля комплексных чисел, на самом деле является лишь оценкой, которая наиболее близка к истинному спектру, если анализируются стационарные сигналы (и еще при ряде обстоятельств). Четко определить, когда сигнал не стационарен, - та еще задачка. Уточнение стабильности статистических характеристик сигнала при переходе от выборки (блока) к выборке упирается в продолжительность выборки. В общем случае речевые и музыкальные сигналы не стационарны. Оценка спектральной плотности мощности может быть выполнена и косвенно: через предварительное вычисление автокорреляционной функции с ее последующим дискретным преобразованием Фурье. Для быстрого вычисления автокорреляционной функции нужны могучие компьютерные ресурсы, поскольку, в отличие от БПФ, в принципе не представляется возможным сократить количество умножений (с накоплением суммы) при данных вычислениях. Результаты прямой и косвенной оценок могут существенно различаться. При косвенной оценке частотные составляющие, обусловленные случайными процессами, будут минимизированы. Как прямые, так и косвенные оценки не состоятельны с точки зрения статистики, то есть их дисперсия не стремится к нулю при увеличении объема выборки. Состоятельная оценка получается усреднением нескольких вычисленных спектров. Поэтому при измерениях абсолютных величин (например, мощности звука на отдельных частотах), выполняемых на основе ДПФ, необходимо быть предельно внимательным.

По большому счету, результат ДПФ (или БПФ) - всего лишь математическая абстракция, заключающаяся в разложении на синусы и косинусы исходного сигнала. Вырезал блок из сигнала - получил спектр, сместил вырезку блока во времени - получил следующий спектр, и т. д. Чтобы оценить величину компонента на любой из частот, необходимо сначала накопить данные в размере одного блока. Следует отметить, что система слуха человека выполняет спектральный анализ совершенно иначе. Ее упрощенная схема - набор узкополосных фильтров с параллельной обработкой во времени выходов этих фильтров (то есть в любой момент времени оценивается вклад на каждой из частот) посредством осцилляторов, связанных между собой если и не по принципу "каждый с каждым", то очень хитрым образом.

Для графического отображения во времени сразу всех вычисленных спектров используется трехмерное представление в виде поверхности или двухмерное - как на физической карте, когда пики и впадины выделяются по контуру цветной закраской (рис. 4б). Наиболее распространенное название такого представления амплитудно-частотных спектров - спектрограмма (не путать с периодограммой).

Не надо удивляться появлению странных пиков в БПФ-спектре (например, боковых лепестков) или расползанию пика, соответствующего заведомо чистому тону. Казалось бы, всего-то параметров: объем блока, форма окна (функция, взвешивающая исходный сигнал в пределах блока) и степень наложения (overlapping, %) окон во времени. Тем не менее, оптимальный выбор сочетания упомянутых параметров зачастую требует удовлетворения взаимоисключающим требованиям. Компромисс на компромиссе, особенно для нестационарных сигналов. Иногда имеет смысл пойти окольными путями. Например, выясняя спектральные нюансы в области низких частот, целесообразно проанализировать некий вспомогательный сигнал, выполнив передискретизацию (resampling) исходного сигнала с применением фильтра НЧ (anti-alias), чтобы устранить нестационарности, вызываемые высокочастотными компонентами (фильтр необходим, чтобы исключить мнимые частотные компоненты, которые могут появиться в низкочастотном вспомогательном сигнале как "отражение" высокочастотных компонентов по аналогии со стробоскопическим эффектом). Данная процедура предусматривается в серьезном программном обеспечении (далее - "софте", от английского software, да простят меня борцы за чистоту русского языка).

SoundForge версии 4.5 (Sonic Foundry, Inc.) - безусловно, софт, хороший во многих отношениях, включая недетский спектральный анализ. Но и цена (499 долларов) хороша. Версия (XP4.5) без спектрального анализа и еще ряда примочек стоит уже 59,95 доллара. Демо-версии, имеющие полный набор функций за исключением сохранения обработанных сигналов на диск, можно скачать с ftp.sonicfoundry.com/demo. Имеется одно досадное ограничение: число спектров, вычисляемых за один подход, не может превышать 64. Поэтому общая продолжительность фрагмента сигнала, для которого вычисляется спектрограмма (здесь она называется sonogram), также ограничена. Зато длительность временного окна (одного блока для вычисления спектра) может меняться по желанию пользователя от 128 до 65536 дискрет, что позволяет, соответственно, варьировать частотное разрешение от ~344,5 до ~0,67 Гц на одну точку спектра в полосе частот до 22050 Гц. По умолчанию блок будет состоять из 2048 дискрет, что дает 1025 точек спектра. Возможности выбора типа функции (аж из шести штук) временного окна также впечатляют, но на практике часто достаточно окна Хэннинга (Hanning), которое представляет собой квадрат полупериода синуса единичной амплитуды. Временные окна могут накладываться друг на друга во времени с большим выбором диапазона наложения (до 99%), что очень полезно при анализе нестационарных сигналов. Таким образом, при отсутствии наложения окон, размере блока 65536 дискрет и частоте дискретизации 44100 Гц максимальная длительность анализируемого фрагмента будет равна (1/44100)х65536х64=~95 с. Естественно, при уменьшении частоты дискретизации длительность анализируемого фрагмента возрастет. Чтобы проанализировать следующий фрагмент, необходимо выставить маркер в нужном месте.

Наконец, три существенных недостатка. Во-первых, сохранение вычисленных спектров на жесткий диск не предусматривается, а для сравнения спектров разных сигналов особо страждущие могут выписать на бумажку значения магнитуды, изящно подсказываемые софтом вслед за движением курсора по спектру. Во-вторых, нет возможности проследить изменение выбранной частотной составляющей во времени. В-третьих, прорисовка более чем 4-8 спектров (не говоря уже о 64) дает на экране такую кашу, что извлечь полезную информацию практически невозможно: остается только пытаться рассмотреть интересующий вас спектр на фоне остальных, что, мягко говоря, не совсем удобно.

WaveLab версии 2.0 (www.steinberg.de) дает профессиональную по качеству обработку звука, но спектральный анализ крайне ограничен, можно сказать, присутствует в качестве красивой игрушки, да и графическое отображение спектрограмм нельзя назвать удачным. Зато в этом софте предусмотрен захват треков с аудиодисков и их сохранение в формате WAV. Правда, в руководстве честно признается, что успешность этой процедуры гарантируется при использовании Plextor (SCSI). Что ж, все строчилы Plextor - убойные аппараты. А вот с HP7200i (IDE) WaveLab работать наотрез отказался, хотя и увидел его. Тогда как прилагаемый к HP в комплекте Adaptec Easy CD Creator добросовестно вылавливал треки, споткнувшись только несколько раз на подозрительно уплотненных сидюках.

В настоящее время каждая уважающая себя компания включает вычисление БПФ спектров в sound-редакторы и прочие программы, работающие со звуком. Сервис-опция получает название спектрального анализа. Однако с точки зрения использования спектрального анализа на практике проку от такого сервиса немного. Что-то мелькает, столбики прыгают, спектрограммки буйно красочным ковром расстилаются. Красота, да и только...

Честно говоря, нахожусь в сильном затруднении по поводу совета, какой софт лучше всего использовать для спектрального анализа, поскольку уже много лет пользуюсь самопальной программулей, полностью удовлетворяющей мои потребности. По долгу службы имею дело с профессиональным софтом, но он слишком дорог и повязан с конкретным "железом", например, с DSP-процессорами. Сильно уважаю датскую фирму Bruel&Kjaer, особенно за великолепную техническую литературу в дополнение к высококачественным спектральным анализаторам. Появляющиеся по миру многочисленные программки не лишены недостатков, отмеченных выше у софта известных производителей. Хорошо продумано графическое отображение спектрограмм в DartPro32 (www.dartpro.com). В реальном времени вычисляет (по входному сигналу или готовому WAV-файлу) и прорисовывает бегущую спектрограмму с высоким разрешением компактный и бесплатный (!) Spectrogram (www.monumental.com/rshorne/gram.html).

Кстати, в математическом программном обеспечении, например, MathCad (очень удобный интерфейс, формулы как будто пишешь на доске) и Matlab, есть прямое и обратное БПФ, циклы и логические операторы, векторно-матричное исчисление и развитая графика. Так что, начитавшись книжек по цифровой обработке, с сигналом можно делать все, что в голову взбредет. Недостаток MathCad: придется сначала конвертировать из формата WAV в ASCII, так что с большими файлами не поработаешь напрямую, без использования собственных DLL, созданных на С или С++.

Дабы по несколько раз "не наступать на грабли", желающие поэкспериментировать со спектральным анализом могли бы воспользоваться следующими нехитрыми правилами.

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

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

- Объем блока в 2048 дискрет минимален для полосы частот до 20 кГц. Увеличением объема блока для улучшения разрешения по частоте следует пользоваться очень осторожно, поскольку резко возрастает вероятность нахватать коварно нестационарных изменений сигнала.

- Применение функции окна типа Хэннинга - обязательно. По мере накопления опыта станет ясно, когда окно можно не применять.

- Для получения спектрограммы начинать следует с 50-процентного наложения окон. То есть, например, при объеме блока 2048 дискрет, каждый последующий блок начинается с 1024-й дискреты предыдущего блока.

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

Рис.

Для спектрального анализа с высоким разрешением коротких импульсных сигналов иногда целесообразно предварительно построить авторегрессионную модель, с помощью которой "продлить" импульс и посчитать спектр уже от искусственного сигнала. Кстати, в DartPro32 это предусмотрено, правда, порядок модели сильно ограничен. Такой подход анализа нестационарных сигналов требует соответствующих вычислительных ресурсов и приемлем, пожалуй, только для подготовленных и экипированных экспертов.

В принципе, ничто не мешает раскладывать анализируемый сигнал не на синусы-косинусы, а используя другие функции, хоть полиномы, - главное, чтобы они были ортогональны. Подбирая определенным образом функции разложения, можно получить такие преобразования сигнала, которые обладают требуемыми свойствами. Основной недостаток ДПФ, как и других дискретно-блоковых преобразований из временной области в частотную, наследуется из неопределенности сродни принципу неопределенности Гейзенберга в квантовой механике. Ну невозможно одновременно точно определить значения частоты и времени, на какие функции ни раскладывай. Чем выше разрешение по частоте, тем хуже разрешение по времени, и наоборот. Тем не менее, оказывается, можно получить оптимальное соотношение разрешения по частоте и по времени.

Бум моды на Wavelet-преобразование (можно перевести как "кусочек волны", но более распространено краткое "вэйвлет") пришелся на середину 90-х годов. Для сокращения введем аббревиатуру WT, подразумевая дискретное вэйвлет-преобразование (DWT). Если где-нибудь встретится CWT, то это также означает вэйвлет-преобразование, но вычисляемое несколько иначе (в частности, гораздо медленнее) и обладающее специфическими свойствами. Следует отметить, что WT не является "полноценным" преобразованием из временной области в частотную. Чтобы было легче представить, что это такое, скажу, на что похож результат WT: на фильтрацию в октавных полосах частот! Связь между частотной полосой осуществляется через параметр масштаба, который влияет на сжатие/растяжение базовой вэйвлет-функции. Количество частотных полос для дискретного WT определяется как log2(N), где N=2m опять размер блока в дискретах. При N=2048 получим 11 частотных полос. Самая высокочастотная полоса простирается от 1/4 до 1/2 от частоты дискретизации. Предыдущая частотная полоса - от 1/4 до 1/8 от частоты дискретизации и т. д. до самых низких частот (для CWT возможен другой расклад). На рис. 2 представлен синтезированный сигнал, частота которого постоянно возрастает во времени, и результат его WT-преобразования с наиболее сжатым (по отношению к базовому) вэйвлетом, то есть для самой высокочастотной составляющей.

Рис.

Из рис. 2 хорошо видно, что низкочастотные компоненты как бы отфильтровываются. При этом в данном WT-результате присутствуют (хоть и ослабленные) высокочастотные компоненты из соседних частотных полос. Тогда как при получении WT для низкочастотных полос проникновение высокочастотных компонентов значительно слабее. Получаем внешнее сходство с обычными полосовыми фильтрами, АЧХ которых имеет слабую крутизну на высоких частотах и крутую - на низких. Однако это сходство очень обманчиво, поскольку WT осуществляется принципиально иначе по сравнению с той же цифровой фильтрацией. Массив "однополосного" результата WT часто называется коэффициентами WT. Здесь следует пояснить, что для фильтрации сигналов с помощью WT просто обнуляются коэффициенты, соответствующие частотной полосе, в которой присутствуют нежелательные компоненты, а затем выполняется обратное преобразование, что дает отфильтрованный сигнал во времени. Для получения спектральной оценки можно вычислять среднеквадратичное значение по каждой группе коэффициентов, в противном случае результаты WT следует интерпретировать как подобие спектрограммы. Таким образом, само по себе WT - скорее масштабно-временное преобразование с разложением в различных частотных полосах, нежели преобразование из временной области полностью в частотную, как в случае БПФ.

Разложение осуществляется обычно с помощью вэйвлет-функций вида второй производной от функции распределения Гаусса. Впрочем, возможны и другие функции, например, обладающие фрактальными свойствами Daubechies. Более подробно об истории WT, типах базовых функций и областях применения можно прочесть в номере "Компьютерры" "Всплеск" (www.computerra.ru/1998/8). Главная особенность WT в том, что длительность окна функций разложения (вэйвлетов) меняется в процессе преобразования в пределах одного блока (в ДПФ длительность временного окна строго постоянна и равна длительности блока). Окно как бы сжимается во времени при переходе к более высоким частотам. Преимущества WT проявляются в спектральном анализе нестационарных сигналов. Используя WT, можно получить импульсные функции отклика на разных частотах. В отличие от традиционных аналоговых и цифровых фильтров, WT не дает "звона" на импульсные воздействия. Математический аппарат WT не так прост, как может показаться на первый взгляд, а настройка параметров еще более трудоемка, чем для БПФ. Изменив тип базового вэйвлета, можно получить иной результат в зависимости от специфики сигнала. Интерпретация результатов WT по зубам эксперту, который пуд соли съел. Поэтому, наверное, и хорошего программного обеспечения с WT, доступного массам, практически нет (может, кто подскажет?) [2]. Главный принципиальный недостаток WT состоит в том, что на низких частотах получаем хорошее разрешение по частоте, но плохое по времени, тогда как на высоких частотах - наоборот. Наглядное руководство по вэйвлет-преобразованию можно найти по адресу www1.iastate.edu/~rpolikar/wavelets.

Для любителей экзотики следует припомнить так называемые полиспектры, например, биспектры и триспектры, называемые также bicoherence и tricoherence. Тут уже частотных осей две и три штуки соответственно. Например, в случае биспектра вычисление БПФ производится от двухмерной матрицы, полученной через свертку исходного и двух задержанных сигналов. Триспектр, врать не буду, самопальным путем еще ни разу не считал. Чтобы отобразить результат графически, мощность частотного компонента приходится связывать с размером неких шариков. Об отображении изменения во времени и говорить нечего. При вычислении триспектра с приличным разрешением по частоте даже очень крутой "персоналке" придется несладко. Интерпретация результатов подобного спектрального анализа - для экспертов не со слабыми нервами. Применение спектров высших порядков позволяет выявить некоторые особенности, которые не удается получить посредством традиционного спектрального анализа.

Рис.

В качестве примера (рис. 3) позвольте привести "нормальные" спектры сигналов, полученных при воспроизведении одного и того же музыкального фрагмента разными устройствами. В обоих случаях аналоговый сигнал с выхода устройств подавался на скромный 16-битный АЦП. На слух различия трудноуловимы, и сразу вспоминаются "термины" типа прозрачно-ясно-воздушное звучание, сбалансированный тон, вялый бас и т. п. Безусловно, по одной паре спектров, извлеченных из одного фрагмента, корректное заключение сделать невозможно. А чтобы выявить все детали спектрального состава сигналов при использовании обычного мультимедийного компьютера, звуковая карта должна иметь соответствующие характеристики (АЧХ, сигнал/шум, динамический диапазон). Тем не менее, из сравнения спектров видно, что можно ожидать от того или иного устройства на разных частотах.

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

Рис.

На десерт приведу пример, как изменяется во времени одна частотная составляющая (~4156 Гц) при звучании так называемой тарелки. Разглядывая один из спектров (рис. 4а, нижняя часть) через zoom по линейной шкале в полосе частот до 11025 Гц, вряд ли обратишь внимание на невзрачный пичок, затерявшийся среди обилия более мощных пиков.

Звук натуральной тарелки вообще очень сложен для спектрального анализа. С одной стороны, обилие составляющих вплоть до самых высоких частот, а с другой, в отличие от того же фортепиано, полная неразбериха с характерными гармониками, которые замаскированы компонентами, близкими к случайному (окрашенному) шуму. Геометрия тарелки приводит к возникновению множества резонансов с близкими частотами. Гармоники "играют" во времени, активизируясь по очереди и взаимодействуя друг с другом. В окошечке, расположенном в верхней части рис. 4а, представлено изменение в течение трех секунд упомянутого выше частотного компонента. Отчетливо видно, как гармоника осциллирует во времени, как бы нехотя повинуясь затуханию. Не звучит, а завораживает, что и говорить!!!

1 (обратно к тексту) - Чистый синус с частотой 30 Гц даст в спектре два близлежащих пика, один на отметке 21,5 Гц (более мощный), а второй на отметке 43 Гц.

2 (обратно к тексту) - Еще раз рекомендую MATLAB 5.2, Wavelet toolbox и исходники на С для сжатия серых изображений на страничке Geoffrey Davies, Darthmouth College (www.cs.dartmouth.edu/~gdavis/wavelet/wavelet.html). Ничего ориентированного именно на звук не знаю, но много информации есть, например, на странице Amara Graps (www.amara.com/current/wavelet.html). - Леонид Левкович-Маслюк.

Список литературы
[1] Рабинер Л., Гоулд Б. Теория и применение цифровой обработки сигналов. М.: "Мир", 1978.