Непрерывное Вейвлет-преобразование

Непрерывное Вейвлет-преобразование



ОГЛАВЛЕНИЕ



ВВЕДЕНИЕ

ПОСТАНОВКА ЗАДАЧИ

	3     ПОДХОДЫ К АНАЛИЗУ НЕСТАЦИОНАРНЫХ СИГНАЛОВ 

3.1  Методы  обработки нестационарных сигналов

3.2  Краткий обзор преобразования Фурье

Основные положения вейвлет-анализа

Методы вычисления непрерывного вейвлет преобразования 

  Во временной области 

  В частотной области

3.3.2   Выбор материнского вейвлета

	4    ОПРЕДЕЛЕНИЕ УЗЛОВЫХ ТОЧЕК ЭКГ НА ОСНОВЕ НЕПРЕРЫВНОГО ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ

4.1  Стандарты описания и обозначения ЭКГ.

4.2  Постановка задачи идентификации.

4.3  Построение модели идеальной ЭКГ.

4.4 Анализ модели ЭКГ

В системе Matlab

С использованием “Vision”

Сравнительный анализ полученных результатов.

ПРОГРАММНАЯ РЕАЛИЗАЦИЯ

 Структура программы

Форматы данных

5.3 Тестирование



ВВЕДЕНИЕ



На сегодняшний день одним из самых распространенных методов диагностики и распознавания сердечно-сосудистых заболеваний является электрокардиография. Сигнал ЭКГ характеризуется набором зубцов, по временным и амплитудным параметрам которых ставится диагноз. До недавнего времени процедуру нахождения характеристик зубцов выполнял врач-кардиолог, использую при этом только чертежные принадлежности. Такая схема достаточно проста и надежна, но требует много времени, и она работала в течении долгого времени из-за отсутствия альтернативных подходов к решению данной задачи.

С развитием компьютеров стали появляться специализированные комплексы, позволяющие выявлять сердечные заболевания, на основе автоматизированного анализа временных параметров ЭКГ. На сегодняшний день известны разработки фирм MedIT, Innomed Medical Co. Ltd. и другие. Кардиографы этих компаний выполняют основные операции, необходимые для работы в реальных условиях. Цены на них колеблются в диапазоне от 1’500$ до 15’000$. 

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

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

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



2. ПОСТАНОВКА ЗАДАЧИ



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

Для этого следует:

изучить форму, стандарты описания и обозначения ЭКГ; 

построить модель идеальной ЭКГ

провести сравнительный анализ эффективности системы Matlab и разработанного модуля

проанализировать временные затраты с учетом параметров съема ЭКГ и аппаратного обеспечения

При проведении исследований использовать систему Matlab 5.0, для разработки модуля использовать среду программирования Delphi 5.0.



3. ПОДХОДЫ К АНАЛИЗУ НЕСТАЦИОНАРНЫХ СИГНАЛОВ



3.1. Методы  обработки нестационарных сигналов



	Большинство медицинских сигналов имеет сложные частотно-временные характеристики. Как правило, такие сигналы состоят из близких по времени, короткоживущих высокочастотных компонент и долговременных, близких по частоте низкочастотных компонент.

	Для анализа таких сигналов нужен метод, способный обеспечить хорошее разрешение и по частоте, и по времени. Первое требуется для локализации низкочастотных составляющих, второе – для разрешения компонент высокой частоты.

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

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

	Существует два подхода к анализу нестационарных сигналов такого типа. Первый – локальное преобразование Фурье (short-time Fourier transform). Следуя по этому пути, мы работаем с нестационарным сигналом, как со стационарным, предварительно разбив его на сегменты (окна), статистика которых не меняется со временем. Второй подход – вейвлет преобразование. В этом случае нестационарный сигнал анализируется путем разложения по базисным функциям, полученным из некоторого прототипа путем сжатий, растяжений и сдвигов. Функция прототип называется материнским, или анализирующим вейвлетом.



3.2  Краткий обзор преобразования Фурье



	Классическим методом частотного анализа сигналов является преобразование Фурье, суть которого можно выразить формулой (1)





	Результат преобразования Фурье – амплитудно-частотный спектр, по которому можно определить присутствие некоторой частоты в исследуемом сигнале.

	В случае, когда не встает вопрос о локализации временного положения частот, метод Фурье дает хорошие результаты. Но при необходимости определить временной интервал присутствия частоты приходится применять другие методы.

	Одним из таких методов является обобщенный метод Фурье (локальное преобразование Фурье). Этот метод состоит из следующих этапов:

1.  в исследуемой функции создается “окно” – временной интервал, для          которого  функция f(x)0, и f(x)=0 для остальных значений;

2.    для этого “окна” вычисляется преобразование Фурье

3. “  окно” сдвигается, и для него также вычисляется преобразование Фурье

	“Пройдя” таким “окном” вдоль всего сигнала, получается некоторая трехмерная функция, зависящая от положения “окна” и частоты.

	Данный подход позволяет определить факт присутствия в сигнале любой частоты, и интервал ее присутствия. Это значительно расширяет возможности метода по сравнению с классическим преобразованием Фурье, но существуют и определенные недостатки. Согласно следствиям принципа неопределенности Гейзенберга в данном случае нельзя утверждать факт наличия частоты 0 в сигнале в момент времени t0  - можно лишь определить, что спектр частот (1,2) присутствует в интервале (t1,t2). Причем разрешение по частоте (по времени) остается постоянным вне зависимости от области частот (времен), в которых производится исследование. Поэтому, если, например, в сигнале существенна только высокочастотная составляющая, то увеличить разрешение можно только изменив параметры метода. В качестве метода, не обладающего подобного рода недостатками, был предложен аппарат вейвлет анализа. [2]



3.3  Основные положения вейвлет-анализа



	Различают дискретный и непрерывный вейвлет анализ, аппарат которых можно   применять    как   для  непрерывных,  так   и   для  дискретных  сигналов.

Cигнал анализируется путем разложения  по  базисным  функциям, полученным из некоторого прототипа путем сжатий,  растяжений и сдвигов (2).         Функция-прототип     называется    анализирующим      (материнским)       вейвлетом. 



     Вейвлет - функция должна удовлетворять 2-м условиям:



1.   Среднее значение (интеграл по всей прямой) равен 0. 

2.   Функция быстро убывает при t  ∞. 



Обычно,  функция-вейвлет обозначается буквой ψ.

В общем случае вейвлет преобразование функции f(t) выглядит так:



          				(2)

где t – ось времени, x – момент времени, s – параметр, обратный частоте, a (*) – означает комплексно-сопряженное.

	Главным элементом в вейвлет анализе является функция-вейвлет. Вообще говоря, вейвлетом является любая функция, отвечающая двум вышеуказанным условиям. Наибольшей популярностью пользуются два изображенных на рисунке 1  вейвлета:

Рис 1. Примеры вейвлетов.



	Сверху изображен  вейвлет “сомбреро” (Mexican Hat), названный так благодаря своему внешнему виду. На нижней части рисунка 1 изображен вейвлет Морле. График любого вейвлета выглядит примерно также, как и вейвлет Морле. Заметим, что вейвлет Морле – комплекснозначный, на рисунке изображены его вещественная и мнимая составляющие.

	Итак, у нас имеется некоторая функция f(t), зависящая от времени. Результатом ее вейвлет-анализа будет некоторая функция W(x,s), которая зависит уже от двух переменных: от времени и от частоты (обратно пропорционально). Для каждой пары x и s рецепт вычисления вейвлет преобразования следующий:

Функция вейвлет растягивается в s раз по горизонтали и в 1/s раз по вертикали.

Далее он сдвигается в точку x. Полученный вейвлет обозначается ψ(x,s).

Производится усреднение в окрестности точки s при помощи ψ(x,s).



В результате “вырисовывается” вполне наглядная картина, иллюстрирующая частотно-временные характеристики сигнала. По оси абсцисс откладывается время, по оси ординат – частота (иногда размерность оси ординат выбирается так: log(1/s), где s-частота), а абсолютное значение вейвлет преобразования для конкретной пары x и s определяет цвет, которым данный результат будет отображен (чем в большей степени та или иная частота присутствует в сигнале в конкретный момент времени, тем темнее будет оттенок).





Рис 2. Вейвлет преобразование стационарного сигнала.



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





Рис 3. Сравнение методов анализа.              

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

В то же время вейвлет преобразование (3b) дает исчерпывающую картину динамики изменения частотных характеристик во времени. Все это указывает на то, что вейвлет преобразование существенно более информативно по сравнению с преобразованием Фурье.



3.3.1 Методы вычисления непрерывного вейвлет-преобразования.



Существует два разных пути проведения вейвлет преобразования. Речь идет о расчетах во временной и частотной областях. При работе во временной области мы имеем дело с функциями, аргументами которых являются временные параметры, а в случае частотной – частотные. В частотной области используется механизм быстрого преобразования Фурье. [5c]



3.3.1.1  Во временной области



	Прежде всего, нам необходимо определить материнский вейвлет. Допустим, мы выбрали некоторую функцию, удовлетворяющую необходимым условиям: ψ0(η), где η – безразмерный период.

Итак, нам дана временная серия X, со значениями xn, в моменты времени n[0,N-1], где N – количество измерений. Каждая величина разделена по времени на постоянную величину dt. Получив основную формулу для материнского вейвлета, необходимо иметь возможность изменять размеры вейвлета. Для этого строится так называемый "масштабированный" вейвлет который будет иметь вид: 



                 			(3)



s – параметр, обратный частоте.



Вычисление вейвлет преобразования является сверткой искомой временной серии с функцией-вейвлетом.  Основная формула имеет вид :



					(4)



в данном случае (*) – означает комплексно-сопряженное.

	Результатом расчета Wn(s) по формуле (4)  будет комплексное число. В качестве конечного результата берется абсолютное значение полученного комплексного числа. [5a]



Блок – схема алгоритма:



Выбор начальной частоты (s=S0)











n =0







Для каждого n'[0,N-1] вычисляются значения функции вейвлета ψ* , где n' – точка в дискретной временной серии













Перемножение полученных результатов на соответствующие значения Xn'











Суммирование полученного числового ряда для всех n'[0,N-1]



Переход к следующей частоте (s=s+1)









n=n+1



								



n=N-1 ? ,?

							нет



						

						да



Получение результата вейвлет преобразования для пары значений n,s







									



Весь спектр частот изучен ?								    нет	





						    да	



Конец





3.3.1.2  В частотной области



	Вейвлет преобразование можно провести в частотной области. Для этого снова в первую очередь необходимо определить материнский вейвлет. Расчет по данной схеме происходит следующим образом: преобразование Фурье самого вейвлета (в данном случае будем рассматривать вейвлет Морле) сконцентрировано вокруг некоторой выделенной частоты 0 ≠ 0. Поэтому преобразование Фурье вейвлета, растянутого в  s  раз, будет сконцентрировано вокруг частоты 0/s  (см рис. 4).

Рис 4. Преобразование Фурье функции вейвлета.



	Так как свертка функций эквивалентна их перемножению в частотной области, “строка” s = const  на изображении вейвлет преобразования показывает эволюцию изучаемой функции на частотах, близких 0/s. То есть умножение Фурье-спектра исходной функции на пик в точке 0/s в частотной области (то есть на Фурье-образ растянутого вейвлета) вырезает из этой функции все то, что дает вклад в ее спектр на частотах, близких 0/s. В результате получается развертка спектрального компонента во времени. [1]

	Основные формулы имеют вид:

					(5)

где (*) – означает комплексно-сопряженное, а знак  (^) – преобразование Фурье.



						(6)												

 							(7)



Блок – схема алгоритма:



Выбор начальной частоты (s=S0)









n =0







Для каждого k[0,N-1] рассчитывается Фурье-спектр исходной временной серии, а также преобразование Фурье функции-вейвлета













Для каждого k[0,N-1] полученные результаты перемножаются.













Суммирование полученного числового ряда для всех k[0,N-1]



Переход к следующей частоте (s=s+1)









n=n+1



								



n=N-1 ? ,?

							нет



						

						да



Получение результата вейвлет преобразования для пары значений n,s







									



Весь спектр частот изучен ?								    нет	





						    да	



Конец





3.3.2  Выбор материнского вейвлета



В качестве материнского вейвлета подходит любая функция, удовлетворяющая двум вышеуказанным условиям. Для реализации алгоритма в качестве анализирующего вейвлета было решено воспользоваться вейвлетом Морле (рис. 5). Это было сделано по трем причинам:

вейвлет Морле один из наиболее популярных [1] и широко применяется

он обладает значительной наглядностью 

он прост в вычислительном плане, что ускоряет работу алгоритма





рис 5. Вейвлет Морле.



Фактически вейвлет Морле является произведением комплексной синусоиды на гауссиан. 

, 				(8)



где  является значением вейвлет функции с безразмерным периодом , а 0 - волновой параметр (при реализации 0=6).

	Необходимо также отметить, что вейвлет Морле является комплекснозначным, то есть имеет действительную и мнимую части.



	4    ОПРЕДЕЛЕНИЕ УЗЛОВЫХ ТОЧЕК ЭКГ НА ОСНОВЕ НЕПРЕРЫВНОГО ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ



4.1  Стандарты описания и обозначения ЭКГ.



Электрокардиограмма (ЭКГ) человека – сигнал, считываемый  в результате распространения волны деполяризации и реполяризации по сердечной мышце. Электрокардиограмма (ЭКГ) представляет собой некоторый сигнал, имеющий пять характерных пиков -  P, Q, R, S и T :

	Обозначенные особенности (пики и интервалы) и являются стандартами описания электрокардиограммы человека.



4.2  Постановка задачи идентификации.



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

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

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

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



4.3.Построение модели идеальной ЭКГ



В медицинских источниках есть сведения о параметрах ЭКГ здорового человека. Обычно эти данные и являются отправной точкой при анализе очередной электрокардиограммы. Для выработки подходов к автоматической идентификации нарушений в работе сердца необходимо построить модель сигнала В результате в ходе выполнения работы были построены две модели идеальной ЭКГ. Первая – для системы Matlab. Вторая – в рамках спецификации компьютерного кардиологического комплекса для анализа ЭКГ человека (для модуля  “Vision”). В обоих случаях модель представляет собой одномерный массив чисел, с частотой дискретизации 225 Гц. Длительность выбрана из расчета 2-3 сердечных сокращений.



4.4.Анализ модели ЭКГ.



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



В системе Matlab.



В системе Matlab с использованием стандартных средств данной системы построена модель двух сокращений сердечной мышцы. Графически она имеет следующий вид:





рис 7. Модель ЭКГ в системе Matlab.

 

	Результаты применения аппарата непрерывного вейвлет-преобразования к данной модели выглядят следующим образом:







рис 7. Вейвлет-преобразование в Matlab.

С использованием “Vision”.



Модель электрокардиограмы в данном случае имеет те же параметры, однако ее длительность увеличена до 3 секунд. Результаты ее обработки в модуле “Vision” выглядят следующим образом:





рис 8. Вейвлет-преобразование в “Vision”.

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



Сравнительный анализ полученных результатов.



Результат, полученный в среде Matlab хорошо локализует особенности ЭКГ, однако возникают проблемы с “чтением” масштаба как по временной оси, так и по частотной. Видно, что продвижение по времени осуществляется не по секундам (или другим временным единицам), а по индексу в массиве, содержащим сигнал. О частотных характеристиках сигнала можно судить лишь приблизительно, ввиду их значительной “размазанности” (см приложение 1.). Все это позволяет сделать вывод о том, что аппарат непрерывного вейвлет-преобразования в среде Matlab не слишком эффективен для решения задачи идентификации.

“Vision” дает наглядную частотно-временную развертку, позволяющую быстро и без дополнительных расчетов определить степень присутствия той или иной частоты в конкретный момент времени. Этому в определенной степени способствует  удобный масштаб (Гц,Сек). Помощь при решении задачи идентификации оказывает механизм выделения максимально присутствующей частоты (см. рис 8). Неудобство составляют небольшие искажения в области малых времен, связанные с особенностями непрерывного вейвлет-преобразования.



5.ПРОГРАММНАЯ РЕАЛИЗАЦИЯ



5.1.Структура программы



Программа содержит несколько модулей выполняющих различные функции: описание формата данных, прорисовка “миллиметровки” на которую непосредственно наносится  ЭКГ, вычисление непрерывного вейвлет-преобразования, визуализация частотно-временной развертки.

“Vision” читает файлы формата “*.1sf”, содержащие  записанную на кардиографе электрокардиограмму. При непосредственном считывании сигнал  автоматически отображается в окне на “миллиметровой бумаге”.  С помощью полосы скроллинга можно просмотреть все двенадцать отведений. На странице “Анализ” представлены необходимые инструменты для построения частотно-временного спектра. Имеется возможность выбора отведения для анализа (1-12), по умолчанию предлагается второе, как наиболее часто используемое. При проведении анализа измеряется время продолжительности процесса.



5.2.Форматы данных

Формат файла “*.1sf” в начале содержит заголовок, который описывает параметры, а также форму записи электрокардиограммы.

	Здесь используется структура имеющая вид:

type

  OtvData     = array [0..MaxArrayOtv] of SmallInt;

  OtvPointer  = ^ OtvData;

  OtvNumber   = array [0..MaxNumberOtv] of OtvPointer;

  OtvYLineP   = array [0..MaxNumberOtv] of Integer;

  PtrOtvYLineP  = ^ OtvYLineP;

  DataPointer   = ^ OtvNumber;

  LocationType  = (Vertical, Horizontal);

  WorkStyleType = (Standart12, Franc, Aculinichev, Reserved);

  ECGType = record

   Frequency    : Word; 

   Time         : Single; 

   OtvQuantity  : Byte;         

   Location     : LocationType;  

   WorkStyle    : WorkStyleType; 

   MultipleCoef : Single;       

  end; 



Frequency – частота съема ЭКГ,  Time – продолжительность съема, OtvQuantity – количество отведений, Location – формат записи (Vertical – предусматривает параллельную запись всех отведений, применяется при считывании  данных с кардиографа, Horizontal – последовательная запись отведений: сначала первое, затем второе и т.д.), MultipleCoef – коэффициент числового ряда сигнала. 



5.3.Оценка временных затрат



Основными критериями, определяющими скорость работы программы являются: частота съема ЭКГ и параметры аппаратной части.

4