УДК (004.352.242)
ИССЛЕДОВАНИЕ МЕТОДА НЕПАРАМЕТРИЧЕСКИХ ОЦЕНОК ДЛИНЫ ПЕРИОДА И ПЕРИОДИЧЕСКОЙ СОСТАВЛЯЮЩЕЙ СИГНАЛА
В.В. Марченко, О.Г. Берестнева1, Д.В. Девятых1, Е.Ф. Суханова1
ЗАО «Элекард Девайсез», г. Томск Томский политехнический университет E-mail: ogb6@yandex.ru
Рассматриваются методы выделения скрытой периодической составляющей сигнала. В качестве основного и наиболее подходящего для использования в данной работе был выбран метод непараметрических оценок длины периода и периодических составляющих, так как именно этот метод позволяет выделять скрытые составляющие из любых периодических сигналов, а не только гармонических. Исследование метода проводилось на модельных сигналах, показана зависимость относительной погрешности метода от количества скрытых периодов и шага дискретизации.
Ключевые слова:
Статистические оценки, показатель разброса, показатель размаха, скрытые периодические составляющие, погрешность. Key words:
Statistical estimates, spread index, range index, hidden periodic components, inaccuracy.
Скрытые периодические процессы протекают практически во всех областях жизни и затрагивают множество сфер человеческой деятельности. Так, например, задачи выделения скрытых периодических составляющих возникают в геологии, сейсмологии, экономике, медицине. На сегодняшний день изучение методов поиска периодичностей, скрытых в наблюдаемых процессах, является по-прежнему актуальной задачей. Статья является продолжением цикла работ по данной тематике [1-3].
Основными параметрами ритмов являются амплитуда, период и фаза. Выход параметров ритмов за пределы нормы или появление их там, где они раньше не обнаруживались, как правило, связано с изменением свойств и возможной неисправностью системы. Существуют различные методы выделения частотных ритмов, такие как метод Стокса, спектральное представление сигнала с помощью ряда Фурье, корреляционный анализ и метод непараметрических оценок длины периода и периодических составляющих [4].
В работе рассмотрен метод непараметрических оценок длины периода и периодических составляющих, так как он позволяет выделять скрытые составляющие из сигналов любой периодической природы, а не только гармонической.
Целью исследования является изучение свойств и ограничений метода статистической оценки скрытой периодической составляющей на модельных сигналах. Объект исследования - сигналы, предполагающие наличие скрытой периодической составляющей, предмет исследования - периодические свойства сигналов, выделенные методом непараметрических оценок.
Кратко опишем суть выбранного метода. Рассмотрим достаточно широкий класс практически полезных непараметрических оценок длины периода и периодической составляющей во временных рядах. Во многих прикладных задачах рассматривают временной ряд (или случайный процесс)
у(0=х(0+е(0, где х(0 - детерминированная периодическая функция от времени ¡, т. е. х(/)=х(/+7) при некотором Т, где Т- длина периода, а е(¡) - «шумы», случайные погрешности, искажающие периодический сигнал. Требуется оценить (минимальную) длину периода Т= Т0 и периодическую составляющую х(г). При этом не предполагается, что функция х(/) входит в какое-либо параметрическое семейство, например конечных сумм синусов и косинусов, т. е. рассматривается задача непараметрического оценивания (минимальной) длины периода и периодической составляющей сигнала [5].
Пусть рассматриваемые функции у(/), х(/), е(¡) определены на отрезке [0; А]. При фиксированном Т рассмотрим «куски» сигнала у(/) на последовательных отрезках длины Т, т. е. на отрезках [0; Т], [Т;2Т], [2Т;3Т],... Удобно ввести последовательность функций на отрезке [0; Т], полученную сдвигами этих кусков к началу координат:
У1<г‘)=у(г), У2(0=у(?+Т), Уз(1)=У(1+2Т),...
Все они определены на отрезке [0;7]. Число этих функций равно числу полных периодов длины Т, укладывающихся на отрезке [0; А], т. е. равно целой части числа А/ Т. В дальнейшем из всех периодов будем рассматривать и оценивать только наименьший [5].
Существуют разные ситуации поведения функции и периодической составляющей:
1. Если Т=Т0 - истинный период (или кратный ему), и погрешности е(0 отсутствуют, то все введенные в предыдущем абзаце функции совпадают между собой и с периодической составляющей:
х(()=уМ=У2(()=Уз(()=... при всех / из [0;7]. При наличии погрешностей полного совпадения не будет. Однако отклонения определяются лишь шумами в различные моменты времени. При этом в качестве оценки периодической составляющей х(/) естественно
взять среднее арифметическое ycp(t) функций
Я(0, й(0, Уз(0,- [5].
2. Если же T отличается от истинного периода T0 (и кратных ему величин), то различия функций y1(t), y2(t), y3(t),... между собой определяются также и различием значений x(t) в точках, отстоящих друг от друга на интервалы, длина которых кратна Т[5].
3. Если число Т/Т0 иррационально, то можно показать, что значения t+mT(modT0), где m - натуральные числа такие, что t+mT<A, асимптотически (при росте A) равномерно заполняют отрезок [0;T0], а потому при выполнении соответствующих условий регулярности, например непрерывности периодической составляющей сигнала, функция y0f(t) приближается к константе - среднему значению периодического сигнала x(t) [5].
4. Если же число T/T0 рационально, то наблюдаем промежуточный случай, в котором уср(0, как можно показать, приближается к периодической функции с периодом T=T0/n при некотором натуральном п. Эта функция получена усреднением п последовательных участков длины T0/n периодического сигнала x(t). Она не является константой, хотя разброс ее значений меньше, чем для исходного периодического сигнала, поскольку T0 - минимальная длина периода [5]. Для оценивания T целесообразно ввести два показателя: показатель разброса F(T;Y) F(T;y1(t),y2(t),y3(t),...) множества функций y1(t), y2(t), y3(t),... на отрезке [0;T и показатель размаха G(T;^=G(T;yср(t)) функции y0f(t) на отрезке [0; T] [5].
Символ Y означает здесь, что показатели разброса и размаха строятся по функции y(t). При этом показатель разброса нацелен на оценку различий в значениях семейства функций при одном и том же значении аргумента. А показатель размаха - на различие значений одной и той же функции при разных значениях аргумента [5].
Для дискретного времени их можно адаптировать двумя способами: либо заменив sup на max, а интеграл на сумму; либо расширив область определения используемых функций на весь отрезок, например, соединив соседние точки отрезками или использовав для заполнения пропусков сплайны более высокого порядка [5].
Ввести показатели разброса F(T;Y)=F(T;y1(t), y2(t),y3(t),...) можно разными способами. Пусть k=[A/T]. Можно использовать различные функционалы супремумного типа (здесь и далее число слагаемых k не будем указывать в обозначении функционалов). Первым рассмотрим максимальный разброс непосредственно между значениями функций:
F,(T, Y ) =
= sup{ | y(t) - y (t)| , i, j = 1,2,..., 0 < t < T} [5].
Второй функционал супремумного типа будет учитывать не произвольные отклонения, а только отклонения от «средней функции», т. е. иметь вид
F2 (7,Y) = sup{ I y.(0 - yp(t) |, i = 1,2,..., k,0 < t < 7}.
Третий функционал показывает, какую зону «заметают» значения функций:
F (7, Y) = sup{ | y.(t), i = 1,2,..., k, 0 < t < 7} --inf{y (t), i -1,2,..., k, 0 < t < 7}.
В качестве показателя разброса представляется полезным рассмотреть то или иное отклонение совокупности функций Y(q), /=1,2,...,£ друг относительно друга. Можно сказать, что эти функции заполняют некую «трубку», которая тоньше всего при истинном значении периода T, а внутри нее проходит периодическая составляющая X(q)=x(t)=x(qT). Естественно рассмотреть различные функционалы интегрального типа. Например, можно проинтегрировать максимум модулей попарных разностей:
1
F,(7, Y) = |max{| Yi(q) - Yy(q) |, i, j =1,2,..., k}dq [5].
0
Вместо максимума можно проинтегрировать сумму:
1 k
F5(7, Y) = JXIY(q) - Yj(q) |dq.
0 j=1
Как и для функционалов супремумного типа, естественно рассмотреть показатели разброса относительно «средней функции»:
1
F6(7, Y) = Jmax{| Y, (q) - YcP (q) |, i=1,2,..., k}dq
0
1 k
F7(7,Y) = J£|Y(q) -Y„(q)dq [5].
0 j=1
Список показателей разброса можно существенно расширить.
Показатели размаха также можно ввести самыми различными способами. Например, можно рассмотреть такой показатель:
G1(7, Y) = sup{|ycp(t) | ,0 < t < 7} --inf{|Уср(t)| ,0 < t < 7}.
При практическом использовании рассматриваемых алгоритмов целесообразно учитывать дополнительные особенности реальных временных рядов. В частности, обратим внимание на неустойчивость супремумов по отношению к выбросам (резко выделяющимся наблюдениям) сравнительно с функционалами интегрального типа. Во многих ситуациях временные ряды дают резкие выбросы (всплески), которые затем, как правило, сглаживаются. Поэтому целесообразно в качестве показателей разброса и размаха использовать функционалы интегрального типа [5].
Выберем для реализации показатели разброса и размаха /¡и G1.
Приступим к исследованию выбранного метода на модельных сигналах и, прежде всего, убедимся в его способности находить ритмические составляющие в различных периодических сигналах.
100
В 50
§ 0
£ -50
-100
Л Л Л . 100 ® 50 § 0 S -50 -100 1 >
і \ 1\ /\ А / \
/ \1 \ /
V/ и 1 С
t, с
а)
50
частоті гц 20 б)
Рис. 1. а) сигнал синуса; б) оценка периодической составляющей сигнала синуса
100
& 50
Е
О
-50
-100
\ Л Л . 100 В 50 І 0 3 -50 -100
\ /\ /\ - \
\/ / \
V V 1
t, с а)
50
частота^ гц 20 Ю
Рис. 2. а) сигнал косинуса; б) оценка периодической составляющей сигнала косинуса
Для этого возьмем гармонические сигналы синуса и косинуса, а также негармонические в виде функций «ступенька», «пила» и тангенса.
Пусть сигнал задан в виде Лі)=А$іп(2ж(охі), где А - амплитуда сигнала, шх - частота, заданная равноотстоящими по аргументу і значениями (рис. 1, а). Выберем следующие значения: А=80, озх=9 и проверим метод на данном сигнале. Оценка частоты сигнала равна 9,6 Гц (рис. 1, б).
Пусть исследуемая функция Лі)=А-со$,(2жсохі), где А - амплитуда сигнала, шх - частота, задана равноотстоящими по аргументу і значениями. Вы-
берем следующие значения: А=80, a)x=9 (рис. 2, а). Оценка частоты сигнала равна 9,6 Гц (рис. 2, б).
Пусть исследуемая функция задана следующим образом
j A, если mod(t, T /2) > 0,5;
f () [A /2, если mod(t, T /2) < 0,5,
где A - амплитуда сигнала, mod (t, T/2) - остаток от деления, заданный равноотстоящими по аргументу t значениями (рис. 3, а). Выберем значения А=10, Т= 1 /9. Оценка частоты сигнала равна 9,9 Гц (рис. 3, б).
Рис. 3. а) сигнал ступенчатой функции; б) оценка периодической составляющей сигнала ступенчатой функции
а)
Рис. 4. а) сигнал функции «пила»; б) оценка периодической составляющей сигнала функции «пила»
Пусть исследуемая функция задана в виде «пилы» f(t)=A-mod(t,7), где А - амплитуда сигнала, mod (t, T) - остаток от деления, заданный равноотстоящими по аргументу t значениями (рис. 4, а). Выберем следующие значения: А=10, T=1/9. Оценка частоты сигнала равна 9,4 Гц (рис. 4, б).
Пусть исследуемая функция f(t)=Atg(2moxt), где А - амплитуда сигнала, а>х - частота, заданная равноотстоящими по аргументу tзначениями (рис. 5, а). Выберем следующие значения: А=10, юх=9. Оценка частоты сигнала равна 9,9 Гц (рис. 5, б).
Проведем исследование зависимости точности обнаружения скрытых параметров от количества периодов, содержащихся в реализации процесса, и шага дискретизации.
Так как гармонические сигналы является наиболее распространенными среди периодических сигналов, встречающихся в повседневной жизни, для дальнейших исследований возьмем функцию синуса.
Исследуем функцию f(t)=A-sin(2nwxt), где А -амплитуда сигнала, wx - частота, заданная равноотстоящими по аргументу t значениями. Выберем следующие значения: А=80, юх=9.
Известно, что для дискретизации сигнала без потери информации согласно теореме Котельникова частота отсчетов должна быть в 2 раза выше верхней граничной частоты спектра сигнала. В качестве оптимального шага, позволяющего, не сохраняя лишних данных, восстановить сигнал с любой заданной точностью, возьмём значение At^O^fmj-1, где fmx - максимальная частота, содержащаяся в спектре сигнала, измеренная в герцах.
Из табл. 1 видно, что при уменьшении шага дискретизации величина относительных погрешностей по частоте также уменьшается. При шаге дискретизации в 5-6 раз меньше оптимального, метод достигает оптимальной точности ~90 %.
Очевидно, что длительность анализируемого сигнала влияет на количество периодов, содержащихся в реализации. Поэтому интересно выяснить, как зависит точность обнаружения скрытых параметров от количества периодов, содержащихся в исследуемой реализации процесса.
При найденном оптимальном шаге дискретизации определим минимальную длину реализации. Результаты исследования представлены в табл. 2.
Таблица 1. Зависимость относительной погрешности метода непараметрических оценок от шага дискретизации
Шаг дискретизации Погрешность по частоте, %
At=At)m 33,33
А^А^т/2 33,33
At=Atau/3 20
А^Аи„,/4 14,28
At=Atопт/5 11,11
At=At)ni/6 9,08
Aí=Atопт/7 7,69
At=Atau/& 6,67
At=At,m/9 5,88
At=Atm/10 5,27
A^Atm/H 4,77
Анализируя полученные результаты можно сказать, что для точного выделения скрытой гармоники в исходной реализации должно содержаться, по крайней мере, два скрытых периода. Таким образом, следуя значениям табл. 2, можно сделать вывод, что погрешность по частоте метода непараметрических оценок составляет 9,09 %.
Таблица 2. Зависимость относительной погрешности метода непараметрических оценок от количества скрытых периодов
Количество периодов Погрешность по частоте, %
1 200
2 9,09
3 9,09
4 9,09
5 9,09
Стоит также отметить, что следует рациональным образом выбирать длину анализируемой зависимости. Помимо общего увеличения объема вычислений, анализ излишне длинного участка невыгоден еще и потому, что при неустойчивом периоде смешиваются в кучу все участки, где период изменяется, что приводит к неубедительным результатам.
Основные результаты, полученные в работе.
1. Выделены основные свойства метода: универсальность в смысле применимости к различным сигналам, а также гибкость за счет использования различных показателей разброса и размаха.
2. Проведено исследование метода на модельном сигнале на зависимость точности обнаружения скрытых параметров:
2.1) шаг дискретизации целесообразно брать в 5-6 раз меньше шага, определяемого по теореме Котельникова, так как дальнейшие уменьшения шага не приводят к существенному уменьшению погрешности;
2.2) для точного выделения скрытой гармоники в исходной реализации должно содержаться, по крайней мере, два скрытых периода, при этом погрешность по частоте метода непараметрических оценок составляет 9,09 %.
СПИСОК ЛИТЕРАТУРЫ
1. Абдулкина Н.Г., Алайцева С.В., Константинова Л.И., Коханов-ская Ю.Г., Кочегуров В.А., Марченко В.В., Степаненко Н.П. Применение геометрического метода анализа фазового портрета для оценки биоэлектрической активности головного мозга у подростков с дисфункцией гипоталамо-гипофизарно-тиреоидной системы // Вестник новых медицинских технологий. - 2009. - Т. 16. - № 1. - С. 14-17.
2. Кочегуров В.А., Константинова Л.И., Марченко В.В. Использование фазового отображения ЭЭГ для контроля состояния здоровья пациентов с тиреопатологией // Известия Томского
3. Наиболее точную оценку периодических сигналов имеет функция «пила». Оценка частоты сигнала данной функции составляет 9,4 Гц. Сигналы синуса и косинуса также имеют достаточно точную оценку, составляющую 9,6 Гц. Оценки частот сигналов тангенса и ступенчатой функции оказались мене точными, их оценки составили 9,9 Гц.
Работа выполнена при финансовой поддержке РГНФ, проект № 12—06—12057в «Создание системы алгоритмических и программных средств обработки, представления и анализа экспериментальных данных в социальных и медицинских исследованиях».
политехнического университета. - 2011 - Т. 319. - № 5. -С. 125-129.
3. Берестнева О.Г., Пеккер Я.С. Выявление скрытых закономерностей в сложных системах // Известия Томского политехнического университета. - 2009. - Т. 315. - № 5. - С. 138-143.
4. Магазинникова А.Л. Основы цифровой обработки сигналов. -Томск: Томский межвузовский центр дистанционного образования, 2002. - 129 с.
5. Орлов А.И. Прикладная статистика. - М.: Изд-во «Экзамен», 2004. - 190 с.