Спектральный анализ сигналов в базисе вейвлетов Текст научной статьи по специальности «Математика»

ISSN G868-5886

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2GGG, том lG, № 3, c. 7G-77

ОРИГИНАЛЬНЫЕ СТАТЬИ

УДК 621.519 © Л. В. Новиков

СПЕКТРАЛЬНЫЙ АНАЛИЗ СИГНАЛОВ В БАЗИСЕ ВЕЙВЛЕТОВ

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

ВВЕДЕНИЕ

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

Один из них — выделение спектральных компонент с помощью устройств с избирательной амплитудно-частотной характеристикой. Это — полосовые фильтры с шириной полосы пропускания значительно меньшей, чем спектр сигнала. Отношение центральной частоты (настройки) фильтра к его ширине полосы носит название добротности и обозначается как Q .

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

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

Второй способ анализа спектра основывается на Фурье-преобразовании реализации сигнала s(t) длительности Т0 по формуле:

„у

5^, ю) = |s(т)e 0dт . (1)

Чтобы избежать так называемых краевых эффектов, реализацию сигнала умножают на гладкую оконную функцию w(t). В результате имеем «оконное» Фурье-преобразование:

5^,ю) =|м(т -1)5(т)е_ 100 dт, (2)

м Т /

где w(t) = 0 при > у2 .

Из (1) и (2) может быть получен текущий фазовый или амплитудный спектр сигнала как соответственно фаза и модуль комплексной функции 5^,ю).

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

ДП =—. Поэтому вычисления по (1) и (2) вы-

Т0

полняются при дискретных значениях частот, равных / ДП (/ = 0,1,2,...). Величина ДП называется еще разрешающей способностью анализа, так как сигнал из двух гармонических колебаний может быть различим, только если их частоты будут отличаться на величину большую, чем ДП. Сам анализ поэтому носит название анализа с постоянным разрешением в отличие от, например, октавнополосной фильтрации, реализующей анализ с множественным разрешением.

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

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

только гармонических функций в качестве базиса преобразования, но и других ортогональных базисов: ортогональных полиномов, функций Уолша, Радемахера, Крестенсена, Карунена—Лоэва и др. Однако трудности интерпретации получаемых результатов или в ряде случаев невозможность построения быстрых алгоритмов явились причиной того, что все эти базисы не получили широкого распространения.

Появление в начале 80-х годов теории вейвлетов, развитой в работах [1, 2, 3], коренным образом изменило отношение к нетрадиционным базисам. Анализ в базисе вейвлетов соединил в себе достоинства оконного преобразования Фурье (точнее, преобразования Габара как его разновидности), в частности локальность анализа с возможностью масштабирования, т.е. растяжения базисной функции или ее сжатия. Это означает, что вейвлеты являются функциями двух переменных — масштаба и сдвига, а сами функции представляют собой кратковременное колебание (иногда говорят — всплеск).

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

где Ж у Ж у+1... — подпространства сигналов

sk (г), занимающих соответствующую полосу частот, У]о — подпространство сигналов (г), включающее константы и у0 = 0, ± 1, ± 2,... — начальный масштаб.

Определим в У}-0 ортонормированный базис

<р]о п (г), представляющий собой множество смещенных во времени функций:

А

Pjo,n (t) = 2 2 p(2Jo t _ n).

(4)

Тогда в соответствии с концепцией кратномасштабного анализа (КМА) [3, 4] функции

Jo+1

p

Jo+1, n

(t) = 2 2 p(2Jo+11 _ n)

образуют ортонормированный базис в подпространстве Vjo +1 = VJo Є Wj, а функции

Pj n (t) = 2/2 p(2 4 _ n), J є Z

образуют ортонормированный базис в подпространстве VУ = ¥]_1 © Жу_1. Величина у определяет

масштаб, а п — смещение базисных функций. Так как последующие подпространства V]■ включают предыдущие, в частности, V, с У1, то справедливо следующее масштабирующее выражение:

Po,o (t) = p(t) = V2 2 h(n)p(2t _ n)

(5)

1. ТЕОРИЯ

Будем рассматривать спектральный анализ множества 8 сигналов л(ґ), принадлежащих пространству Ь2 квадратично интегрируемых функций. В результате такого анализа мы разбиваем (разделяем) сигнал на отдельные компоненты, такие что

S(t) = SJ0 (t) + 2 Sk (t)

k

(3)

k =1

где s]o (г) — низкочастотная компонента, sk (г) —

высокочастотные компоненты сигнала.

Представление (3) продуктивно рассматривать как разделение Ь 2 на неперекрывающиеся подпространства Жк такие, что

L 2 = V, Є WЄ W, +1 Є....

2 Jo Jo Jo+1

где p(t) называется масштабной функцией, h(n) — коэффициентами масштабной функции (иногда — масштабным фильтром или масштабным вектором).

Определим в W0 ортонормированный базис

V j,n (t)) называемый вейвлетным базисом. Тогда, учитывая, что W0 с Ух, получим для функции ¥oo(t) = V(t) соотношение, аналогичное (5):

W(t) = V2 2 g (n)P(2t _ n),

(б)

где из условия ортогональности подпространств ...Vo 1W01W1 ... коэффициенты

g (n) = (_1) nh(1 _ n).

(7)

Базисные функции для подпространств Жу образуются смещением и масштабированием функции щ(г):

¥у,п (г) = 2 У^(2 Ч _ п).

Из (6) следует, что вейвлеты полностью определяются масштабными функциями (5).

Теперь сигнал э(г) может быть представлен в виде линейной комбинации:

^) = 2су0(к)^^0,k(г) + 22'й](к)¥у к({).

к=1 к =1 1=]

Это означает, что

^ = 2СУ0(к)(] (() (8а)

к

и

^ (г) = 2(к)¥у,к (г), (8б)

к

где коэффициенты определяются из условия ортогональности базисных функций:

су0 (к) =|л(0] (г)Ж, (9а)

йу (к) =|э(г)Щ] к (г)йг. (9б)

Преобразование сигнала по формулам (9) называется анализом, а восстановление по формулам (8) — синтезом, который приводит к вычислению спектральных компонент в выражении (3). Коэффициенты с]0(к) и йу (к) называются спектром

сигнала s(t) соответственно в базисах масштабных функций и вейвлетов.

Подставляя выражения (9) в (8), получим:

э^) =| 5(Т)Ф у0(г _ т)йг (10а)

и

^ (г) =| э(т)Ту (г _ т)йт, (10б)

где функции

фу-0(г) = E^]0(m>P]0(m _ t), (11а)

т

V] (г) = 2^ ] (т)¥ ] (т _ г) (11б)

т

называются ядрами Дирихле интегральных преобразований (10). Они представляют собой импульсные отклики соответствующих фильтров: низкой

частоты, «вырезающего» из сигнала э(г) компоненту Эу (г), и полосовых фильтров, «вырезающих» из сигнала высокочастотные компоненты.

Выполняя преобразование Фурье этих ядер, получим:

Ф ]0 (®) “МЮ)Р

и

'V у (®) “К у (®)|2.

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

В наибольшей степени для этой цели подходят базисные функции, заданные на ограниченном множестве точек (базисы с компактным носителем) и многократно дифференцируемые. При этом компактность базисов минимизирует число вычислительных операций, а к-кратное, например, дифференцирование гарантирует затухание частотной характеристики как к . Здесь мы рассмотрели только базисы Добеши [2, 5], так как исследование всего множества вейвлетных базисов с точки зрения их применения для спектрального анализа требует отдельной работы. Для иллюстрации частотных свойств базисных функций Добеши на рис. 1 и 2 показаны Фурье-образы ядер Дирихле двух функций Добеши йЬ4 и йЬ20. Ниже будет введено понятие порядка вейвлета. Здесь же лишь заметим, что йЬ20 имеет более высокий порядок, чем йЬ4.

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

2. АЛГОРИТМЫ

Уравнения (5) и (6) позволяют получить быстрые алгоритмы для вычислений с у (к) и й у (к).

Выражение для масштабной функции, полученное из (5) для произвольного значения у, имеет вид:

ф],к (г) = 2 2 2Мп)ф(2]+1 г_2к _п) =

п

]+1

= 2 2 2к(п _ 2к)ф(2]+ г_ т).

т

Аналогично из (6) имеем:

]+!

¥],к (О = 2 2 2ё(п _ 2к)ф(2]+' г _ т).

т

Подставляя эти выражения в (9) и выполняя преобразования, получим следующие рекуррентные формулы для вычисления спектров с у (к) и

йу (к), которые получили название анализ от

«тонкого» к «грубому» масштабу:

с у (к) = 2 ^(т _ 2к) с у+Дт), (12а)

т

й] (к) = 2 & (т _ 2к) с у+Дт). (12б)

т

В качестве коэффициентов су (к) на самом

«тонком» значении масштаба у = ]т принимаются /

временные отсчеты сигнала, т.е.

с1т (к) = э(к).

Итерационная процедура (12) обычно завершается при некотором значении у = у0, которое выбирается исходя из априорной информации о сигнале, а именно, о его продолжительности. Напомним, что величина у0 определяет ширину полосы пропускания низкочастотного фильтра и самого первого полосового фильтра вейвлет-разложения. Одна итерация этой процедуры может быть представлена в виде блок-схемы, показанной на рис. 3.

В первых блоках вычисляется свертка (цифровая фильтрация) по формуле:

N-1

У (к) = 2 х(п)Ь(к - п),

п=0

где N — число отсчетов в окне данных.

Так как соответствующие суммы в выражениях (12) не являются свертками, то для вычисления коэффициентов су (к) и йу (к) по (12) методом

цифровой фильтрации аргументы весовых коэффициентов фильтров к и & должны быть взяты с обратным знаком. Во втором блоке выполняется операция прореживания в два раза (из-за множителя 2 при аргументе к в формуле (12)).

Для последующей итерации используются значения су предыдущей и т.д. до у = у0.

Рис. 1. Базисные функции Добеши йЬ4 — а) и Фурье-образы ядер Дирихле — б)

Рис. 2. Базисные функции Добеши йЬ20 — а) и Фурье-образы ядер Дирихле — б)

Рис. 3. Один шаг итерационной процедуры вейвлет-анализа — а) и его условное обозначение — б)

з(к)-

а)

N

А

Л

Ы/ /л

Я/

СІ і - і J■m 1

^т-2

N

с]0+1 А /2 у т у 0 1

Чо

I К/

/2 іт-10 "І

%

з(к)

А

5)

I Г“

^т-2 <1

%

с

1 т -2

Рис. 4. Многошаговая схема вейвлет-анализа — а), условное обозначение — б)

Схема многошаговой итерационной процедуры анализа с использованием обозначения рис. 3, б показана на рис. 4.

В первом блоке производится обработка временных отсчетов сигнала, а в каждом последующем — соответствующих коэффициентов Су . В

результате будем иметь ]т - ]0 -1 векторов вейвлет-коэффициентов ёу и один вектор коэффициентов разложения по масштабной функции — Су .

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

Для вычисления спектральных компонент по формулам (8) с учетом выражений для ф]к (V) и

Рис. 5. Один шаг итерационной процедуры синтеза сигналов по их коэффициентам разложения в вейвлетном базисе — а), условное обозначение — б)

фу к (ґ) может быть получена формула восстановления отсчетов сигнала для каждого масштаба, которая имеет вид [5]:

с у+і (к) = 2 с у (т)к(к - 2т) + ^ (т) g (к - 2т).

т т

Процесс восстановления идет от «грубого» масштаба к «тонкому». Схема блоков фильтров, осуществляющих эту процедуру на одной итерации, показана на рис. 5.

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

Схема многошаговой итерационной процедуры синтеза с использованием обозначения на рис. 5, б приведена на рис. 6, причем для удобства «стыковки» со схемой анализа рис. 4 процесс восстановления изображен справа налево.

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

Аналогично вычисляются все другие спектральные компоненты из формулы (4).

Число операций умножения, необходимое для вычисления всех DWT-коэффициентов для массива данных N и длины векторов к и g, равной Ь, будет 2ЬК [5]. Столько же операций нужно

а)

сі? _? J^n * і I $0 1

5 %-1 5 [о+1 я

и

СЇ

']т-1

!І, _п СІ

2

I

%

сі ЛТП с >1

О

ф)-

б)

Рис. 6. Многошаговая схема вейвлет-синтеза —а), условное обозначение — б)

Чо

Это означает, что вейвлеты позволяют выявлять не только частотные особенности любого сигнала, но и определять время его появления. Кроме того, вследствие масштабирования длительность базисной функции всегда согласована с продолжительностью сигнала. При выборе для выполнения анализа вейвлет-функции нужно учитывать два параметра: длину Ь вектора к и порядок функции К . Для вейвлетов Добеши они связаны соотношением

К = Ь — 1. Чем больше порядок вейвлетов, тем

большее число раз они дифференцируемы. При Ь < 4 они недифференцируемы. При Ь = 6 они дифференцируемы только один раз; при Ь = 14 — дифференцируемы дважды [5]. Степень дифференцируемости определяет крутизну фронтов полосовых фильтров спектрального анализа (см. рис. 1 и 2).

Так как порядок вейвлета К определяет число нулевых моментов, то он определяет также максимальную степень полинома, аппроксимирующего сигнал, для которого коэффициенты вейвлет-разложения будут иметь нулевые значения. В этом случае вся энергия сигнала будет сосредоточена в коэффициентах разложения по масштабной функции, т.е. в су0 (к) (и соответственно в спектральной компоненте Sj0(t) ).

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

Рис. 7. Схема вычисления спектральной компоненты 5- 2(ґ)

-т -2 ' '

выполнить, чтобы восстановить сигнал или вычислить все его спектральные компоненты. Следовательно, для анализа-синтеза сигнала 5(0 в вейвлетном базисе необходимо выполнить 4ZN операций. Напомним, что число операций комплексного умножения для БПФ равно N log2 N, что сравнимо или даже больше, чем в случае БІТ.

3. МОДЕЛИРОВАНИЕ

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

Рис. 8. Две последовательные синусоиды — а) и их разложение на компоненты в вейвлетном базисе ёЬ16 — б)

Рис. 9. Смесь всплеска гармонического сигнала и коррелированного шума — а) и разложение по спектральным компонентам в вейвлетном базисе db8 — б)

коэффициентам: 70, 70 + 1, 70 + 2 и 70 + 3. Четко видна граница раздела двух синусоид, выявленных фильтрами соседних частот настройки (в масштабах 70 + 1 и 70 + 2).

Рис. 9 иллюстрирует эффект выделения сигнала из шума (отношение сигнал/шум = 0.5) при разложении по спектральным компонентам при восьми значениях масштаба (показаны последние четыре).

ВЫВОДЫ

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

СПИСОК ЛИТЕРАТУРЫ

1. Grossman A., Morlet J. Decomposition of Hardy into square integrable wavelets of constant shape // Siam journ. of Math. Analm. 1984. V. 15. P.723-736.

2. Daubechies I. Orthonormal bases of compactly supportid wavelets. // Communications on pure and Applied Mathematics. November. 1988. V. 41. P.909-926.

3. Mallat S.G. Multiresolution approximations and Wavelet of ortonormal bases of L2(R) // Transactions of the American Mathematical Society. September. 1989. V. 315. P. 69-87.

4. Петухов А.П. Введение в теорию базисов всплесков. СПб.: Изд. СПбГТУ, 1999. 131 с.

5. Burrus C.S., Gopinath R.A. and Haital Guo. Introduction to Wavelets and Wavelet Transform.: New Jersey: Prentice Hall., 1998. 268 p.

Институт аналитического приборостроения РАН, Санкт-Петербург

Материал поступил в редакцию 20.04.2000.

WAVELET-BASED SPECTRAL ANALYSIS OF SIGNALS

L. V. Novikov

Institute for Analytical Instrumentation RAS, Saint-Petersburg

Octave-band spectral analysis of signals based on the discrete wavelet transform (DWT) is studied. It is shown that the slope of band-pass filter response fronts increases and the inter-band interference (crosstalk) decreases with increasing wavelet function order. Fast computational algorithms and circuits are presented.