Содержание отчета

1. Название, цель и задачи работы.

2. Тексты программ и изображения графических окон, которые требовалось сохранять в ходе выполнения работы.

3. Объяснение полученных результатов и выводы.

Лабораторная работа 4. Спектральный анализ экг

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

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

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

,

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

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

а

б

в

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

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

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

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

Порядок выполнения работы

Расчет спектральных оценок для тестового сигнала. Значения частоты дискретизации (F_d), длительности (t_max), амплитуд (A₁, A₂) и частот (F₁, F₂) синусоид, а также постоянной составляющей (c) для тестового сигнала указаны в табл. П.4 приложения.

1. Запустите систему MATLAB.

2. В главном окне установите путь к рабочей папке.

3. Создайте новый М-файл и сохраните его под каким-нибудь именем.

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

5. Рассчитайте и выведите в верхнем поле первой колонки тестовый сигнал (сумма двух синусоид и константы, см. табл. П.4). Нанесите на график координатную сетку, пользуясь командой

grid on % Нанесение координатной сетки на график

6. Получите разложение тестового сигнала, пользуясь функцией fft, позволяющей вычислить БПФ:

ft=fft(s); % Расчет БПФ для сигнала “s”

7. Рассчитайте амплитудный спектр:

N=tmax*Fd; % Число отсчетов в тестовом сигнале

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

end

8. Постройте в верхнем поле второй колонки график амплитудного спектра для частот от 0 Гц до 0,5F_d :

df=Fd/N; % Шаг по оси частот

for j=1:(N/2)

f(j)=df*(j-1); % Расчет аргумента для построения графика

end

stem(f,as(1:(N/2),'.') % График амплитудного спектра

9. Убедитесь, что линии спектра совпадают по значениям с заданными для соответствующих частот амплитудами синусоид (постоянной составляющей соответствует частота 0 Гц).

10. Рассчитайте СПМ и постройте ее график в правом верхнем поле:

for j=1:N % Цикл для перебора всех элементов разложения

if (j==1) % Случай нулевой частоты:

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

else % Для всех частот, кроме нулевой:

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

end;

end;

axes(hAxes3) % Активизация 3-го поля для вывода графика

plot (f,psd(1:(N/2)))% График СПМ

Рассчитанные с использованием БПФ спектральные оценки имеют два существенных недостатка: 1) в них может присутствовать сильная составляющая на нулевой частоте; 2) из-за необходимости предположения о периодичности сигнала могут возникнуть ложные высокочастотные составляющие, вызванные скачками на краях анализируемого фрагмента. Для устранения этих нежелательных эффектов из сигнала перед БПФ вычитают среднее значение и далее умножают сигнал на так называемую “оконную” функцию, подавляющую скачки на краях.

11. Вычтите из сигнала его среднее значение:

m=mean(s); % Вычисление среднего для сигнала “s”

s0=s-m; % Удаление из сигнала среднего значения

12. В полях второго ряда постройте графики полученного сигнала (с нанесением координатной сетки), а также амплитудного спектра и СПМ, рассчитанных для этого сигнала.

В случае правильного выполнения последнего пункта, спектральная составляющая на нулевой частоте должна быть равна нулю.

13. Умножьте сигнал с удаленным средним на оконную функцию, называемую “окном Хэмминга”:

w=hamming(N); % Создание массива значений оконной функции “w”

for j=1:N % Умножение сигнала на оконную функцию:

sw(j)=s0(j)*w(j);

end

14. Постройте в левом нижнем окне одновременно два графика: сигнала, умноженного на оконную функцию, и самой этой функции. Нанесите на график координатную сетку.

15. В остальных полях третьего ряда постройте графики рассчитанных для полученного сигнала амплитудного спектра и СПМ.

16. Сохраните графическое окно и текст программы.

17. Установите дробное значение частоты для одной из синусоид.

На двух верхних спектрах должен наблюдаться эффект “утечки” мощности. На нижнем спектре этот эффект должен быть менее заметным.

18. Сохраните графическое окно.

Расчет спектральных оценок для фрагментов ЭКГ. Имя файла ЭКГ указано в табл. П.4. Частота дискретизации F_d = 250 Гц.

19. Перепишите в рабочую папку файл с реальной записью ЭКГ.

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

20. Создайте новый М-файл и сохраните его под каким-нибудь именем.

21. Наберите в тексте программы команду чтения заданного файла ЭКГ и выполните получившуюся программу.

В результате последнего действия в рабочем пространстве MATLAB (Workspace) образуется переменная, содержащая прочитанные значения сигнала. С помощью специального компонента MATLAB – программы SPTool можно просмотреть этот сигнал и выбрать участки для анализа.

22. Через кнопку Start, расположенную в нижнем левом углу главного окна MATLAB, вызовите программу SPTool (Start->Toolboxes->Signal Processing->SPTool). В окне SPTool выберите пункт меню File->Import.

23. В поле Workspace Contents окна Import to SPTool отметьте переменную, в которую программа прочла содержимое файла, и кнопкой “>” перенесите ее в поле Data. Нажмите OK для закрытия окна Import.

В результате последнего действия в левом поле окна SPTool появится имя массива для просмотра сигнала.

24. Нажмите кнопку View для вызова окна просмотра Signal Browser.

25. С помощью кнопки растяните изображение сигнала так, чтобы на экране помещалось несколько циклов сердечных сокращений.

26. Просмотрите запись сигнала, пользуясь рамкой прокрутки, расположенной в нижней части окна. Выберите три участка длительностью по 4 с, которые по виду подходят под определения трех описанных нарушений ритма. Для каждого участка с помощью маркерной линии (которая инициализируются кнопкой “1” над правым верхним углом поля просмотра сигнала) отметьте начало участка и запишите номер соответствующего индекса, отображаемый в поле Marker 1.

27. После выполнения последнего пункта воспользуйтесь значениями уже определенных индексов для перезаписи выбранных фрагментов сигнала в отдельные массивы. Например:

ECG=load('R4_00.txt'); % Чтение сигнала в переменную “ECG”

Fd=250; % Частота дискретизации

tmax=4; % Длительность фрагмента

T=1/Fd; % Интервал дискретизации

t=0:T:tmax-T; % Аргумент для оси времени

N=tmax*Fd; % Число отсчетов в фрагменте ЭКГ

ind=13500; % Индекс, определенный с помощью SPTool

Norma=ECG(ind:ind+N-1); % Чтение фрагмента в отдельный массив

plot(t,Norma) % Контрольный вывод графика

28. В качестве заготовки используйте программу, подготовленную в первой части работы, сохранив ее под другим именем. Переделайте эту программу так, чтобы в левых полях она выводила по два графика для каждого из трех фрагментов ЭКГ: в исходном виде, а также после удаления среднего и умножения на окно. В средних и правых полях должны выводиться графики соответствующих амплитудных спектров и СПМ для частот от 0 до 20 Гц, т. е. первые 80 значений.

29. Сохраните графическое окно и текст программы.