Анализ медико-биологических сигналов ЭЭГ средствами MATLAB

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Целью курсовой работы является получение знаний о программном обеспечении для решения практических задач по анализу медико-биологических сигналов и изображений, использование современных пакетов прикладных программ для анализа медико-биологических сигналов и изображений.                                                                                                                         В данной работе проводится анализ результатов ЭЭГ одного отведения средствами MATLAB, а именно производится статистический и спектральный анализ, а так же цифровая фильтрация сигнала.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 Основные сведения о сигналах ЭЭГ

       Электроэнцефалограмма (ЭЭГ) (электро- + др.-греч. ἐγκέφαλος — головной мозг + …грамма) Г - метод регистрации электрической активности (биопотенциалов) головного мозга через неповрежденные покровы головы, позволяющий судить о его  функциональном состоянии, наличии очаговых поражений, мозговых расстройств и их характере. Одной из основных характеристик ЭЭГ является частота. Однако из-за ограниченных возможностей человека при визуальном анализе ЭЭГ, применяемом в клинической электроэнцефалографии, целый ряд частот не может быть достаточно точно охарактеризован оператором (рис.1 )

 

 

 

Рисунок-1. Сигнал ЭЭГ, записанный в течение 1 с.

 

Ниже приведены основные ритмы, наблюдаемые при записи ЭЭГ при различных состояниях человека. Спонтанные колебания биопотенциалов различаются по частоте: Могут меняться также их амплитуда и форма.

Рисунок - 2. Альфа ритм: частота 8-13 Гц, амплитуда 30-70 мкВ

 

 

 

 

 

          Рисунок -3. Бета ритм - частота 14-30 Гц,  амплитуда 5-30 мкВ

 

             Рисунок - 4. Гамма ритм: частота 30-100 Гц, амплитуда < 10 мкВ

 

Рисунок - 5. Дельта ритм: частота до 4 Гц,  амплитуда 20-30 мкмВ

              Рисунок - 6. Тета ритм: частота 4-8 Гц,  амплитуда 100-150 мкВ

 

 

 

 

 

 

 

 

 

 

 

                        Регистрация ЭЭГ

Для проведения ЭЭГ на голове фиксируются электроды, которые соединяются проводами с электроэнцефалографом. Аппарат усиливает потенциалы, полученные с датчиков в сотни тысяч раз и записывает их в память компьютера.                                                                                              Важное значение при регистрации ЭЭГ имеет расположение электродов, при этом электрическая активность одновременно регистрируемая с различных точек головы может сильно различаться.                                                                       При записи ЭЭГ используют два основных метода (монтажа): биполярный и монополярный. В первом случае оба электрода помещаются в электрически активные точки скальпа, во втором один из электродов располагается в точке, которая условно считается электрически нейтральной (мочка уха, переносица).При биполярной записи регистрируется ЭЭГ, представляющая результат взаимодействия двух электрически активных точек (например, лобного и затылочного отведений), при монополярной записи — активность какого-то одного отведения относительно электрически нейтральной точки (например, лобного или затылочного отведения относительно мочки уха). Выбор того или иного варианта записи зависит от целей исследования. В исследовательской практике шире используется монополярный вариант регистрации, поскольку он позволяет изучать изолированный вклад той или иной зоны мозга в изучаемый процесс. Международная федерация обществ электроэнцефалографии приняла так называемую систему "10-20", позволяющую точно указывать расположение электродов. В соответствии с этой системой у каждого испытуемого точно измеряют расстояние между серединой переносицы (назионом) и твердым костным бугорком на затылке (инионом), а такжемежду левой и правой ушными ямками.                     Возможные точки расположения электродов разделены интервалами, составляющими 10% или 20% этих расстояний на черепе.

 При этом для удобства  регистрации весь череп разбит  на области, обозначенные буквами: F — лобная, О — затылочная область, Р — теменная, Т — височная, С — область центральной борозды. Нечетные номера мест отведения относятся к левому, а четные — к правому полушарию. Буквой Z — обозначается отведение от верхушки черепа. Это место называется вертексом.

 

 

Рисунок- 7. Регистрация  ЭЭГ.

 

 

 

 

 

Рисунок  -8. Условная схема регистрации ЭЭГ при монополярном отведении с референтным электродом на мочке (слева) уха и биполярных отведениях(справа).

В обычной практике ЭЭГ отводят с помощью электродов, расположенных на интактных покровах головы. Электрические потенциалы усиливают и регистрируют. В электроэнцефалографах предусмотрено 16-24 и более идентичных усилительно-регистрирующих блоков (каналов) , позволяющих одномоментно записывать электрическую активность от соответствующего количества пар электродов, установленных на голове пациента. Современные электроэнцефалографы создают на базе компьютеров. Усиленные потенциалы преобразуют в цифровую форму; непрерывная регистрация ЭЭГ отображается на мониторе и одновременно записывается на диск.

После обработки ЭЭГ может быть распечатана на бумаге. Электроды, отводящие потенциалы, представляют собой металлические пластины или стержни различной формы с диаметром контактной поверхности 0,5-1 см. Электрические потенциалы подаются на входную коробку электроэнцефалографа, имеющую 20-40 и более пронумерованных контактных гнёзд, с помощью которых к аппарату можно подсоединить соответствующее количество электродов.

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

 

Рисунок-9 Восходящая ретикуло-кортикальная неспецифическая система регуляции уровня функциональной активности мозга: Д1 и Д2 - десинхронизующие активирующие системы среднего и переднего мозга соответственно; С1 и С2 - синхронизующие тормозящие сомногенные системы продолговатого мозга и моста и неспецифических ядер промежуточного мозга соответственно.

 

 

 

 

Статистический анализ сигналов ЭЭГ

Для выполнения статистического анализа сигналов ЭЭГ использовали следующие функции:

Mean- определение среднего арифметического значения выборки;

Std- расчет среднего квадратичного отклонения выборки;

Var- определение дисперсии выборки;

Min- определение наименьшего значения выборки;

Max- определение наибольшего значения выборки;

Sort- сортировка выборки по возрастанию.

Статистический анализ можно проводить, используя различные математические пакеты или разрабатывать специальные приложения. Для исследовательских работ предпочтение следует отдать математическим системам, таким как электронные таблицы Excel, MATLAB, MATHCAD и др.

При статистическом анализе ВСР определяются следующие основные параметры:

среднее арифметическое значение MeanRR,

 

;        (1)

 

выборочная несмещенная дисперсия VarRR,

;

или

; (2)

выборочное среднее квадратическое отклонение (стандартное отклонение) StdRR,

;        (3)

 

коэффициент эксцесса,

    (4)

коэффициент асимметрии,

;   (5)

functionvarargout=Stat(GG);

load('eeg01_16.txt')

GG=eeg01_16

%Минимальное число  GG интервалов равно 50

%StatGG определяет статистические параметры ЭЭГ сигнала

%Контроль числа  отсчетов

iflength(GG)<50

error('число должно быть больше 50')

end

%Определение среднего  арифметического

MeanGG=mean(GG)

%Определение СКО(стандартное отклонение)

StdGG=std(GG)

%Определение коэффициента  вариации

VarGG=(StdGG/MeanGG)*100

%расчет ассиметри

SkGG=skewness(GG)

%расчет эксцесса

KurtGG=kurtosis(GG)

%Определяем минимальное  значение ЭЭГ

MinGG=min(GG)

%Определяем наибольшее  значение ЭЭГ 

MaxGG=max(GG)

%медиана

MeGG=median(GG)

%размах

RGG=range(GG)

Return

Используя функции SkGG и KurtGG, дополнили текст функции StatGG с целью расчета асимметрии и эксцесса.                                                                  Командой [MeanGG, StdGG, VarGG, SkGG, KurtGG]=StatGG определили статистические характеристики.

MeanGG = 2.4302

StdGG =38.1228

VarGG = 1.5687e+003

SkGG = 0.1640

KurtGG = 3.4411

MinGG = -115

MaxGG = 128

MeGG = 2

RGG =243

Выполнили построение графиков ЭЭГ сигналов:

functionvarargout=PlotGG(GG)

GG=load ('eeg01_16.txt')

%закрытие всех  окон 

closeall

% Построение графика  ЭЭГ сигнала

figure(1)

plot(GG),grid

title('графикЭЭГсигнала')

ylabel('ЭЭГ ,мкВ')

xlabel('НомеротсчетасигналаЭЭГ')

pause

%Построениеритмограммы

figure(2)

bar(GG),grid

title('диаграмма')

ylabel('ЭЭГ,мкВ')

xlabel('НомеротсчетасигналаЭЭГ')

pause

%Построениегистограммы

figure(3)

histfit(GG,10),grid

title('Гистограмма')

 

 

ylabel('частота')

xlabel('ЭЭГ,мкВ')

pause

%Построение скатерограммы

EGnumber=length(GG);

GGx=GG(2:2:GGnumber);

GGy=GG(1:2:GGnumber-1);

figure(4)

plot(GGx,GGy,'r*'),grid

title('Скатерограмма')

ylabel('ЭЭГ,мкВ')

xlabel('ЭЭГ,мкВ')

pause

 

 

 

 

 

 

 

 

 

 

 

 

Построили диаграмму:

 

Построили гистограмму:

 

Построили скатерограмму:

Проанализируем полученную гистограмму:

Определим моду:

X0-нижняя граница модального интервала: -20

(верхняя граница модального  интервала: 8)

Fm0-частота в модальном интервале:345

Fm0-1-частота в предыдущем интервале:225

Fm0+1- частота в интервале следующем за модальным:305

∆-величина интервала  24,3

M0= -20+((345-225)/(345-225+345-305))*24,3 =1,775мкВ

Полученное значение попадает в модальный интервал.

 

 

 

 

eeg=load('eeg03_16.txt');

eegmin=-115;eegmax=128;eegdelta=24.3;

% Определение моды RR-интервалов

% Ввод вектора RR интервалов

Edg=eegmin:eegdelta:eegmax;

nhist=histc(eeg, Edg);

[AModa nMo]=max(nhist);

fmo=nhist(nMo);

fmo1=nhist(nMo-1);

fmo2=nhist(nMo+1);

x0=min(Edg)+(nMo-1)*eegdelta;

% Мода

Moda=((fmo-fmo1)/(2*fmo-fmo1-fmo2))*eegdelta+x0

% Амплитуда моды

AmpMod=AModa/length(eeg)

 

Moda =

    0.4250

 

AmpMod =

 

    0.2630

 

 

 

 

 

 

 

 

 

 

 

Определили  зависимость среднего арифметического значения от числа интервалов:

GG=load('eeg03_16.txt')

%Закрытие всех окон 

close all

%определение среднего  значения,дисперсии и СКО

n=1:length(GG);

for i=1:length(n);

varr(i)=var(GG(1:n(i)));

stdr(i)=std(GG(1:n(i)));

meanr(i)=mean(GG(1:n(i)));

end

%построение графика зависимости  дисперсии,от числа отсчета сигнала ЭЭГ

figure(1)

subplot(3,1,1),plot(varr),grid

ylabel('Дисперсия')

xlabel('Число отсчетов сигнала ЭЭГ','FontName','CourierMew')

%построение графика зависимости  среднего квадратичного отклонения  от числа

%отсчетов сигнала ЭЭГ

subplot(3,1,2),plot(stdr),grid

ylabel('СКО','FontName','CourierMew')

xlabel('Число отсчетов сигнала ЭЭГ ','FontName','CourierMew')

%построение графика зависимости  среднего значения сигнала ЭЭГ

subplot(3,1,3),plot(meanr),grid

ylabel('Среднее значение GG,мкВ','FontName','CourierMew')

xlabel('Число отсчетов сигнала ЭЭГ','FontName','CourierMew')

Рассмотрели спектральный анализ сигнала ЭЭГ. Выполнили спектральный анализ ЭЭГ сигнала отведения F7-A1.

%Выделим отведения F7-А1.

eegF7A1=eeg03_16(:,1);

%Сформируем вектор времени  записи сигнала ЭЭГ

fs=250;dt=1/fs;ns=length(eegF7A1);t=0:dt:dt*(ns-1);

%Выведем график сигнала  отведения О2-А2

figure(1);plot(t,eegF7A1);grid

title('EEG F7-A1 Signal');

xlabel('Time,s');ylabel('Magnitude')

%ВыполнимБПФсигнала.

Ye=fft(eegF7A1);Pe=Ye.*conj(Ye);

f=fs*(0:ns/2)/ns;np=round(ns/5);

%Выведем график спектра  мощности сигнала F7-A1

figure(2);plot(f(1:np),Pe(1:np));grid

title('EEG Signal Power Spektrum');xlabel('Frequency,Hz')

%Расчитаем ВКФ (Взаимную корреляционную функцию)соседних отведений F8-А2 &F3-А1.

eeg F8A2=load('eeg03_15.txt');

nccf=xcorr(eegF7A1,eegF8A2,'coeff');

%Выведемграфикнаэкран

figure(3);ncf=length(nccf);plot(t,nccf(ns:ncf));

title('Normalized Cross-Correlation Sequence')

xlabel('Time,s');grid

 

Получили графическое отображение: Сигнал ЭЭГ: отведение F7-A1:

Рисунок -10. Сигнал ЭЭГ: отведение F7-A1:

Спектральная мощность сигнала ЭЭГ, по которой видно что присутствуют частоты от 0 до 20, но доминирует частота=10 Гц. Следовательно, преобладает альфа ритм.

 

График корреляционной функции соседних отведений  F7-A1 и F8-А2

Коэффициент корреляции находится в пределах от -0.2 до 0.2, следовательно, влияние  отведений F7-A1 и F8-А2 друг на друга слабое.

               Рассмотрим методику удаления  помех из сигнала ЭЭГ

%Удалим высокочастотные  помехи из сигнала ЭЭГ

%Эти помехи могут быть  вызваны различными причинами.

% Для удаления помех  выполним цифровую фильтрацию  сигнала.

% Применим ФНЧ Баттерворта 8 порядка.

clear all%Очистка рабочей области

close all

fs=250;% частота дискретизации сигнала 250 ГЦ.

n=8;% порядок фильтра

fc=35;%частота среза ФНЧ

Wp=fc/(fs/2);%нормированный частотный диапазон.

% коэффициенты полинома  цифрового фильтра

[B,A]=butter(n,Wp);

nl=128;%длина вектора комплексных значений АЧХ

% Вектор комплексных значений  АЧХ цифрового фильтра

h=abs(freqz(B,A,nl));ff=fs/(2*nl)*(0:nl-1);

%АЧХ фильтра

figure(1);subplot(3,1,1);plot(ff,h)

title('Frequency response');xlabel('Frequency(Hz)')

ylabel('Magnitude');axis([0 2*fc 0 1.4]);grid

% Введем данные сигнала ЭЭГ

eeg=load('eeg03_16.txt');

eegn=length(eeg);

t=(1:eegn)/fs;

s=100*sin(2*pi*50*t);% синусоидальный сигнал с частотой 50 Гц

% Добавим к сигналу  ЭЭГ высокочастотную помеху 50 Гц

slen=eeg+s';

subplot(3,1,2);plot(t,slen);

xlabel('Time in seconds');ylabel('EEG');

axis tight; grid

% Отфильтруем сигнал ЭЭГ

eegf=filter(B,A,slen);

subplot(3,1,3);plot(t,eegf);axis tight;

title('The Noised EEG Signal');xlabel('Time in Seconds');

ylabel('EEG');grid

%Спектральный анализ сигнала  до и после фильтрации

nf=256;Y=fft(slen,nf);Yf=fft(eegf,nf);

%Определим спектральную  мощность сигнала

P=Y.*conj(Y)/nf;Pf=Yf.*conj(Yf)/nf;

 

 

%значения частотной оси

f=fs*(0:nf/2)/nf;

figure(2);subplot(2,1,1);plot(f(2:nl),P(2:nl))    grid;title('Original Signal Spectrum');       ylabel('FFT Magnitude');                     subplot(2,1,2); plot(f(2:nl), Pf(2:nl));grid           title ('After filtering Spectrum');

Получили графики: