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

Интерполяция сигнала с нерегулярной дискретизацией

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

Аджемов С.С.,

начальник НИО МТУСИ, adjemov@srd.mtuci.ru Курахтенков Л.В.,

младший научный сотрудник МТУСИ,

kurakht@srd.mtuci.ru

Романов Э.Ю.

заведующий лабораторией МТУСИ, romanov@srd.mtuci.ru

Постановка задачи

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

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

Для простоты будем считать: Г :[0ДЬ К |/() < А при некотором A > 0, яирр / с [-а, а]. Имеется разбиение отрезка [0,1]: 0 < ^ <...

< г < 1, Ч;) = у,. При аппроксимации предполагается, что функция А периодически продолжается за пределы отрезка [0,1], так что

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

ї ()= X (апЯІП (2пп)+ ъ с°8 (2пп))= X спе2ш"

п=0 п=-О

Такая аппроксимация сразу дает набор значений для спектра сигнала.

Интерполяция тригонометрического

многочлена

Условие, при котором такая аппроксимация позволяет восстановить функцию f (т.е. она дает интерполяцию тригонометрического многочлена), хорошо известно. Здесь предполагается, что значения тригонометрического многочлена /^) в точках ^, і = 1,..., г заданы абсолютно точно.

Лемма 1. Тригонометрический многочлен

ї а) = X спе1піп

п=-О

однозначно определяется своими значениями в различных точках ^,..., .интервала длины 1 в том и только в том случае, когда выполняется следующее неравенство Котельникова:

г > 2О +1.

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

Следствие. Неравенство остается справедливым при рассмотрении сигнала на временном интервале любой длины І.

Интерполяция в виде sinc-ряда

Если неравенство Котельникова выполняется, то функцию f можно однозначно восстановить не только в виде тригонометрического многочлена, но и в виде так называемого sinc-ряда. Это следует из того факта, что в этом случае последовательность функций

gi(t) = slnc(2Q(t — ti)), где sinc(x) = sln(px) со- Px

ставляет фрейм [2] на подпространстве функций со спектром из [-Q,Q]. Например, если tn = nTдля некоторого T> 0 и всех n є Z (равномерная дискретизация), то теорема Котельникова дает представление сигнала в виде sinc-ряда

f (t )= 2TQ X f (tn )gn (t)

n=—^

Опишем аналогичный метод построения представления функции f(t) в виде sinc-ряда в нерегулярном конечном случае, т.е. метод вычисления коэффициентов dj представления

f (t) = 2TQ X dj sinc(2Q(t — tj))

j=—~

в случае конечной суммы. Пусть Ь = [Ь(] — вектор-столбец с компонентами bj = f (tj), j = 1,..., r. Составим Тёплицеву

матрицу

R = (tj — ti) j=i.r,;=i,..,r.

Тогда искомый вектор d = Ц] находится, как решение системы линейных уравнений Rd = Ь.

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

Аппроксимация функции тригонометрическим многочленом

Задача состоит в аппроксимации функции { тригонометрическим многочленом с заданной шириной спектра О. Напомним, что функция считается заданной на отрезке и при аппроксимации периодически продолжается на всю ось вне этого отрезка с периодом 1. Искомый тригонометрический многочлен Р будет, по методу наименьших квадратов, решением задачи минимизации суммы

XI )- Р(11 )|2 ^ тп

]

Теорема 1. Предположим, что для данного набора точек {}с[0,1] ассоциированная Тёплицева матрица

тк,1 = Xе ->У,к,1 = -а, -, а

;=1

обратима в С. Зафиксируем некоторый набор положительных весов м>1 > 0,1 = 1, ..., г. Тогда решение задачи минимизации суммы

XIP (!,)- f (', )2

j=1

может быть получено следующим образом:

1. Вычислить элементы Тёплицевой матрицы

-2п/(Л-IУJ

(TQ )k ,l =X’

j=1

k, l = —Q,..., Q

и вектор у є C(2Q+l): уk =

= X wjf (;j )

e iniktj, k = —Q,..., Q.

J=1

ческого многочлена

, от-

личного от

Pp (t),

выполняется

XI p, (o)-f (о-j <X| p (j )-f (j )2 -

J=l J=l

X| p (tj)—f t )■ -j=X p fc) -j —

J=1 J=1

—X(f (o яти- тер fc )>Xi f )2 >-

Теперь преобразуем каждое из трех слагаемых. Первое слагаемое является представлением мощности отсчетов:

Xp (tj )2 -j = X (X a^jX

j=l (k=-Q l=-Q

2ni (k-l ) j

k,l=-Q

= X aAX-je J = X aA(Tq\j =(a,TQa).

k ,l=-Q j=l k,l=-Q

Среднее слагаемое преобразуется следующим образом:

2. Вычислить решение аорІ є С(2О+1) системы линейных уравнений Тааор1 = у.

3. Составить искомый тригонометрический многочлен

РоР1 (*)= X а^е2^'.

к=-О

Доказательство

Докажем, что для любого^ ригонометри-

X P ('j у h )w = X ak X f (j>~2,ni4 = (у,a )■

J=1 k =-Q j = 1

Таким образом:

X|P (tJ)- y (tJ )2 w=<a ’r«a>-

J=1

-(a, у) -( у a)+X|y (tJ)2 wj ■

J=1

При Tq aopt = У имеем:

(P(fJ ) - f ('j )|2 -| Pp ('j ) - f ('j )|2 )j =

j=1

= (a,Tna) - (a,y)-(У,a) - (aop, ,Tnaop.) +

+ (aopt > У) + (у> Opt ) = (a - Op (a - aopt ^ ^ 0.

Полученное выражение всегда неотрицательно и обращается в 0 только при a = aopt, что и доказывает теорему.

С помощью теоремы 1 можно построить эффективный алгоритм аппроксимации.

Алгоритм 1. В предположениях теоремы

1 решение Popt может быть вычислено за O(QlogQ) операций по следующему алгоритму:

1. Вычислить элементы Тёплицевой матрицы Tq и вектора у с помощью быстрого преобразования Фурье. Заметим, что для вычисления Тёплицевой матрицы необходимо вычислить 4Q +1 ее элементов. На это потребуется O (Q log Q) арифметических действий.

2. Решить систему линейных уравнений TQa = у с помощью итеративного метода сопряженных градиентов, который гарантированно сходится за 2Q +1 итераций, но во многих случаях дает хорошее приближение значительно быстрее. Вначале определяются векторы

Op q» r0 е C(2Q+1): a0 = 0, r0 = q0 = У ■

Для каждого n > 1 определяются:

a„ = a„ , +

\r„-l, qn—l)

тр--------------г q—l,

T qn-l, qn—l)

(rn-l, qn—l) T

tq q,-l,

(TQ qn—l, qn—l)

q = r

in n

(rn, tq qn—^

'tf---------------\q—l

(tq qn-^ qn—l)

и искомая аппроксимация ищется в виде

Pn (t)= 'X an,ken■

k=-Q

После очередной итерации вычисляется расстояние

=1 |p- (:,)- f ('j ) I, w,\p. (<-)- f ('j )2

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

Самая сложная часть в этом шаге — умножение TQq^ Тёплицевой матрицы на вектор. Как будет показано далее, это возможно за O (Q log Q ) арифметических операций, благодаря вложению матрицы Tq в циркулянтную и выполнению нескольких быстрых дискретных преобразований Фурье.

Таким образом, данный шаг осуществим за O (k Q log Q) арифметических операций, где к — количество необходимых итераций.

3. В качестве искомой аппроксимации выбираем теперь действительную часть тригонометрического многочлена Pn(t). Здесь опять понадобится осуществить обратное дискретное быстрое преобразование Фурье, на которое понадобится O (Q log Q) операций.

Следствие: Если в предположениях алгоритма 1 матрица Tq необратима, то метод сопряженных градиентов дает решение в смысле наименьших квадратов.

Оптимизация алгоритма аппроксимации

Подбор весов

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

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

k=—Q

нести никакой информации об основной концентрации энергии. Соответственно энергия X|Р (^ ) , информацию про которую дают известные отсчеты, будет слабо связана с общей энергией сигнала

l = j|р (t)2 d;.

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

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

Идея подбора весов заключается в компенсации различных плотностей отсчетов во временной области.

Итак, пусть имеется разбиение:

0 < х1 <... < хг < 1.

Для сохранения периодичности введем Х0 = хг - 1 и х, + ] = х, + 1. Тогда определим веса следующим образом:

Теорема 2. Предположим, что

тах (х.+, - х ) = 8 < —^.

і=і.Л '+1 '> 2О

Тогда спектр матрицы То содержится в следующем интервале

у (Гд)с[(1 - 2Од)2, (1 + 2Од)2 ],

а число обусловленности оценивается как (1 + 2Од)2

(1 — 20д)2

Умножение тёплицевой матрицы на вектор

Определение. Матрица

C = (С,.,), 0 < i, j < m — 1 называется цир-кулянтной, если существуют такие числа

c0,..., cm-1, что Ci,j = Ci-jmodm. В частности, любая циркулянтная матрица является Тёплицевой, и, обратно, каждую Тёплицеву матрицу T размера nxn можно рассматривать как подматрицу циркулянтной матрицы порядка 2n

/ T и А

U T

CT =

' 2

То есть величина м. равна длине интервала, содержащего точку х;, с концами в средних точках + X и х±± +^1 .

2 2

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

В [3] была сформулирована и доказана следующая теорема.

где и — другая Тёплицева матрица порядка п, которая строится по Т следующим образом: если Т.= для всех 0 < < п, то и,.=

V где и0 = °, ик = к - п при к > 0 и ик = П- к

при к < 0. Тогда для вычисления произведения матрицы Т на произвольный п-вектор у, достаточно вычислить произведение циркулянтной матрицы

C

поскольку

у T'

0 иу _

Таким образом, при оценке асимптотической сложности умножения Тёплицевой матрицы на вектор, матрицу можно считать циркулянтной.

Как было показано в [6], известен алгоритм, основанный на быстром дискретном преобразовании Фурье, позволяющий умножать циркулянтную матрицу на вектор.

Теорема З. Умножение циркулянтной матрицы порядка m на m-вектор осуществимо за O (m log m ) операций.

Литература

1. Chang R.H., Ng M.K. Conjugate gradient methods for Toeplitz systems, SIAM Review, 38 (1996), 3. -pp. 427-482.

2. Duffin R. and Schaeffer A A class of nonharmonic Fourier series, Trans. Amer. Math. Soc., 72 (1952). -pp. 341-366.

3. Bertrand C. and Sehier P. A Novel Approach for Full Digital Modems Implementing Asynchronous Sampling Techniques. Proc. IEEE Global Telecommun. Conf., 2: 1320-1324, London, UK, November 1996.

4. Farokh A. Marvasti (editor), Nonuniform Sampling: Theory and Practice, Springer, 2001.

5. Feichtinger H.G., Grochenig K., Strohmer T. Efficient numerical methods in non-uniform sampling theory, Numerical Mathematics, 69 (1995). — pp. 423-440.

6. Strang G. A proposal for Toeplitz matrix calculation, Stud. Appl. Mat., 74 (1986). — pp. 171-176.