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

УДК 681.518.22

Булатов В.Н., Худорожков О.В., Семёнов А.П.

Оренбургский государственный университет E-mail: promel@unpk.osu.ru

МЕТОДИКА ЦИФРОВОГО СПЕКТРАЛЬНОГО АНАЛИЗА СИГНАЛА С НЕРАВНОМЕРНЫМИ ПО ВРЕМЕНИ ВЫБОРКАМИ

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

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

В области современных информационно-измерительных систем (ИИС) источниками измерительной информации в подавляющих случаях являются датчики с выходными аналоговыми сигналами х(Ь), и вся последующая их обработка производится с помощью специализированных контроллеров, которые оперируют с «цифрой». Для «оцифровки» аналоговых сигналов применяется в основном аналого-цифровое преобразование (АЦП) с равномерной дискретизацией. Для цифровой обработки в частотной области (фильтрация, частотная коррекция и тому подобное) массивов выборок значений аналоговых сигналов, полученных в результате АЦП с равномерной дискретизацией, хорошо отработаны математические методы и их алгоритмизация, например, в виде быстрого преобразования Фурье (ПБФ). При этом массив выборок сигнала разбивается на так называемые окна с одинаковым временным интервалом Т и с помощью БПФ формируется линейчатый спектр «оконного» фрагмента сигнала в предположении, что данный фрагмент бесконечно повторяется с периодом Т.

Подобное АЦП сопровождается предельной методической погрешностью Дкв - погрешностью квантования А Кв =х+1-х=х.+2-х-1=... (рисунок 1), так как моменты выборок определяются частотой считывания (с постоянным интервалом дискретизации) цифрового кода N которые с вероятностью, близкой к единице, не соответствуют моментам перехода сигнала х(Ь) через ¿-уровень квантования.

В данной работе рассматривается цифровой спектральный анализ по совокупности выборок [Щ] цифровых значений сигнала х(Ь), которые зарегистрированы по моментам ^ перехода сигнала х(р) через уровни квантования х.=х(1;), при этом формируется одномерный массив пар чисел [Щ]. Подобный подход позволяет исключить методическую погрешность в виде А кв, но при этом интервал дискретизации для общего случая становится неравномерным, и спектральный анализ -по аналогии с БПФ - требует разработки специальной методики. Причем, данная методика может быть и полезна и для случаев, указанных в [2].

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

N(t о) = N о ;...;N(tI) = N,; N(tI+1) = N,+1; N(ti+2) = Nit2

(1)

Рисунок 1

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

^ -10 Ф ... Ф 11+1 -Ф 11+2 -11+1 Ф .... Ф ^ - ^ (2) - в том числе интерполяционные формулы Лагранжа и Ньютона [1].

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

N„(1) = N о +АП (1 -1 о) + Д 21 (1 -11 )(1 -1 о) +... +

( к-1

"А „1 П(1 - О = N0 + £ А к1 П(1 -1,)

(3)

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

- для разделенной разности А11 :

А„ =

N(1.) - N(1 о)

(4)

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

А12 =

А13 =

N(12) - N(11) . 12 - 11 ' N(13) - N(12),

А _ N(0 - Ж1т-1) . А1т _ 1 -1 '

т т-1

Ащ _

N(1n) - N(1n-l) .

1 - 1

1п 1п-1

(5)

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

А _ А12 -А11 .

21 . , ;

(6)

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

А _ А13 - А12

А_

13 -11

А _ А14 -А13

А_

14 - 12

А 2т _-

А

1(т+1)

-А

1 -1

т+1 т-1

А

2(п-1)

А1п -А1(п-1) 1 -1

(7)

для разделенной разности А31 : А -А

а _ ^ 22 21 . ¿А,. ■

13 - 1о

(8)

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

А _ А 23 -А 22 .

А32 _ 1 -1 .

14 11

А _А 24 -А 23 . А33 _ .

15 - 12

А3т _

3т

А - А

2(т+1) 2т

1 -1

т+ 2 т-1

А

3(п-2 )

А - А

2(п-1) 2(п-2)

1п 1п-3

(9)

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

А ш _

А - А

(п-1)2 (п-1)1

1п - 1„

(10)

Из анализа процесса составления выражений (3)-(10) для разделенных разностей следует, что для вычисления их численных значений можно составить блок-схему алгоритма с рекуррентной формулой вычисления двухмерного массива с элементами Акт и выделение из него массива разделенных разностей Акд (рисунок 2). В качестве примера реализации данного алгоритма на рисунке 3 приведена программа в среде МаШсаё для вычисления конечных разностей Акд при шести парах выборок {ЛСД.

_о

к_1

11 - 1о

13 - 12

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

n-1

П (t - tj ), входящего в выражение (3) и форми-

i=0

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

F(t) = ¿a itj = (1 + A ot) • (1 + Ait) •... • (1 + An-it) =

i=0

n-1

=П(1+Ait), (ii)

i=0

где a i являются коэффициентами производящей функции (11), содержащими информацию

о сочетаниях числом

/ \ n

i

V /

без повторений из n

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

I o

I n '

II

I n '

I 2

i n '

a o = 1;

a1 = A 0 + A1 + A2 +... + An-1;

a 2 = Ao -A^ + .A. o -A. 2 + -A^-A. 2 + ••• + -A-n—2 ■—•n-1 ;

a3 = Ao -A^-A. 2 + Ao •A-1-A-3 + -A-o •A'2 ■•'З + ••• + ■—■n-3 ■—■n-2 ■—■n-1 ; a n = o —2 ••'З • •••• • ••'n-2 1 •

(12)

Если в выражении (11) произвести замену Ai = -1/ti, то это выражение примет вид:

F(t) = П(1 - Ait) = П(1 - t/ti) =

i=o t-t

i=o

(-1)n

n+1 n-1

t ttt • • t хП(t-ti)- (13)

ti tot1t2 ... tn-1 i=o

П (-1)

i=o

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

число, а второй сомножитель в виде

П (t - ti)

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

F1(t) = П(1 - О = (-1Г1 о 1112 • ...• 1п-1 х

хЦ (1 - t/ti).

i=o

(14)

После подстановки в (12) А; = -1 и формирования сумм произведений с общим знаменателем ^^2 •... • 1п-1 (который впоследствии сокращается) окончательно получим выражение для производящей функции И (Ь):

F1(t) = Xв= (1 - 1о )(1 - 11 )(1 - 12 ) • ... • (1 - 1п-1 ) ,

(15)

где для соответствующих сочетаний коэффициенты в к=п-;, полученные из коэффициентов а;

' 0 1 " 0 "

-1.45 1.1

-3 3.03

t :=

2.33 4.43

3 6.74

, 126 j , S

Рисунок 2

n := 5 k := 1..

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

i := I..11 - 1 Д0,1 := Nt

1-1

for к Е 1.. п for m s 1.. n - k + 1

^k-l.m-l ~ ^k-1.n

0 0 -1.45 -3 133 3

0 -1.31* -ОИБ 3.S07 029 -1JS

0 0.17 1J84 4).MS -0.46S 0

0 027- -0.414 0.097 0 0

0 -0.102 0.074 0 0 0

0.022 0 0 0 0

126^

^k.i =

0 -1.313

0 0.17

0 D.Z74

0 -D.1Q2

0.022

о J

Рисунок 3. Вид программы вычисления A k1 в среде Mathcad

3

=0

=0

в результате замены А, _ -1 /1,, принимают значения:

: воп]_ 1;

: Р1п1_-(1о +11 +12 +... + 1п-1);

2 |: Р2п1_ 1о11 + 1о12 + 1112 + 1о1э + 1113 + 1213 + ■■■ + 1п-21п-1.

I 3 I: Р^ _-(1о1112 + 1о1Л + 1о1213 + 111213 + ... + 1п-э1п-21п-1 ).

: Рпп1_ (-1)п1о111213 ■■■■■ ■ 1п-21п-1.

(16)

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

Теперь выражение (3) с учетом (16) примет вид:

^(1) _ N о + £ А к1 ¿р^1

(17)

Введя в (17) формальную величину разделенной разности А о1 _ N о и систематизируя множители при окончательно получим:

Л_А оДо] 1о +А11 (Р11]+во1] 1)+

Nn(t)Ак1 ¿р^1

+ А 21 (Р [22] +Р12]1 + Р[о2112 ) +

+ А31 (Р331 +Р[2311 + Р[[3112 +во3113 ) + ■■■ +

+а ш (рпп1+рпп-11+рпп1212+■■■+Р1п11п-1+воп11п) _

Г

о

¿А цР'11 Iо + ¿А цр1!11 I1 + ¿А в

\ ( +

12 +... +

Существует несколько алгоритмов формирования сочетаний на основе рекуррентных методов, аналогичных приведенному в [3], необходимых для вычисления (16). Авторский алгоритм от традиционных методов принципиально не отличается, а его особенности связаны только с решением конкретной задачи вычисления коэффициентов рк. Поэтому из-за большого объема описания этого алгоритма он в данной работе не приводится.

Методика цифрового спектрального преобразования, составляющая второй этап работы, основана на использовании общих теоретических положений спектрального анализа для аппроксимированных сигналов [4] и частного решения для спектральной плотности «оконного» сигнала в(Ь), аппроксимированного степенным полино-

мом вида

¿мк,

представленного на интервале

[ 1 о, 1 о + т ] п+1 выборками [5]:

- |Ю

Iа11! I(-1)к

1_о к _о

'¿М !(-1)к

11-

(-|ю)к(1 - к)!

(-|ю)к(1 - к)!

! а.1!

(I®У

ехр( - |Ю 1 о )

- _|ю

(20)

Для получения массива комплексных чисел - выборок из (20), необходимого для технологии цифровой обработки спектра, делается предположение, что данный фрагмент сигнала в(Ь) (как и в случае с БПФ) повторяется с периодом т. В этом случае в частотной области формируется линейчатый спектр с т с номерами гармоник т и интервалом между ними, равным 2п / т, который можно получить из (20) в виде следующего преобразования: . _ 8(ю_ 2пт/т) _

¿А пР!

[ 1

-(п-1)

1п-1 +

¿А ири

¿А иР1 -1п

1_п

_]!ак1к

1п _

(18)

где коэффициенты степенного полинома N (1) имеют вид

_ ¿А 11Р[-к.

(19)

_±_ Г

2пт I

¿а.т11! ¿

У а^ т11!—

' (г

1

(|2ят)"

С|2лт)к(1 - к)!

ехр(-|2пт1о / т).

(21)

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

о

1

к_о

1 о + т

о

т

1

_о

2

т

+

1

_п-1

к_о

_о

! а1 т11!'

ьо

к_о

1

В случае необходимости определения фазового спектра фазовый множитель ехр(-|2пт1 о / т) можно учитывать уже с использованием массива чисел, полученных по (21):

1_о

(21)

Ст _ Ст (о) ■ еХР(-|2пт1 о / т).

11.04.2014

Список литературы:

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

2. Булатов, В.Н. Спектрально-временной метод определения частоты Доплера на основе целенаправленного изменения времени / В.Н.Булатов, Н.А.Косарев, О.В.Худорожков // Вестник Оренбургского государственного университета. -Оренбург: ОГУ, 2011. - №1. - С. 192-196.

3. Мамонтов, Д.В. Алгоритм формирования комбинаций при расчете перестановок, размещений и сочетаний [Электронный ресурс] / Д.В. Мамонтов, С.Б. Волошин. - Режим доступа: http://www.voloshin-sb.rU/PortaIs/0/ Рошпк^/^гкз/ combinator.pdf. - Дата обращения 12.04.2014.

4. Булатов, В.Н. Метод оценки погрешности определения фазового спектра кусочно-аппроксимированного сигнала / В.Н.Булатов // Вестник Оренбургского государственного университета. - Оренбург: ОГУ, 2009. - №2. - С. 84-88.

5. Булатов, В.Н. Спектральный анализ цифровых сигналов с неравномерной дискретизацией / В.Н. Булатов, Е.С. Тимонов, Д.А. Даминов // Вестник Оренбургского государственного университета. - Оренбург: ОГУ, 2006. - №6. - Т.2. - С. 185-190.

Сведения об авторах: Булатов Виталий Николаевич, профессор кафедры промышленной электроники и информационно-измерительной техники Оренбургского государственного университета,

доктор технических наук, профессор 460018, г. Оренбург, Шарлыкское шоссе, 5, ауд. 15319, е-тай: bulat@o56.ru

Худорожков Олег Викторович, заведующий кафедрой промышленной электроники и информационно-измерительной техники Оренбургского государственного университета, кандидат технических наук, доцент 460018, г. Оренбург, Шарлыкское шоссе, 5, ауд. 15237, е-таП: oren.khov@mail.ru Семёнов Алексей Петрович, лаборант, студент кафедры промышленной электроники и информационно-измерительной техники Оренбургского государственного университета 460018, г. Оренбург, Шарлыкское шоссе, 5, ауд. 15315, е-mail: alexeiSP@ya.ru