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

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

Загрузка и построение графика формы волны ЭКГ, где peaks 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), чтобы увеличить peaks R в форме волны ECG. MODWT является неопределенным вейвлет, которое обрабатывает произвольные размеры выборки.

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

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

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

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

Если у вас есть Signal Processing Toolbox™, можно использовать findpeaks для определения местоположения peaks. Постройте график формы волны 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 мс от истинного пика (msec).

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

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

Если вы пытаетесь работать с квадратными величинами исходных данных, вы находите возможность вейвлета преобразования для изоляции R, peaks делает задачу обнаружения намного легче. Работа с необработанными данными может вызвать неправильные идентификации, такие как когда квадратный пик 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);

В дополнение к переключателям в полярности peaks 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, чтобы изолировать peaks 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')

Вы видите, что peaks в деталях уровня 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. Данные отбираются с частотой дискретизации 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')

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

Другой способ изоляции излучения - это использование обратной 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 для получения частотно-временного анализа данных OAE и icwt для получения частотно-локализованного приближения к сигналу.

Ссылки

Кюи, Х., Брайант, Д. М. и Рейсс. Гиперсканирование на основе NIRS обнаруживает повышенную межличностную когерентность в верхней лобной коре во время сотрудничества ", Нейроизображение, 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

Mallat, S. «A Wavelet Tour of Signal Processing: The Sparse Way», Academic Press, 2009.

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

Moody GB, Mark RG. «влияние базы данных аритмии MIT-BIH». IEEE Eng in Med and Biol 20 (3): 45-50 (май-июнь 2001).

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

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