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

Булатов В.Н., Тимонов Е.С., Даминов Д.А.

Оренбургский государственный университет

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

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

Цифровая обработка сигналов является основополагающей технологией преобразования сигналов в современных информационно-измерительных и управляющих системах (ИИУС). Существенная часть информационного преобразования, как правило, производится над образом сигнала в частотной области. Преимущество такого подхода особенно сказывается при выделении информативного колебания с конечной шириной спектра, представленного в смеси с широкополосным шумом.

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

Однако при решении ряда задач по цифровой обработке сигналов приходится сталкиваться с неравномерной дискретизацией сигналов, обусловленной, например, использованием АЦП с генерацией выходных кодов по моментам сравнения входного сигнала с уровнями квантования или нелинейными преобразованиями во времени цифрового сигнала с равномерной дискретизацией. Для решения ряда подобных задач авторами были разработаны частные методики спектральных преобразований с использованием ранее полученных теоретических результатов [1]. Но возрастающие требования к разрешающей способности обозначенного выше класса ИИУС требуют создания универсального аппарата преобразования сигнала в частотной области, разработки спектрально-временного метода преобразования

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

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

1-й этап - аппроксимация сигнала с неравномерной дискретизацией аналитическим выражением (интерполяция по множеству значений сигнала) с установленной погрешностью;

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

1 Выбор аппроксимации сигнала с неравномерной дискретизацией

Пусть на отрезке ?п] заданы п+1 точек

?0,?р?2,...,?п , которые называются узлами интерполяции, и имеются значения некоторого сигнала е($ в этих точках (рисунок 1):

е({с)=е0, е({1)=е1, е({2)=е2 >■■■■> е(0=еп- (1)

При этом в общем случае допускается, что

¡1 - ¿0 ф ¿2 - ¿1 ф ¿3 - ¿2 ф ■■■■ ф ¿п - ¿п-1 ■

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

т=е„т=е,т=*г.■.т=е, (2)

В такой общей постановке задача может иметь бесчисленное множество решений. Однако эта задача становится однозначной, если вместо произвольной функции /(^ искать полином Р () степени не выше п, удовлетворяющий условиям (2), то есть такой, чтобы:

РЛ)=*о> РА)=е, Р^=е2рп(0=еп (3)

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

Из известных интерполяционных формул (Бесселя, Стирлинга, Стеффенсена, Эверетта, Г аусса, Чебышева и других) для решения поставленной задачи подходят только интерполяционные формулы Лагранжа и Ньютона, так как остальные известные авторам интерполяционные формулы рассчитаны на равномерную или специальную дискретизацию [2].

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

Полином Ньютона строится на основе разделенных разностей, что при изменении значения п приводит к вычислению или исключению высшего порядка разделенных разностей [3].

Обозначим разделенные разности следующим образом.

Разделенная разность первого порядка функции е(0:

е(?0; ¡1) = А11 =

Є(І1) - є(І0 ).

Є^2 ) ~ е(?1) .

(4)

е(1п-1,1п) = Аь =

е(*п ) - е(?п-1)

Разделенная разность второго порядка функции е(ґ):

Рисунок 1

е(?0 ;?1;?2) = А21 = А*2 А11

А13 - А,,

Є(?1; ?2 ; ¡3 ) = А22 = 1 '

(5)

А1п А1(п—1)

Разделенная разность третьего порядка функции е(ґ):

е(?0; ¡1,12; ¿3) = А31 = А22 А21;

е(ч; ¡2 'А; ¡4) = А32 = А23 А22

¡4 ¡1

(6)

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

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

ч * А(п-1)2 А(п-1)1

е(?0; Ц;...;іп) = А п1 =. (7)

¡п ¡0

Интерполяционный многочлен Ньютона для неравных промежутков с учетом (4)-

(7) будет определяться алгебраическим многочленом п-й степени:

Рп (О = е(^)+(? - ¡0 )Ац + (г - ¡0)(? - ¡1) А21 +

+ ... + (1 - - Iп )Ап1. (8)

Это есть ни что иное, как степенной полином вида

Рп (0 = ап(" + ап-1(П 1 + ап-2(П 2 + ... + а0. (9)

1 -1 о

п п-2

I - I -7

п п-3

¡1 ¡0

¡2 ¡1

Из выражения (8) очевидно определение только двух коэффициентов:

an = А n1;

i-1

Ч = e(to) + В(-1У 4'1 П tj ] i=1 j=0

(10)

Формула для вычисления остальных коэффициентов а. для (9) при п>3 была получена в результате достаточно сложной систематизации сумм произведений из выражения

(8) для весовых множителей перед степенным аргументом ї , которая выглядит следующим образом:

= £[(-1)m+‘ Ki(m-i)Ая]

П=1

(11)

где коэффициенты К., принадлежат неполной двухмерной матрице значений, определяемые как:

Ki) = 1 при і є [1, n];

і

Ki1 = П tm при і є [1,(n -1)];

m=1

j-l

Kj = K(j-1)tj +ntk при j є [2,(n-1)];

(E)

4j-^i( j-1)‘ i+( j-1)^^(i-1) j

Из анализа выражений (4)-(7) и (10)-(12) следует, что процесс вычисления коэффициентов а. для степенного полинома вида (9) хорошо алгоритмизируется и, с точки зрения цифровой обработки сигнала е(^ по его выборкам е является эффективной основой для получения аналитического выражения аппроксимированного сигнала е(^ на интервале [^п].

Существенной проблемой при этом остается определение погрешности аппроксимации (9) функции сигнала е(^, представленного выборками с неравномерной дискретизацией. Методику оценки этой погрешности лучше разрабатывать для частотной области, учитывая, что как распознавание, так и определение информативных параметров в современных ИИУС производится в основном в частотной области - в достаточно узкой области информативной части спектра сигнала е() Для этого необходимо получить корректное решение спектрального преобразования для выражения вида (9).

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

По условиям интегрируемости степенных рядов [2] спектральную характеристику сигнала, представленного (9), можно записать как Б(ю) = Б0 (ю) + Б1 (ю) + Б2 (ю) + Б4 (ю) +

+... + S n (w) = £ Si (w)

(13)

где каждое слагаемое в общем случае можно представить интегральным преобразованием в виде первообразной (без учета пределов интегрирования):

S0 (w)= fa0 exp (- jwt )dt = exp( ) a0;

J - jw

Sj (w)=J a jt exp( - jw t )dt =

exp(- jwt)

jw

t-

jw

52 (w)= Ja2 t2 exp (- jw t)dt =

exp (- jw t) a f 12 21 + 2-1

a2 t +

- jw - jw (-jw)2

53 (w)=Ja313 exp (- jwt )dt =

exp (- jw t

jw

3 312 3 - 2t

t3-------------+ -

3-2-1

jw (- jw)2 (- jw)3

S4 (w)=Ja4t4 exp (- jwt)dt=

exp (- jwt)

- jw

a4 (t -

4t3

4 - 3t2

- jw (- jw)2

- 4 - 3 - 21 4 - 3 - 2-1

(- jw)3 + (- jw)4 И так далее. Учитывая очевидную закономерность формирования S;- (w), можно составить выражение спектральной плотности для последнего слагаемого atn:

Sn (w)= Jan tn exp (- jwt )dt =

exp

(-jw t]

jw

tn --

ntn 1 n

(n - 1)t

n-2

jw (- jw j2

n(n - 1)(n - 2)tn-3 n(n - 1)(n - 2)(n - 3)t' (- j w)3 (- j w)4

-4

-... + (-1)n

Подставив в выражение (13) полученные выше выражения для 8г- (ю) и произведя сис-

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

Б(ю) = £ Б , (ю) =

г =0

- jwt

I ati! I (-1)k

i-k

j— i=o k=o (-j—) (i - k)!

(14)

Полученное выражение (14) замечательно тем, что оно для функций вида (9):

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

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

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

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

т = гп - ?0 , (15)

на основе первообразной (14) получим выражение спектральной плотности для любого отрезка [0, х ] функции е() представленного п+1 выборками и отстоящего от начала координат на величину t0 (рисунок 1):

Б(ю) = ]Г Б,- (ю) =

„- jWt

i-k

—— I ai

- j— i=0

I (-1)k-----------------~k-Г-

k=o (- j—)k (i - k)!

I ai i=0

I (-1)k k=0

_i-k

I ai i=0

1

( j— )'

(- j—)k (i - k)! exp( - j— t0 )

- j— .

(16)

Программирование процесса вычисления по (16) вплоть до п = 31 показало эффективность вычислений в виде высокой производительности и точности для функций, представ-

ленных степенным многочленом со степенью меньшей, чем n (при таких соотношениях погрешность аппроксимации отсутствует).

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

3 Оценка погрешности аппроксимации гармонического колебания полиномом Ньютона

При исследовании характера погрешности аппроксимации была взята норма в виде отношения числа узлов интерполяции (моментов отсчетов), приходящейся на период синусоидального колебания, при этом крайние узлы совпадали соответственно с началом и окончанием периода колебаний. В результате исследования было получено семейство функций ошибки 5n (t) = e(t) - Pn (t), общий вид которой соответствует частному случаю на рисунке 2.

Значения амплитуд функции ошибки в зависимости от числа узлов интерполяции на периоде колебаний представлены в таблице 1.

Анализ полученных характеристик 5 n (t) показал, что, во-первых, функция ошибки по амплитуде резко убывает с увеличением числа n узлов интерполяции на одном и том же интервале интерполяции. Во-вторых, функция ошибки по амплитуде растет по мере приближения дискретизации к равномерной и достигает максимума при равномерной дискретизации, при этом функция ошибки становиться нечетной функцией с быстро затухающими (примерно -40 дБ/(21 /n)) от краев сечения колебаниями с периодом 21 /n.

Наибольший интерес для анализа влияния ошибки аппроксимации на процесс цифровой обработки сигнала e(t) представляет аналитическое выражение этой ошибки. Для нечетных значений n=7,9,11, ,17,19, которые являются наибо-

t0 + t

Рисунок 2

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

, ^0,0165«3 - 0,113га2 - 0,242и +1,93 5(0 = (-1) 2 -----

10

n-5

sin

exp

exp

- np t)

(17)

На рисунке 2 для n=15 пунктирной линией показан график зависимости (17).

Первый множитель в выражении (17) является весовым и зависит от числа узлов интерполяции (то есть, от числа отсчетов), который можно обозначить как

, 0,0165n3 - 0,113n2 - 0,242n +1,93 ,10Ч

gn = (-1) 2 --------------------------------^-■(18)

С учетом (18), а также учитывая быстрое затухание экспонент в выражении (17) в пределах интервала аппроксимации t, можно вывести выражение для спектральной плотности функции ошибки: t

Ser (®) = ¡gn sin

( np t (-np(t- t)N

sin exp

t V / _ t V /

exp

r - np t)Л

Ser (w) = gn ¡

np t

exp

e jwtdt»

np(t - t) t

e -j wtdt -

-gn JS

n t (- np t) 4

exp

t V / _ t V > _

e-jwtdt =

„ 2 (Wt ,

2ю cosí — | + j4 np v 2

np (wt

w cosí

t V2

. (wt siní —

V2

t

J np , 4

4| — | + ю

.(19)

Выражение (19) было проверено на функции sin(0,2pt) при различных значениях n на

основе сравнения с выражением

t

Ser(ю) = J(e(t) - Pn(t)) - j w‘dt (20)

0

с использованием численных методов интегрирования комплексных функций. При этом было установлено, что при равномерной дискретизации в области несущей спектры (19) и (20) полностью совпадают. На рисунке 3 представлен амплитудный спектр | Serl( ю )| при равномерной дискретизации.

При случайном распределении интервалов дискретизации с равномерным законом распределения амплитудный спектр ошибки достигает минимума. На рисунке 3 он представлен как |Ser2( ю )|. Это объясняется тем, что спектр функции ошибки при неравномерной дискретизации приближается к спектру «белого» шума, то есть становиться более равномерным, а не сосредотачивается в области частот, кратных частоте несущей.

Проведенное исследование по сравнению этих двух видов амплитудных спектров функции ошибки интерполяции полиномом Ньютона гармонической функции на интервале, равном ее периоду повторения, при различных значениях n показало, что существует устойчивое соотношение между их максимальными значениями в пределах 20-и процентной погрешности:

t

t

t

t

t

Таблица 1

n 7 9 11 13 15 17 19

1 d(n) 1 max 1,7 X 10'3 7,2 X 10'5 2,2 X 10'6 5 X 10'8 9 X 10'10 1,3 X 10'11 2,1 X 10'13

Рисунок 3

|Se,2 (М)! m

|Serl(®)l m

(21)

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

cos| -^2- |+ j4

©г

w— cos

- I 2

.(22)

Полученное выражение (22), по сути, представляет собою предельную функцию

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

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

Se,O Sn

2np

-

+ 4w4

+m4

(23)

Данная зависимость представлена на рисунке 3. Отличительное свойство этой функции - ее независимость от четности п.

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

w-

-

-

+w

Список использованной литературы:

1. Булатов В.Н., Дегтярев С.В. Метод выделения информативной части спектра зашумленных доплеровских сигналов с использованием нелинейной системы времени // Вестник ОГУ, 2004. - №2. С. 163 - 167.

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

3. Волков Е.А. Численные методы: Учебн. пособие для вузов. - 2-е изд., испр. - М.:Наука. Гл. ред. физ. мат. лит., 1987. - 248 с.