CC BY
131
УДК 681.3.06: 621.372
А.В. КОЛОМЕЙЦЕВА, Г.В. МИШУГОВА, А.П. МУЛ, Г.Ю. РЯБЫХ ПРИМЕНЕНИЕ ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ
И МЕТОДА ПРОНИ ДЛЯ ИДЕНТИФИКАЦИИ БИОГЕННЫХ СИГНАЛОВ
Рассмотрены нестационарные сигналы биогенного характера (кардиограмма, электроэнцефалограмма и др.). Применен анализ сигналов с использованием вейвлет-преобразования и метода Прони, которые в отличие от классического метода Фурье дают улучшенную спектральную оценку биогенных сигналов. Ключевые слова: дискретное косинус- и синус-преобразование, оконное преобразование Фурье, вейвлет-преобразование Добеши, метод Прони, тёплицевы структуры.
Введение. Технологии цифровой обработки сигналов и их изображений находят все более широкое применение в различных областях, в частности при определении характеристик и идентификации биогенных сигналов.
Разработка алгоритма для анализа биогенного сигнала, однако, является непростой задачей; довольно часто это даже не целенаправленный процесс. Инженер или компьютерный аналитик часто бывает поражен изменчивостью и разнообразием признаков в биогенных сигналах и системах, где эти факторы проявляются в большей степени, чем в физических системах или наблюдениях. Учет всех возможностей и степеней свободы в биогенных системах является наиболее сложной проблемой для большинства применений. Методы, показавшие свою пригодность для работы с определенными системами или наборами сигналов, могут оказаться несостоятельными в других, на первый взгляд похожих, ситуациях.
Так, для получения правдивого результата биогенного сигнала следует иметь расслабленное и спокойное состояние при отсутствии движений. Кашель, напряжение мышц, движения глаз, конечностей вызывают соответствующие сигналы, играющие роль нежелательных артефактов. В результате рекомендацией при снятии сигнала служит задержка дыхания на несколько секунд. Но это предложение неприемлемо для пациентов в критическом состоянии или при проведении диагностики у младенцев.
В таких и других случаях для удаления артефактов используются методы цифровой обработки сигналов.
Обработка сигналов, связанная с анализом их спектров, называется спектральным анализом. Спектральный анализ используется при распознавании, обнаружении и сжатии сигналов. В основе спектрального анализа сигналов лежит интегральное преобразование и ряды Фурье.
Идентификация особенностей биогенных сигналов - весьма важный этап, так как допущенные ошибки сказываются на правильности врачебного заключения [1].
Постановка задачи. Исследовать методы идентификации для устранения артефактов, способных искажать биогенные сигналы, без ухудшения качества исследуемых сигналов. Применить хорошо разработанные методы Фурье с помощью дискретного преобразования Фурье (ДПФ), используя алгоритм быстрого преобразования Фурье (БПФ). Применить вейвлет-преобразование и метод Прони для анализа и параметризации нестационарных сигналов. Провести оценку погрешности и сравнительный анализ полученных результатов.
Методы обработки нестационарных сигналов. Большинство биогенных сигналов имеет сложные частотно-временные характеристики. На рис.1 представлен нестационарный кардиосигнал биогенного характера, на примере которого рассматривается идентификация сигналов. Как правило, такие сигналы состоят из близких по времени, короткоживущих высокочастотных и долговременных, близких по частоте низкочастотных компонент.
Рис.1. Кардиосигнал (ЭКГ), содержащий шумы
Для анализа таких сигналов нужен метод, способный обеспечить хорошее разрешение и по частоте, и по времени. Первое условие требуется для локализации низкочастотных составляющих, второе - для разрешения компонент высокой частоты.
Преобразование Фурье - сигнал, заданный во временной области, в виде разложения по ортогональным базисным функциям (синусам и косинусам), при этом выделяются частотные компоненты [2]. Недостаток преобразования Фурье заключается в том, что частотные компоненты не могут быть локализованы во времени, что накладывает ограничения на применимость данного метода к ряду задач (например, в случае изучения динамики изменения частотных параметров сигнала на временном интервале).
Вейвлеты представляют собой семейства математических функций определенной формы, которые локальны во времени и по частоте, и в которых все функции получаются из одной базовой (порождающей) посредством ее сдвигов и растяжений по оси времени. Вейвлетные функции базиса позволяют сконцентрировать внимание на тех или иных локальных особенностях анализируемых процессов, которые не могут быть выявлены с помощью традиционных преобразований Фурье.
Метод Прони - это метод анализа коротких отрезков сигнала, основанный на аппроксимации сигнала конечной суммой комплексных экспонент (рис. 2-3), т.е. на подгонке экспоненциальной модели к измеренным эквидистантным (равноотстоящим, равноудаленным) значениям и последующем вычислении дополнительных значений посредством оценки параметров этой экспоненциальной модели в промежуточных точках.
Рис. 2. Короткий отрезок кардиосигнала, обработанный методом Прони (сумма комплексных экспонент)
Рис. 3. Спектр короткого отрезка сигнала и фаза после анализа методом Прони
Краткий обзор преобразования Фурье. Классическим методом частотного анализа сигналов является прямое преобразование Фурье:
ГО
£ (Ш) = | 5 ^ )е ~ Jatdt . (1)
— ГО
Обратное преобразование Фурье:
1 ГО
5(t) =----- Г5(ю)е]юtdю . (2)
2 л ■>
— ГО
Результат преобразования Фурье - амплитудно-частотный спектр, по которому можно определить присутствие некоторой частоты в исследуемом сигнале (рис.4).
Рис. 4. Фурье-спектр сигнала
Преобразование Фурье имеет ряд недостатков:
- исходный сигнал заменяется на периодический с периодом, равным длительности исследуемого сигнала;
- со временем плохо работает при изменении параметров процесса (нестационарности), поскольку дает усредненные коэффициенты для всего исследуемого образца;
- не дает представления о динамике изменения спектрального состава сигнала. Рассмотрим три подхода к анализу нестационарных сигналов биогенного происхождения -
метод оконного преобразования Фурье, метод вейвлет-преобразования и метод Прони.
Оконное преобразования Фурье. Классическое преобразование Фурье связано со спектром сигнала, взятым во всем диапазоне существования переменной. При анализе нестационарных сигналов интерес представляет только локальное распределение частот, в то время как требуется сохранить изначальную переменную (обычно время). В этом случае используется оконное преобразование Фурье, представляющее собой сегментирование сигнала на фрагменты («окна»). Для начала необходимо выбрать некоторую оконную функцию:
(3)
где F(t, ю) - распределение частот (несколько искаженное) части оригинального сигнала /^) в окрестности времени t.
Недостатки оконного преобразования Фурье. Когда нестационарный сигнал имеет высокочастотные компоненты в течение короткого промежутка времени и низкочастотные колебания при рассмотрении больших временных масштабов, оконные преобразования позволяют проанализировать либо высокие частоты в коротком окне времени, либо низкочастотную компоненту, но не оба колебания одновременно (рис.5).
Рис.5. Кардиосигнал после обработки методом Фурье
В результате был разработан подход, при котором для различных диапазонов частот использовались временные окна различной длительности. Оконные функции получались после растяжения-сжатия и смещения по времени гауссиана. Эти базисные вейвлет-функции называются компактными волнами.
Функция у (0 может называться вейвлетом, если она удовлетворяет следующим свойствам:
- область определения функции - вся числовая ось;
- интеграл от функции у (^ по области определения равен нулю;
- функция у (^ быстро убывает при ^±±о .
На рис.6 изображены примеры простейших базисных функций, посредством которых могут быть представлены произвольные колебания. Вейвлет-разложение по ним проводится путем вычисления сверток сигнала/ с компактной волной у при различных масштабах и сдвигах аргумента.
Рис.6. Компактные волны, полученные дифференцированием функции Гаусса
Рассмотрим свертку произвольного сигнала f (t) с функцией:
ад
w (т, а ) = 1/а | f (t)у (t - т/а )dt, (4)
—ад
где т - дрейфовая переменная; а - масштабная переменная.
Результатом свертки является функция, определенная в частотно-временной области, в отличие от преобразования Фурье, которое определено только в частотной области. Таким образом, вейвлет-преобразование позволяет судить не только о частотном спектре сигнала, но также о том, в какой момент времени появилась та или иная гармоника.
Параметр 1/а - частотная характеристика сигнала, подвергшегося вейвлет-преобразованию, привязана к «временной точке» - т. Это означает, что с помощью вейвлет-анализа мы можем изучать сигналы, частотные характеристики которых локально зависят от времени. Фурье-анализ такой возможности не дает.
Как и Фурье, вейвлет-преобразование имеет обратное преобразование, которое записывается формулой:
ад ад
f ( * ) = g 11 у (* — t/a) w (t, а) d^/а2 dt, (5)
—ад 0
где g - некоторая константа [5].
Вейвлет-преобразование Добеши (рис.7).
Рис.7. Структура вейвлета Добеши
fQ С С2 Сз
с, - С2 С1 - Со
Со С1 С2 Сз
Сз - С2 С1 - Со
. • • Со С1 С2 зС
• • • С, -С2 С1 -Со
С2 Сз • • Со С1
1С1 - Со • • Сз -С2
- матрица Добеши.
Преобразование Добеши от некоторого вектора соответствует умножению его на матрицу.
Каждая строка матрицы соответствует свертке вектора сигнала с локализованным вейвлет-фильтром. При этом нечетные строки сглаживают, а соответствующие четные фильтры вычисляют поправки к сглаживанию, обеспечивая обратимость преобразования. Коэффициенты С0...С3 могут быть выбраны так, чтобы поправки обращались в нуль для гладких функций. За-
нуление свертки с постоянной и линейной функциями дает два условия на неизвестные коэффициенты:
с3 - с2 + с - с0 = о
ос3 - 1С2 + 2С, - 3С0 = 0.
-2
'0
Далее потребуем, чтобы преобразование было ортогональным - обратное преобразование получалось простым транспонированием. Из условия равенства произведения матриц прямого и обратного преобразования единичной матрице получаем еще два условия [4]:
^2 . . /"»2
С 20 + С 21 + С \ + С 23 = 1
С2 С0 + С3С1 = 0
Откуда
С0 =(1 ^л/з)мТ2 С =( з + 7з) /^л/2
с2 =( з -л/з) мТ2 С3 =(1 -л/з)/4>/2 .
Полное вейвлет-преобразование получается путем умножения вектора сигнала на полученную матрицу Добеши, далее нечетные компоненты сигнала (сглаженная версия сигнала) снова умножаются на матрицу Добеши (вдвое меньшей размерности), и так до тех пор, пока не будет получено полностью сглаженное представление, т.е. среднее значение сигнала (рис.8-9).
Рис.8. Кардиосигнал, обработанный вейвлет-методом
Рис.9. График кардиосигнала и спектрограмма
Исходный подход Прони. Если число используемых отсчетов данных равно числу экспоненциальных параметров, то возможна точная подгонка экспонент под имеющиеся данные. Рассмотрим функцию дискретного времени, представляющую собой сумму /»-экспонент:
х[п\ = 2\4~1 .
(6)
к=1
Заметим, что в этом выражении используется х[п\, а не ~ [п], поскольку точно 2р-ком-плексных отсчетов х[1\, х[2р\ используется для точной подгонки к экспоненциальной модели с 2р-комплексными параметрами hh ...Мр Zl ..., zp. Входящие в р-уравнения (6) можно записать в матричной форме:
( 0 0 0 Л
г1 г 2 гр 'hl Л (х[1\ ^
7171 71 212 2 р 1^2 = х[2\
7р—1 7р—1 7р—1 ^ 1 ^2 ...¿р V / hp V р V х[ р\,
(7)
Матрица с временными индексами элементов z имеет структуру матрицы Вандермонда. Если может быть найден метод для раздельного определения элементов z, то уравнение (7) можно рассматривать как систему уравнений, решив которую, определяют неизвестный вектор комплексных амплитуд.
Ключ к разделению основан на том факте, что уравнение (7) является решением некоторого однородного линейного разностного уравнения с постоянными коэффициентами. Для того чтобы определить вид этого разностного уравнения, определим сначала полином ф (г), корнями
которого являются экспоненты гк :
Ф(г) = П (г _ гк) .
(8)
к=1
Если произведения в (13) выразить в виде степенной последовательности, то полином можно представить в следующем виде:
Ф( г) = 2 а[т\г
р—т
(9)
т=0
с комплексными коэффициентами а[т], для которых а[0]=1.
Осуществляя в уравнении (6) сдвиг индекса от п к (п-т) и умножая обе его части на параметр а[т], получаем
(10)
а[т\х[п — т\ = а[т\2 hkzkn~т—1 .
к=1
Записывая аналогичные произведения а[0\х[п\, ...,а[т-1\х[п-т+1\ и осуществляя суммирование, получаем:
р р р
2 а[т\х[п — т\ = 2 2 а[т\ггп—т—1 , (11)
т=0 I=0 т=0
которое справедливо при р+1 < п < 2р. Получаем:
]Га[т\х[п — т\ = 2^-7п—р ]Га[т\ггр—т—1 = 0 . (12)
т=0 г=0 т=0
Сумму в правой части (12) можно рассматривать как полином, определяемый уравнением
(9), который записан через свои корни, что и обеспечивает в уравнении (12) равенство нулю. Уравнение (12) - это линейное разностное уравнение, однородное решение которого выражается
формулой (6). Полином (9), ассоциированный с этим линейным разностным уравнением, называ-
ется характеристическим.
Можно записать в виде следующего рхр-матричного уравнения р-уравнения, представляющие истинные значения коэффициентов а[п], удовлетворяющих (12):
'х[р]х[р -1]......х[1] Vа[1] ^ (х[ Р л 1] Л
х[ р + 1]х[ р].
.х[2]
х[2р - 1]х[2р - 2].. .х[р]Да[р]
а[2]
х[ р +1]
х[ р + 2]
х[2 р]
(13)
Из уравнения (13) следует, что, имея 2р-отсчетов комплексных данных, возможно разделение множеств параметров ^ и zk. Комплексные полиномиальные коэффициенты а[1], ..., а[р], которые являются функциями только зависящих от времени компонентов zk экспоненциальной модели, позволяют по временным отсчетам сформировать соотношения для линейного предсказания. Матрица в уравнении (13) имеет тёплицеву структуру [3].
Процедуру Прони для подгонки р-экспонент к 2р-отсчетам данных можно теперь представить в виде следующих трех этапов.
На первом этапе получается решение уравнения (13) для коэффициентов полинома.
На втором этапе вычисляются корни полинома, определяемого уравнением (9). Используя корень zi можно определить коэффициент затухания а и частоту синусоиды f с помощью соотношений:
аі = 1п Ьг-
|/Т ( С -!) ,
(14)
fi = агС£ [1т {zi }/ Re {zi }] / 2%Т (Гц) . (15)
Для завершения процедуры Прони корни полинома, вычисленные на втором этапе, используются далее для формирования элементов матрицы уравнения (7), которое затем решается относительно р комплексных параметров ЭД1], ...., h[р]. Каждый параметр Ь используется далее для определения амплитуды А{ и начальной фазы , которые вычисляются с помощью ва-
риаций
А; = \К\ , (16)
0г = аг^[1т {к1 }/Re {hг }] (рад). (17)
Результаты вычислений показаны на рис.3, 10-11.
Рис.10. Кардиосигнал, обработанный методом Прони (сумма комплексных экспонент)
Рис.11. Спектр сигнала и фаза после анализа методом Прони
Метод наименьших квадратов Прони. На практике число отсчетов данных Ы, как правило, превышает то минимальное их количество, которое необходимо для подгонки модели из р-экспонент, т.е. N > 2р. В этом переопределенном случае последовательность отсчетов данных может быть аппроксимирована лишь как экспоненциальная последовательность:
р
~[п] = 2 ^Г1 , где 1 < П < N. (18)
к=1
Ошибка аппроксимации в данном случае распределяется выражением е[п] = х[п] - ~[п].
Одновременное нахождение порядка р и параметров {кк, zk}, где 1 < к <р, которые минимизируют
сумму квадратов ошибки и представляет собой трудную нелинейную задачу:
N
Р = £|в[п]|2 . (19)
п=1
Используя вариант метода Прони можно определить субоптимальное решение, которое обеспечивает удовлетворительные результаты. Используя на первом и втором этапе трехэтапного метода Прони соответствующие линейные процедуры наименьших квадратов, получим процедуру экспоненциального моделирования, которую иногда называют обобщенным методом Прони [3]. Выводы. Применены хорошо разработанные методы Фурье с помощью дискретного преобразования Фурье (ДПФ), использовался алгоритм быстрого преобразования Фурье (БПФ), вейвлет-преобразование и метод Прони для анализа и параметризации нестационарных сигналов. Проведена оценка погрешности и сравнительный анализ полученных результатов.
На основе результатов анализа ЭКГ человека показано, что вейвлетные базисы могут быть хорошо локализованными как по частоте, так и по времени. При выделении в сигналах хорошо локализованных разномасштабных процессов можно рассматривать только те масштабные уровни разложения, которые представляют интерес.
Вейвлетные базисы, в отличие от преобразования Фурье, имеют достаточно много разнообразных базовых функций, свойства которых ориентированы на решение различных задач. Базисные вейвлеты могут иметь и конечные, и бесконечные носители, реализуемые функциями различной гладкости.
Использование метода Прони позволяет получить более точные спектральные оценки для коротких последовательностей данных. Указанные возможности позволяют получить новую информацию о динамике изменений ЭКГ в процессе деятельности человека.
Библиографический список
1. Рангайан Р.М. Анализ биомедицинских сигналов. Практический подход / Р.М. Рангайан; пер. с англ.; под ред. А.П. Немирко. - М.: ФИЗМАТЛИТ, 2007. - 440 с.
2. Основы цифровой обработки сигналов: курс лекций / А.И. Солонина, Д.А. Улахович, С.М. Арбузов, Е.Б. Соловьева. - 2-е изд., испр. и перераб. - СПб.: БХВ-Петербург, 2005. - 768 с.
3. Марпл-мл. С.Л. Цифровой спектральный анализ и его приложения / С.Л. Марпл-мл.
- М.: Мир, 1990. - 584 с.
4. Добеши И. 10 лекций по вейвлетам / И. Добеши. - Ижевск: НИЦ «Регулярная и хаотическая динамика», 2001. - 464 с.
5. Терехов С.А. Вейвлеты и нейронные сети: лекции для школы-семинара «Современные проблемы нейроинформатики» / С.А. Терехов. - М.: МИФИ, 2010.
References
1. Rangaian R.M. Analiz biomedicinskih signalov. Prakticheskii podhod / R.M. Rangaian; per. s angl.; pod red. A.P. Nemirko. - M.: FIZMATLIT, 2007. - 440 s. - in Russian.
2. Osnovy cifrovoi obrabotki signalov: kurs lekcii / A.I. Solonina, D.A. Ulahovich, S.M. Arbuzov, E.B. Solov'eva. - 2-e izd., ispr. i pererab. - SPb.: BHV-Peterburg, 2005. - 768 s. - in Russian.
3. Marpl-ml. S.L. Cifrovoi spektral'nyi analiz i ego prilojeniya / S.L. Marpl-ml. - M.: Mir, 1990.
- 584 s. - in Russian.
4. Dobeshi I. 10 lekcii po veivletam / I. Dobeshi. - Ijevsk: NIC «Regulyarnaya i haoticheskaya dinamika», 2001. - 464 s. - in Russian.
5. Terehov S.A. Veivlety i neironnye seti: lekcii dlya shkoly-seminara «Sovremennye problemy neiroinformatiki» / S.A. Terehov. - M.: MIFI, 2010. - in Russian.
Материал поступил в редакцию 08.06.10.
A.V. KOLOMEYTSEVA, G.V. MISHUGOVA, A.P. MUL, G.Y. RYABYKH
APPLICATION OF WAVELET-TRANSFORM AND PRONY METHOD FOR THE BIOGENIC SIGNALS IDENTIFICATION
Non-stationary signals of biogenic character (cardiogram, eleсtroencephalogram, etc.) are considered. Analysis of signals with the use of wavelet-transform and Prony method - which unlike classical Fourier method make an improved spectral evaluation - is applied.
Key words: discrete cosine- and sine- transform, Fourier window transform, Daubechies wavelet-transform, Prony method, Toeplitz structures.
КОЛОМЕЙЦЕВА Анна Васильевна окончила Донской государственный технический университет (2010).
Область научных интересов: математическое моделирование, идентификация биогенных сигналов.
Автор 1 публикации. annakolomeytceva@mail.ru
МИШУГОВА Галина Васильевна, окончила Донской государственный технический университет (2010).
Область научных интересов: математическое моделирование, идентификация биогенных
сигналов.
Автор 1 публикации. galinamishugova@mail.ru
МУЛ Александра Павловна, старший преподаватель кафедры «Математика» Донского государственного технического университета. Окончила механико-математический факультет Ростовского государственного университета (1972).
Область научных интересов: теория функции, педагогика высшей школы.
Автор 25 публикаций.
РЯБЫХ Галина Юрьевна, кандидат физико-математических наук (1994), профессор (1996) кафедры «Математика» Донского государственного технического университета. Окончила механико-математический факультет Ростовского государственного университета (1976).
Область научных интересов: операторы свертки, алгебры сингулярных интегральных операторов со сдвигами и разрывными коэффициентами в пространствах суммируемых и гёльдеровских функций с весом, проблемы фредгольмовости, краевые задачи.
Автор более 32 публикаций.
ryabich@aaanet.ru
Anna V. KOLOMEYTSEVA graduated from the Maths Department, Don State Technical University (2010).
Research interests: mathematical simulation, biogene signals identification.
Author of 1 scientific publication.
annakolomeytceva@mail.ru
Galina V. MISHUGOVA graduated from the Maths Department, Don State Technical University (2010). Research interests: mathematical simulation, biogene signals identification.
Author of 1 scientific publication.
galinamishugova@mail.ru
Alexandra P. MUL, Senior Lecturer of the Mathematics Department, Don State Technical University. She graduated from the Faculty of Mechanics and Mathematics, Rostov State University (1972). Research interests: function theory, higher school pedagogics.
Author of 25 scientific publications.
Galina Y. RYABYKH, Professor of the Mathematics Department, Don State Technical University. Candidate in Physics and Maths (1994), Professor (1996). She graduated from the Mechanics and Mathematics Faculty, Rostov State University (1976).
Research interests: operators of convolution; аlgebras of singular integrated operators with shifts and explosive factors in spaces of summarised and Holder functions with weight; problems of Fredholm property; regional problems.
Author of more than 32 scientific publications.
ryabich@mail.ru
