CC BY
38
УДК 519.233.2
Вестник СПбГУ. Сер. 1. Т. 1 (59). 2014. Вып. 1
ТОЧНОСТЬ ОЦЕНИВАНИЯ
ПЕРЕДАТОЧНОЙ ФУНКЦИИ ЛИНЕЙНОГО ФИЛЬТРА
А. Ю. Груша
Санкт-Петербургский государственный университет,
Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7/9
Рассматривается задача определения параметров линейной регрессии. Неизвестная им-пульная характеристика оценивается по методу наименьших квадратов и методу максимального правдоподобия. Для обоих методов приводятся способы быстрого расчета точности оценок импульсной характеристики и передаточной функции с помощью многочленов Сегё и Шура. Определяется зависимость точности оценивания от значений параметра регуляризации. Биб-лиогр. 8 назв. Ил. 3.
Ключевые слова: идентификация параметров, точность оценивания, вычислительная сложность, сверхбыстрый алгоритм Шура, многочлены Сегё, формула Гохбегра—Семенцула.
1. Введение. Рассматривается задача идентификации параметров линейной регрессии. Предполагается, что известен вектор х = (хь)^=1 и наблюдается случайный вектор у = (уг)Т=н, для которого
N-1
Уь = Е хкхь-к + V, N < * < Т, (1)
к=0
где V = (уь)т=1 —белый шум с дисперсией
Вектор коэффициентов х0 = (хк )й=-/ оценивается по методу наименьших квадратов (МНК),
т
J (х) = е
t=N
N-1
yt - Е Xkxt-k
^ mm.
х
k=0
Известно, что оптимальная оценка может быть записана в виде
тт х = R-1Y, R Y =E$tyt, (2)
t=N t=N
где для вектора-столбца введено обозначение Фt = (xt-k)д=)1-
Ставится задача быстрого расчета оценки импульсной характеристики объекта X, соответствующей оценки значений передаточной функции, а также точностей этих оценок.
При N ^ ж и T > N сложность решения системы линейных уравнений RX = Y в общем виде имеет порядок O(N3). Для теплицевых систем алгоритм Левинсона— Дурбина [1] позволяет понизить количество операций до O(N2). При помощи сверхбыстрого алгоритма Шура [2] для поиска многочленов Сегё [3] и формулы Гохберга— Семенцула [4] сложность понижается до O(N log2 N). В данной статье представлен алгоритм расчета оценки точности передаточной функции по многочленам Сегё, который имеет сложность порядка O(N log N).
Как известно, в предположении о гауссовском распределении шумов и вектора X0 метод максимума правдоподобия (ММП) сводится к МНК с дополнительным параметром регуляризации. В приложениях точность оценки передаточной функции
2
существенно зависит от этого параметра. В разделе 4 представлен алгоритм расчета точности оценивания передаточной функции по регуляризованному МНК, имеющий сложность O(N log N).
2. Сверхбыстрый алгоритм решения теплицевых систем. Предположим, что T ^ N и что сигнал x начинается и заканчивается паузами: xt =0 при 1 < t < N — 1 и при T — N +1 < t < T. Тогда информационная матрица R = (Rij )fj=0 является теплицевой:
T-N
Rij ri-j ^ ^ xtxt-|i-j| t=N
0 < i,j < N — 1,
где г = (гк) —корреляционная функция последовательности х.
Многочленом Сегё [3] степени N — 1, соответствующим корреляционной функции г, называется многочлен
N1
-(z) = 1+^ фкzk
к= 1
с корнями вне единичного круга, для которого разложение в ряд Лорана функции
S(z)
d
-(z)-(z-1) ^
Ж
^ г к zk
при некотором числе d > 0 имеет те же N начальных коэффициентов: г к = Гк при 0 < к < N — 1.
Такой многочлен -(z) и число d существуют, единственны и являются решением уравнения Юла—Уокера [5, 6]. Для расчета коэффициентов фк обычно применяется алгоритм Левинсона—Дурбина [1], имеющий сложность O(N2). В задачах акустической эхокомпенсации это слишком много для системы реального времени. Г. Амма-ром [2] был разработан сверхбыстрый алгоритм расчета вектора (фк)£=-/ сложности O(N log2 N) на основе теории многочленов Шура [7, 8]. Далее будем предполагать, что по вектору (гк)N=-1 построен многочлен -(z) и число d, называемое остаточной дисперсией.
При заданном многочлене Сегё обращение теплицевой самосопряженной матрицы с первым столбцом (гк)N—i1 осуществляется при помощи формулы Гохберга—Се-менцула [4]
R
1
= -(MiMi - M0M0J),
где
M1
1
ф1
V-N-1 -N-2
0
0 1
M0
0
-N -
0 0
ф1 ф2
-N-1 0/
Вычисление R-1Y сводится к умножению на треугольные теплицевы матрицы Mo, M1 и транспонированные, что выполняется при помощи ДПФ. Вычисления содержат всего 7 ДПФ на N точек. Их сложность пропорциональна O(N log N). Таким образом, расчет оценки х по МНК содержит следующие операции:
1
• расчет корреляционной функции r = (rk)^=01 • Сложность O(log N) операций на отсчет;
• сверхбыстрый алгоритм Шура [2] по вектору r вычисляет пару многочленов Шура и многочлен Сегё. Сложность O(N log2 N);
• реализация формулы Гохберга—Семенцула в частотной области и расчет оценки МНК £ Сложность O(N log N).
3. Точность оценки передаточной функции. Оцениваемый вектор х0 можно рассматривать как импульсную характеристику линейного блока (1). Тогда х = (Xfc )fc=o1 —оценка импульсной характеристики, и соотвутствующую передаточную функцию и ее оценку обозначим
N-1 N-1
W0(z)= £ xkz-k, W(z)= £ Xkz-k. fc=0 fc=0
В следующем утверждении представлен способ расчета точности оценок передаточной функции в N равноотстоящих значений, требующий O(N log2 N) вычислительных операций.
Теорема 1. Пусть ^(z) —многочлен Сегё порядка N — 1, построенный по корреляционной функции r = (rfc )N=o1 последовательности x, и d — остаточная дисперсия. Они определяют спектральную плотность S(z).
Представим комплексное значение ^(z) в показательной форме:
Ф(^ф) = р(ФКа(ф), 0 < ф < п, р(ф) > 0, (3)
причем функция а(ф) непрерывна на [0, п]. Тогда для любого ф G [0, п]
¥
- ' = ( 1 - ^Р- ) . (4)
2
2а'(ф)
S (в*ф) V N
Доказательство. Пусть Дх = х — х0 — погрешность оценки импульсной характеристики. Матрица ковариаций оценок импульсной характеристики есть
21
е(Дх)(ДхГ = <R
где R — теплицева информационная матрица МНК с первым столбцом r = (rfc )^=-)1-
z-V+1
Для любого комплексного числа z определим строку q(z) = (l,z
).
Тогда передаточная функция и ее оценка есть W 0(z ) = q(z )х, W (z) = q(z)X Поэтому
1
е
W 0(е4ф ) — W (е4ф )
2 a,2
Здесь и далее символ * обозначает транспонирование и комплексное сопряжение.
По формуле Сегё [3] теплицева матрица R и соответствующий многочлен Сегё ф(г) связаны следующим условием: для любых комплексных чисел А и ц на единичной окружности
q(A)R-1 q(M)*
l
d(1 — A-V)
^(A-1 ЖМ) — (A-1M)V ^(A)^(M)
1
Рассмотрим предельный случай. Пусть ц = г и Л ^ г. По правилу Лопиталя
д(г)Е-1д(г)* =
Л-г
Л ф(\-1)ф(г) - {\-1г)мф{\)ф{г)*
Л—d
Ит —
Л—г а
1
N
и
ф(Л-1 )ф(г) - (Л-1 г)*ф(Л)ф(г)* - 2Re[zф'(z)ф(г)*] + N|ф(г)|2 |ф(г )|2
Л=г
1
Пусть 2 = е1ф и ф(г) = р(ф)е1а(ф). Тогда
ф(е<ф)ф = ггф'(г), ^(ф(егф )))ф =
¿0)))/ =
. гф'(г)
у -
и поэтому
Ие
гф'(г)
ф(г)
что завершает доказательство.
= 1т (^(ф(е^)))ф = 1т ((^(р(ф)))ф + га'(ф)) = а '(ф),
(5)
□
4. Регуляризация. В задачах эхокомпенсации, рассмотренных в разделе 5, расчет передаточной функции по МНК дает слишком большие погрешности. Если спектр входного сигнала мал в определенной полосе частот, Б(егф) « 0 при ф £ [ф-, ф+], то из теоремы 1 следует, что точность оценки очень плохая. Этого удается избежать при помощи регуляризации, к которой сводится метод максимума правдоподобия.
Предположим, что оцениваемый вектор х0 является случайным и гауссовским со средним 0 и матрицей ковариаций аXI*, где I* — единичная матрица. Пусть последовательность (^) гауссовская, с независимыми значениями, имеющими распределение
N (0,а2).
Метод максимума правдоподобия сводится к минимизации функционала
л (х) = j (х) + 7|х|2 = е
г=*
N-1
Уг - Хкхг-к
к=0
+ 7 |х|2
при выборе параметра 7 = а^/аХ = 70.
Предположим, что для расчетов выбран параметр 7, который может отличаться от 70 = а"22/аХ. Минимум функционала JY достигается при
Х7
^ФгФ* + 71*, У ^^ФгУг.
г=N
г=N
Матрица Д7 теплицева. Ее первый столбец совпадает с эмпирической корреляционной функцией г = (гк)^=о1 последовательности х, кроме первого элемента. Первый элемент, стоящий на главной диагонали матрицы Д7, содержит дополнительное слагаемое 7.
Л
2
Таким образом, ММП сводится к той же вычислительной задаче, что и в разделе 1, но соответствующая спектральная плотность Б7(г) отделена от нуля. Параметр 7 регуляризирующий, но точность оценок чувствительна к нему. Обычно отсутствуют надежные данные о дисперсиях а2 и аХ. Поэтому необходимо оценить зависимость точности оценок передаточной функции от выбора числа 7.
Следующие два утверждения определяют значение ковариации ошибки оценивания импульсной характеристики и быстрый способ расчета ковариации ошибки оценивания передаточной функции, сложность которого не превосходит log(N)).
Лемма 1. Пусть х7 — точка минимума функционала JY и Дх7 = Х7 — Х° — ошибка оценивания импульсной характеристики. Тогда матрица ковариаций ошибки оценивания есть
ЕДх7 (Дх7 )т = а2й-1 + вй-2, (6)
где
в = 7(7ах - а
2
Доказательство. Запишем все уравнения измерения в векторной форме у = Фх0 + V. Тогда по определению оценки х7 ,
Дх7 = (Д + )-1(Дх0 + Ф«) - Х0 = Д-1Ф« - 7Д-1Х°.
Из независимости V и х0 следует, что
еДх7(Дх7)т = а2Д-1ДД-1 + 72а£Д-2 = а2Д-1 + (72аХ - 7а2)Д-2, (7)
что совпадает с утверждением леммы. □
Теорема 2. Пусть 7 > 0 — произвольный регуляризирующий параметр. По теп-лицевой матрице Д7 построим многочлен Сёге ф(г) и остаточную дисперсию а. Определим спектральную плотность Б (г) = а/1 ф (г) |2.
Тогда дисперсия ошибки оценивания, передаточной функции в точке г = егф равна
а2 в е|ДИ^)|2 = ^.(М-2а'(ф)) + ^(В(г)-С(г)), (8)
£
где а(ф) — непрерывная ветвь комплексного аргумента функции ф(егф) и
(1 2*-2 \ 1 N-1
^ , од = щ^ур Е *<*(** + о,
где вещественные коэффициенты (Ьк)к=сТ2 и (ск)д=/ определяются из разложения
2N-2 N-1
ф(г)2 = £ Ькгк, ф(г)ф(г-1) = ]Т С|к|гк.
к=0 k=-N+1
Доказательство. Определим строку д(г) = (1, г-1,..., z-N+1). По аналогии с доказательством теоремы 1 из леммы 1 следует, что
Е|Д^ (г)|2 = а2^(г )Д-19(г )* + в|«(г )Д-1|2.
Первое слагаемое вычислено в теореме 1. Второе слагаемое без множителя в представим в виде
(9)
где ¿ш — нормированная мера Лебега на единичной окружности. Это верно, поскольку из определения д(^) следует, что
/ д(м)*д(м) ¿ш(м) = /« «'1м1 = 1
— единичная матрица.
По формуле Сегё для любых комплексных чисел А и м
д(А)Д-1д(м)* а(1 - Ам*) = ^(АЖм)* - Ам*<Д(АЖм)*,
где ^ — многочлен в обратных степенях, = -1).
Отсюда подынтегральное выражение в (9), в котором |м| = 1, А — произвольное комплексное число, не равное м, можно записать в виде
к(А)Д-1 ?(м)*|2
^(АЖм)* - Ам*<МАЖм)*
¿(1 - Ам*)
^(АЖм-1) - А«м-«^(А-1Жм)
¿(1 - Ам-1)
№(м)|
12
а211 - Ам-112
^(м)2м-«
—2 И.е
а211 - Ам-1|2
А« ^(А)>(А-1)
Оба слагаемых как функции от А имеют полюсы на единичной окружности. Поэтому выберем А = г^ при | = 1 и 0 < г < 1, вычислим интегралы от каждого слагаемого и затем вычислим предел суммы при г ^ 1 слева.
Воспользуемся следующим тождеством. Пусть |А| < 1. Тогда по теореме о вычетах
мк
=1 |1 - Ам-1|2
мк
=1 (1 - А*м)(м-А)
Ак (1 -1А|2) 1, к > 0,
(А*)|к|(1 -1А|2) 1, к < 0.
Подставим разложения в квазимногочлены функций ^(м)2 и |^(м)|2 = ^(м)^(м 1), введенные в формулировке теоремы. Получим, что при | А| < 1
|д(А)д-1д(м)* |2 аш(м) =
1
■У|м| = 1
х(|^(А)|2 + |А«^(А-1)|2) - 2Ие 28
а2 (1 -1А |2)
т-1
1
]Г ейАк +Е (А*)к
к=1
т-1
2«—2
][>(А*)|к-№| + Е ЬкА
к-« 1
А« ^(А)>(А-1)
2
2
2
Пусть Л = гг, где |г | = 1. Выражение в правой части есть отношение двух бесконечно малых при г ^ 1. Действительно, первое выражение в фигурных скобках при Л = г равно 2|-0(г)|4, а выражение в квадратных скобках равно г-N^(г)2 • г^(г-1).
Вычислим предел по правилу Лопиталя. Дифференцирование по г в точке г = 1 при |г | = 1 дает
^(rz)^(rz 1) + r2N^(r 1 z
„2N
-1.
11
- r r = 1
(rz)N ^(rz-1)^(r-1 z-1)
2N |^(z)
Nz N V(z-1)2.
r=1
Эти слагаемые в общей сумме сокращаются и остается
Г -1 Г N-1
/ ЫА^УмЛ^тЫ = — |VW|2 fecfc(zfc +z-fc) -
•/|м| = 1 d I fc=1
2N-2 -| s
^(z-1)2 ^ |k—N|6fczk
U — П -I S
-2 Re
Отсюда заключение теоремы получается подстановкой S(z) = d/|^(z)
□
2 N 1 n) i2 ) n —n вычисля-
Замечание 1 (о сложности расчета точности оценки передаточной функции). Поскольку длина импульсной характеристики равна N, передаточная функция полностью определяется по набору своих значений в N точках. Этими точками принято выбирать zn = exp(2nm/N), 0 < n < N — 1. Значения (W(zn))N—П1 вычисляются при помощи одного БПФ длины N от оценок (xk)N—П1, и поэтому сложность составляет O(N log N).
Точность оценок определяется теоремой 2. Значения e|AW(zn)|2)n=n ются по значениям нескольких многочленов, что также выполняется при помощи БПФ. Поэтому расчет точности оценок передаточной функции также имеет сложность O(N log N).
Замечание 2 (о выборе параметра регуляризации). Оценка импульсной характеристики х = xY существенно зависит от параметра регуляризации y. Оптимальный выбор y = ст^/^Х требует каких-либо оценок заранее неизвестных величин и о"Х. Эти величины можно оценить по усредненным значениям нормы xY и показателя качества
D.
y -
t=N
|yt — 12
2
r
2
Лемма 2. Пусть 7 > 0 — произвольный параметр и по нему построена оценка Х7. Тогда
ЕБ1 = а2(Т - 2Ж +1 + 7Ш2) + ^72(«1 - Ш2), е|х7|2 = - 7ш2) + (Ж - 27Ш1 + 72Ш2),
где
т1 = tr й-1, т2 = tr Д-2.
Доказательство. Поскольку у = Ф*х° + V, Д7 = ФФ* + , случайные векторы
V и х°независимы и х7 = Д- 1Фу
Е|х7 |2 = а^2 ^(Д—1ФФ*Д—1) + аХ МД—1ФФ*ФФ*Д—1) =
= а2 МД—1(Д7 - 7/*)Д—1] + аХ МД—1(Д7 - 7/*)(Д7 - 7/*)Д-1] =
= а^(Ш1 - 7Ш2) + аХ(Я - 27Ш1 + 72Ш2),
что соответствует второму уравнению в лемме.
По определению минимальное значение функционала качества ^у равно + 7 |х712. Вычислим
е7ътш = е[|у|2 - у*Ф*Д—1Фу] =
= а2 [Т - N +1 - ^ Д—1(Д7 - 7/*)] +
+ аХ tr[(Д7 - 7/*) - (Д7 - 7/*)Д—1(Д7 - 7/*)] =
= а^(Т - 2Я +1+ 7Ш1)+ аХ(7Ж - 72Ш1).
Вычитанием получаем утверждение леммы. □
В практических расчетах в уравнениях леммы можно заменить е|х712 на сумму квадратов оцененных коэффициентов х7, а е|Д71 —на сумму квадратов невязок. Величины Ш1 и Ш2 определяются из уравнений
1 N —1 1 N —1
= — Е ч(^п)Щ1д(гп)*, то2 = — ^ ч(^„)Щ2д(гп)
п=° п=°
вычисленных в теоремах 1 и 2, где ) = (1, г71,..., г7*+1). Затем решается система
22
двух линейных уравнений для определения оценок а^ и аХ.
5. Эхокомпенсация. Примером применения данной теории на практике может служить задача компенсации акустического эха. Уравнение (1) может быть интерпретировано как линейный блок, где х — входной сигнал, поступающий из сети, влияние которого нужно подавить; у — сигнал, записываемый микрофоном, содержащий как полезный сигнал, так и акустическое эхо, которое может быть представлено в виде
ХиХ(—и; х° —импульсная характеристика эхо-канала, описывающая акустические свойства помещения, в котором распространяется эхо, и оборудования, которе воспроизводит и записывает сигнал.
Информационная матрица Д, рассчитанная по входному сигналу х, почти всегда оказывается плохо обусловленной, поэтому нахождение обратной матрицы зачастую вызывает определенные трудности и приводит к значительным погрешностям. Регу-ляризирующий параметр 7 призван бороться с плохой обусловленностью матрицы Д и, тем самым, повышать точность оценки. На рис. 1 изображена передаточная функция эхо-канала Ш°(егф), подлежащая оцениванию. На рис. 2 приведено две ее оценки. Синим цветом изображена оценка, сделанная без учета регуляризирующего параметра. Как можно видеть, в области высоких частот погрешность оценивания становится значительной. Зеленым цветом изображена оценка, проведенная с учетом регуляризирующего параметра. Она лишена подобного недостатка и гораздо точнее оценивает идеальную передаточную функцию Ш°(егф).
4000 6000 8000 10000 12000
Рис. 1. Передаточная функция эхо-канала.
Рис. 2. Оценки передаточной функции: синим цветом изображена оценка, полученная без регуляризации, зеленым — с регуляризирующим параметром 7 = 0.01 • го.
Рис. 3. Среднеквадратичное отклонение передаточной функции: синим цветом изображено значение отклонения при 7 = 0, зеленым — при 7 = /а^ , красным — при 7 = 2а2 /а2.
Значения регуляризирующего параметра 7 тоже непосредственно влияют на точность оценивания передаточной функции. На рис. 3 показано среднеквадратичное отклонение передаточной функции из теоремы 2.
В первом эксперименте параметр 7 был принят равным нулю, что соответствует самой большой погрешности на графике (синий цвет). Во втором эксперименте использовалось оптимальное расчетное значение 7; полученное отклонениие изображено на графике зеленым цветом. В третьем эксперименте намерено было использовано значение в два раза превосходящее расчетный оптимум, чтобы оценить размер погрешности при допущении такой ошибки. Полученная погрешность изображена красным цветом на графике. Таким образом показано, что увеличение параметра регуляризации приводит к росту погрешности.
Литература
1. Durbin J. The fitting of time series models // Review of International Statistical Institute. 1960. Vol.28, N3. P. 233-244.
2. Ammar G., Gragg W. Superfast solution of real positive definite Toeplitz systems // SIAM J. Matrix Anal. Appl. 1988. Vol. 9. P. 61-76.
3. Szego G. Orthogonal polynomials. American Mathematical Society. Providence, RI. 1939.
4. Гохберг И. ЦСеменцул А. А. Об обращении конечных теплицевых матриц и их континуальных аналогов // Математические исследования. 1972. Кишинёв. Т. 7. Вып. 2, №24. С. 201-223.
5. Kay S. M. Modern spectral estimation: Theory and application. Englewood Cliffs, NJ: Prentice-Hall, 1988.
6. Marple S. L. Digital spectral analysis with Applications. Englewood cliffs, NJ: Prentice-Hall, 1987.
7. Fulton W. Young tableaux: with applications to representation theory and geometry. Cambridge University Press, 1997.
8. Sagan B. E. The symmetric group: representations, combinatorial algorithms, and symmetric functions. 2nd ed. Springer, 2001.
Статья поступила в редакцию 24 октября 2013 г.
Сведения об авторе: Груша Андрей Юрьевич —аспирант; [email protected]
ACCURACY OF THE ESTIMATION OF
THE TRANSFER FUNCTION OF THE LINEAR FILTER
Andrei Yu. Grusha
St.Petersburg State University, Universitetskaya nab., 7/9, St. Petersburg, 199034, Russian Federation; [email protected]
The identification problem of linear regression parameters is studied. Unknown impulse response is estimated by Least Squares and Maximum Likelihood methods. The main application considered is the acoustic echo cancellation problem with about two thousand parameters to be estimated. The number of measurements is ten times more only and the problem is ill conditioned due to specific speech spectrum. The superfast Toeplitz inversion numerical algorithm is applied for estimation of the echo impulse response and echo transfer function. As a rule, speech signal has a poor spectrum that leads to singularity in matrix inversion. Accuracy of the transfer function estimate is explicitly expressed by initial data. Complexity of calculation of the accuracy is shown to be the same as complexity of FFT for the length of impulse response. Regularization of the matrix inversion problem is studied. An explicit formula for the loss of the cost function after regularization is presented. The underlying mathematical technique is based on theory of Szego and Schur polynomials. Application results of echo cancellation are shown. Refs 8. Figs 3.
Keywords: parameters identification, computational complexity, superfast Schur algorithm, Szego polynomials.
