Обнаружение волны R в электрокардиограмме с MODWT

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

Нагрузка и график формы сигнала ЭКГ, где пики R комплекса QRS были аннотированы двумя или более кардиологами. Данные ЭКГ и аннотации взяты из базы данных аритмии MIT-BIH. Данные дискретизируются на частоте 360 Гц.

load mit200
figure
plot(tm,ecgsig)
hold on
plot(tm(ann),ecgsig(ann),'ro')
xlabel('Seconds')
ylabel('Amplitude')
title('Subject - MIT-BIH 200')

Можно использовать вейвлеты для построения автоматического QRS-детектора для использования в таких приложениях, как оценка интервала R-R.

Существует две клавиши для использования вейвлетов в качестве детекторов общих признаков:

Вейвлет «sym4» напоминает комплекс QRS, что делает его хорошим выбором для обнаружения QRS. Чтобы проиллюстрировать это более четко, извлеките комплекс QRS и постройте график результата с расширенным и переведенным вейвлетом «sym4» для сравнения.

qrsEx = ecgsig(4560:4810);
[mpdict,~,~,longs] = wmpdictionary(numel(qrsEx),'lstcpt',{{'sym4',3}});
figure
plot(qrsEx)
hold on
plot(2*circshift(mpdict(:,11),[-2 0]),'r')
axis tight
legend('QRS Complex','Sym4 Wavelet')
title('Comparison of Sym4 Wavelet and QRS Complex')

Используйте максимальное перекрывающееся дискретное вейвлет-преобразование (MODWT) для увеличения пиков R в форме сигнала ЭКГ. MODWT является недекимированным вейвлет-преобразованием, которое обрабатывает произвольные размеры выборки.

Во-первых, разложить сигнал ЭКГ до уровня 5, используя вейвлет «sym4» по умолчанию. Затем реконструируют локализованную по частоте версию формы сигнала ЭКГ, используя только вейвлет-коэффициенты на шкалах 4 и 5. Шкалы соответствуют следующим приблизительным полосам частот.

Это охватывает полосу пропускания, показанную для максимизации энергии QRS.

wt = modwt(ecgsig,5);
wtrec = zeros(size(wt));
wtrec(4:5,:) = wt(4:5,:);
y = imodwt(wtrec,'sym4');

Использовать возведенные в квадрат абсолютные значения аппроксимации сигнала, построенной из вейвлет-коэффициентов, и использовать алгоритм поиска пиков для идентификации R пиков.

При наличии Toolbox™ обработки сигналов можно использовать findpeaks для определения местоположения пиков. Постройте график сигнала R-пика, полученного с помощью вейвлет-преобразования, аннотированного автоматически обнаруженными местоположениями пиков.

y = abs(y).^2;
[qrspeaks,locs] = findpeaks(y,tm,'MinPeakHeight',0.35,...
    'MinPeakDistance',0.150);
figure
plot(tm,y)
hold on
plot(locs,qrspeaks,'ro')
xlabel('Seconds')
title('R Peaks Localized by Wavelet Transform with Automatic Annotations')

Добавьте аннотации экспертов в форму сигнала R-пика. Время автоматического обнаружения пика считается точным, если оно находится в пределах 150 мс от истинного пика (мс).

plot(tm(ann),y(ann),'k*')
title('R peaks Localized by Wavelet Transform with Expert Annotations')

В командной строке можно сравнить значения tm(ann) и locs, которые являются временем эксперта и временем автоматического пикового обнаружения соответственно. Увеличение пиков R с помощью вейвлет-преобразования приводит к скорости попадания 100% и отсутствию ложных срабатываний. Вычисленная частота сердечных сокращений с использованием вейвлет-преобразования составляет 88,60 ударов в минуту по сравнению с 88,72 ударов в минуту для аннотированной формы сигнала.

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

figure
plot(tm,ecgsig,'k--')
hold on
plot(tm,y,'r','linewidth',1.5)
plot(tm,abs(ecgsig).^2,'b')
plot(tm(ann),ecgsig(ann),'ro','markerfacecolor',[1 0 0])
set(gca,'xlim',[10.2 12])
legend('Raw Data','Wavelet Reconstruction','Raw Data Squared', ...
    'Location','SouthEast')
xlabel('Seconds')

Используя findpeaks по квадратам необработанных данных получается двенадцать ложных срабатываний.

[qrspeaks,locs] = findpeaks(ecgsig.^2,tm,'MinPeakHeight',0.35,...
    'MinPeakDistance',0.150);

Помимо переключателей в полярности пиков R, ЭКГ часто искажается шумом.

load mit203
figure
plot(tm,ecgsig)
hold on
plot(tm(ann),ecgsig(ann),'ro')
xlabel('Seconds')
ylabel('Amplitude')
title('Subject - MIT-BIH 203 with Expert Annotations')

Используйте MODWT для выделения пиков R. Использовать findpeaks для определения местоположения пиков. Постройте график R-пиковой формы сигнала вместе с экспертными и автоматическими аннотациями.

wt = modwt(ecgsig,5);
wtrec = zeros(size(wt));
wtrec(4:5,:) = wt(4:5,:);
y = imodwt(wtrec,'sym4');
y = abs(y).^2;
[qrspeaks,locs] = findpeaks(y,tm,'MinPeakHeight',0.1,...
    'MinPeakDistance',0.150);
figure
plot(tm,y)
title('R-Waves Localized by Wavelet Transform')
hold on
hwav = plot(locs,qrspeaks,'ro');
hexp = plot(tm(ann),y(ann),'k*');
xlabel('Seconds')
legend([hwav hexp],'Automatic','Expert','Location','NorthEast')

Частота попаданий снова составляет 100% при нулевых ложных тревогах.

В предыдущих примерах использовался очень простой вейвлет QRS-детектор, основанный на приближении сигнала, построенном из modwt. Цель состояла в том, чтобы продемонстрировать способность вейвлет-преобразования изолировать компоненты сигнала, а не создавать наиболее надежный QRS-детектор на основе вейвлет-преобразования. Можно, например, использовать тот факт, что вейвлет-преобразование обеспечивает многомасштабный анализ сигнала для улучшения обнаружения пиков. Изучите шкалу 4 и 5 в квадрате длин волн, нанесенную на график вместе с временем пика R, как аннотировано экспертами. Детали уровня 4 сдвинуты для визуализации.

ecgmra = modwtmra(wt);
figure
plot(tm,ecgmra(5,:).^2,'b')
hold on
plot(tm,ecgmra(4,:).^2+0.6,'b')
set(gca,'xlim',[14.3 25.5])
timemarks = repelem(tm(ann),2);
N = numel(timemarks);
markerlines = reshape(repmat([0;1],1,N/2),N,1);
h = stem(timemarks,markerlines,'k--');
h.Marker = 'none';
set(gca,'ytick',[0.1 0.6]);
set(gca,'yticklabels',{'D5','D4'})
xlabel('Seconds')
title('Magnitude-Squared Level 4 and 5 Details')

Вы видите, что пики в деталях уровня 4 и уровня 5 имеют тенденцию совпадать. Более совершенный алгоритм поиска пиков вейвлет может использовать это, используя информацию из нескольких масштабов одновременно.

Изменяющийся во времени вейвлет-анализ когерентности динамики мозга

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

Чтобы проиллюстрировать это, изучите данные спектроскопии ближнего инфракрасного диапазона (NIRS), полученные у двух людей. NIRS измеряет активность мозга, используя различные характеристики поглощения кислородсодержащего и дезоксигенированного гемоглобина. Данные взяты из Cui, Bryant, & Reiss (2012) и были любезно предоставлены авторами для этого примера. Место записи было верхней лобной корой для обоих субъектов. Данные дискретизируются при частоте 10 Гц.

В эксперименте испытуемые альтернативно сотрудничали и соревновались в выполнении задачи. Период выполнения задания составил семь секунд.

load NIRSData
figure
plot(tm,[NIRSData(:,1) NIRSData(:,2)])
legend('Subject 1','Subject 2','Location','NorthWest')
xlabel('Seconds')
title('NIRS Data')
grid on

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

cwt(NIRSData(:,1),10,'bump')
figure
cwt(NIRSData(:,2),10,'bump')

Анализ CWT показывает сильные частотно-модулированные колебания в обоих наборах данных около 1 Гц. Это связано с сердечными циклами у двух субъектов. Кроме того, появляется более слабое колебание в обоих наборах данных около 0,15 Гц. Эта деятельность является более сильной и более последовательной в предмете 1, чем предмет 2. Вейвлет-когерентность может улучшить обнаружение слабых колебаний, которые совместно присутствуют в двух временных рядах.

[wcoh,~,F] = wcoherence(NIRSData(:,1),NIRSData(:,2),10);
figure
surf(tm,F,abs(wcoh).^2); view(0,90)
shading interp
axis tight
hc = colorbar;
hc.Label.String = 'Coherence';
title('Wavelet Coherence')
xlabel('Seconds')
ylabel('Hz')
ylim([0 2.5])
set(gca,'ytick',[0.15 1.2 2])

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

taskbd = [245 1702 2065 3474];
tvec = repelem(tm(taskbd),2);
yvec = [0 max(F)]';
yvec = reshape(repmat(yvec,1,4),8,1);
hold on
stemPlot = stem(tvec,yvec,'w--','linewidth',2);
stemPlot.Marker = 'none';

Этот пример используется cwt для получения и построения графика частотно-временного анализа отдельных временных рядов NIRS. Также используется пример wcoherence для получения вейвлет-когерентности двух временных рядов. Использование вейвлет-когерентности часто позволяет обнаружить когерентное колебательное поведение в двух временных рядах, которые могут быть довольно слабыми в каждой отдельной серии. Обратитесь к Cui, Bryant и Reiss (2012) за более подробным анализом вейвлет-когерентности этих данных.

Частотно-временной анализ данных об отоакустических выбросах с использованием CWT

Отоакустические излучения (OAE) представляют собой узкополосные колебательные сигналы, излучаемые улиткой (внутренним ухом), и их наличие свидетельствует о нормальном слухе. Загрузите и постройте график некоторых примеров данных OAE. Данные дискретизируют при частоте 20 кГц.

load dpoae
figure
plot(t.*1000,dpoaets)
xlabel('Milliseconds')
ylabel('Amplitude')

Эмиссия вызвана стимулом, начинающимся через 25 миллисекунд и заканчивающимся через 175 миллисекунд. Исходя из экспериментальных параметров, частота излучения должна составлять 1230 Гц. Получение и построение графика CWT как функции времени и частоты. Используйте аналитический вейвлет Морса по умолчанию с 16 голосами на октаву.

[dpoaeCWT,f] = cwt(dpoaets,2e4,'VoicesPerOctave',16);
helperCWTTimeFreqPlot(dpoaeCWT,t.*1000,f,...
    'surf','CWT of OAE','milliseconds','Hz')

Можно исследовать эволюцию времени OAE, найдя коэффициенты CWT, наиболее близкие по частоте к 1230 Гц, и рассматривая их величины как функцию времени. Постройте график величин вместе с маркерами времени, обозначающими начало и конец вызывающего стимула.

[~,idx1230] = min(abs(f-1230));
cfsOAE = dpoaeCWT(idx1230,:);
plot(t.*1000,abs(cfsOAE))
hold on
AX = gca;
plot([25 25],[AX.YLim(1) AX.YLim(2)],'r')
plot([175 175],[AX.YLim(1) AX.YLim(2)],'r')
xlabel('msec')
title('CWT Coefficient Magnitudes')

Существует некоторая задержка между началом вызывающего стимула и ОАЭ. Как только возбуждающий стимул прекращается, ОАЭ немедленно начинает распадаться по величине.

Другим способом изолирования излучения является использование обратного CWT для восстановления частотно-локализованного приближения во временной области.

Построение частотно-локализованной аппроксимации излучения путем извлечения коэффициентов CWT, соответствующих частотам между 1150 и 1350 Гц. Используйте эти коэффициенты и инвертируйте CWT. Постройте график исходных данных вместе с реконструкцией и маркерами, указывающими начало и конец вызова стимула.

frange = [1150 1350];
xrec = icwt(dpoaeCWT,f,frange);
figure
plot(t.*1000,dpoaets)
hold on
plot(t.*1000,xrec,'r')
AX = gca;
ylimits = AX.YLim;
plot([25 25],ylimits,'k')
plot([175 175],ylimits,'k')
grid on
xlabel('Milliseconds')
ylabel('Amplitude')
title('Frequency-Localized Reconstruction of Emission')

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

Важно отметить, что даже если для реконструкции был выбран диапазон частот, аналитическое вейвлет-преобразование фактически кодирует точную частоту излучения. Чтобы продемонстрировать это, возьмем преобразование Фурье аппроксимации излучения, реконструированное из аналитического CWT.

xdft = fft(xrec);
freq = 0:2e4/numel(xrec):1e4;
xdft = xdft(1:numel(xrec)/2+1);
figure
plot(freq,abs(xdft))
xlabel('Hz')
ylabel('Magnitude')
title('Fourier Transform of CWT-Based Signal Approximation')
[~,maxidx] = max(abs(xdft));
fprintf('The frequency is %4.2f Hz\n',freq(maxidx))

The frequency is 1230.00 Hz

Этот пример используется cwt для получения частотно-временного анализа данных ОАЭ и icwt получение частотно-локализованного приближения к сигналу.

Ссылки

Кюи, Х., Брайант, Д.М. и Рейсс. A.L. «Гиперсканирование на основе NIRS выявляет повышенную межличностную когерентность в верхней лобной коре во время сотрудничества», Neuroimage, 59 (3), 2430-2437, 2012.

Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG, Mietus JE, Moody GB, Peng C-K, Stanley HE. «PhysioBank, PhysioToolkit и PhysioNet: компоненты нового исследовательского ресурса для сложных физиологических сигналов». Циркуляционный 101 (23): e215-e220, 2000.http://circ.ahajournals.org/cgi/content/full/101/23/e215

Маллат, С. «Вейвлет-тур обработки сигналов: разреженный путь», Академическая пресса, 2009.

Муди, Г.Б. «Оценка анализаторов ЭКГ». http://www.physionet.org/physiotools/wfdb/doc/wag-src/eval0.tex

Муди, Марк РГ. «Влияние базы данных аритмии MIT-BIH». IEEE Eng in Med and Biol 20 (3): 45-50 (май-июнь 2001).

Приложение - Вспомогательные функции

В этом примере используются следующие вспомогательные функции.