Wavelet packet decomposition EEG on the basic frequency rhythms.pdf Имеется обширная литература, посвященная методам анализа сигналов электроэнцефалограммы (см., например, [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 Гц. Вейвлет-разложение сигнала ЭЭГ на отдельные частотные диапазоны позволяет рассматривать их независимо друг от друга и дает возможность изучать частотные и другие свойства каждой компоненты. В данной работе мы предлагаем новый подход для получения числовых характеристик ЭЭГ на основе многоуровневого вейвлет-разложения и применения преобразования Гильберта, что позволяет выделить новые частотные диапазоны, получить их визуализацию и новые числовые характеристики ЭЭГ. 1. Многоуровневое вейвлет-разложение Основная операция вейвлет-анализа [9] представляет собой разложение S ^ {Dj, A1} изучаемого сигнала S = {Sn} на две компоненты Di = {D1k} и Ai = {Ai,k} при помощи некоторых фильтров. Массив Ai представляет сглаженную часть сигнала и называется массивом коэффициентов аппроксимации. Массив D1 представляет детали, которыми исходный сигнал S отличается от его сглаженной части. С точки зрения анализа сигналов ортогональные вейвлеты представляют собой четыре цифровых фильтра {hn}, {gn}, {h*n} , {g*} [9]. Фильтры {h*}, {g*n} используются для разложения сигнала по формулам Результат действия фильтра {И*} представляет низкочастотную аппроксимацию сигнала. Результат действия фильтра {g*n} представляет высокочастотную часть сигнала. Фильтры {hn} и {gn} используются для восстановления сигнала S = {Sn} по формуле Sn = ^k (hn -2 kA1,k + gn -2 k Dk). (2) При многоуровневом вейвлет-анализе процедура вейвлет-разложения (1) применяется многократно к массивам коэффициентов аппроксимации. Это может быть изображено схематично следующим образом (рис. 1): S ^ {Di, Ai} ^ {Di, D2, A2} ^ •■■ ^ {Di, D2, DN-i, Dn, AN}. S-* Aj -> A2 ---- -*~AW Рис. 1. Многоуровневое вейвлет-разложение сигнала S Восстановление сигнала производится поэтапно в обратном порядке. Полученные массивы вейвлет-коэффициентов D1, D2, ..., DN, AN представляют интерес при анализе сигналов. Однако важнее получить компоненты сигналов RecD1, RecD2, ..., RecDNи RecAN, которые получаются, если применить процедуру восстановления сигнала только по одному набору коэффициентов, считая, что остальные коэффициенты равны нулю. При этом сумма всех таких компонент будет равна исходному сигналу: S = RecD1 + RecD2 + ... + RecDN + RecAN. При пакетном вейвлет-анализе раскладываются не только коэффициенты аппроксимации A1, A2, ..., An, но и детализирующие коэффициенты D1, D2, ..., DNпо тем же формулам (1). В результате получается так называемое дерево разложения (рис. 2), где вершина (0,0) - это исходный сигнал S. При пакетном разложении получается более широкий набор из 2N+1-1 массивов коэффициентов, соответствующих узлам дерева разложения. Узлы дерева обозначаются либо парами (n, 0), (n, 1), ... , (n, 2n - 1), как на рис. 2, либо обычными цифрами по порядку сверху вниз, слева направо, когда исходный сигнал обозначается цифрой 0, коэффициенты первого уровня - 1, 2, а последний ряд коэффициентов нумеру- 2N л ^N+1 ^ - 1 до 2 - 2. f3,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: {RecDuk} ^ {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. Числовые характеристики компонент сигнала Для всех полученных компонент RecD\, RecD2, ..., RecD6 и RecA6 представляют интерес следующие числовые характеристики: спектр мощности, стандартные отклонения, дисперсия, l2- и ^-нормы и энтропия. Все они легко вычисляются в MATLAB. Дополнительно определим еще три числовые характеристики компоненты сигнала. 1. Относительная энергия компоненты сигнала - это отношение Дэнергии компоненты к ^-энергии всего сигнала (^-энергией сигнала S = {Sn} мы называем сумму квадратов модулей элементов Sn). Поскольку компоненты сигнала локализованы по частоте, то вполне естественно рассматривать некоторое усредненное, «основное» значение частоты каждой компоненты. Мы определим такую среднюю частоту компоненты сигнала двумя способами. 2. Средняя статистическая частота. Определяется как такая частота, которая делит пополам мощность спектра. Для вычисления этого параметра делается дискретное преобразование Фурье компоненты сигнала RecD, затем вычисляется мощность спектра Фурье Si =2П=12\ cn \2 при суммировании до половины длины сигнала (поскольку вторая половина спектра мощности симметрична первой) и находится первое значение n, когда сумма первых квадратов становится больше или равной половине всей > S; / 2. Значение частоты, соответствующее n, назовем средней статистической ча 2 суммы: стотой i-ой компоненты сигнала RecDi. 3. Средняя мгновенная частота. Определяется на основе преобразования Гильберта. Напомним, что преобразование Гильберта y(t) = H(x(t)) функции x(t) определяется формулой , . 1 ^ x(x) , y(t) =- J- dx, Л t -X если, конечно, этот интеграл существует в смысле главного значения. Одно из основных свойств преобразования Гильберта заключается в том, что H(H(x)) = - x. Тогда комплексная функция z(t) = x(t) + i y(t) является собственной для преобразования Гильберта: H(z) = -i z. Запишем комплексную функцию z(t) в форме z(t) = A(t)e1 e(t). Тогда амплитуда A(t) определяется как модуль функции z(t), а мгновенная частота колебаний - формулой ю = d9/dt, где 0(t) = arctg(y/x). В MATLAB есть функция z = hilbert(x) для выполнения преобразования Гильберта дискретного сигнала x, которая вычисляет комплексную функцию z. Тогда формула instfreq = Fs/(2*pi)*diff(unwrap(angle(z))) дает нам мгновенную частоту сигнала. Применим дискретное преобразование Гильберта z = hilbert(x) к каждой компоненте разложения и вычислим мгновенную частоту в каждый момент. Вычисления и графики мгновенных частот показывают, что мгновенная частота колебаний меняется незначительно около некоторого среднего значения F0. Будем называть это значение F0 средней мгновенной частотой компоненты. Преобразование Гильберта позволяет также найти модулирующую функцию A(t) = |z|, которая определяет амплитуду колебаний компоненты сигнала. Теперь наша компонента имеет основную гармонику в виде X = A(t)•cos(2л•F0 -t). Указанные выше числовые характеристики вычислены для 256 фрагментов сигнала ЭЭГ (по 4 фрагмента длительностью 8 с для каждого из 64 каналов) для пациента с открытыми / закрытыми глазами. Вычисления показывают, что эти числовые характеристики существенно отличаются для разных каналов ЭЭГ, что говорит о том, что они регистрируют разную электрическую активность участков мозга. На рис. 4 в качестве примера представлены значения относительных энергий и средних частот для Альфа-диапазона каждого из 64 каналов. Вычисления для случая ЭЭГ, записанной с закрытыми глазами, показывают, что относительная энергия Альфа-диапазона стала почти в два раза меньше, но средние частоты увеличились на 1-2 Гц. При записи ЭЭГ с закрытыми глазами основными являются Бета- и Дельта-ритмы, которые в сумме дают более 60% энергии сигнала. Рис. 4. Относительная энергия и средние частоты (в Гц) Альфа-диапазона всех 64 каналов ЭЭГ (сплошная линия - с открытыми глазами, штрих-пунктир - с закрытыми глазами) 2.4. Разложение компоненты Гамма-1 (RecD3) Результаты вычислений показывают, что средняя мгновенная и средняя статистическая частоты для высокочастотных компонент RecD2 и RecD3 ведут себя неустойчиво при небольших сдвигах фрагмента сигнала. Это объясняется тем, что существенная часть их частотных спектров мощности распределена достаточно широко: от 28 до 50 Гц и от 57 до 75 Гц соответственно. Эти компоненты требуют дальнейшего разложения на более узкие частотные диапазоны. Для решения этой задачи воспользуемся пакетным вейвлет-разложением. Функция T = wpdec (X,N,'wname') MATLAB производит пакетное разложение уровня N сигнала X, а функция X = wprcoef (T,p) производит прямое восстановление только по одному набору коэффициентов в узле p, т.е. создает компоненту сигнала, соответствующую выбранному узлу пакетного дерева [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, поскольку сумма их относительных /2-энергий составляет 0,98 от l 2-энергии RecD3. Получаем разложение первой высокочастотной компоненты на следующие диапазоны: Gamma-1 = RecD46, 10) + RecD46, 12) + RecD46, 13) + RecD46, 14) + RecD46, 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,(6,12) и RecD3,(ey13). Результаты вычислений относительной энергии и средних частот новых диапазонов показывают их зависимость от выбора канала ЭЭГ, что говорит о том, что они регистрируют разную электрическую активность участков мозга. Наблюдается также зависимость параметров от состояния пациента (с открытыми и закрытыми глазами). Пакетное разложение компоненты Гамма-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, сумма их относительных 12-энергий составляет 0,96 от 12-энергии всей компоненты RecD2. В табл. 2 представлены относительная энергия и средние частоты компонент RecD2,(5, 13), RecD2,(7, 48), RecD2,(7, 49), RecD2,(7, 50) и RecD2,(7, 5i) для первого канала ЭЭГ. Более детальные результаты получаются при дальнейшем разложении в узле (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 каналов ЭЭГ) для пациента как с открытыми, как и с закрытыми глазами. Вычисления показывают, что эти параметры существенно отличаются для разных каналов ЭЭГ, что говорит о том, что они регистрируют разную электрическую активность участков мозга. Показано также, что они зависят от условий регистрации ЭЭГ (с открытыми или закрытыми глазами). Введенные параметры, относительная энергия, средняя статистическая частота и средняя мгновенная частота компонент сигнала могут использоваться при автоматизированной обработке ЭЭГ.

Скачать электронную версию публикации