CC BY
104
ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
2016 Управление, вычислительная техника и информатика № 2 (35)
УДК 51-76
БО!: 10.17223/19988605/35/6
П.Н. Подкур, Н.К. Смоленцев
ВЕЙВЛЕТ-ПАКЕТНОЕ РАЗЛОЖЕНИЕ ЭЭГ НА ОСНОВНЫЕ ЧАСТОТНЫЕ РИТМЫ
Методами вейвлет-анализа производится разложение электроэнцефалограммы на основные частотные ритмы и вычисляются их числовые характеристики. Более подробно изучается высокочастотная часть ЭЭГ. Она разбивается на несколько компонент с хорошей частотной локализацией, для которых находятся их основные частоты, амплитуда колебаний, относительная энергия и другие характеристики.
Ключевые слова: вейвлет-анализ; пакетный вейвлет-анализ; электроэнцефалограмма; ритмы ЭЭГ; высокочастотные компоненты ЭЭГ.
Имеется обширная литература, посвященная методам анализа сигналов электроэнцефалограммы (см., например, [1-3] и библиографии этих работ). При анализе электроэнцефалограммы (ЭЭГ) используются как традиционные статистические методы и методы анализа Фурье, так и более современные методы вейвлет-анализа [1, 2] и преобразования Гильберта-Хуанга [4]. В последнем случае сигнал ЭЭГ раскладывается на эмпирические моды нелинейных и нестационарных колебаний, а затем применяется преобразование Гильберта. В случае вейвлет-анализа распространено использование методов как непрерывного [1, 2, 5], так и дискретного вейвлет-анализа [6, 7]. В работах [8, 9] дискретный вейвлет-анализ успешно использовался для анализа электрокардиограмм. Традиционно выделяют и исследуют следующие основные диапазоны (ритмы) ЭЭГ [3]: Дельта-ритм (0-4 Гц), Тета-ритм (4-8 Гц), Альфа-ритм (8-16 Гц), Бета-ритм (16-31 Гц) и Гамма-ритм - это частоты выше 30 Гц. В последнее время при использовании электроэнцефалографов высокого разрешения (с частотой оцифровки до 20 тыс. Гц) наблюдается интерес и к высокочастотным диапазонам [10, 11]. Например, в работе [10] показано, что повышенная частотная активность некоторых участков мозга в диапазоне от 60 до 100 Гц является предвестником эпилептического приступа и наблюдается также во время медленного сна. Использование вейвлет-анализа сигнала позволяет выделить указанные выше ритмы в виде отдельных компонент сигнала. При этом получается еще несколько высокочастотных компонент, соответствующих частотам от 30 до 50 Гц, от 50 до 75 Гц и от 75 до 150 Гц. Вейвлет-разложение сигнала ЭЭГ на отдельные частотные диапазоны позволяет рассматривать их независимо друг от друга и дает возможность изучать частотные и другие свойства каждой компоненты. В данной работе мы предлагаем новый подход для получения числовых характеристик ЭЭГ на основе многоуровневого вейвлет-разложения и применения преобразования Гильберта, что позволяет выделить новые частотные диапазоны, получить их визуализацию и новые числовые характеристики ЭЭГ.
Основная операция вейвлет-анализа [9] представляет собой разложение £ ^ {Д, А1} изучаемого сигнала £ = {£п} на две компоненты Д = {А^} и А1 = {А^} при помощи некоторых фильтров. Массив А1 представляет сглаженную часть сигнала и называется массивом коэффициентов аппроксимации. Массив Д представляет детали, которыми исходный сигнал £ отличается от его сглаженной части. С точки зрения анализа сигналов ортогональные вейвлеты представляют собой четыре цифровых фильтра {И„}, {£„}, {И*п} , {£*} [9]. Фильтры {^П*}, {g*п} используются для разложения сигнала по формулам
1. Многоуровневое вейвлет-разложение
Результат действия фильтра {И*} представляет низкочастотную аппроксимацию сигнала. Результат
действия фильтра {ё*п} представляет высокочастотную часть сигнала. Фильтры {к„} и {ёп} используются для восстановления сигнала £ = {£п} по формуле
-2 кА1,к + ёп -2 к А*). (2)
При многоуровневом вейвлет-анализе процедура вейвлет-разложения (1) применяется многократно к массивам коэффициентов аппроксимации. Это может быть изображено схематично следующим образом (рис. 1):
£ ^ {Д, А} ^ {Д, ¿>2, Л} ^ - ^ {А, А, • • Ат-1, ¿м, Ам}. Я-* А1 -> А2 ---- -*>АУ
Рис. 1. Многоуровневое вейвлет-разложение сигнала £
Восстановление сигнала производится поэтапно в обратном порядке. Полученные массивы вейвлет-коэффициентов Вь В2, ..., Бм, Ам представляют интерес при анализе сигналов. Однако важнее получить компоненты сигналов ЯесВ\, ЯесВ2, ..., ЯесВ^ и ЯесАм, которые получаются, если применить процедуру восстановления сигнала только по одному набору коэффициентов, считая, что остальные коэффициенты равны нулю. При этом сумма всех таких компонент будет равна исходному сигналу:
£ = ЯесВ1 + ЯесВ2 + ... + ЯесБы + ЯесАы .
При пакетном вейвлет-анализе раскладываются не только коэффициенты аппроксимации А1, А2, ..., Ам, но и детализирующие коэффициенты В1, В2, ..., по тем же формулам (1). В результате получается так называемое дерево разложения (рис. 2), где вершина (0,0) - это исходный сигнал При пакетном разложении получается более широкий набор из массивов коэффициентов, соответ-
ствующих узлам дерева разложения. Узлы дерева обозначаются либо парами (п, 0), (п, 1), ... , (п, 2п- 1), как на рис. 2, либо обычными цифрами по порядку сверху вниз, слева направо, когда исходный сигнал обозначается цифрой 0, коэффициенты первого уровня - 1, 2, а последний ряд коэффициентов нумеру- 1 до 2 - 2.
(3,0) (3,1) (3,2) (3,3) (3,4) (3,5) (3,6) (3,7)
Рис. 2. Схема пакетного разложения
Полное дерево пакетного вейвлет-разложения содержит много коэффициентов, что затрудняет их изучение. Кроме того, некоторые из коэффициентов могут быть малоинформативными. Поэтому на самом деле важно получить не все дерево, а только некоторое поддерево оптимальной величины в смысле числа коэффициентов и их информативности.
2. Вейвлет-анализ сигнала ЭЭГ
Для изучения мы используем оцифрованный 64-канальный сигнал ЭЭГ длительностью 129 с, снятый на электроэнцефалографе высокого разрешения (500 отсчетов в секунду). Сигнал ЭЭГ регистриру-
ется по 64 общепринятым каналам по системе «10^20», рекомендованной Международной федерацией обществ электроэнцефалографии и клинической нейрофизиологии [3]. Часть сигнала ЭЭГ записана с закрытыми глазами пациента. Для каждого из 64 каналов рассматриваются фрагменты сигнала, соответствующие записи сигнала как с открытыми, так и с закрытыми глазами пациента. Для выбранных фрагментов производятся вейвлет-разложение и вычисление характеристик. Вычисления производятся в системе MATLAB [12] с использованием пакета вейвлет-анализа MATLAB Wavelet Toolbox [9]. Функции этого пакета вейвлет-анализа предусматривают корректную обработку граничных значений при действии фильтров путем симметричного продолжения сигнала.
2.1. Выбор вейвлета
Мы будем использовать ортогональный вейвлет Мейера dmey, который получается из вейвлета Мейера [9] бесконечной импульсной характеристики усечением его фильтра до 102 членов. Он имеет носитель на промежутке [0,101] и центральную частоту Fr = 0,6634 Гц. Выбор этого вейвлета объясняется хорошей локализацией частотных спектров компонент сигнала. Дело в том, что данный вейвлет имеет самый широкий частотный спектр среди ортогональных вейвлетов с компактным носителем. В нем в равной степени представлены частоты, находящиеся в достаточно большой окрестности его центральной частоты 0,6634 Гц. Именно поэтому он дает хорошее разложение сигнала на слагаемые, соответствующие определенным полосам частот.
Поскольку частота дискретизации составляет 500 отсчетов в секунду, то максимальная регистрируемая частота сигнала равна 250 Гц. Поэтому при первом уровне разложения вейвлет Мейера будет выделять элементы сигнала с частотами, близкими к центральной частоте первого уровня разложения, равной Fr1 = 0,6634 250 = 165,85 Гц. При втором уровне разложения вейвлет Мейера будет замечать частоты, близкие частоте в два раза меньше Fr2 = 82,93 Гц, а при третьем уровня разложения вейвлет Мейера будет замечать частоты, близкие к 41,46 Гц. В коэффициентах детализации 4-го уровня разложения будут отражаться элементы сигнала с частотами, близкими к 20,73 Гц, для 5-го уровня разложения - 10,37 Гц, для 6-го - 5,18 Гц и для 7-го уровня - 2,59 Гц.
2.2. Вейвлет-разложение
Достаточно сделать разложение сигнала ЭЭГ до 6-го уровня: S ^ {Db D2, ..., D6, A6}. При этом элементы сигнала с частотами в пределах от 0 до 4 Гц (дельта-ритм) будут представлены коэффициентами аппроксимации A6. Для вейвлет-разложения фрагмента сигнала до 6-го уровня используем следующую команду MATLAB: [c,l] = wavedec(Fragment,6,'dmey'). В результате получаем структуру [c,l], которая содержит набор вейвлет-коэффициентов {Db D2, ..., D6, A6}, где Db D2, ..., D6 - коэффициенты деталей и A6 - аппроксимирующие коэффициенты. Восстановление первоначального сигнала производится последовательно в обратном порядке. Если мы применим процедуру восстановления только к одному набору коэффициентов, когда все остальные коэффициенты состоят из нулей, то получим часть сигнала, соответствующую одному набору коэффициентов. Будем называть такую часть компонентой сигнала. Компоненты сигнала, восстановленные только по коэффициентам деталей Db D2, ., D6, будем называть высокочастотными и обозначать RecD\, RecD2, ., RecD6 соответственно. Например, RecD2 -это компонента сигнала, восстановленная по следующему набору вейвлет-коэффициентов: {0, D2, 0, 0, 0, 0, 0}, где 0 означает массив из нулей. Аналогично низкочастотные компоненты RecA1, RecA2, ., RecA6 получаются восстановлением только по одному набору аппроксимирующих коэффициентов. В MATLAB Wavelet Toolbox есть функция wrcoef, которая позволяет восстановить как высокочастотные, так и низкочастотные компоненты сигнала по полученной ранее структуре [c,l] вейвлет-коэффициентов {Db D2, ..., Dn,An } [9]. Тогда фрагмент нашего сигнала раскладывается в сумму следующих компонент:
Fragment = RecDx + RecD2 + RecD3 + RecD4 + RecD5 + RecD6 + RecA6. На рис. 3 показаны графики исходного сигнала (фрагмент первого канала ЭЭГ) и нескольких его компонент. Частотный спектр мощности P = {Pn} каждой компоненты RecD,k находится в MATLAB обыч-
ным образом. Выполняется дискретное преобразование Фурье fft: {RecD^} ^ {cn}, вычисляется квадрат модуля каждого коэффициента Фурье cn и делится на длину L сигнала: Pn =| cn |2 /L .
Проведенные расчеты показали, что частотные спектры компонент сигналов достаточно хорошо локализованы и полученные компоненты представляют основные диапазоны [3] ЭЭГ: RecA6 - это Дельта-ритм (0-4 Гц), RecD6 - Тета-ритм (4-8 Гц), RecD5 - Альфа-ритм (8-16 Гц), RecD4 - Бета-ритм (1631 Гц). Часть сигнала, спектр которого содержит частоты выше 30 Гц, обычно называется Гамма-диапазоном. В полученном разложении Гамма-диапазон полностью представлен тремя компонентами RecD1, RecD2 и RecD3. Частотные спектры компонент RecD2 и RecD3 локализованы в следующих пределах: RecD3 - от 28 до 50 Гц и RecD2 - от 57 до 75 Гц. Назовем эти диапазоны Gamma-2 и Gamma-1 соответственно. Первая компонента разложения RecD1 имеет относительную /2-энергию порядка 10-7, а ее частотный спектр распределен от 75 до 250 Гц. Поэтому она не рассматривается нами как отдельный диапазон сигнала ЭЭГ. На рис. 3 приведены графики исходного сигнала и нескольких его компонент. В результате получаем разложение сигнала в классических терминах диапазонов:
Fragment = Gamma-2 + Gamma-1 + Beta + Alpha + Theta + Delta.
Рис. 3. Графики исходного сигнала и нескольких его компонент
2.3. Числовые характеристики компонент сигнала
Для всех полученных компонент RecD1, RecD2, ..., RecD6 и RecAв представляют интерес следующие числовые характеристики: спектр мощности, стандартные отклонения, дисперсия, ¡2- и ^-нормы и энтропия. Все они легко вычисляются в МЛТЬЛБ. Дополнительно определим еще три числовые характеристики компоненты сигнала.
1. Относительная энергия компоненты сигнала - это отношение Дэнергии компоненты к ^-энергии всего сигнала (^-энергией сигнала £= {£„} мы называем сумму квадратов модулей элементов S„).
Поскольку компоненты сигнала локализованы по частоте, то вполне естественно рассматривать некоторое усредненное, «основное» значение частоты каждой компоненты. Мы определим такую среднюю частоту компоненты сигнала двумя способами.
2. Средняя статистическая частота. Определяется как такая частота, которая делит пополам мощность спектра. Для вычисления этого параметра делается дискретное преобразование Фурье компоненты сигнала RecDi, затем вычисляется мощность спектра Фурье Si = 2П=12| cn |2 при суммировании до половины длины сигнала (поскольку вторая половина спектра мощности симметрична первой) и находится первое значение п, когда сумма первых квадратов становится больше или равной половине всей
суммы:
> Sj / 2. Значение частоты, соответствующее п, назовем средней статистической ча-
стотой /-ой компоненты сигнала ЯгеВ.
3. Средняя мгновенная частота. Определяется на основе преобразования Гильберта. Напомним, что преобразование Гильберта у(?) = Н(х(?)) функции х(?) определяется формулой
, . 1 ^ х(т) , У(?) =-1— Л
Л ? ~Т
если, конечно, этот интеграл существует в смысле главного значения. Одно из основных свойств преобразования Гильберта заключается в том, что Н(Н(х)) = - х. Тогда комплексная функция г(?) = х(?) + ' у(?) является собственной для преобразования Гильберта: Н(г) = -/ г. Запишем комплексную функцию г(?) в
форме г(?) = А(?)г'е(?). Тогда амплитуда А(?) определяется как модуль функции г(?), а мгновенная частота колебаний - формулой ю = где 0(?) = arctg(y/x). В МЛТЬЛБ есть функция г = Ы1ЬеЛ(х) для выполнения преобразования Гильберта дискретного сигнала х, которая вычисляет комплексную функцию г. Тогда формула
instfreq = Fs/(2*pi)*diff(unwrap(angle(z))) дает нам мгновенную частоту сигнала. Применим дискретное преобразование Гильберта г = Ы1ЬеЛ(х) к каждой компоненте разложения и вычислим мгновенную частоту в каждый момент. Вычисления и графики мгновенных частот показывают, что мгновенная частота колебаний меняется незначительно около некоторого среднего значения Будем называть это значение средней мгновенной частотой компоненты. Преобразование Гильберта позволяет также найти модулирующую функцию А(?) = |г|, которая определяет амплитуду колебаний компоненты сигнала. Теперь наша компонента имеет основную гармонику в виде X = А(?)-со8(2л-^0 •?).
Указанные выше числовые характеристики вычислены для 256 фрагментов сигнала ЭЭГ (по 4 фрагмента длительностью 8 с для каждого из 64 каналов) для пациента с открытыми / закрытыми глазами. Вычисления показывают, что эти числовые характеристики существенно отличаются для разных каналов ЭЭГ, что говорит о том, что они регистрируют разную электрическую активность участков мозга. На рис. 4 в качестве примера представлены значения относительных энергий и средних частот для Альфа-диапазона каждого из 64 каналов. Вычисления для случая ЭЭГ, записанной с закрытыми глазами, показывают, что относительная энергия Альфа-диапазона стала почти в два раза меньше, но средние частоты увеличились на 1-2 Гц. При записи ЭЭГ с закрытыми глазами основными являются Бета- и Дельта-ритмы, которые в сумме дают более 60% энергии сигнала.
Рис. 4. Относительная энергия и средние частоты (в Гц) Альфа-диапазона всех 64 каналов ЭЭГ (сплошная линия - с открытыми глазами, штрих-пунктир - с закрытыми глазами)
2
2.4. Разложение компоненты Гамма-1 (ЕееВ3)
Результаты вычислений показывают, что средняя мгновенная и средняя статистическая частоты для высокочастотных компонент RecD2 и RecD3 ведут себя неустойчиво при небольших сдвигах фрагмента сигнала. Это объясняется тем, что существенная часть их частотных спектров мощности распределена достаточно широко: от 28 до 50 Гц и от 57 до 75 Гц соответственно. Эти компоненты требуют дальнейшего разложения на более узкие частотные диапазоны. Для решения этой задачи воспользуемся пакетным вейвлет-разложением. Функция Т = wpdec (Х^,^пате') МЛТЬЛБ производит пакетное разложение уровня N сигнала X, а функция X = wprcoef (Т,р) производит прямое восстановление только по одному набору коэффициентов в узле р, т.е. создает компоненту сигнала, соответствующую выбранному узлу пакетного дерева [9].
Оптимальное пакетное дерево разложения компоненты RecD3 строится по принципу разложения только тех узлов, для которых относительная энергия соответствующих компонент больше или равна 1% от энергии всего сигнала. В частности, на четвертом уровне разложения имеется всего два узла (4, 2) и (4, 3), в которых относительная энергия компонент будет составлять 2,5 и 96,8% от энергии всего RecD3. На пятом уровне - это узлы (5, 5), (5, 6) и (5, 7). На шестом уровне получается разложение RecD3 на пять компонент, соответствующих узлам (6, 10), (6, 12), (6, 13), (6, 14) и (6, 15). Соответствующие этим узлам компоненты разложения RecD3 имеют достаточно хорошо локализованный спектр: узел (6, 10) - частоты от 46 до 50 Гц, узел (6, 12) - частоты от 28 до 35 Гц, узел (6, 13) - частоты от 35 до 39 Гц, узел (6, 14) - частоты от 43 до 48 Гц и узел (6, 15) - частоты от 38 до 43 Гц. Соответствующие этим узлам компоненты хорошо представляют весь сигнал RecD3, поскольку сумма их относительных ^-энергий составляет 0,98 от l 2-энергии RecD3. Получаем разложение первой высокочастотной компоненты на следующие диапазоны:
Gamma-1 = RecDЧ6, 10) + RecDЗXв, 12) + RecDш 13) + RecDЗX6, 14) + RecDЗXв, 15). Относительная энергия этих диапазонов и средние частоты (для первого канала) указаны в табл. 1. При седьмом уровне разложения получается еще более «тонкое» разложение на 9 компонент, соответствующих узлам (7, 20), (7, 24), (7, 25), (7, 26), (7, 27), (7, 28), (7, 29), (7, 30), (7, 31).
Т а б л и ц а 1
Значения относительных энергий и средних частот компонент Гамма-1
Параметры / узлы (6, 10) (6, 12) (6, 13) (6, 14) (6, 15)
Относительная энергия 0,0129 0,4018 0,3020 0,0836 0,1792
Мгновенная частота 47,38 33,64 36,72 44,31 40,49
Средняя статистическая частота 46,875 32,875 37,75 44,25 40,0
Мы видим, что основную часть составляют компоненты RecD3,(в,l2) и RecD3I(6д3). Результаты вычислений относительной энергии и средних частот новых диапазонов показывают их зависимость от выбора канала ЭЭГ, что говорит о том, что они регистрируют разную электрическую активность участков мозга. Наблюдается также зависимость параметров от состояния пациента (с открытыми и закрытыми глазами).
Пакетное разложение компоненты Гамма-2 (RecD2) производится совершенно аналогично. Локализованные по частоте компоненты сигнала RecD2 находятся в узлах (5, 13), (7, 48), (7, 49), (7, 50), (7, 51). Соответствующие диапазоны имеют частоты: диапазон узла (7, 48) - от 60 до 65 Гц, диапазон узла (7, 49) - от 64 до 67 Гц, диапазон узла (7, 51) - от 66 до 68 Гц, диапазон узла (7, 50) - от 68 до 71 Гц, диапазон узла (5, 13) - от 70 до 75 Гц. Найденные компоненты хорошо представляют весь сигнал RecD2, сумма их относительных ^-энергий составляет 0,96 от ^-энергии всей компоненты RecD2. В табл. 2 представлены относительная энергия и средние частоты компонент RecD2(5, 13), RecD2y(^y 48), RecD2,(^, 49), RecD2^^y 50) и RecD2^^y для первого канала ЭЭГ. Более детальные результаты получаются при дальнейшем разложении в узле (5, 13). Тогда дополнительно к указанным выше компонентам добавляются (вместо (5, 13)) следующие: диапазон узла (7, 52) - от 76 до 80 Гц, диапазон узла (7, 53) - от 74 до 76 Гц, диапазон узла (7, 54) - от 69 до 73 Гц и диапазон узла (7, 55) - от 72 до 75 Гц.
Т а б л и ц а 2
Значения относительных энергий и средних частот компонент Гамма-2
Параметры / узлы (5, 13) (7, 48) (7, 49) (7, 50) (7,51)
Относительная энергия 0,1224 0,1764 0,2433 0,1710 0,2501
Мгновенная частота 72,2201 63,5380 65,2604 69,3311 67,4535
Средняя статистическая частота 71,875 62,875 65,125 69,125 67,625
Заключение
В данной работе методами вейвлет-анализа производится разложение ЭЭГ на частотные диапазоны, которые включают как классические ритмы ЭЭГ, так и ряд новых высокочастотных ритмов, которые получаются пакетным вейвлет-разложением. Для полученных вейвлет-компонент определены новые числовые характеристики, такие как относительная энергия и средние частоты двух типов. Вейвлет-разложения и вычисления указанных параметров проведены для 256 фрагментов сигнала ЭЭГ (по 4 фрагмента длительностью 8 с для каждого из 64 каналов ЭЭГ) для пациента как с открытыми, как и с закрытыми глазами. Вычисления показывают, что эти параметры существенно отличаются для разных каналов ЭЭГ, что говорит о том, что они регистрируют разную электрическую активность участков мозга. Показано также, что они зависят от условий регистрации ЭЭГ (с открытыми или закрытыми глазами). Введенные параметры, относительная энергия, средняя статистическая частота и средняя мгновенная частота компонент сигнала могут использоваться при автоматизированной обработке ЭЭГ.
ЛИТЕРАТУРА
1. Aldroubi A., Unser M. Wavelets in Medicine and Biology. CRC Press, 1996. 640 p.
2. Павлов А.Н. и др. Вейвлет-анализ в нейродинамике // УФН. 2012. Т. 182, № 9. С. 905-939.
3. Гнездицкий В.В. Обратная задача ЭЭГ и клиническая электроэнцефалография (картирование и локализация источников
электрической активности мозга). М. : МЕДпресс-информ, 2004. 624 с.
4. Swarnalatha R., Prasad D.V. Detection of Sleep Bruxism Based on EEG Hilbert Huang Transform // 5th International Conference on
Biomedical Engineering and Technology (ICBET). IPCBEE, 2015. V. 81. P. 33-39.
5. Габова А.В. и др. Использование вейвлет-преобразований для анализа электрической активности мозга при болезни Пар-
кинсона // Нервные болезни. 2012. № 3. С. 1-6. URL: http://www.atmosphere-ph.ru/modules/Magazines/articles/ner-vo/an_ 3_2012_02.pdf
6. Siddiqi A.H. et al. Spectral Analysis of Eeg Signals by using Wavelet and Harmonic Transforms // Istanbul Aydin Universitesi
Dergisi. 2014. V. 3(9). P. 1-20.
7. Gajic D., Djurovic Z., Di Gennaro S. and Gustafsson F. Classification of EEG signals for detection of epileptic seizures based on
wavelets and statistical pattern recognition // Biomedical Engineering: Applications, Basis and Communications. 2014. V. 26, No. 2. 1450021.
8. Подкур П.Н. О высокочастотных компонентах кардиосигнала // Седьмая Всероссийская научно-практическая конференция
«Новые достижения в развитии электрокардиографии» ГУ НИИ кардиологии Томского научного центра СО РАМН Тюменский кардиологический центр. Тюмень, 2005. С. 123-126.
9. Смоленцев Н.К. Основы теории вейвлетов. Вейвлеты в MATLAB. М. : ДМК Пресс, 2013. 628 с.
10. Worrell G.A. et al. High frequency oscillations and seizure generation in neocortical epilepsy // Brain. 2004. Vol. 127. P. 14961506. URL: http://brain.oxfordjournals.org/content/127/7/1496
11. Grenier F., Timofeev .I, Steriade M. Focal synchronization of ripples (80-200 Hz) in neocortex and their neuronal correlates // J. Neurophysiology. 2001. V. 86, No. 4. P. 1884-1898.
12. Смоленцев Н.К. MATLAB: программирование на Visual C#, Borland JBuilder, VBA: учебный курс. М. : ДМК Пресс ; СПб. : Питер, 2009. 456 с.
Подкур Полина Николаевна, канд. физ.-мат. наук, доцент. E-mail: paulina.podkur@gmail.com Российский экономический университет им. Г.В. Плеханова, Кемеровский институт (филиал) Смоленцев Николай Константинович, д-р. физ.-мат. наук, профессор. E-mail: smolennk@mail.ru
Кемеровский государственный университет Поступила в редакцию 29 января 2016 г.
Podkur Polina N., Smolentsev Nikolai K. (Kemerovo Institute (branch) of Plekhanov Russian University of Economics, Kemerovo State University, Russian Federation).
Wavelet packet decomposition EEG on the basic frequency rhythms.
Keywords: wavelet analysis; wavelet packet analysis; electroencephalogram; rhythms of EEG; high-frequency components EEG. DOI: 10.17223/19988605/35/6
This work is devoted to the study of signal electroencephalogram (EEG) by wavelet analysis methods. Traditionally for study of EEG, statistical methods and methods of the analysis Fourier are used. It is considered that the upper boundary frequency of normal EEG, noticeably influencing on its form does not exceed 30 Hz. The fact that higher frequencies are not visible, and in the Fourier spectrum of all EEG signal, the power spectrum of low-frequency part of EEG signal far exceeds and overrides the power spectrum of high frequency oscillations therefore the high frequency ranges are also practically closed for study. Wavelet analysis allows decomposing the entire EEG signal into separate components, which include both traditional rhythms: Delta, Theta, Alpha, Beta and Gamma and series of new frequency localized rhythms. This allows to consider them independently of each other and makes it possible to study separately the frequency and other properties of each component, including high frequency which poorly known at present. In this paper, we propose a new approach for obtaining numerical characteristics of the EEG based on the multilevel wavelet package decomposition and applying the Hilbert transformation, which allows to clearly allocate new frequency bands, to obtain visualization, and new numerical characteristics of the EEG. It is shown that by the wavelet decomposition up to the 6th level with the Meyer wavelet, the EEG signal is decomposed into wavelet components
EEG = RecDi + RecD2 + RecD3 + RecD4 + RecD5 + RecD6 + RecA6, where RecD4, RecD5, RecD6, and RecA6 are classic ranges of Beta, Alpha, Theta and Delta, respectively, and the other RecD1, RecD2 , and RecD3 represent the part of the signal, which is usually called Gamma range. In this case, the range RecD1 has a relative l2-energy of about 10-7, and its frequency spectrum is allocated from 75 to 250 Hz, therefore RecDi not considered in this work as a separate range of EEG signal. Frequency spectra of the components RecD2 and RecD3 are localized: RecD3 - from 28 Hz to 50 Hz and RecD2 - from 57 Hz to 75 Hz. In the work these ranges are named Gamma-2 and Gamma-1, respectively. In addition to the common numeric characteristics of all the bands, in this work three characteristics are defined: the relative energy of component, the average instantaneous frequency and the average statistical frequency of component. The calculation of these parameters for 256 fragments of EEG signal are conducted (in 4 fragments with a duration of 8 seconds, for each of the 64 channels of EEG) for the patient as open as with eyes closed. Calculations show that these numerical characteristics differ significantly for different EEG channels, this suggests that these parameters record the electrical activity of different brain regions. It is also shown that they depend on the recording conditions of EEG (with open or closed eyes). It turned out that the average instantaneous and average statistical frequency for the high-frequency components RecD2 and RecD3 behave unstably in case of small shifts of the fragment of the signal. The reason is that a significant part of their frequency spectra of power although is localized, but is distributed widely: from 28 Hz to 50 Hz or from 57 Hz to 75 Hz, respectively. These EEG components need further decomposition into narrower frequency ranges. For better frequency localization, these components are decomposed into several parts using the wavelet packet decomposition. In each case, the packet nodes of the tree decomposition are found, which allow to decompose RecD2 and RecD3 on new bands with good frequency localization:
RecD3 = RecD3,(6, 10) + RecD3(6, 12) + RecD3(6, 13) + RecD3(6, 14) + RecD3(6, 15), RecD2 = RecD2,(5, 13) + RecD2,a, 48) + RecD2,(7, 49) + RecD2,(7, 50) + RecD2,(7, 51) . The relative energy and average frequencies are calculated for obtained ranges. Thus, in this work, additionally to four classical EEG rhythms 10 bands with good frequency localization are found, which allow to determine their frequency characteristics.
REFERENCES
1. Aldroubi, A. & Unser, M. (1996) Wavelets in Medicine and Biology. CRC Press.
2. Pavlov, A.N. Hramov, A.E., Koronovskii, A.A., Sitnikova, E.Yu., Makarov, V.A. & A.A. Ovchinnikov (2012) Wavelet analysis in
neurodynamics. Uspekhi fizicheskikh nauk - Advances in Physical Sciences. 182(9). pp. 905-939. (In Russian).
3. Gnezditskiy, V.V. (2004) Obratnaya zadacha EEG i klinicheskaya elektroentsefalografiya (kartirovanie i lokalizatsiya istochnikov
elektricheskoy aktivnosti mozga) [Converse task of EEG and clinical electroencephalography (mapping and localization of sources of brain electrical activity)]. Moscow: MEDpress-inform.
4. Swarnalatha, R. & Prasad, D.V. (2015) Detection of Sleep Bruxism Based on EEG Hilbert Huang Transform. 5th International Conference on Biomedical Engineering and Technology (ICBET). 81. pp. 33-39. DOI: 10.7763/IPCBEE. 2015. V81. 7
5. Gabova, A.V., Gnezditskiy, V.V., Karabanov, A.V. et al. (2012) Ispol'zovanie veyvlet-preobrazovaniy dlya analiza elektricheskoy aktivnosti
mozga pri bolezni Parkinsona [The use of wavelet transforms for analysis of electrical brain activity in Parkinson's disease]. Nervnye bolezni. 3. pp. 1-6. [Online] Available from: http://www.atmosphere-ph.ru/modules/Magazines/aiticles/nervo/an_3_2012_02.pdf!
6. Siddiqi, A.H. et al. (2014) Spectral Analysis of Eeg Signals by using Wavelet and Harmonic Transforms. istanbul Aydin Universitesi
Dergisi. 3(9). pp. 1-20.
7. Gajic, D., Djurovic, Z., Di Gennaro, S. & Gustafsson, F. (2014) Classification of EEG signals for detection of epileptic seizures based
on wavelets and statistical pattern recognition. Biomedical Engineering: Applications, Basis and Communications. 26(2). DOI: 10.4015/S1016237214500215.
8. Podkur, P.N. (2005) [About high frequency components of ECG]. Novye dostizheniya v razvitii elektrokardiografii [New achieve-
ments in the development of electrocardiography]. The Seventh All-Russian Scientific-Practical Conference, State Scientific Research Institute of Cardiology Tomsk Scientific Centre SB RAMS, Tyumen Cardiology Center. Tyumen. pp. 123-126. (In Russian).
9. Smolentsev, N.K. (2013) Osnovy teorii veyvletov. Veyvlety v MATLAB [Fundamentals of the theory of wavelets]. Moscow: DMK Press.
10. Worrell, G.A., Parish, L., Cranstoun, S. & Litt, B. (2004) High frequency oscillations and seizure generation in neocortical epilepsy. Brain. 127. P. 1496-1506. DOI: 10.1093/brain/awh149
11. Grenier, F., Timofeev, I. & Steriade, M. (2001) Focal synchronization of ripples (80-200 Hz) in neocortex and their neuronal correlates. Journal of Neurophysiology. 86(4). pp. 1884-1898.
12. Smolentsev, N.K. (2009) MATLAB: programmirovanie na Visual C#, Borland JBuilder, VBA [MATLAB: programming in Visual C#, Borland JBuilder, VBA]. Moscow: DMK Press; St. Petersburg: Piter.
