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

УДК 534.1

ВЫЯВЛЕНИЕ СКРЫТЫХ ПЕРИОДИЧНОСТЕЙ ЗАШУМЛЕННЫХ СИГНАЛОВ С ПОМОЩЬЮ МОДЕЛИ СТРУКТУРНОЙ ФУНКЦИИ В ФОРМЕ НЕПРЕРЫВНОЙ ДРОБИ

М. А. Новосельцева

DETECTION OF HIDDEN PERIODICITIES OF NOISY SIGNALS BY MEANS OF MODEL OF STRUCTURAL FUNCTION IN THE FORM OF CONTINUED FRACTION

M. A. Novoseltseva

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

In article the algorithm of detection of hidden periodicities in a signal containing an additive noise component is offered. The algorithm is based on the theory of continued fractions and structure analysis and allows to define automatically an amount of harmonics and their parameters without presence of the aprioristic information.

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

Keywords: hidden periodicity, noisy signal, continued fraction, structural function.

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

Развитие теории случайных процессов позволило дать проблеме выявления скрытых периодичностей более общую формулировку и дало возможность объединить в рамках единой концепции целый ряд ранее известных приёмов и одновременно выдвинуть новые методы.

Задача о выявлении скрытой периодичности в зашумленном сигнале может быть сформулирована следующим образом. Имеются измерения сигнала в равноотстоящие промежутки времени с шагом дискретизации Д1 - х(кАЬ) (к=0,1,2,...), представляющие результат наблюдения некоторого физического процесса хф. На основании общих соображений о существе изучаемого физического процесса может быть высказана гипотеза о том, что функция х(.() содержит слагаемое, представляющее собой периодическую функцию времени s(t) [9]: х(к АЬ) = в(кАЬ) + п(кАЬ) =

п (1)

= ^ А АЬ + ф) + г/(кАЬ),

где - круговая частота /-ой гармоники, А1 - амплитуда /-ой гармоники, ф - начальная фаза /-ой гармоники, п - число гармоник, п(кАЬ) - случайная компонента, представляющая собой непериодическую функцию и являющаяся помехой. При практическом анализе процессов фазовые углы ф обычно игнорируются [1], тогда

c(k At) = s(k At) + ц(к At) =

(2)

= ^ А АЬ) + п(кАЬ).

г :=1

Задачей методов выявления скрытых периодичностей является определение числа гармонических компонент и их параметров (круговая частота, амплитуда).

Необходимо отметить, что задача выявления скрытых периодичностей принципиально отлична от задачи разложения в ряд Фурье функции, заданной на интервале [0, Ь]. Хорошо известно, что практически любую функцию можно представить на интервале [0, Ь] рядом Фурье, считая её периодической с периодом равным длине интервала Ь. Для этого требуется лишь выполнение условий Дирихле, всегда справедливых для функций, являющихся записью реального физического процесса. Так как величина Ь, представляющая собой продолжительность обрабатываемой записи процесса, обычно ни в коей мере не связана с существом самого процесса, то такая обработка записи не позволяет, вообще говоря, полностью выявить характерные черты исследуемого явления. В методах выявления скрытых периодичностей период искомой компоненты не навязывается заранее, а определяется в процессе самого исследования.

Обзор основных разработок по данной проблеме был сделан М. Г. Серебренниковым и А. А. Перво-званским [9] ещё в 1965 г. Позднее ряд авторов в своих трудах [1, 2, 3] использовали приводимые ими методы, не внося особых изменений и модификаций. Данные методы выявления скрытых периодичностей обладают рядом недостатков. Одним из них является необходимость иметь априорную информацию о процессе х(р. Для большинства методов так же необходимо знать число периодических компонент. Сами алгоритмы являются достаточно сложными и громоздкими. Стоит отметить, что при изучении всех методов авторы уделяют внимание

П

вопросам точности и достоверности проводимых теоретических исследования. Однако практическому применению отводится меньшая роль. Также возникают сложности с определением параметров почти периодических сигналов [1, 9]. На практике почти периодические процессы встречаются достаточно часто и порождаются физическими явлениями, в которых одновременно действуют гармонические процессы, не связанные между собой. Примером такого процесса является вибрация многомоторного винтового самолета, в котором двигатели не синхронизированы [1].

Использование цепных дробей для выявления скрытых периодичностей впервые предложил Лагранж в 1772 году. Предложенный им метод решает эту задачу несколько окольным путём, при помощи цепных дробей и рядов специальных типов. Этот довольно громоздкий метод был малоизвестен и считался малопригодным на практике. Позднее метод был модифицирован, однако по-прежнему имел ряд недостатков: итерационная процедура определения числа содержащихся периодических компонент на основе перебора непрерывных цепных дробей требует выполнения значительного числа операций; наличие областей неопределенности при нахождении параметров гармоник приводит к невозможности получения модели сигнала.

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

Используем правильные С-дроби и модифицированный алгоритм В. Висковатова [4, 5, 6], а также структурный анализ для выявления скрытых периодичностей сигнала. Этот аппарат позволит автоматически определить количество гармоник и их параметры без наличия априорной информации. Данный метод может использоваться для широкого класса сигналов s(t): гармонических, полигармонических, почти периодических. Кроме того, существует возможность определить, относится ли сигнал к какому-либо из перечисленных классов и вообще, обладает ли он свойством периодичности как таковой.

Раздел теории случайных процессов, связанный с изучением случайных процессов на базе исследования структурных функций, носит название структурного анализа [8]. Структурная функция [6, 8], предложенная А. Н. Колмогоровым, имеет вид:

Сх (ь, Ь + т) = М{ х (Ь) - х (Ь + т)} , (3)

где х^) - некоторый случайный процесс. Очевидно, что функция всегда неотрицательна, четна и удовлетворяет условию Сх (Ь,0) = 0. Структурные функции отражают в своём поведении наличие осциллирующих компонент исследуемого случайного процесса, что может быть использовано для оператив-

ного определения параметров скрытых периодических составляющих процесса.

Практическое построение структурной функции более надежно по сравнению с корреляционной, поскольку на нее не влияют ошибки определения среднего значения процесса х^). Оно осуществляется по формуле:

1 М1 -к 2

Сх(к ) = ТТЛ £(х(гх(г + к)) , (4)

N - к г=1

где N1 - число измерений процесса х^). В дальнейшем будем называть операцию нахождения структурной функции по формуле (4) первым структурным преобразованием и обозначать:

С1 (к ) = Сх (к ).

В (4) в качестве исходных данных при нахождении структурной функции выступает сам процесс х(.1). Тогда по аналогии можно ввести второе структурное преобразование:

с2 (к) =

і

Ж,-к

N - к

Т.(с1 (г (-сі (г + к)),

'2 'ь г=1

N <^ (5)

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

Третье структурное преобразование имеет вид:

с3 (к ) =

1 N.-* ,2

1 £ ((і (-сх (г + к)),

N3 - к г = 1

N3 <N2.

(6)

Следовательно, т-ое структурное преобразование представляется в виде:

ст (к ) =

-к

і

Е(ст-((г (- ст-1 (г + к))

г =1

(7)

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

Получим модель третьего структурного преобразования сигнала на основании его значений, используя при этом теорию непрерывных дробей и в частности модифицированный алгоритм В. Висковатова [4, 5, 6]. По результатам значений третьего структурного преобразования (6) строится идентифицирующая матрица:

(—1) — строка 1 1 1 ... 1

(0) — строка С3(0) С1(ДЬ) СК2Д1) ... СКпДЬ) ...

1 — строка а1(0) а1 (ДЬ) а1 (2ДЬ) ... а1 (пДЬ) ...

а2 (0 ) а2 (ДЬ) а2 (2ДЬ) а ьо ... п > , (8)

т — строка ат (°; ат (ДЬ ) ат (2ДЬ) а 3 .. п Ь

в которой (-1)-строка содержит значения единичной риодическим. В противном случае, сигнал не отно-

функции 1(ь), а (0)-строка - значения третьего сится к вышеперечисленным классам и не имеет пе-

риода.

структурного преобразования С3 (пДЬ), Дг - шаг В случае наличия комплексных полюсов полу-

дискретизации, а элементы ат (пДЬ) последовательно определяются с помощью соотношения:

М) =

т —2

(П + ат—1 (П +

,—2(0)

,—1(0)

(9)

где а—1(п) = 1(пД£).

а0(п) = С®(nД^), т=1,2,3,..., п=0,1,2,...

Элементы первого столбца матрицы (8) порождают частные числители правильной С-дроби [5], что позволяет получить модель структурной функции сигнала в форме дискретной передаточной функции (ДПФ) виртуального объекта [6, 7]:

ао(0)а—1(0) ai (0)

1 ; Г-

С3 (1)*—1

1 +

а1 (0 ^ г 1

(10)

1 +

а2 (0)г

-1

1 +

а2(0)г

-1

При аппроксимации дробно-рациональной

функции конечного порядка в матрице (8) наблюдается появление нулевой строки, номер которой позволяет идентифицировать порядок функции. При изучаемом процессе вида (2) в матрице (8) нулевой является 4-я строка, отвечающая за порядок знаменателя. Если в некоторой /-той строке (1=0,1,2,...)

матрицы (8) конечное число к первых элементов

равны нулю, а последующие элементы отличны от

нуля, то необходимо осуществить сдвиг влево на к

элементов до появления в нулевом столбце ненулевого элемента и далее продолжить определение других элементов матрицы (8) по правилу (9). Для 1-той строки при восстановлении модели третьего структурного преобразования элемент а,1 (0) умножается

на г г .

Если ДПФ содержит только комплексные полюса, то сигнал является периодическим или почти пе-

ченная модель ДПФ (10) позволяет определить параметры гармоники с максимальной амплитудой. Далее приступают к нахождению круговой частоты ^ гармоники с максимальной амплитудой

^ =Д Ш^) , (11)

где г1 = и ± т - полюса ДПФ (10).

Амплитуду гармоники находим по формуле:

А1 = ^шах(3 (к )) / 2 . (12)

Тогда выявленная гармоника с максимальной амплитудой имеет вид:

81(кД Ь) = А1 вт(-ш 1кД Ь) . (13)

Поскольку число гармоник исходного сигнала неизвестно, необходимо провести следующую процедуру. Используя полученную модель гармоники (13), вычтем её из исходного сигнала х(кДг), тогда получим:

ж1(к Д Ь) = х (к Д Ь) — в1(к Д Ь), (14)

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

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

Первый этап

Ввод исходных данных, то есть измерений случайного процесса х(кД1), к = 0, И1 . Ввод счетчика гармоник 1=0.

а

Второй этап

Получение первого, второго, третьего структурных преобразований процесса х(кД1) по формулам (4), (5), (6).

Третий этап

Составление идентифицирующей матрицы В. Висковатова (8) на основе третьего структурного преобразования (6). Если 4-ая строка не является нулевой, то сигнал имеет і гармоник, представляется в виде (2), процедура выявления гармоник в х(кД1) закончена. Иначе осуществляется переход на следующий этап.

Четвёртый этап

Получение ДПФ в виде (10) на основе полученной матрицы В. Висковатова. Нахождение полюсов полученной ДПФ. Если ДПФ содержит только комплексные полюса, то сигнал обладает периодичностью, счетчик гармоник становится равным і=і+1, осуществляется переход на следующий этап. В противном случае, сигнал не обладает периодичностью и не может быть представлен в виде (2).

Пятый этап

Определение параметров гармоники: круговой частоты:

ї(0

2

1

0

-1

-2

-зН

ю- = ахг(гг), (15)

! ДЬ к 1

амплитуды

Аг = .^шах( (к )) / 2 . (16)

Тогда выявленная гармоника имеет вид:

811 (кДЬ) = А вт(югкДЬ). (17)

Вычитаем полученную гармонику (17) из исходного сигнала х(кДг) при 1=1:

х1(кДЬ) = х(кДЬ) — вг (кДЬ), (18)

а при 1>1 из х1(кДг):

хг+1(кДЬ) = хг (кДЬ) — вг (кДЬ) (19)

и возвращаемся ко второму этапу.

Пример

Рассмотрим случай, когда исходный процесс является почти периодическим

х (Ь) = ) + 2sin(^/2п^) + п(Ь), (20)

где п(г) - представляет собой белый шум с нулевым средним и дисперсией Бл=0.25. Выберем шаг Дг=0.2 исходя из условия 8Р-идентифицируемости [4, 5], а также зададим число измерений процесса N1=400. График процесса (15) приведён на рис. 1.

<1411

»

Рис. 1. График процесса (20)

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

1 1 1 1 1 1

102.5791 332.5411 518.7493 526.1131 348.4081 115.4487

-2.2418

-4.0571

-4.1289 -2.3965 -0.1254

1.4321 3.2153

-0.4355 -0.9932

« 0 « 0

4.0599

-1.2637

3.3405

ДПФ, рассчитанная по элементам первого столбца идентифицирующей матрицы, имеет вид: г + 0.9966

х2 - 1.2452* + 0.9763

Параметры гармоники с наибольшей амплитудой, рассчитанные по формулам (15), (16) равны:

ю, = 4.4454 , ^=2.0068.

Таким образом, определили гармонику с наибольшей амплитудой

в1(кД Ь) = 2.0068 sin(4.4454kД Ь).

Вычитаем из исходного сигнала х(к) выявленную гармонику 8а(к) по формуле (18): х1(к)=х(к)-81(к).

С Х](к) проводим те же вычисления, что и с х(к). Идентифицирующая матрица, вычисленная по значениям значений третьего структурного преобразования (6), имеет вид:

1 1 1 1 1 1

1.3459 1.8667 0.1983 0.7081 2.0612 0.7169

-0.3869 -0.8527 0.4739 -0.5315 0.4673

3.5906 1.3720 -0.8474 2.7393

-2.5858 -0.9887 0.6107

« 0 « 0

1 1 1

0.3470 10-6 0.6454 ■ 10-6 0.1047 ■ 10-

-0.8599 -2.0182 -2.3176

-0.4872 0.3229 1.5538

3.0100 5.8850 5.4544

— 2.6181 — 5.0019

Тогда ДПФ будет иметь вид:

г + 1.048

С(г) = —

г2 + 0.6179г + 1.0005

Параметры второй гармоники равны:

ю2 = 9.4240, А2=1.0038. Определили гармонику, у которой амплитуда меньше, чем у первой гармоники: в2(кДЬ) = 1.0038 вт(9.4240кДЬ).

Вычитая из сигнала х^к) найденную гармонику по формуле (19), получим: х2(к)=х1(к)-82(к).

Рассчитав первое, второе и третье структурные преобразования для процесса х2(к), вычислим идентифицирующую матрицу:

-1.5166

1.7978

— 0.6181

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

вт (кД Ь) = в1(кД Ь) + в2(кД Ь) =

= 2.0068 вт(4.4454кД Ь) + (21)

+1.0038 вт(9.4240кД Ь).

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

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

Литература

1. Бендат, Дж. Прикладной анализ случайных данных / Дж. Бендат, А. Пирсол. - М.: Мир, 1989. -540 с.

2. Бокс, Дж. Анализ временных рядов. Прогноз и управление / Дж. Бокс, Г. Дженкинс. - М.: Мир, 1974. - Вып. 1. - 406 с.

3. Волгин, В. В. Оценка корреляционных функций в промышленных системах управления / В. В. Волгин, Р. Н. Каримов. - М.: Энергия, 1979. -80 с.

4. Карташов, В. Я. Анализ и исследование ап-проксимационных свойств непрерывных дробей при решении задачи структурно-параметрической идентификации динамических объектов // Препринт № 22 / В. Я. Карташов. - Барнаул: Изд-во Алтайского госуниверситета, 1996. - 40 с.

5. Карташов, В. Я. Непрерывные дроби / В. Я. Карташов. - Кемерово: Изд-во Кемеровского госуниверситета, 1999. - 88 с.

6. Карташов, В. Я. Идентификация стохастических объектов: учеб. пособие / В. Я. Карташов, М. А. Новосельцева; ГОУ ВПО «Кемеровский государственный университет». - Кемерово, 2008. -104 с.

7. Патент 2399060 Российская Федерация,

МПК в01Я 23/16 (2006.01). Способ анализа многочастотных сигналов, содержащих скрытые периодичности / В. Я. Карташов, М. А. Новосельцева; заявитель и патентообладатель КемГУ. -

№ 2009114193/28; заявл. 14.04.2009; опубл.

10.09.2010, Бюл. N 25. - 12 с.

8. Романенко, А. Ф. Вопросы прикладного анализа случайных процессов / А. Ф. Романенко, Г. А. Сергеев. - М.: Советское радио, 1968. - 247 с.

9. Серебренников, М. Г. Выявление скрытых периодичностей / М. Г. Серебренников, А. А. Пер-возванский. - М.: Наука, 1965. - 244 с.