Диссертация (Система скрининговой ранней диагностики паркинсонизма на основе анализа параметров мониторинга сигналов многоканальных электроэнцефалограмм, электромиограмм и мышечного тремора), страница 12

Соответствующая функция butter вMATLAB возвращает описание фильтра в виде векторов-строк b и a, имеющихдлину n+1 и содержащих коэффициенты полиномов числителя и знаменателяфункции передачи в порядке убывания степеней переменной z:()(1) + (2) −1 + . . . +( + 1) −() ==.()1 + (2) −1 + . . . +( + 1) −(3.2)Частотой среза фильтра Баттерворта называется частота, на котороймодуль коэффициента передачи равен √1 / 2. Функция butter в MATLABиспользует значение частоты среза Wn, нормированное к частоте Найквиста.Таким образом, нормированное значение частоты среза Wn должно лежать впределах от 0 до 1, при этом 1 соответствует частоте Найквиста (то естьполовине частоты дискретизации).

Если параметр Wn является двухэлементным79вектором, Wn = [w1, w2], функция butter производит синтез дискретногополосового фильтра порядка 2*n, полоса пропускания которого лежит впределах от w1 до w2 [121].4.ДлявыделенияогибающейсигналаЭМГприменялосьпреобразование Гильберта, рассмотренное в разделе 2.1.2.Блок № 2.Осуществлялось вейвлет-преобразование с применением материнскойфункцией Морле, с помощью которого строялись вейвлет-спектрограммысигналов многоканальной ЭЭГ (см. Рисунок 2.3), ОЭМГ (см.

Рисунок 2.11) имышечного тремора (см. Рисунок 2.16). Далее вычислялись локальныемаксимумы этих вейвлет-спектрограмм (см. Рисунок 2.4, 2.12, 2.17). Алгоритмпоиска локальных максимумов будет описан ниже.Блок № 3.Вычислялись частотно-временные (см. Рисунок 2.7), динамическиечастотные и интегральные гистограммы (см. Рисунок 2.5, слева, 2.6, слева). Спомощью интегральных гистограмм осуществлялась визуальная частотнаясинхронизация сигналов многоканальной ЭЭГ, ОЭМГ и мышечного тремора(см. Рисунок 2.20, 2.21).Блок № 4.Проводился расчет количественных электрофизиологических признаковPi(n) (см. формулу 2.25) согласно модели, описанной во второй главедиссертационной работы. Также осуществлялся расчет эвклидового расстоянияR с учетом мышечного тремора и без учета мышечного тремора для групппациентов и контроля (см.

формулу 2.26). Отметим, что эвклидово расстояниеR – это количественная мера в пространстве этих электрофизиологическихпризнаков паркинсонизма на ранней стадии.На Рисунке 3.1 показано, что, если эвклидовое расстояние R меньшезаданного порогового значения, то считается, что испытуемый здоров, если жеR больше заданного порогового значения, то предполагается, что у80испытуемого БП, он попадает в группу риска и направляется на болеедетальное обследование к врачам.Алгоритм поиска локальных максимумов на вейвлет-спектрограммезаключается в следующем. Рассмотрим пример поиска локального максимумана вейвлет-спектрограмме (Рисунок 3.3) [122, 123]. Локальным максимумом(всплеском) мы будем называть пик M на вейвлет-спектрограмме, полуширина(на полувысоте) по времени t и полуширина по частоте f которого превышаютнекоторые пороги TH и FH соответственно.

При этом величина порога повремени (TH) зависит от частоты f этого пика и определяется следующимобразом:TH = NP / (2f) ,(3.3)где NP – это коэффициент полуширины локального максимума по времени,опытным путем были выбраны NP = 2, а FH = 1 Гц.Рисунок 3.3. Пример локального максимума ЭЭГ в частотно-временнойобласти. Диаграмма слева – срез вейвлет-спектрограммы по времени (осьабсцисс – время, ось ординат – мощность), диаграмма справа – по частоте (осьабсцисс – частота, ось ординат – мощность). Обозначения: WH – полуширинапика на полувысоте, HH – полувысота пика, TH - величина порога по времени,FH - величина порога по частоте.81Характерной особенностью спектрограммы ЭЭГ является то, чтомощность всплесков в одних частотных диапазонах может значительнопревышать мощность всплесков в других частотных диапазонах.

Более того,всплески могут иметь продолжения (хвосты) в соседних частотных диапазонах.Вследствие этого, классические методы анализа ЭЭГ, основанные на вейвлетах,могут ошибочно относить хвосты всплесков в одних частотных диапазонах кэлектрической активности мозга в соседних частотных диапазонах.Рассмотрим участок вейвлет-спектрограммы пациента на 1-й стадии БП вчастотных диапазонах тета (4-8 Гц) и альфа (8-12 Гц). На спектрограмме (см.Рисунок 3.4) показано, что на рассматриваемом участке есть нескольковсплесков (A1, А2, А3, А4), которые имеют чётко выраженную частотновременную локализацию в альфа-диапазоне.Рисунок 3.4.

Спектрограмма всплесков (A1, А2, А3, А4) в диапазоне альфа навременном отрезке от 12 до 18 секунды. Между всплесками A1, А2, А3, А4наблюдается длинный голубой «хвост» некоторого другого локальногомаксимума В, вершина которого находится в диапазоне тета. Ось абсцисс –время (сек), ось ординат – частота (Гц).82Вместе с тем, между ними наблюдается длинный голубой «хвост»некоторого другого локального максимума В, вершина которого находится вдиапазонетета.Разработанныйалгоритмпредотвращаетошибочноераспознавание хвоста указанного локального максимума В в качествеэлектрической активности в частотном диапазоне альфа. То есть, алгоритмпозволяет выделить всплески с малыми амплитудами даже при наличии рядомс ними всплесков с большими амплитудами.Данныйкомплекспрограммбылапробированприобработкеэкспериментальных данных, полученных при исследованиях пациентов наранней стадии болезни Паркинсона и контрольной группы испытуемых.Получено свидетельство о регистрации пакета программ в ФИПС [123].3.2.

Программная реализация алгоритма анализа мониторинговыхизмеренийсигналовэлектроэнцефалограмм,электромиограммимышечного тремораИспользуя, набор инструментов, встроенных в систему MATLAB, былирешены задачи обработки, анализа особенностей сигналов многоканальнойЭЭГ, ЭМГ и мышечного тремора, фильтрации шумов, удалении выбросов всигналах, использования режекторных фильтров для удаления сетевой наводки.В ходе выполнения работы был построен алгоритм обработки данныхэлектроэнцефалограмм, электромиограмм и мышечного тремора во временнойобласти (Рисунок 3.1), который позволил создать комплекс программ в средеMATLAB, осуществляющий анализ ЭЭГ, ЭМГ и мышечного тремора [122,123]. Данная программа является актуальной с точки зрения автоматизацииобработки и анализа биомедицинской информации.

Ниже будет описаналгоритмглавныхисполняемыхфайловeegemg_all_metrics2mat.mиacs_all_metrics2mat.m программы Peaks_1, в которых вызывается в видефункции основная подпрограмма group2mat.m. Подпрограмма group2mat.m83вызывает в виде функции подпрограмму txt2mat.m, которая в свою очередьвызывает следующие подпрограммы:-plot_signal.m;-preprocess_signal.m;-calc_and_plot_spectrum.m;-calc_and_plot_wavelet.m;-plot_and_save_flashes.m.Висполняемыхфайлахeegemg_all_metrics2mat.mиacs_all_metrics2mat.m задаются все основные параметры программы, выбранныеоптимальным образом, такие как:выбор метода обработки (либо преобразование Фурье, либонепрерывный вейвлет Морле);нижнее значение частоты для тета-ритма, равное 4 Гц;верхнее значение частоты для тета-ритма, равное 8 Гц;нижнее значение частоты для альфа-ритма, равное 8 Гц;верхнее значение частоты для альфа-ритма, равное 12 Гц;нижнее значение частоты для бета-ритма, равное 12 Гц;верхнее значение частоты для бета-ритма, равное 24 Гц;нижнее значение частоты для гамма-ритма, равное 24 Гц;верхнее значение частоты для гамма-ритма, равное 30 Гц;частота дискретизации, равная 500 Гц для файла eegemg_all_metrics2mat.m и 1378 Гц для файла acs_ all_metrics2mat.m;длина окна в минутах для метода удаления выбросов, равная 1минуте;фактор децимации, равный 8;порог нижнего значения частоты для вычисления вейвлет-спектрограмм, равный 1 Гц;шаг по частоте для вычисления вейвлетов, равный 0.1;84порог верхнего значения частоты для вычисления вейвлет-спектрограмм, равный 26 Гц;полуширина частотного окна, равная 1 Гц;число периодов, равное 2;длина окна при вычислении спектра, равная 3;параметр сглаживания спектра, равный 100;размер перекрытия окна, равный 7 / 8;минимальная частота при построении спектра, равная 2 Гц;максимальная частота при построении спектра, равная 25 Гц;параметры для вычисления вейвлетов, равные 1;параметр режекторного фильтра (добротность), равный 35;количество столбцов в гистограммах, равное 25;папка с входными файлами, равная папке либо с правостороннимипациентами, либо с левосторонними пациентами, либо папке с практическиздоровыми испытуемыми.

В каждой из папок хранятся записи оцифрованныхЭЭГ, ЭМГ и мышечного тремора;имена каналов (fp1, fp2, f3, f4, c3, c4, p3, p4, o1, o2, f7, f8, t3, t4, t5,t6, cz, emg1, emg2, emg3, emg4 – для файла eegemg_ all_metrics2mat.m и lh, rh –для файла acs_ all_metrics2mat.m);имена каналов, выбранных пользователем.Подпрограмма group2mat.m выполняет следующие действия:вычисляетномеравыбранныхпользователемканалов(ChannelNamesSelectedByUser) в исходном списке каналов (ChannelNames);выбирает для обработки только те файлы, которые имеютрасширение .txt;вызывает подпрограмму txt2mat.m;формирует таблицу значений для тета (4-8 Гц) диапазона, альфа (8-12 Гц) диапазона, бета (12-24 Гц) диапазона, гамма (24-30 Гц) диапазона,сформированных в подпрограмме txt2mat.m, и затем записывает эти значения в85соответствующие 4 файла с указанием выбранного пользователем методаобработки.Подпрограмма txt2mat.m выполняет следующие действия:для каждого канала, выбранного пользователем, данная функциястроит и выводит на экран сигнал (Рисунок 2.2, 2.10, 2.15), с помощьюподпрограммы plot_signal.m;для каждого канала, выбранного пользователем, происходитпредобработка (Рисунок 3.1, блок № 1), содержащегося в нем сигнала спомощью функции preprocess_signal.m, которая включает себя:o удаление выбросов в сигналах методом Хьюбера Х84;o фильтрацию режекторным фильтром на частотах 50, 100, 150,200 Гц для удаления сетевой наводки на этих частотах;o фильтрацию фильтром Баттерворта 4-го порядка;o преобразование Гильберта для выделения огибающей сигналаЭМГ.для каждого канала, выбранного пользователем, данная функциястроит и выводит на экран спектр сигнала (Рисунок 2.5, справа, 2.6, справа,2.13, справа, 2.14, справа, 2.18, справа, 2.19, справа), с помощью подпрограммыcalc_and_plot_spectrum.m, которая вычисляет значения для диапазонов тета (4-8Гц), альфа (8-12 Гц), бета (12-24 Гц), гамма (24-30 Гц);для каждого канала, выбранного пользователем, данная функциявыполняетсявейвлет-преобразованиесматеринскимвейвлетомМорле(Рисунок 3.1, блок № 2) и выводит на экран вейвлет-спектрограмму (Рисунок2.3, 2.11, 2.16), с помощью подпрограммы calc_and_plot_wavelet.m, котораявычисляет значения всплесков на вейвлет-спектрограммах для диапазонов тета(4-8 Гц), альфа (8-12 Гц), бета (12-24 Гц), гамма (24-30 Гц).для каждого канала, выбранного пользователем, вычисляются,выводятся на экран и сохраняются в файл значения локальных максимумов(всплесков) на частотно-временной плоскости для частоты, времени иамплитуды с помощью подпрограммы plot_and_save_flashes.m.