Цифровая обработка сигналов на основе интерполяции по Ньютону элементов базиса

4. Карташева, О. В. Использование адаптивной системы тестирования АСТ-Тест для контроля знаний при дистанционном изучении темы «Базы данных» // Мат-лы конф. «Информационные технологии в образовании». - 2001.

5. Компьютерные технологии в высшем образовании / ред. кол.: А. Н. Тихонов и др. - М.: Изд-во Мос-ков. гос. ун-та, 2004. - С. 272.

6. Мамиконов, А. Г. Принятие решений и информация / А. Г. Мамиконов. - М.: Наука, 1993. - С. 184.

7. Петрушин, В. А. Обучающие системы: архитектура и методы реализации (обзор) // Известия РАН. Техническая кибернетика. - 1993. - № 2. - С. 164-190.

8. Петрушин, В. А. Экспертно-обучающие системы / В. А. Петрушин. - Киев: Наукова Думка, 1992. -С. 196.

9. Приобретение знаний : пер.с япон. / под ред. С. Осуги, Ю. Саэки. - М.: Мир, 1990. - С. 304.

10. Рыбина, Г. В. Обучающие интегрированные экспертные системы: некоторые итоги и перспективы // Искусственный интеллект и принятие решений. - 2008. - № 1. - С. 22-46.

11. Стефанюк, В. Л. Теоретические аспекты разработки компьютерных систем обучения : учеб. пособие / В. Л. Стефанюк. - Саратов: СГУ, 1995. - С. 98.

УДК 519.6:681.3 ББК 74.262.21

Я. Е. Ромм, Л. Н. Аксайская

ЦИФРОВАЯ ОБРАБОТКА СИГНАЛОВ НА ОСНОВЕ ИНТЕРПОЛЯЦИИ ПО НЬЮТОНУ ЭЛЕМЕНТОВ БАЗИСА

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

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

Y. E. Romm, L. N. Aksayskaya

DIGITAL PROCESSING OF SIGNALS ON THE BASIS OF INTERPOLATION ON NEWTON OF ELEMENTS OF BASIS

Abstract. The computer scheme piecewise polynominal is stated to approximation of functions, derivatives and certain integrals on the basis of Newton's interpolyatsionny polynom with the minimum temporary complexity for any way set border of an error. On this basis the basis of discrete transformations of Fourier with a logarithmic assessment of temporary complexity is in parallel calculated at any number of counting, thus from computing algorithm multiplication operations are excluded. Besides, the full complex of algorithms of digital processing of signals with logarithmic temporary complexity is in parallel carried out.

Key words. Functions evaluation, parallel schemes, digital processing of a signalov.operatsiya of multiplication.

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

Обсуждаются аспекты построения последовательных и параллельных алгоритмов минимизации временной сложности компьютерного вычисления функций, производных и определенных интегралов на рассматриваемой основе. Синтезированные алгоритмы прилагаются к вычислению элементов базиса, коэффициентов, а также непосредственно к выполнению дискретных прямых и обратных преобразований Фурье. Помимо того, рассматривается их приложение к другим схемам цифровой обработки сигналов. Временная сложность T (R), где R - число процессоров, ниже оценивается на модели неветвящихся параллельных программ [7] и измеряется количеством последовательных шагов алгоритма без учета обмена.

Для функции одной действительной переменной вида

7 =

выбирается система непересекающихся подынтервалов равной длины:

г = 0

хг+1 -хг- = {¡-а^Р, 1=0,1, ...,Р-1, Р = 2к, А:е{ОД,...}.

При априори заданной границе £ погрешности аппроксимации функции (1) для каждого отдельного подынтервала из (2), (3) строится интерполяционный полином Ньютона степени п, где п выбирается минимальным для заданной точности приближения одновременно на всех подынтервалах. Полином преобразуется по дистрибутивности так, что принимает канонический вид для каждого подынтервала:

+аИ/х + а21/х2 +--- + ап1/х" ^ П = СОШ, / = О, ,Р - 1 .

Построение (4) выполняется для всех подынтервалов при условии:

|/0Рп^е, / = 0Д...,Р-1 .

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

(4). Если х е ^ . хг+1 , то, в обозначении V = хг+1 — хг, дешифрация индекса осуществляется по

формуле г = int | —— J. где int - целая часть числа, а - из (1). Если для рассматриваемой функции

вычислить и хранить коэффициенты для всех 2k подынтервалов, то время вычисления функции зависит только от степени полинома (4). По схеме Горнера значение этого полинома вычисляется с временной сложностью t Сj= п {¿у + 1С где ty,tc - время бинарного умножения и сложения.

За счет сужения подынтервала степень n можно сделать «сколь угодно» малой при соответственном возрастании P в (3). Схема минимизации степени полинома Ньютона конкретно осуществляется следующим образом. Если границы i -го подынтервала из (2), (3) обозначить как ai 0 , b 10,

bi0 ~ai0

шаг интерполяции - w. =-, то равноотстоящие узлы интерполяции на текущем шаге при-

п

мут вид

XIJ =аг0+ Jwi , J = п - 1 •

Полином Ньютона на i -м подынтервале записывается в виде

^rn О f

j =1 j k=0

где А'' у, Q - конечная разность у-го порядка в точке х, () .Пусть

тогда [2]

0

х — xi0

t =-

n Ajy-n 1

^rnC=f

j=1 j  k=0

w

Процесс приведения (8) к виду (4) влечет значения коэффициентов

ao,f = /О,о), aitf = Yjb4dij

¿Уи

где bjj =—, dij - коэффициенты полиномов вида Pnj\j=d0j+dljt + d2jt2+... + dnjtl

„ lJ

с натуральными корнями, входящими в состав полинома (8).

Значения степеней полинома Ньютона, интерполирующего функцию у = Их, ie [/2,1

Таблица 1

£

I

7

8

Во входном столбце таблицы г = 10 4 ДО 5, ... обозначает априори задаваемую границу

погрешности. Во входной строке k задает показатель степени Р = 2к количества подынтервалов из (2), (3). На пересечении строки, содержащей е, и столбца, содержащего к , указывается минимальное значение степени n , при которой функция аппроксимируется полиномом данной степени

с точностью до s на каждом из 2'',; подынтервалов. Пустующая клетка означает, что в используемой версии языка программирования s оказалось недостижимым для данного числа подынтервалов ни при одном значении n (в работе использован Delphi).

Данная схема минимизации степени полинома, интерполирующего функцию (1), позволяет построить схему вычисления определенного интеграла и производной [5]. Пусть имеют место соотношения (1) - (3). Вычисление интеграла от функции (1) по промежутку |a. b | по аддитивности

сводится к сумме интегралов по подынтервалам , хг+1 из (2), (3):

Р-1*7+1

jf \f (x)dx.

г= 1

Для г -го подынтервала справедливо соотношение

/(х)*Рт(х), V хе

где Рп I

b P-1x!+1

« £ jpm (x)dx .

a 7=1 Xj

Используется замена переменной (7). При ее выполнении значениям х = х1 и х = х(+| соответствуют 1 = 0 и / = и. С учетом 1'п, = у\'п, К. получится РпЛ (x)dx = wjx¥nj(t)dt. В результате

4+1

\pm(x)dx =wl\^m(t)dt.

Интеграл в правой части вычисляется непосредственно:

w, ^nl(t)dt = M>14>(n+l)l(t)\

где

(и+1) 1

Ос-

а0 if a\if 2 a2if 3 anif И+1

-—t+--Г +-—Г +... +-.

1 2 3 и+1

b

a

x

n

0

x

n

0

0

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

Подставляя верхние и нижние пределы интегрирования в правую часть (14), получаем

wг"*{п + \)

Учитывая (16), равенство (14) можно записать в виде

(п + 1) / (и) •

Значения ^(и+1)г («) вычисляются по схеме Горнера:

( fin „ .. Л

(и + 1)г

anif a(n-l)if П + —--1—

V W

п+1

п + -

(n-2)if П- 1

А

и + •

ог f

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

а0I/ аИ/ а2I/ ап1/

1 ' 2 ' 3 ' ... ' п + 1 в (15) можно рассчитать и занести в память компьютера для постоянного хранения, при этом вычисление интеграла в дальнейшем будет осуществляться путем выборки коэффициентов, соответствующих заданной функции и заданному подынтервалу, по аналогии с таблично-алгоритмической схемой аппроксимации функций. Более того, в данном частном случае можно хранить заранее вычисленные значения (18) для каждого номера подынтервала.

В общем случае, используя аддитивность интеграла по промежутку, с учетом (13) - (18), получаем значение интеграла на всем промежутке [а, Ь]:

р-1

\f О w i ^ ( Я +1) i (n) •

a i=1

Существенно, что (19) строится для минимизирующего значение показателя п алгоритма кусочно-полиномиальной аппроксимации. Этого значения, как показывает эксперимент, оказывается достаточным для достижения наперед заданной, сравнительно высокой, точности вычисления определенного интеграла.

Точность вычисления на основе (19) сопоставлялась с точностью формулы трапеций и метода Симпсона.

В обозначениях

- априори заданная граница погрешности

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

Ь - а

лах длины h =

2

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

Погрешности интегрирования для функции

Таблица 2

Дх) = sin(x) , 11, b =

II

2

г = 10"

n k 8 для метода на основе вариации полинома Ньютона степени п и числа подынтервалов 3 для метода трапеции для метода Симпсона

5 7 5,42*10-20 2,81*10-1° 7,04*10-18

4 13 1,08*10-19 2,81*10-1° 7,04*10-18

3 14 6,78*10-19 2,81*10-10 7,04*10-18

4 17 1,08*10-19 2,81 *1°-10 7,04*10-18

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

w

о

п

1

n

b

=о

19

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

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

а01/ аУ/ а2 i f anif

из (15), (18) являются хранимыми.

1 2 3 и + 1

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

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

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

/(х)*Рт(х), Ухе|г,хг+1^, где Р„ I (х) из (5). При интерполяции по Ньютону РП1 (х) = у¥П1 (/) в этом случае

Дх) «¥„,.

Взятие производной по х от обеих частей (20) влечет приближение (/(х))'х «(Тиг (/))'х • С учетом (7) получится

где 0¥ш (0)'г = «1г /• + 2^2; + З«3г Г+ ••• + « «иг , ?х = — .

Значение п 1 (I) вычисляется по схеме Горнера:

Ш 4 > > > 1

^ „,

Соотношение (21) дает формулу кусочно-полиномиального приближения производных. С целью экспериментальной оценки погрешности данного приближения на каждом подынтервале из (2), (3) в проверочных точках (7) выполнялось сравнение (21) с аналитическим значением производной. Результаты сравнения сведены в таблицу:

Таблица 3

Погрешность кусочно-полиномиальной аппроксимации производной функции

е = Ю-19

/{х) = $,т(х), |з, Ь 3=

1Д 2

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

11 0 9,048*10-17

10 1 -5,421 *10-20

10 2 0

8 3 9,926*10-1/

7 4 9,758*10-19

7 5 0

10 6 -5,963*10-19

5 7 7,589*10-19

5 8 -2,602*10-18

4 9 -2,711*10-19

4 10 5,421*10-19

4 11 -2,873*10-18

4 12 -7,589*10-19

4 13 -3,524*10-19

4 14 -9,487*10-19

4 15 -2,927*10-18

Здесь п - степень полинома Ньютона, 2к - число подынтервалов, при этом сама функция /(х) = 5т(х) вычислялась с точностью е = 10 Из табл. 3 видно, что с помощью предложенной

схемы производная функции вычисляется практически с точностью аппроксимации функции. Существенно, что функция и производная могут вычисляться одновременно. При построении библиотеки стандартных подпрограмм, для случая наперед известной функции, коэффициенты аппроксимирующего полинома могут быть хранимыми для каждого подынтервала. В этом случае производная, как и функция, будет вычисляться за время n сложений и умножений (в табл. 3 достигается п = 4 для £ = 10~19).

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

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

, ЛГ-1

При выполнении прямого дискретного преобразования Фурье Х

N п=0

N-1

и обратного ДПФ х(пТ)= ^X(k)ei2nk!N (ОДПФ) изложенный метод можно использовать для

£=0

динамического вычисления действительной и мнимой части элементов тригонометрического базиса el2nk/N = cos(2nk / N) + i sin(2/?/i / Y). По рассмотренной схеме кусочно-полиномиальной аппроксимации элементы базиса ДПФ и ОДПФ можно представить в виде алгебраического полинома с заданными коэффициентами

N

новременно для всех номеров подынтервалов г, все элементы базиса в форме (22), (23) могут быть вычислены параллельно и синхронно. Для наперед заданной точности аппроксимации степень каждого интерполирующего полинома (23) и (24) может быть выбрана сколь угодно малым натуральным числом, поэтому временная сложность параллельного вычисления всех элементов базиса оценивается из соотношения 1(11) = (9(1), где число процессоров К совпадает с количеством одновременно вычисляемых элементов базиса, Я - N .

Далее, вычисленные значения (22), (23) на тех же N процессорах одновременно умножаются на соответствующие коэффициенты х(пТ) ДПФ.

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

оценкой временной сложности: т

2 =[Гоё2Ж tc.

Отсюда время параллельного выполнения по рассматриваемой схеме ДПФ с произвольным числом отсчетов составит

Т{Н) = о

Аналогичная оценка имеет место для вычисления ОДПФ.

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

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

^ N-1

сху

где хп - выборка х(1), уп - выборка у(1), >„+/_- - последовательность, полученная сдвигом последовательности уп относительно х„ на к шагов, ох, <т г - среднеквадратичные значения х„

и Уп

л N-1 , N-1

п~0 п~0

В случае комплексных х(/) и у(/) корреляция равна

, N-1 1 *

--\хпУп+к •

*

М ^ (п

х У п=и

Если комплексные амплитудные коэффициенты величин х* и уп обозначить X т и Ут

* *

сп и уп обозначить X* то конечные фурье-преобразования Сх у (к) и X* У т можно представить в виде

ху у ' т т

N-1

Е

к=о \ " )

ы-1

у Сху (к) = Xх™ Тт ехР т =0

Ил т к

N

Функция ковариации представляет собой корреляционную функцию, умноженную на среднеквадратичные значения хп и уп :

со ч[х,у

Для двух конечных причинных последовательностей к(п) и у (к) длиной соответственно N и К свертка определяется выражением:

N

s

п= о

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

входе системы, к(п) - импульсный отклик системы, функция s(k) - выходной сигнал системы.

Для вычисления корреляции (25) целесообразно предварительно по схемам сдваивания параллельно вычислить значения левых частей (26). При этом квадратный корень аппроксимируется по

1

кусочно-полиномиальной схеме применительно к функции у = х2. После этого с сохранением оценки (24) по схеме сдваивания выполняется параллельное вычисление (25). Порядок оценки (24) сохраняется также в случае комплексных значений х(/) и у(/) при параллельном вычислении корреляции с использованием схемы сдваивания. В выражениях (28), (29) элементы базиса можно определить по предложенной параллельной схеме кусочно-полиномиальной аппроксимации. Очевидно, порядок рассматриваемой оценки временной сложности не изменится при вычислении ко-вариации (30). С распространением схемы сдваивания на операцию свертки (31) весь представленный набор вычислительных операций цифровой обработки оказывается выполнимым параллельно с временной сложностью порядка

т «г У о <[о82 ж ;

где N - произвольно задаваемое значение переменной (количество отсчетов).

Таким образом, фактически полный набор операций одномерной цифровой обработки сигналов с использованием предложенных схем можно выполнить параллельно с логарифмиче-

ской оценкой временной сложности. При этом из асимптотики схемы сдваивания остаются исключенными операции умножения.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Айфичер, Э. Цифровая обработка сигналов: практический подход / Э. Айфичер, Б. Джервис. - М. : Вильямс, 2004. - 989 с.

2. Березин, И. С. , Методы вычислений / И. С. Березин, Н. Г. Жидков. - М. : Наука, 1970. - 464 с.

3. Люк, Ю. Специальные математические функции и их аппроксимации / Ю. Люк; - М. : Мир, 1980. - 608 с.

4. Ромм, Я. Е. Бесконфликтные и устойчивые методы детерминированной параллельной обработки : дис. ... д-ра техн. наук. 05.990.001006 / Ромм Яков Евсеевич. - Таганрог : ТРТУ, 1998. - 546 с. ; ВНТИ Центр.

5. Ромм, Я. Е. Кусочно-полиномиальное вычисление функций, производных и определенных интегралов на основе интерполяции по Ньютону / Я. Е. Ромм, Л. Н. Аксайская. ;Таганрог, 2007. - 36 с. - Деп. в ВИНИТИ 02.05.07, № 487.

6. Ромм, Я. Е. Параллельная схема дискретного и быстрого преобразований Фурье на основе полиномиального представления базиса / Я. Е. Ромм, С. А. Фирсова // Известия РАН. Математическое моделирование. - 2006. - Т. 18. - № 11. - С. 3-12.

7. Солодовников, В. И. Верхние оценки сложности решения систем линейных уравнений // Теория сложности вычислений. - ЛОМИ АН СССР, 1982. - Т. 118. - С. 159-187. - I: Записки научных семинаров.

УДК 681.3.06:681.323(519.6) ББК 27.41

Я. Е. Ромм, Н. А. Брайцара

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

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

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

Y. E. Romm, N. A. Braytsara

PARALLEL SORTING ON THE BASIS OF THE BINARY TREE AND MERGE ON THE MATRICES COMPARISONS

Abstract. The article is devoted to the synthesis of parallel sorting on the basis of the tree select and merge. Building combines the properties and purpose of schemes of search on the basis of the tree and on the basis of merger for a single, parallel algorithm for sorting. The proposed algorithm is different resistance, a way of addressing the links to the input elements. Presented logarithmic assessments of the time complexity of the sort of the model is not branching of parallel programs.

Key words: parallel sorting, binary tree, parallel merger, the upper estimate of the time complexity.

Постановка вопроса. Ставится задача синтеза алгоритма устойчивой параллельной сортировки на основе бинарного дерева с использованием параллельного слияния по матрицам сравнений. Цель такой постановки - совместить свойства и назначение схем поиска на основе дерева и на основе слияния для получения единого параллельного алгоритма сортировки. Требуется, помимо построения, представить способ адресации по ссылкам к входным элементам, а также дать оценку временной сложности полученной сортировки.

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

А = {7,1,-2,-7,0,-3,5,17,8,-1,9,-§11,13,-5,6}.