Айфичер Э., Джервис Б. Цифровая обработка сигналов, практический подход (2-е изд., 2004) (1095888), страница 167
Текст из файла (страница 167)
Очевидно, что обычная схема наименьших квадратов не подходит для оперативной (в реальном времени) оценки, тогда как предпочтителен оперативный алгоритм, включающий ограниченное фиксированное число арифметических операций и не содержащий прямого обращения матриц. При использовании этих оперативных алгоритмов д так обновляется с поступлением каждой новой выборки, что изменения окулярных движений могут отражаться в д, и устранение окулярных артефактов в этом случае можно рассматривать как адаптивную фильтрацию окулярных артефактов с ЭЭГ [38).
Существует два метода оценки параметров, подходящих для оперативной оценки д, — метод наименьших квадратов [38) и рекурсивный метод наименьших квадратов [7, 30]. С точки зрения требований к вычислениям и памяти алгоритм, построенный на схеме наименьших квадратов, более эффективен. Кроме того, для него несущественна проблема численной неустойчивости, присущая алгоритму, построенному на рекурсивной схеме (см.
ниже), В то же время рекурсивная схема наименьших квадратов дает лучшую сходимость и по этой причине предпочтительнее. Обычно используется следующая форма схемы наименьших квадратов: д(т + 1) = д(т) + 2т»х(лт + 1) [у(т + 1) — хт(т + 1)д(т)), (14.7) Глава 14. Приложения и разобранные примеры где д(т) и 6(лг + 1) — оценки д в пз-й и (пз + 1)-й точке выборки соответственно, уэ — константа, отвечающая за скорость сходимости и устойчивость алгоритма. Для обеспечения сходимости уз должно входить в диапазон 0<л < 1УЛ., где Л вЂ” максимальное собственное значение матрицы (ХтХ ) в формуле (14.6). Впрочем, время сходимости алгоритма прямо пропорционально отношению максимального и минимального собственных значений матрицы (ХтХ ), которое может быть очень большим, если входные переменные близки к коллинеарным, что обычно для ЭОП В то же время схема наименьших квадратов широко используется в биомедицинских приложениях для снижения шума артефактов в основном из-за своей простоты.
Подходящий рекурсивный алгоритм, построенный на схеме наименьших квадратов, получается путем экспоненциального взвешивания данных с целью постепенного устранения влияния на оценку старых данных. Таким образом, .У=~~~ у 'е (1), О< у<1. (14.8) Минимизация Т относительно значений 9 приводит к следующему рекурсивному алгоритму: В(т+ 1) = д(т) + с~у(пз+ 1) — х~(т+ 1)д(пз)] (14,9, а) Р(т+1) = — ~Р(т) — — Р(т)х(т+1)х (го+1)Р(гл), (14.9,б) т у~ а где и = у+ х~(т+ 1)Р(т)х(т+ 1), хт = ]хг(т+ 1) хз(пз+ 1) ...
х„(т+ 1)], С = Р(т+ 1)х(т+ 1) = Р(т)х(гп+ 1)/а. Здесь аргумент пз используется, чтобы подчеркнуть тот факт, что величины получаются в каждой точке выборки, у называется коэффициентам забывания и не дает матрице Р(гп + 1) устремиться к нулю (а д(т + 1) — к константе) с увеличением гл, таким образом становится возможным отслеживание медленно меняющегося параметра.
Обычно э выбирается между 0,98 и 1. При меньших значениях самым свежим данным присваиваются слишком большие весовые коэффициенты, что приводит к сильно флуктуирующим оценкам. Существует, впрочем, две основные проблемы, которые могут возникнуть при непосредственной реализации рекурсивного алгоритма наименьших квадратов.
Первая, именуемая "раздуванием" ("Ыоъ-пр") возникает, когда сигнал не инициируется непрерыв- 14.2. Приложения ЦОС но, например, когда нет окулярных движений, что приводит к экспоненциальному росту элементов Р в формуле (14.9, б): йш [Рп(т+1)[ = йш ~ — "~ — оо. ГРп(т)) (14.10) Р(т+ 1) = 0(т+ 1)0(т+ 1)Ю~(т+ 1), (14.11) где Щт + 1) — единичная верхняя треугольная матрица, П (т + 1) — обратная к ией, Р(т + 1) — диагональная матрица.
Таким образом, вместо обновления матрицы Р обновляются ее факторы П и О. Используя уравнение (14.11), формулу (14.9, б) можно переписать следующим образом: Р(т+1) = — П(т) ~0(т) — — хх ~ ц (т), ,г (14.12) х = 0(т) ц (т) х(т + 1). Если выражение в квадратных скобках факторизовать на верхнюю треугольную и диа- гональную матрицы, получим 6(пх)6(т)В (т) = 0(т) — -ххт, ( 14.13) Впрочем, из-за появления миниатюрных окулярных движений (которые есть всегда) и других сигналов, обычно встречающихся в каналах ЗОГ, данная проблема не настолько серьезна. Вторая проблема — это чувствительность алгоритма к компьютерным ошибкам округления, которые дают отрицательно определенную матрицу Р и в конечном итоге— неустойчивость.
Для успешных расчетов необходимо, чтобы матрица Р была положительно квазиопределенной, что эквивалентно следующему требованию: в автономном алгоритме матрица ХгХ имеет обратную, но из-за наличия в уравнении (14.9, б) разности членов положительную определенность Р гарантировать можно не всегда [4, 7, 30]. Данная проблема острее в многопараметрических моделях, особенно если переменные (в данном случае — электроокулограммы) линейно зависимы [30), н если алгоритм реализован в малой системе с конечной длиной слова [7). Проблему численной неустойчивости можно решить путем подходящей факторизации матрицы Р, чтобы избежать вычитания в формуле (14.9, б). Подобные алгоритмы факторизации лучше подходят для численных расчетов, а их точность сравнима с точностью алгоритмов, построенных на рекурсивной схеме наименьших квадратов, в которых используется двойная точность [4„5).
Два таких алгоритма — алгоритм квадратного корня и алгоритм 1313-факторизации. С точки зрения памяти и сложности вычислений 1ЛЭ-алгоритм эффективнее. Фактически И3-алгоритм — это обычный алгоритм квадратного корня„в котором не используется квадратный корень, а значит, свойства этих алгоритмов аналогичны. В 1Л3-алгоритме матрица Р(т + 1) факторизуется следующим образом: Глава 14. Приложения и разобранные примеры где черта используется, чтобы отличить факторы 11 и 0 матрицы 0(т) — (1/сз)нтт от факторов матрицы Р. Далее получаем р(т+1) = -П(т)Р(т)П(т)В (т)0 (т). (14.14) 7 Сравнивая уравнения (14.11) и (14.14) и заметив, что произведение верхних треугольных матриц также является верхней треугольной матрицей, с помощью симметрии уравнения (14.14) получаем и(т+ 1) = 0(тп)0(т) (14.15, а) 0(т+1) = — 5(т).
1— 7 Таким образом, для обновления 11(т+1) и Р(т+1) необходимы подходяшие рекурсивные формулы для 0 и Р(т). В работе [4) приводится алгоритм рекурсивного обновления Щт+1) и Р(т+1) для фильтра Калмана, в котором используетсядисперсия члена ошибки е(!), но ие 7. Данный алгоритм тривиально приводится к виду, удобному для использования в задаче устранения окулярных артефактов. Данный модифицированный алгоритм, в котором фигурирует величина 7, приведен в приложении 14А. Вектор усиления С, получаемый на этапе 10 указанного алгоритма, используется для обновления оценок параметров, как следует из формулы (14.9, а). Так, хотя Р(т + 1) можно получить, обновляя 1)0-элементы согласно формуле (14.11), непосредственно вычислять Р(т + 1) необязательно.
(14.15, б) 14.2.2.3. Некоторые свойства 00-алгоритма Чтобы лучше понять 1)0-алгоритм, полезно записать его явно. Таким образом, для двухпараметрической модели алгоритм, приведенный в приложении, выглядит тато ° этап 1: о! = х,(т+ 1); оз — — хз(т+ 1) + с!зз(т)хз(т+ 1); ° этап 2: 6, = дз(т)х!(т+ 1); Ьз — — !(з(т)нз, ° этап 3; !з! = 7+ 6 о, = 7+ !1з(т)хзз(т + 1); ° этап 4: !(з(т+ 1) = г(!(т)/аз, ° этап 5: ез = сз! + Ьзоз = о! + !1з(т)оз' ° этап б: Ьгзз(т+ 1) = Ьгзз(т+ 1) — Ь| из~с!!; ° этап 7: Ь, = Ь, + ЬЭУзз(т) = !1з(т)хз(т+ 1) + !)з(т)оз(Гзз(т); ° этап 8: !(з(т + 1) = г(зглзГ7глз.
Из этапа 3 видно, что, если начальные значения диагональных элементов (те. дз(0) и Ы, (0)) положительны„сз,, а следовательно, г1! (пз+ 1) всегда будет больше нуля. То же справедливо для !Эз (а в формуле (!4.9)) и !(з(т + 1). 1. Положительная определенность Р. Матрица Р положительно определена тогда и только тогда, когда х Рх ) О, исключая случай х, = хз — — ... —— х„= 0 (2).
Из формул (14.9) и (14.11) получаем, что !з = 7 + хтРх = 7 + к~11011~к, откуда для двухпараметрической модели вычисляем (см. выше этапы 3 и 5): !."2 ~ "1 + !12оз г + ~~! (т)х! (т + 1) + ~~з(™)о2 14.2. Приложения ЦОС вы так что в этом случае *тР =,( (т)хг(го+ ц+( (ж)ог Следовательно, знак хг Рх зависит от знаков диагональных элементов (т.е. Н, и Иг), юторые, как отмечалось ранее, всегда положительны, поэтому матрица Р является положительно определенной.
Значит, алгоритм ~Л)-факторизации гарантирует положительную определенность Р. 2. Проблема раздувания. При отсутствии данных (т.е. югда хг = кг — — ... — — х„= О) сг, = аг — — т (этапы 3 и 5) и обе величины 4, и 4 постоянно масштабируются на т, которое меньше единицы, в результате диагональные элементы экспоненциально растут. Таким образом, похоже, что проблема "раздувания" не исчезает при факторизации, как утверждают некоторые авторы. Для уменьшения эффекта "раздувания" можно использовать другие рекурсивные схемы наименьших квааратов (см., например, (14)), но в задаче устранения окулярных артефактов минимальные окулярные движения и другие присущие системе шумы гарантируют, что значения х никогда ие станут равными нулю. Оперативный алгоритм устранения окулярных артефактов с помощью (Л)-факторизации и алгоритма квадратного корня моделировался на компьютере.
Было показано, что полученные результаты подобны результатам, которые позволяют получить автономные эквивалентные средства. 14.2.2.4. Аппаратное обеспечение дпя оперативной системы устранения окупярных артефактов В данном и следующем разделах описывается оперативная система устранения окулярных артефактов, в юторой используется Ш)-алгоритм, Вначале устанавливается целевая спецификация системы.