В этом примере показано, как использовать вейвлеты, чтобы анализировать физиологические сигналы.
Физиологические сигналы являются часто неустановившимся подразумевать, что их содержимое частоты изменяется в зависимости от времени. Во многих приложениях эти изменения являются мероприятиями.
Вейвлеты разлагают сигналы на изменяющуюся во времени частоту (шкала) компоненты. Поскольку функции сигнала часто локализуются вовремя и частота, анализ и оценка легче при работе с более разреженными (уменьшаемыми) представлениями.
Этот пример представляет несколько иллюстративных случаев, где способность вейвлетов обеспечить локальный частотно-временной анализ сигналов выгодна.
Комплекс 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. Шкалы соответствуют следующим аппроксимированным диапазонам частот.
Масштабируйтесь 4 - [11.25, 22.5), Гц
Масштабируйтесь 5 - [5.625, 11.25), Гц.
Это покрывает полосу пропускания, которая, как показывают, максимизировала энергию 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 (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).
Следующие функции помощника используются в этом примере.