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

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

Загрузите и постройте форму волны ECG, где peaks R комплекса QRS был аннотирован двумя или больше кардиологами. Данные о ECG и аннотации взяты из Базы данных Аритмии 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 является неподкошенным вейвлетом, преобразовывают, который обрабатывает произвольные объемы выборки.

Во-первых, анализируйте форму волны ECG вниз к уровню 5 с помощью значения по умолчанию 'sym4' вейвлет. Затем восстановите локализованную версию частоты формы волны ECG с помощью только коэффициенты вейвлета в шкалах 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 ударами/минута для аннотируемой формы волны.

При попытке работать над квадратными величинами исходных данных, вы находите, что возможность вейвлета преобразовывает, чтобы изолировать peaks 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);

В дополнение к переключателям в полярности peaks R ECG часто повреждается шумом.

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 стремится к co-occur. Более усовершенствованный алгоритм нахождения пика вейвлета мог использовать это при помощи информации от нескольких шкал одновременно.

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

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

Чтобы проиллюстрировать это, исследуйте спектроскопию в ближней инфракрасной области (NIRS) данные, полученные в двух испытуемых людях. NIRS измеряет мозговое действие путем использования различных характеристик поглощения окисленного и deoxygenated гемоглобина. Данные взяты из 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 получить когерентность вейвлета двух временных рядов. Использование когерентности вейвлета часто позволяет вам обнаружить когерентное колебательное поведение в двух временных рядах, которые могут быть довольно слабыми в каждом отдельном ряду. Consult Cui, Bryant, & Reiss (2012) для более подробного анализа когерентности вейвлета этих данных.

Частотно-временной анализ данных об эмиссии Otoacoustic с CWT

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

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

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

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

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

[~,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. Если вызывающий стимул отключен, OAE сразу начинает затухать в величине.

Другой способ изолировать эмиссию состоит в том, чтобы использовать обратный CWT, чтобы восстановить локализованное частотой приближение во временном интервале.

Создайте локализованное частотой приближение эмиссии путем извлечения коэффициентов CWT, соответствующих частотам между 1150 и 1 350 Гц. Используйте эти коэффициенты и инвертируйте 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 получить локализованное частотой приближение к сигналу.

Ссылки

Цуй, X., Брайант, D.M., и Reiss. A.L. "ОСНОВАННОЕ НА БИК-СПЕКТРОСКОПИИ гиперсканирование показывает увеличенную межабонентскую когерентность в превосходящей лобной коре во время сотрудничества", Нейроизображение, 59 (3), 2430-2437, 2012.

Голдбергер ЭЛ, LAN Amaral, Стекло L, Гаусдорф ДЖМ, Иванов PCh, Марк РГ, Mietus JE, Капризный Гбайт, Пенг C-K, Стэнли ХЭ. "PhysioBank, PhysioToolkit и PhysioNet: Компоненты Нового Ресурса Исследования для Комплексных Физиологических Сигналов". Циркуляция 101 (23): e215-e220, 2000. http://circ.ahajournals.org/cgi/content/full/101/23/e215

Mallat, S. "Тур вейвлета по обработке сигналов: разреженный путь", Academic Press, 2009.

Капризный, G.B. "Оценивая ECG Анализаторы". http://www.physionet.org/physiotools/wfdb/doc/wag-src/eval0.tex

Капризный Гбайт, Марк РГ. "Удар Базы данных Аритмии MIT-BIH". Инженер IEEE в Медиане и Biol 20 (3):45-50 (мочь-июнь 2001).

Приложение - поддерживание функций

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