Добавил:

artemtvi Опубликованный материал нарушает ваши авторские права? Сообщите нам.

Вуз:

Санкт-Петербургский государственный электротехнический университет "ЛЭТИ"

Предмет:

Методы обработки и анализа биомедицинских сигналов

Файл:

Лаба 4 / laba4

.docx

Скачиваний:

Добавлен:

08.04.2022

Размер:

410.36 Кб

Скачать

☆

МИНОБРНАУКИ РОССИИ

Санкт-Петербургский государственный

электротехнический университет

«ЛЭТИ» им. В.И. Ульянова (Ленина)

Кафедра БТС

отчет

по лабораторной работе №4

по дисциплине «Методы обработки и анализа биомедицинских сигналов и данных»

Тема: Спектральный анализ ЭКГ

Студентка гр. 7502		Екимова Н.В.
Студент гр.7502		Дегилевич А.А.
Студент гр.7502		Звягинцев Г.Р.
Преподаватель		Калиниченко А.Н.

Санкт-Петербург

2019

ЛАБОРАТОРНАЯ РАБОТА №4

«СПЕКТРАЛЬНЫЙ АНАЛИЗ ЭКГ»

Цель работы: исследование разложение сигнала по методу быстрого преобразования Фурье (БПФ), изучение спектров сигнала ЭКГ в норме и при различных патологиях.

Основные положения

При исследовании медико-биологических сигналов часто используются методы спектрального анализа, позволяющие получить численные оценки частотного состава сигнала. Наиболее распространен спектральный анализ, основанный на дискретном преобразовании Фурье:

где x (n)- отсчеты дискретного сигнала,N - число отсчетов, а k= 0, 1, ..., N-1 - номера частотных составляющих разложения. Для ускорения расчета применяют алгоритмы быстрого преобразования Фурье (БПФ).

В данной работе предлагается воспользоваться этим методом для исследования ЭКГ при трех различных видах сердечного ритма: нормальном ритме (часть a рисунка), желудочковой тахикардии (часть б) и фибрилляции желудочков (часть в).

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

Задание на выполнение работы

Рассчитать и исследовать спектральные оценки тестового сигнала. Построить графики сигнала, а также амплитудного спектра и спектральной плотности мощности (СПМ).
Рассчитать спектральные оценки для фрагментов ЭКГ, соответствующих трем различным видам сердечного ритма.

РАСЧЕТ ФУНКЦИЙ И ВЫВОД ГРАФИКОВ

Расчет спектральных оценок для тестового сигнала.

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

close;

clear;

X0=250;

Y0=100;

W0=1200;

H0=800;

figure('Position',[X0,Y0,W0,H0])

x1=60;

y1=60;

dx=50;

dy=50;

w=250;

h=200;

hAxes1=axes('Units','pixels','Position',[x1,y1+2*h+2*dy,2*w,h]);

hAxes2=axes('Units','pixels','Position',[x1+dx+2*w,y1+2*h+2*dy,w,h]);

hAxes3=axes('Units','pixels','Position',[x1+2*dx+3*w,y1+2*h+2*dy,w,h]);

hAxes4=axes('Units','pixels','Position',[x1,y1+h+dy,2*w,h]);

hAxes5=axes('Units','pixels','Position',[x1+dx+2*w,y1+h+dy,w,h]);

hAxes6=axes('Units','pixels','Position',[x1+2*dx+3*w,y1+h+dy,w,h]);

hAxes7=axes('Units','pixels','Position',[x1,y1,2*w,h]);

hAxes8=axes('Units','pixels','Position',[x1+dx+2*w,y1,w,h]);

hAxes9=axes('Units','pixels','Position',[x1+2*dx+3*w,y1,w,h]);

Fs=100;

tmax=1;

A1=1.5;

F1=5;

% F1=6.969;

A2=0.8;

F2=13;

C=0.2;

T=1/Fs;

t=0:T:tmax-T;

S=C+A1*sin(2*pi*F1*t)+A2*sin(2*pi*F2*t);

Asum=A1+A2;

ft=fft(S);

N=tmax*Fs;

for j=1:N

if (j==1)

as(j)=sqrt(real(ft(j))^2+imag(ft(j))^2)/N;

else

as(j)=sqrt(real(ft(j))^2+imag(ft(j))^2)/N*2;

end;

df=Fs/N;

for j=1:(N/2)

f(j)=df*(j-1);

end;

for j=1:N

if (j==1)

psd(j)=1/N*(real(ft(j))^2+imag(ft(j))^2)/Fs;

else

psd(j)=2/N*(real(ft(j))^2+imag(ft(j))^2)/Fs;

end;

m=mean(S);

S0=S-m;

ft0=fft(S0);

for j=1:N

if (j==1)

as0(j)=sqrt(real(ft0(j))^2+imag(ft0(j))^2)/N;

else

as0(j)=sqrt(real(ft0(j))^2+imag(ft0(j))^2)/N*2;

end;

for j=1:N

if (j==1)

psd0(j)=1/N*(real(ft0(j))^2+imag(ft0(j))^2)/Fs;

else

psd0(j)=2/N*(real(ft0(j))^2+imag(ft0(j))^2)/Fs;

end;

w=hamming(N);

for j=1:N

Sw(j)=S0(j)*w(j);

end;

ftw=fft(Sw);

for j=1:N

if (j==1)

asw(j)=sqrt(real(ftw(j))^2+imag(ftw(j))^2)/N;

else

asw(j)=sqrt(real(ftw(j))^2+imag(ftw(j))^2)/N*2;

end;

for j=1:N

if (j==1)

psdw(j)=1/N*(real(ftw(j))^2+imag(ftw(j))^2)/Fs;

else

psdw(j)=2/N*(real(ftw(j))^2+imag(ftw(j))^2)/Fs;

end;

axes(hAxes1)

plot(t,S)

grid on

set(hAxes1,'Ylim',[C-Asum C+Asum])

axes(hAxes2)

stem(f,as(1:N/2),'.')

axes(hAxes3)

plot(f,psd(1:(N/2)))

axes(hAxes4)

plot(t,S0)

grid on

set(hAxes4,'Ylim',[-Asum Asum])

axes(hAxes5)

stem(f,as0(1:N/2),'.')

axes(hAxes6)

plot(f,psd0(1:(N/2)))

axes(hAxes7)

plot(t,Sw)

hold on

plot(t,w)

grid on

set(hAxes1,'Ylim',[-Asum Asum])

axes(hAxes8)

stem(f,asw(1:N/2),'.')

axes(hAxes9)

plot(f,psdw(1:(N/2)))

Вид программы отображен на Рисунке 1.

Рис.1 Графическое окно тестовой программы

Рассчет спектральных оценок для спектров ЭКГ.

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

close;

clear;

X0=250;

Y0=100;

W0=1200;

H0=800;

figure('Position',[X0,Y0,W0,H0])

x1=60;

y1=60;

dx=50;

dy=50;

w=250;

h=200;

hAxes1=axes('Units','pixels','Position',[x1,y1+2*h+2*dy,2*w,h]);

hAxes2=axes('Units','pixels','Position',[x1+dx+2*w,y1+2*h+2*dy,w,h]);

hAxes3=axes('Units','pixels','Position',[x1+2*dx+3*w,y1+2*h+2*dy,w,h]);

hAxes4=axes('Units','pixels','Position',[x1,y1+h+dy,2*w,h]);

hAxes5=axes('Units','pixels','Position',[x1+dx+2*w,y1+h+dy,w,h]);

hAxes6=axes('Units','pixels','Position',[x1+2*dx+3*w,y1+h+dy,w,h]);

hAxes7=axes('Units','pixels','Position',[x1,y1,2*w,h]);

hAxes8=axes('Units','pixels','Position',[x1+dx+2*w,y1,w,h]);

hAxes9=axes('Units','pixels','Position',[x1+2*dx+3*w,y1,w,h]);

ECG=load('R4_01.txt');

Fs=250;

tmax=4;

T=1/Fs;

t=0:T:tmax-T;

Norma=ECG(:,1);

VTach=ECG(:,2);

VFibr=ECG(:,3);

ft=fft(Norma);

N=tmax*Fs;

for j=1:N

if (j==1)

as(j)=sqrt(real(ft(j))^2+imag(ft(j))^2)/N;

else

as(j)=sqrt(real(ft(j))^2+imag(ft(j))^2)/N*2;

end;

df=Fs/N;

for j=1:(N/2)

f(j)=df*(j-1);

end;

for j=1:N

if (j==1)

psd(j)=1/N*(real(ft(j))^2+imag(ft(j))^2)/Fs;

else

psd(j)=2/N*(real(ft(j))^2+imag(ft(j))^2)/Fs;

end;

m=mean(Norma);

Norma0=Norma-m;

ft0=fft(Norma0);

for j=1:N

if (j==1)

as0(j)=sqrt(real(ft0(j))^2+imag(ft0(j))^2)/N;

else

as0(j)=sqrt(real(ft0(j))^2+imag(ft0(j))^2)/N*2;

end;

for j=1:N

if (j==1)

psd0(j)=1/N*(real(ft0(j))^2+imag(ft0(j))^2)/Fs;

else

psd0(j)=2/N*(real(ft0(j))^2+imag(ft0(j))^2)/Fs;

end;

w=hamming(N);

for j=1:N

Normaw(j)=Norma0(j)*w(j);

end;

ftw=fft(Normaw);

for j=1:N

if (j==1)

asw(j)=sqrt(real(ftw(j))^2+imag(ftw(j))^2)/N;

else

asw(j)=sqrt(real(ftw(j))^2+imag(ftw(j))^2)/N*2;

end;

for j=1:N

if (j==1)

psdw(j)=1/N*(real(ftw(j))^2+imag(ftw(j))^2)/Fs;

else

psdw(j)=2/N*(real(ftw(j))^2+imag(ftw(j))^2)/Fs;

end;

ftT=fft(VTach);

N=tmax*Fs;

for j=1:N

if (j==1)

asT(j)=sqrt(real(ftT(j))^2+imag(ftT(j))^2)/N;

else

asT(j)=sqrt(real(ftT(j))^2+imag(ftT(j))^2)/N*2;

end;

for j=1:N

if (j==1)

psdT(j)=1/N*(real(ftT(j))^2+imag(ftT(j))^2)/Fs;

else

psdT(j)=2/N*(real(ftT(j))^2+imag(ftT(j))^2)/Fs;

end;

mT=mean(VTach);

VTach0=VTach-mT;

ftT0=fft(VTach0);

for j=1:N

VTachw(j)=VTach0(j)*w(j);

end;

ftTw=fft(VTachw);

for j=1:N

if (j==1)

asTw(j)=sqrt(real(ftTw(j))^2+imag(ftTw(j))^2)/N;

else

asTw(j)=sqrt(real(ftTw(j))^2+imag(ftTw(j))^2)/N*2;

end;

for j=1:N

if (j==1)

psdTw(j)=1/N*(real(ftTw(j))^2+imag(ftTw(j))^2)/Fs;

else

psdTw(j)=2/N*(real(ftTw(j))^2+imag(ftTw(j))^2)/Fs;

end;

ftF=fft(VFibr);

for j=1:N

if (j==1)

asF(j)=sqrt(real(ftF(j))^2+imag(ftF(j))^2)/N;

else

asF(j)=sqrt(real(ftF(j))^2+imag(ftF(j))^2)/N*2;

end;

for j=1:N

if (j==1)

psdF(j)=1/N*(real(ftF(j))^2+imag(ftF(j))^2)/Fs;

else

psdF(j)=2/N*(real(ftF(j))^2+imag(ftF(j))^2)/Fs;

end;

mF=mean(VFibr);

VFibr0=VFibr-mF;

ftF0=fft(VFibr0);

for j=1:N

VFibrw(j)=VFibr0(j)*w(j);

end;

ftFw=fft(VFibrw);

for j=1:N

if (j==1)

asFw(j)=sqrt(real(ftFw(j))^2+imag(ftFw(j))^2)/N;

else

asFw(j)=sqrt(real(ftFw(j))^2+imag(ftFw(j))^2)/N*2;

end;

for j=1:N

if (j==1)

psdFw(j)=1/N*(real(ftFw(j))^2+imag(ftFw(j))^2)/Fs;

else

psdFw(j)=2/N*(real(ftFw(j))^2+imag(ftFw(j))^2)/Fs;

end;

axes(hAxes1)

plot(t,Norma)

hold on

plot(t,Normaw)

grid on

axes(hAxes2)

stem(f,as(1:N/2),'.')

hold on

stem(f,asw(1:N/2),'.')

grid on

set(hAxes2,'Xlim',[0 20])

axes(hAxes3)

plot(f,psd(1:(N/2)))

hold on

plot(f,psdw(1:(N/2)))

grid on

set(hAxes3,'Xlim',[0 20])

axes(hAxes4)

plot(t,VTach)

hold on

plot(t,VTachw)

grid on

axes(hAxes5)

stem(f,asT(1:N/2),'.')

hold on

stem(f,asTw(1:N/2),'.')

grid on

set(hAxes5,'Xlim',[0 20])

axes(hAxes6)

plot(f,psdT(1:(N/2)))

hold on

plot(f,psdTw(1:(N/2)))

grid on

set(hAxes6,'Xlim',[0 20])

axes(hAxes7)

plot(t,VFibr)

hold on

plot(t,VFibrw)

grid on

axes(hAxes8)

stem(f,asF(1:N/2),'.')

hold on

stem(f,asFw(1:N/2),'.')

grid on

set(hAxes8,'Xlim',[0 20])

axes(hAxes9)

plot(f,psdF(1:(N/2)))

hold on

plot(f,psdFw(1:(N/2)))

grid on

set(hAxes9,'Xlim',[0 20])

Вид программы отображен на Рисунке 2.

Рис.2 Графическое отображение программы для расчета спектральных оценок ЭКГ

Вывод

В ходе лабораторной работы были изучены разложения сигнала по методу быстрого преобразования Фурье (БПФ), спектры сигналов ЭКГ в норме и при различных патологиях; были исследованы спектральные оценки тестового сигнала; построены графики сигнала, а также амплитудного спектра и спектральной плотности мощности (СПМ); были рассчитаны спектральные оценки для фрагментов ЭКГ, соответствующим трем различным видам сердечного ритма.

Соседние файлы в папке Лаба 4

#
25.12.20192.59 Кб57Lab4_1.m
#
25.12.20193.34 Кб63Lab4_2.m
#
08.04.2022410.36 Кб39laba4.docx
#
25.12.201919 Кб50R4_07.txt
#
08.04.2022360.13 Кб42Отчет 4.docx