УДК 681.51.015
БЫСТРЫЙ АЛГОРИТМ ОПРЕДЕЛЕНИЯ ОРДИНАТ ИМПУЛЬСНОЙ ПЕРЕХОДНОЙ ФУНКЦИИ ПРИ ВОЗБУЖДЕНИИ ДИНАМИЧЕСКОГО ОБЪЕКТА ТЕСТ-СИГНАЛОМ НА ОСНОВЕ ДВОИЧНОЙ М-ПОСЛЕДОВАТЕЛЬНОСТИ
© 2012 В.Ф. Яковлев
Самарский государственный технический университет
Поступила в редакцию 14.12.2012
Рассмотрен алгоритм вычисления ординат импульсной переходной функции линейного динамического объекта при его возбуждении тест-сигналом на основе двоичной М-последовательности с помощью быстрого преобразования Уолша-Адамара.
Ключевые слова: импульсная переходная функция, тест-сигнал, двоичная М-последовательность, бы-строе преобразование Уолша.
Для идентификации линейных динамических объектов используются их различные модели, в том числе и импульсная переходная функция (ИПФ). В этом случае объект может быть описан следующим уравнением:
N -1
y[i] = h0 + At ■ 2 h[j] ■ x[i - j]. (1)
j =0
Здесь y[i] - реакция объекта, At - шаг дискретизации, x[i] - входной тест-сигнал, h[j] - ординаты импульсной переходной функции, h0 -реакция объекта на постоянный рабочий сигнал, N - время памяти объекта [1].
На практике для независимой оценки ординат ИПФ при идентификации применяют ортогональные к сдвигу кусочно-постоянные тест-сигналы небольшой амплитуды, не нарушающие нормальное функционирование объекта, тогда при измерении реакции объекта один раз на такте тест-сигнала:
p
2 y[i] ■x[i - j] h[j ] = ^-. (2)
2 x2[i - j]
i=0
В качестве тест-сигнала удобно использовать легко реализуемые псевдослучайные двоичные М-последовательности [1, 2]. Определение взаимокорреляционной функции для различных j по (2) требует большого объема вычислений.
Известен алгоритм, уменьшающий объем вычислений для определения ординат ИПФ при возбуждении динамического объекта тест-сигналом в виде двоичной М-последовательности, на
Яковлев Вадим Фридрихович, кандидат технических наук, доцент кафедры "Теоретическая и общая электротехника". E-mail: vf [email protected]
основе быстрого преобразования Адамара [3, 4, 5]. Целью этой работы является модернизация быстрого алгоритма для вычисления ординат ИПФ.
Двоичная М-последовательность является упорядоченным с помощью сопровождающей матрицы (характеристического полинома Р(х)), множеством компонент Б. вектора координат элементов поля Галуа СР(2п) в степенном базисе [6]. М-последовательность, генерируемая с помощью характеристического полинома степени п, имеет период Р = (2П - 1) тактов.
При генерации тест-сигналов эти компоненты заменяются реальными сигналами с нормированными значениями:
х(0) = +1; х(1) = -1. (3)
Умножение для реальных сигналов оказывается эквивалентным сложению по модулю 2 для компонент Б.:
1Ф 0 = 1 ^ (-1) • (+1) = -1 1Ф1 = 0 ^ (-1) • (-1) = +1 0 Ф 0 = 0 ^ (+1) • (+1) = +1 (4) 0 Ф1 = 1 ^ (+1) • (-1) = -1
В тест-сигнал на основе двоичной М-последовательности вводятся дополнительные такты для получения несмещенной оценки ЬЦ] ("нулевая строка"), а также для устранения погрешности от неверного задания исходного состояния объекта перед началом тестирования [1, 2].
Для упрощения изложения алгоритм вычисления ординат ИПФ далее рассматривается для небольшого количества Ь[]]. На рис. 1 приведен тест-сигнал на основе двоичной М-последова-тельности с характеристическим полиномом Р(х) = х3+х+1 и периодом 7 тактов, точками указаны моменты измерения реакции объекта.
На практике используются М-последова-тельности с характеристическими полиномами степенью не менее 8.
x(t)
У[1] * У[7]
у[0]
Рис. 1. Тест-сигнал на базе двоичной
М-последовательности. а - формирование начальных условий, б - период М-последовательности, формирование "нулевой строки".
В соответствии с (2):
(h > "0 ' + + + + + + + +^ (
h[0] +---+ - + + y[1]
h[1] = 1 + +---+ - + y[2]
h[2] + + +---+ - y[3]
h[3] = 8 + - + +---+ y[4]
h[4] + + - + +--- y[5]
h[5] +-+-++-- y[6]
1 h[6]j 1 +-- + - + + -у 1 y[7]J
(5)
H =
P +1
-■ X • Y
(6)
Здесь H и Y - векторы-столбцы оценок ординат ИПФ и реакции объекта, X - задержанные реплики входного сигнала, матрица плана эксперимента. С целью упрощения записи в матрице Х приведены только знаки нормированных значений тест-сигнала, дополнительная "нулевая строка" размещена в левом столбце Х.
В аппаратном или программном генераторе непосредственно доступны М-последователь-ность и ее реплики с запаздыванием от 1 до (n-1) такта. М-последовательности с запаздыванием j > n являются линейными комбинациями последовательностей с запаздываниями меньше n того же характеристического полинома:
Si_j = «1 ■ Si ©a2 ■ S_1 ©■......Фa„ ■ S_„+1. (7)
Для тест-сигналов выражение (7) записывается так:
x[i - j] = (x[i])a ■ (x[i -1])«2 ■ ■(x[i - n + 1])a" (8)
Коэффициенты а в (7, 8) совпадают с коэффициентами полинома-остатка R(x) в GF(2n) [6]:
R (x) = x} mod F(x) . (9)
Выражение (9) иногда называют алгоритмом Дэвиса, его применяют для генерации задержанных М-последовательностей при идентификации
линейных динамических объектов. Например, в пятой строке матрицы Х приведен тест-сигнал с запаздыванием на 4 такта по отношению к входному сигналу генератора М-последовательности. Вычислим коэффициенты полинома-остатка по (9) для этого случая:
x4 mod( x3 + x +1) = (x2 + x), т.к.
R4
x
x3 + x +1
x+
x + x x3 + x + 1
(10)
Из (10) следует, что а1 = 1, а2 = 1, а3 = 0 и х[/ - 3] = х[/] • х[/ -1].
Двоичная запись матрицы Х2 имеет вид:
(00000000^ 01110100 00111010 00011101 01001110 . (11) 00100111 01010011 01101001
X 2
В (12) представлены матрица Уолша-Адама-ра [7] W и ее двоичный аналог W2 размерностью 8 ■ 8.
W
( + + + + + + + +^ +-+-+-+-++--++--+-+--+-+
++++----
+--+-++-+--++-- + ++----++
W2
(00000000^ 01010101 00110011 01011010 00001111 01101001 01100110 00111100
. (12)
При умножении вектора-столбца на матрицу Уолша-Адамара используют эффективное быстрое преобразование Уолша-Адамара (БПУ), значительно уменьшающее объем вычислений [7].
Согласно (7) любая строка матрицы Х2 является линейной комбинацией по модулю 2 п строк, соответствующих последовательной смене п-разрядных двоичных чисел в генераторе М-последовательности. В Х2 это вторая, третья и четвертая строки. Для тест-сигналов в Х этот алгоритм трансформируется в (8).
В матрице Уолша-Адамара размерностью 2" • 2" строки также определяются линейной комбинацией по модулю 2 п строк, соответствующих последовательному возрастанию двоичных п-разрядных чисел, образованными пересечением столбца и этими строками. Верхний (второй) элемент столбца соответствует младшему разряду.
1
В [3, 4] предложена методика применения БПУ для вычисления оценок ординат ИПФ при возбуждении динамического объекта тест-сигналом на основе М-последовательности.
В матрице Х2 столбцы упорядочены в порядке появления двоичных п-разрядных чисел в генераторе М-последовательности. После перестановки столбцов в порядке возрастания двоичных п-разрядных чисел матрица Х2 совпадет с (соответственно Х с Ш). Чтобы результат умножения матрицы на вектор-столбец не изменился, следует переставлять и компоненты у[1] вектора У и спектральные коэффициенты Уолша, порядок следования которых отличается от расположения элементов Ь^] в векторе Н. Таким образом:
н = а2 ■ ж ■ а1 ■ у. (13)
Здесь А1 и А2 - матрицы перестановок, определяемые характеристическим полиномом М-последовательности.
Для перестановки элементов в векторе У достаточно записывать замеры реакции объекта у[1] в массив не в естественном порядке, а по адресам, задаваемым двоичными числами в генераторе М-последовательности. Например, из (5) и (11) следует, что в векторе У элементы должны быть переписаны в следующем порядке: у[0], У[1], у[6], у[2], у[7], у[5], у[4], у[3], матрица А, тогда имеет вид:
(10000000^ 01000000 00000010 00100000 00000001 00000100 00001000 00010000
а =
(14)
У
ве факторизации матрицы преобразования Уол-ша-Адамара на произведение п одинаковых матриц. В отличие от алгоритма типа "бабочка" [7] в этом случае в каждом из п этапов БПУ выполняется одна и та же процедура, что упрощает программную реализацию алгоритма.
На рис. 2 алгоритм БПУ изображен в виде графа для рассматриваемого примера. Спектральные коэффициенты, вычисленные с помощью БПУ упорядочены по возрастанию номеров функций Уолша.
В [3, 4] показано, что адреса в двоичном формате, по которым расположены оценки ординат ИПФ ЬЦ] в векторе спектральных коэффициентов Уолша, могут быть определены из матрицы К размерностью п ■ 2п. В матрицу К отбираются п столбцов из Х2, так чтобы верхние строки матрицы К, начиная со второй, образовали единичную матрицу размерностью п ■ п . Младший бит адреса располагается в левом столбце.
Для случая Р(х) = х3+х+1 матрица К и матрица перестановок А2 на ее основе имеют вид:
К =
(000> (10000000>
100 01000000
010 00100000
001 110 , А = 00001000 00010000
011 00000010
111 00000001
I101) 00000100
(15)
Разумно применить алгоритм БПУ на осно-
Соответствие спектральных коэффициентов Уолша и ординат ИПФ указано и на рис. 2.
Недостатком рассмотренного быстрого алгоритма является необходимость генерации в процессе вычислений всей матрицы Х2 размерностью 2п ■ 2п, например, при п = 8 это уже 65536 элементов.
Рис. 2. Алгоритм БПУ
В этой статье предлагается другой способ сопоставления ординат ИПФ ЬЦ] со спектральными коэффициентами Уолша в БПУ.
Спектральные коэффициенты, полученные в результате применения БПУ (Рис 2), расположены по возрастанию номеров функций Уолша при их естественном упорядочивании (по Ада-мару) [7].
Полная ортонормированная система 2П функций Уолша (Рис.3) определяется на основе п функций Радамахера К(1) [7]. Функции Уолша определяются через произведение функций Радамахера. Наличие или отсутствие ^й функции Радамахера при получении 1-й функции Уолша определяется ;|-м разрядом двоичного представления числа 1. Например, = ,№001(1) = И^), = *гш(1;) = Я^) • Я3() и т.д.
Матрицы Х (5) и Ш (12) отличаются только перестановкой столбцов, следовательно строка х[1] в (5) соответствует И1(1) в (12), строка х[1] • х[1 -1] в (5) соответствует Я^^) • Я2(г) = №оц(1) = №з(1) в (12) и т.д.
Таким образом, двоичная запись номера спектрального коэффициента Уолша с1 указывает сигналы с каких разрядов генератора М-пос-ледовательности использовались при его вычислении. Адреса, по которым находятся оценки ординат ИПФ ЬЦ] в векторе спектральных коэффициентов Уолша определяются так: Ь0 = с0 и ЬШ = ск , к = при ^п, к определяется по (9) при ] > п. Необходимости генерировать всю матрицу Х при таком алгоритме не возникает.
Процесс идентификации связан с автоматическим управлением объектом исследования, оборудованием для подачи тест-сигнала и сбора данных, обработкой данных. Это удобно осуществлять в специализированной среде программирования LabVIEW [8]. Разумно и подготовительную работу делать в той же среде. Автором был разработан несложный виртуальный прибор для определения ординат ИПФ ЬЦ] линейного динамического объекта при его возбуждении тест-
JE Генератор М-гюсл.vi Front Panel *
File Edit View Project Operate lools Window Help
Ro(t), Wooo(t) Ri(t), Wooi(t) R2(t), Woio(t) Won(t) R3(t), wioo(t)
Wioi(t) Wno(t) Win(t)
"1 п г 1 п г.
и и 1 г J и ~~1 г
1_1 1_1 г.
~1 п 1 г
U L J и
=Kf d=tp~
Рис. 3. Функции Уолша, упорядоченные по Адамару
сигналом на основе двоичной М-последователь-ности. На рис. 4 приведен фрагмент лицевой панели виртуального прибора, на рис. 5 часть блок-схемы.
В разработанном виртуальном приборе динамический объект моделировался набором ординат ИПФ (100, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 и 12). На рисунке 4 представлено состояние лицевой панели виртуального прибора после проведения идентификации. Использовалась М-пос-ледовательность с характеристическим полиномом х8 + х6 + х5 + х4 + 1, который задан в двоичной форме. Ординаты ИПФ полностью восстановлены, указаны их адреса в массиве спектральных коэффициентов Уолша.
ВЫВОДЫ
При вычислении оценок ординат ИПФ линейных динамических объектов при возбужде-
I 15pt Application Font - 11 ||±i-11 I
Характеристический полином М-по£лцдов«тел»пости
Степень полинома
[шшооГ
Ординаты импульсной переходной функции
kJlioo и "Г H ьп3 F- 10 111 |l2 |o lô
|o I< И.................... ............>|
Адреса ординат импульсной переходной функции в спектре Уолша
2 1 3 16 32 64 123 113 226 131 27 5-1 103 2
Ш . |гГГГГГГГГГГf
1 1-1 ц и 1-11-1и и
TT
Рис. 4. Лицевая панель виртуального прибора
Рис. 5. Часть блок-схемы виртуального прибора
нии их тест-сигналом на основе двоичных М-последовательностей используется эффективный алгоритм быстрого преобразования Уолша-Адамара. При этом необходимо в зависимости от вида характеристического полинома М-последо-вати произвести перестановки компонентов в векторах реакции объекта и спектральных коэффициентов Уолша. В статье предложен более простой алгоритм определения соответствия между оценками ординат ИПФ и спектральными коэффициентами, полученными в результате БПУ, не требующий генерации всего плана эксперимента размерностью 2" • 2".
СПИСОК ЛИТЕРАТУРЫ
1. Ikonen E. Advanced process identification and control.
New York: Marcel Dekker Inc., 2002. 316 p.
2. Яковлев В.Ф. Выбор характеристического полино-
ма двоичной М-последовательности для идентификации нелинейного динамического объекта // Известия Самарского научного центра РАН. 2011. Т.13. №4. С.133-135.
3. Cohn M, Lempel A. On fast M-sequences transforms // IEEE Trans Inform. Theory. - 1977. IT - 23. C.135-137.
4. Jens Hee Impulse response measurements using MLS // Bruel & Kj®r, Denmark. - 2003. 16 pp. URL: http:// jenshee.dk (дата обращения 18.11.2011).
5. Perrett M. Implementation of a M-sequence pseudo random binary sequence audio measurement system based on the Hadamard transform // University College London. -2010. 4 pp. URL: http://www.ee.ucl.ac.uk/lcs/ previous/LCS2010/lens2010_submission_25.pdf (дата обращения 18.11.2011).
6. Davies W.D.T. System identification for self-adaptive control. New York: Wiley-Interscience, 1970. 290 р.
7. ЗалманзонЛ.А. Преобразования Фурье, Уолша, Хаа-ра и их применение в управлении, связи и других областях. М.: Наука, 1989. 496 с.
8. Тревис Дж. LabVIEW для всех. М.: ДМК Пресс, 2005. 540 с.
FAST ALGORITHM FOR COMPUTING THE IMPULSE RESPONSE OF A LINEAR DYNAMIC OBJECT USING BINARY M-SEQUENCE AS INPUT
© 2012 V.F. Yakovlev
Samara State Technical University
The paper considers algorithm of computing the impulse response of a linear dynamic object using binary M-sequence as input via fast Walsh-Hadamard transform.
Keywords: impulse response, test-signal, binary M-sequence, fast Walsh-Hadamard transform.
Vadim Yakovlev, Candidate of Technics, Associate Professor at the of Theoretic Electro Technology Department. E-mail: [email protected]