АНАЛИТИЧЕСКИЙ МЕТОД ОПРЕДЕЛЕНИЯ СПЕКТРА СИГНАЛА С НЕРАВНОМЕРНЫМИ ПО ВРЕМЕНИ ВЫБОРКАМИ Текст научной статьи по специальности «Математика»

УДК 519.216

A.С. Боровский, доктор технических наук, доцент кафедры управления и информатики в технических системах, ФГБОУ ВПО «Оренбургский государственный университет»

e-mail: borovski@mail.ru

B.Н. Булатов, доктор технических наук, профессор кафедры промышленной электроники и информационно-измерительной техники, ФГБОУ ВПО «Оренбургский государственный университет»

e-mail: vita.bulatov@yandex.ru

АНАЛИТИЧЕСКИЙ МЕТОД ОПРЕДЕЛЕНИЯ СПЕКТРА СИГНАЛА С НЕРАВНОМЕРНЫМИ ПО ВРЕМЕНИ ВЫБОРКАМИ

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

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

Вводная часть. В области современных информационно-измерительных систем (ИИС) перед обработкой вся аналоговая измерительная информация переводится в «цифру». Для «оцифровки» аналоговых сигналов применяется в основном аналого-цифровое преобразование (АЦП) с равномерной дискретизацией.

Для цифровой обработки в частотной области (фильтрация, частотная коррекция и тому подобное) массивов выборок значений аналоговых сигналов, полученных в результате АЦП с равномерной дискретизацией, хорошо отработаны математические методы и их алгоритмизация, например, в виде быстрого преобразования Фурье (БПФ), основанного на методе дискретного преобразования Фурье (ДПФ) [6]. При этом массив выборок сигнала разбивается на окна [7] (фрагменты сигнала) с временным интервалом Т, и с помощью БПФ формируется линейчатый спектр оконного сигнала в предположении, что данный фрагмент бесконечно повторяется с периодом Т.

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

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

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

Задача вычисления спектра фрагмента сигнала с неравномерной дискретизацией аналитическим методом решается в три этапа:

1-й этап - полиномиальная аппроксимация сигнала с неравномерной дискретизацией подходящим аналитическим выражением;

2-й этап - получение выражения спектральной характеристики в комплексном виде для аппроксимированного сигнала;

3-й этап - получение решения спектрального преобразования для вычисления массива выборок из спектра аппроксимированного фрагмента сигнала (получение линейчатого спектра).

1. Полиномиальная аппроксимация фрагмента сигнала с неравномерной дискретизацией

Для того чтобы использовать интегральное преобразование Фурье для функции, представленной одномерным массивом пар чисел на первом этапе формируется полином, аппроксимирующий функцию N(^1, для которой выполняется условие (1):

= ЛТ0; = К,; Щм) = Км; Щч2) = ЛТ,.+2;... (1)

Из области численных методов известны несколько методов интерполяции многочленами N(1) п-ой степени, удовлетворяющих условию (1) в п+1 узлах интерполяции - временных отсчетах где i = 0, 1, 2, ... , п, которые справедливы для случаев, когда

¿1 - 'о * - * ',+1 * ',+2 * •••• * К -К-г

(2)

= Ы0 + д„(? - г0) + Д2|(г - ?,)(?-г0)+...

(3)

¡=0 4=1 V ¡=0 )

где Ди, к = 1, 2, 3, ... , п - разделенные разности к-го порядка, которые определяются по известным формулам [5]:

- для основной разделенной разности первого порядка Дп:

Ап =

(4)

а последующие разделенные разности первого порядка, необходимые для вычисления разделенных разностей второго порядка и выше, вычисляются следующим образом:

Д12 =-

Аи =

1-4 Щг )~Щ2).

А1т =

3

1

Д.„ =

(5)

■ для разделенной разности Д21:

л _А12-А„. 21 ~ ^ ^ >

(6)

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

А _ А13~А12.

ч ч

К ним относятся интерполяционные формулы Лагранжа и Ньютона [5]. Анализ составления этих формул показывает, что более подходящим для решения обозначенной в начале статьи задачи является полином Ньютона. Это связано с тем, что при необходимости изменения значения п коэффициенты полинома Лагранжа надо рассчитывать опять сначала. Полином Ньютона составляется с использованием разделенных разностей, что при изменении значения п вызывает дополнительное вычисление или исключение разделенных разностей более высокой степени [5].

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

_ Аи-А13 . 23 " * '

4 2

А1(и+1) А1т 1 _ ¡т-Х

А1„ -А1(Л-1)

*2(п-1)

г.

п-2

(7)

■ для разделенной разности Д31:

_Д22-Д21.

5

(8)

Таким образом, алгоритм, заложенный в составление формул (3) - (9), позволяет составить выражение для вычисления разделенной разности и-порядка:

Д _ А(Я-1)2 А(«-1)1 "1 ~~ ^ ^

(10)

Для вычисления разделенных разностей по выражениям (3) - (10) можно составить блок-схему алгоритма с рекуррентной формулой вычисления двухмерного массива с элементами Д* и выделение из него массива разделенных разностей Акд (рисунок 1).

В качестве примера реализации данного алгоритма на рисунке 2 приведена программа в среде MathCAD (для данной статьи среда программирования Mathcad выбрана для наглядности) с использованием встроенных операторов [4] для вычисления конечных разностей Д^, при семи парах выборок {Ж;/} с неравномерной дискретизацией из сигнала ЛТ(0 = со8(0,2-;г-г) на интервале Т=10.

п := 6

F(0 = £a/=0+4,0-0+40'-•0+4-i0 = fl0+40,

(11)

N

( 1 1 ГоЛ

0.309 2

-0.309 3

-1 t := 5

-0.809 6

0.809 9

v 1 V v10y

к := 1 ..II Формальные конечные разности для нулевой строки двухмерного массива:

I := 1 ..п -ь 1

Дп ■ := N. 1

0,1 1-1

Д :=

for k е 1 ..п for mel..n k+l

Ai- < к, m

Л , , - ,

к-l,m+l к-l,m I к 1 ^т—1

Д

к 0..п

к, 1

-0.346

-0.091

0.036

-2.389-10-3

-3.265-10-4

6.435-10-5

Рисунок 2. Пример реализации алгоритма с рекуррентной формулой вычисления двухмерного массива в среде МаШСАБ

Рисунок 1. Блок-схема алгоритма с рекуррентной формулой вычисления двухмерного массива Д^

л—1

Способ вычисления произведения Wit-t,), яв-

i=0

ляющегося составной частью (3), реализован авторами на основе свойств одного из видов производящих функций [5]:

: а0= 1;

: a1=4,+4+4+-+4-i;

: а2 = AaAl + AqA2 + AlA1 +...+Ап_гАп_j; : a3 = A0AlA1 + A0A1A3 + A0A2A3 +... + Ап_гАп_2Ап_х;

(12)

• an — 4444 '•••• '4-г4-1-

При замене в выражении (11) Al =-l/tl выражение (11) примет вид:

^(0 = n0-40 = n0-i/i,) =

где а1 являются коэффициентами производящей функции (11), содержащими информацию о сочетаниях числом ^ без повторений из п объектов

А0, А1,..., Ап-1 в виде сумм произведений всех сочетаний:

¡=0 \л+1

Л-1 t f ( 1\Л+1 и-1

-1) ХП(,_,Д

i=0 Ч 0 12 ¡=0

(13)

Левый сомножитель в (13) - это просто число.

л-1

Второй сомножитель в виде П повторяет

произведение в (3), и согласно полученному выражению (13) оно является основой производящей функции F1(t) сочетаний для объектов {^ при соответствующей замене в (12) А, = -1/гг:

= = ИГ'МЛ ---^-1 кПО-^)- (14)

¡=0 1=0

После подстановки в (12) А1 =-Шг и формирования сумм произведений с общим знаменателем t0tlt2 который, как видно из (14), далее со-

кратится, получаем выражение для искомой производящей функции F1(t):

^(0=ЕЛ./(15)

г=о

где коэффициенты полученные из коэффициентов а, в результате замены Д- = -1/*„ для соответствующих сочетаний принимают значения:

: /Г=1;

: А1"1

■ Р™ = ад+ад+'А+ад+'а++ -+; : А1"1 = "(УА + УА + '0^3 +'А'з +

(16)

V 1=о

(17)

ч<=о

=1

ЕМ5

5А/Й =2>/,

Г=

где коэффициенты степенного полинома Рф имеют вид:

(19)

Для вычисления (16) был разработан алгоритм, который отличается от традиционных особенностями вычисления коэффициентов вк из-за плавающего п. Если сформировать матрицу из слагаемых (16), где номер строки будет соответствовать номеру столбца матрицы произведений, а слагаемые в строке (16) «слева-направо» будут расположены в соответствующем столбце в матрице произведений «сверху-вниз», начиная с нулевой строки, то на основе этой матрицы можно получить искомую матрицу для вычисления (16). Ниже (рисунок 3) приведен пример формирования этой матрицы в среде Mathcad в продолжение примера на рисунке 2 (без учета /6 - для упрощения).

р:=

Особенностью левой части выражения (16) является то, что в нее введен верхний индекс (в квадратных скобках) для указания плавающего числа п объектов, участвующих в соответствующих сочетаниях с числом объектов в сочетании, указанных в нижнем индексе, причем первое сочетание всегда начинается с объекта /0.

С учетом (16) выражение (3) будет выглядеть следующим образом:

к

1 -:

1 -

1 -

1 -

1 -

"П0Д П0,2 "П0,3 П0,4 "П0!5

1 2 3 4 5

2>,2

! — 0 1=0 1 = 0 1 = 0 ¡ = 0

2 5 9 14

2Х2 -2Хз 2Х,4 0

¡ = 0 1= 0 1 = 0 1 = 0

3 9 19

-2>1,1 0 0

1 — 0 ¡=0 1 — 0

4 14

0 0 0

1=0 1= 0

5

1=0 0 0 0 0

Введя в (17) для постоянной составляющей формальную величину разделенной разности Д01 = и систематизируя множители при /', окончательно получим полином Ньютона в виде степенного полинома:

ы>ч ¡=о ,

+...+Дж +№+/№ +0-Ч-) = (18)

Рисунок 3. Пример реализация вычисления (16) для среды Mathcad

Вычислив в среде Mathcad коэффициенты (19), заменив дискретную величину t на непрерывную /1, получим искомый полином Р(/1) (рисунок 4).

2. Решение спектрального преобразования сигнала, аппроксимированного степенным полиномом

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

Методика спектрального преобразования осно-

+ 84(®) + ... + 8,(в») = £8,(®),

(20)

У®

8[(<») = -

хр(-уд?*).

-у«

1

-]со

2

г2--+

-у® 2-1

-}со {-]а>)2)

83(ю)=|а,г3 ехр(-у'юг )<Й= ех^-уаи)

-У®

"«5

ЗГ

3-2? 3-2-1

-У ю (-Ую)2 (-у®)3

(ю)=| а 4 ехр (-7<м/)<Л= _ехр(->г) 4 4?3 _ 4-3?2

-7®

У® (-У®)2

4 -3 -2 Г 4-3-2-1 • +

(-У®)5 (-У®)4

8, (®)= /е. ехр(- ](о1 = _ ехр(-]<Ш) [ „_т"~1 и(и-ф"~2 _

-у®

(-У® )2

Рисунок 4. Вид полинома P(t1) и его оригинала N^1)

вана на использовании общих теоретических положений спектрального анализа для непрерывных функций [3] и частного решения для спектральной плотности оконного сигнала N(1), представленного на интервале [¿„,*0 + Г] п+1 выборками и аппроксимированного степенным полиномом

и

вида = По условиям интегрируемости

4-0

степенных рядов [1] функцию спектральной плотЯ

ности сигнала Р© можно записать как

4=0

8(й>) = в0 (е>) + в! (о) + в2 (а?) +

т[п-\)[п-2)Г3 г{п-\){п-2){п-Ъ)^ п!

' " "' ("У'®)"

Тя

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

ной функции вида Р(г) = ^ак1к::

п >

ы -У®« м (-М(1-к)\

(21)

где каждое слагаемое в общем случае можно получить в результате взятия интеграла Фурье, то есть получить решение в виде первообразной от интеграла соответствующего слагаемого а/:

80(ю)= \а0 ехр(-уйЯ )<И = 6Хр^ а0;

* — 1Л1

Полученное выражение (21) полезно тем, что оно для функций, представленных степенным полиномом:

- получено в виде первообразной для неопределенного интеграла;

- не содержит погрешности спектрального преобразования;

- прозрачно для алгоритмизации и программирования вычисления значений S(ю) для различных частот т.

При вычислении спектральной плотности «оконного» фрагмента сигнала, представленного функцией Р(^), с целью исключения накопления ошибок вычислений, связанных с возведением в степень больших значений кроме нормирования желательно привязывать окно к началу координат, используя теорему о смещении. Применяя понятие окна спектрального преобразования [7] длительностью

(22)

на основе (21) получаем выражение спектральной плотности для любого отрезка [0,Т] сигнала РО) представленного п+1 выборками и отстоящего от начала координат на величину

п л /

ы -7®« ы {-т(1-к)%

п 1 Уаг'!——

£ <М

(23)

-У®

3. Вычисление массива выборок из спектра аппроксимированного полиномом фрагмента сигнала

Для получения массива комплексных чисел -выборок из (23), необходимого для цифровой обработки спектра в ИИС, делается предположение, что фрагмент сигнала P(t) (как и в случае с ДПФ) повторяется с периодом T. В этом случае в частотной области формируется линейчатый спектр с с номерами гармоник т и интервалом между ними, равным 2л/Т, который можно получить из (23) в виде следующего преобразования [3]:

^ = 80'о = 2лт!Т) = ^р^-п^ш ^ (24)

т т ч

В этом случае временную функцию (17) на интервале Т можно представить в виде ряда Фурье в экспоненциальной форме:

Р( 0=

т=-оо

При определении постоянной составляющей спектра фрагмента оконного сигнала можно использовать выражение (24) для формальной гармоники с номером т=0 (в системах связи ИИС обычно равняется нулю): ,0+г

|>(0* у N.

и + 1

•(25)

Для гармоник с номерами тф-0 выражение спектра (24) преобразуется следующим образом:

_ 8(й) = 2ят/Т) ,0+г_

Спф0 - т 1<„ -

- ¡т1м!Т

- }2тт

-]т 2Я/Т

Цт

Т(-}т2л П)

Т<-к

1м 1(-1)

(~]т2я /ТУ (1-к)\ 1>~к

(—ут 2л / Т) (1 — к)1 ехр( -у'т 2ш0 /Т) -

1«. -

1

- )2тш1 1

]2лт)

±"4 К-1)

(-/т2я-/7')*(1-£)! Г

= 3

- ¡2лт) 1

2 пт

. 1 2 пт

К"1*

(-]т2пу{1-к)\ Г

Е а>Т'щ Е

?{г]т2м1ту{1-ку.

1 1 ^ а,ТЧ

ехр( -2л1а /Т) = ехр(-уот2Л0/Г)-

ехр(-]т2тП„ /Т) =

(26)

5>.г'я Е

_1_

■й (7'2ята)' ехр( -ут 2я-<0 / Г).

ехр( -уш 2гг<„ !Т) =

В частности, если начало окна отнести к началу временной оси, то выражение (26) приобретет для практики цифрового спектрального анализа простую и удобную для программирования форму:

Г

с0т — 7

1

2ят

I«*™ 1т7

1

(27)

.„'с/гяио'о-*)!,

В случае необходимости определения фазового спектра при т^0 фазовый множитель ехр(-у'2л7и?0 / Г) можно учитывать уже с использованием массива чисел, полученных по (27):

¿т=с0т-ехр(-./-2^0/Г). (28)

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

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

с0(т)

п

а.Т1 Н-г

1

Е-т-

к = 0 0-2-тг-т) -(1-к)!

2-тг-т

Амплитуда основного колебания : А1:=2|с0(1)| А! = 0.998 Коэффициент гармоник с учетом 2-5 гармоник:

1/(2|С0(2)|)2 + (2|с0(3)|)2 + (2|С0(4)|)2+(2|С0(5)|)2 : А1

- 2.941 х 10" 3

Рисунок 5. Фрагмент программы вычисления коэффициента гармоник ^ в среде Mathcad

Литература

1. Бронштейн, И.Н. Справочник по математике для инженеров и учащихся втузов / И.Н. Бронштейн, К.А. Семендяев. - Москва: Наука, 1986. - 544 с.

2. Булатов, В.Н. Применение спектрально-временных методов в доплеровских системах: монография / В.Н.Булатов, Н.А. Косарев, О.В. Худорожков. - Оренбург: ООО ИПК «Университет», 2012. - 196 с.

3. Гоноровский, И.С. Радиотехнические цепи и сигналы: учебник для вузов / И.С. Гоноровский. -Москва: Сов. радио, 1977. - 672 с.

4. Дьяконов, В.П. MathCAD-8.0 в математике, физике и в Internet / В.П. Дьяконов, И.В. Абраменкова. -Москва: Нолидж, 1998. - 352 с.

5. Корн, Г. Справочник по математике / Г. Корн, Т. Корн; под ред. И.Г. Арамановича. - Москва: Наука, 1978. - 832 с.

6. Нуссбаумер, Г. Быстрое преобразование Фурье и алгоритмы вычисления сверток / Г. Нуссбаумер; пер. Ю.Ф. Касимов - Москва: Радио и связь, 1985. - 246 с.

7. Херрис, С.М. Использование окон при гармоническом анализе методом дискретного преобразования Фурье / С.М. Херрис // ТИИРЭ. - 1968. - Т.61. - №1. - С. 95-96.