Фильтр сигнала ошибки

Основные понятия адаптивной обработки сигналов

Важнейшим признаком классификации является наличие или отсутствие образцового, или опорного сигнала (desired signal, reference signal). При наличии образцового сигнала процесс адаптации называется обучением с учителем (supervised learning). В данном случае адаптивный фильтр стремится сделать свой выходной сигнал максимально близким к образцовому сигналу. О том, что именно подразумевается под близостью сигналов, сказано далее при описании конкретных алгоритмов. Адаптация без использования образцового сигнала называется слепой адаптацией (blind adaptation) или обучением без учителя (unsupervised learning). Разумеется, в этом случае требуется некоторая информация о структуре полезного входного сигнала (например, знание типа и параметров используемой модуляции). Очевидно, что слепая адаптация является более сложной вычислительной задачей, чем адаптация с использованием образцового сигнала.

ЗАМЕЧАНИЕ

У внимательного читателя должен возникнуть вопрос — какой практический смысл могут иметь алгоритмы с использованием образцового сигнала, если при этом выходной сигнал должен быть заранее известен? Однако есть целый ряд практических задач, при решении которых образцовый сигнал оказывается доступен. Подробнее об этом далее в разд. «Применение адаптивных фильтров» этой главы. Забегая вперед, заметим, что в ряде случаев при этом полезным сигналом является не выходной сигнал фильтра, а сигнал ошибки, т. е. разность между образцовым сигналом и выходным сигналом адаптивного фильтра.

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

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

К нелинейным адаптивным системам относятся, в частности, нейронные сети (neural networks), в определенной степени моделирующие работу нервной системы живых организмов. Рассмотрение таких систем выходит за рамки тематики данной книги. За дополнительной информацией можно обратиться, например, к книгам [43, 44].

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

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

Рис. 9.1. Общая структура адаптивного фильтра

Входной дискретный сигнал x(k) обрабатывается дискретным фильтром, в результате чего получается выходной сигнал y(k). Этот выходной сигнал сравнивается с образцовым сигналом d(k), разность между ними образует сигнал ошибки e(k). Задача адаптивного фильтра — минимизировать ошибку воспроизведения образцового сигнала. С этой целью блок адаптации после обработки каждого отсчета анализирует сигнал ошибки и дополнительные данные, поступающие из фильтра, используя результаты этого анализа для подстройки параметров (коэффициентов) фильтра.

В данной главе рассмотрены два адаптивных алгоритма, часто применяемых на практике в различных радиотехнических системах. Это алгоритмы LMS (Least Mean Square, метод наименьшего квадрата, в отечественных источниках иногда используется аббревиатура МНК), и RLS (Recursive Least Squares, рекурсивный метод наименьших квадратов, в отечественных источниках иногда используется аббревиатура РНК).

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

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

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

Оптимальный фильтр Винера

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

Итак, пусть входной дискретный случайный сигнал {x(k)} обрабатывается дискретным фильтром порядка N с коэффициентами {wn} , n = 0, 1, …, N (рис. 9.2).

Отсюда получаем равенство, называемое уравнением Винера — Хопфа (Wiener — Hopf equation):

Умножив слева обе части равенства на обратную корреляционную матрицу R−1 , получаем искомое решение для оптимальных коэффициентов фильтра:

Такой фильтр называется фильтром Винера (Wiener filter). Подстановка (9.8) в (9.5) дает минимально достижимую дисперсию сигнала ошибки:

для фильтра Винера некоррелирован с входным и выходным сигналами фильтра.

В качестве примера с помощью MATLAB рассчитаем оптимальный фильтр для коррекции искажений, вносимых в передаваемый сигнал x0 (k) каналом связи, имеющим четырехэлементную импульсную характеристику следующего вида:

{hk } = {−2, − 4, 6, 3}, k = 0, 1, 2, 3.

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

Корреляционная матрица входного сигнала строится как матрица Теплица на основе данной корреляционной функции:

Восстановление переданного сигнала неизбежно требует внесения некоторой временной задержки, поэтому образцовый сигнал должен представлять собой задержанную копию переданного:

d (k) = x0 (k − k) .

В линии задержки фильтра на k-м шаге находятся отсчеты искаженного сигнала с номерами k, k – 1, k – 2, …, k – N, где N — порядок фильтра. Каждый из этих отсчетов представляет собой линейную комбинацию отсчетов переданного сигнала:

m=−∞

Поскольку отсчеты исходного сигнала считаются статистически независимыми, при вычислении n-го элемента вектора p результат усреднения будет отличен от нуля только для одного слагаемого (9.10). Сразу же учтем тот факт, что средний квадрат сигнала x0 (k) равен единице:

∞

pn = x(k − n)d(k) = ∑ x0 (m) hk−n−m x0 (k − k) =

m=−∞

∞

= ∑ hk−n−m x0 (m) x0 (k − k) = h k−m.

m=−∞

Таким образом, вектор p содержит перевернутую импульсную характеристику канала, при необходимости обрезанную или дополненную нулями с одной или двух сторон:

h k

k−1

Далее приводится код MATLAB-программы, реализующей расчет оптимального фильтра. Код составлен так, чтобы можно было легко варьировать вид импульсной характеристики канала и порядок рассчитываемого фильтра (единственное ограничение — предполагается, что импульсная характеристика канала связи целиком укладывается в пределах вектора p):

h = [-2 -4 6 3]; % импульсная характеристика канала связи N = 32; % порядок рассчитываемого фильтра

r = xcorr(h, N); % двусторонняя КФ импульсной характеристики канала r = r(N+1:end); % односторонняя КФ импульсной характеристики канала R = toeplitz(r); % корреляционная матрица искаженного сигнала

p = zeros(N+1, 1);

delay = 16; % задержка фильтрации

p(delay-length(h)+2:delay+1) = fliplr(h); % вектор взаимных корреляций w = Rp; % коэффициенты оптимального фильтра

subplot(1,2,1)

Рис. 9.3. Импульсные характеристики корректирующего фильтра (слева) и скорректированного канала связи (справа)

ЗАМЕЧАНИЕ

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

На рис. 9.3 приведены графики импульсной характеристики корректирующего фильтра и сквозной импульсной характеристики скорректированного канала связи. Видно, что сквозная характеристика близка к единичному импульсу, что говорит о хорошем качестве коррекции. Столь качественная коррекция оказалась возможна потому, что частотная характеристика канала не имеет глубоких провалов (с помощью функции freqz легко убедиться в том, что неравномерность АЧХ канала составляет около 11 дБ).

Градиентный поиск оптимального решения

Один из наиболее распространенных адаптивных алгоритмов основан на поиске минимума целевой функции (9.4) методом наискорейшего спуска. При использовании данного способа оптимизации вектор коэффициентов фильтра зависит от номера итерации k: w(k). На каждой итерации вектор коэффициентов смещается на величину, пропорциональную градиенту целевой функции в данной точке (см. формулу (9.6)):

2

где μ — положительный коэффициент, называемый размером шага.

Анализ сходимости данного процесса приведен, например, в [13]. Показано, что алгоритм сходится, если

где λmax — максимальное собственное число корреляционной матрицы R. Скорость сходимости при этом зависит от разброса собственных чисел корреляционной матрицы R — чем меньше отношение λmax/λmin, тем быстрее сходится итерационный процесс.

Адаптивный алгоритм LMS

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

pˆ (k) = d (k) x(k) .

При использовании данных оценок формула (9.11) принимает следующий вид:

w(k +1) = w(k) + µ d (k)x(k) − µ x(k)xT (k) w(k) =

(9.13)

= w(k) + µ x(k)(d (k) − xT (k) w(k)).

Выражение, стоящее в скобках, согласно (9.3), представляет собой разность между образцовым сигналом и выходным сигналом фильтра на k-м шаге, т. е. ошибку фильтрации e(k). С учетом этого выражение для рекурсивного обновления коэффициентов фильтра оказывается очень простым:

Алгоритм адаптивной фильтрации, основанный на формуле (9.14), получил название LMS (Least Mean Square, метод наименьшего квадрата). Можно получить ту же формулу и несколько иным образом: использовав вместо градиента статисти-

чески усредненного квадрата ошибки e2 (k) градиент его мгновенного значения e2 (k) .

Несмотря на простоту алгоритма LMS, анализ его поведения оказывается чрезвычайно сложной задачей, для которой не существует точного аналитического решения. Частные результаты, полученные при использовании различных приближений, можно найти в большом количестве книг и статей. Достаточно подробный анализ выполнен в упоминавшейся в начале главы книге С. Хайкина, при этом получены следующие результаты. Алгоритм LMS сходится в среднем (converges in the mean values; это означает, что математические ожидания коэффициентов фильтра при

k → ∞ стремятся к оптимальному решению (9.8)), если коэффициент μ удовлетворяет тому же условию (9.12), которое требуется для сходимости детерминированного градиентного спуска. Однако это еще не гарантирует того, что средние квадраты коэффициентов будут стремиться к фиксированным значениям. Алгоритм

сходится в среднем квадрате (converges in the mean square) в более узком диапазоне значений μ:

n=0

где λk — собственные числа корреляционной матрицы R, а σ2x — средний квадрат входного сигнала фильтра.

На формуле (9.15) основан нормированный (normalized) LMS-алгоритм, в котором коэффициент μ на каждом шаге рассчитывается исходя из энергии сигнала, содержащегося в линии задержки:

где μ0 — нормированное значение μ, лежащее в диапазоне от 0 до 2, а ε — малая

положительная константа, назначение которой — ограничить рост μ при нулевом сигнале на входе фильтра (из формулы (9.16) видно, что максимально возможное значение μ составляет μ0 ε ).

Даже если LMS-алгоритм сходится, дисперсии коэффициентов фильтра при k → ∞ не стремятся к нулю — коэффициенты флуктуируют вокруг оптимальных значений. Из-за этого ошибка фильтрации в установившемся режиме оказывается боль-

ше, чем ошибка для винеровского фильтра e2 (k)min (см. ранее формулу (9.9)):

e2 (k)LMS = e2 (k)min + Eex ,

где Eex — средний квадрат избыточной ошибки (excess error) алгоритма LMS. В той же книге С. Хайкина приводится следующая приближенная формула для так называемого коэффициента расстройки (misalignment), равного отношению средних квадратов избыточной и винеровской ошибок:

n=0

Значение коэффициента μ влияет на два главных параметра LMS-фильтра — скорость сходимости и коэффициент расстройки. Чем больше μ, тем быстрее сходится алгоритм, но тем больше, согласно (9.17), становится коэффициент расстройки, и наоборот.

Основным достоинством алгоритма LMS является предельная вычислительная простота — для подстройки коэффициентов фильтра на каждом шаге нужно выполнить N + 1 пар операций «умножение-сложение». Платой за простоту являются медленная сходимость и повышенная (по сравнению с минимально достижимым значением (9.9)) дисперсия ошибки в установившемся режиме.

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

Детерминированная задача оптимальной фильтрации

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

Итак, пусть обработке подвергается последовательность отсчетов {x(k)}, коэффициенты нерекурсивного фильтра порядка N образуют набор {wn} , n = 0, 1, …, N, а отсчеты образцового сигнала равны {d(k)}. Выходной сигнал фильтра определяется формулой (9.1), а ошибка воспроизведения образцового сигнала — формулой (9.2) или, в векторном виде, (9.3).

Теперь сформулируем детерминированную оптимизационную задачу: мы хотим отыскать такие коэффициенты фильтра {wn} , чтобы суммарная квадратичная ошибка воспроизведения образцового сигнала была минимальной:

k =0

Для решения задачи необходимо перейти в формуле (9.3) к матричной записи, так сказать, вдоль координаты k, получив формулы для векторов-столбцов выходного сигнала y и ошибки воспроизведения входного сигнала e:

Здесь d — вектор-столбец отсчетов образцового сигнала, а X — матрица, столбцы которой представляют собой содержимое линии задержки фильтра на разных тактах:

Выражение (9.18) для суммарной квадратичной ошибки в матричном виде можно записать следующим образом:

Подставив (9.19) в (9.20), имеем

J (w) = (d − XT w)T (d − XTw) =

=dTd − (XTw)T d − dT (XTw) + (XTw)T (XTw) =

=dTd − wT Xd − dT XTw + wT XXTw.

Для нахождения минимума необходимо вычислить градиент данного функционала и приравнять его к нулю:

J (w) = −2Xd + 2XXTw = 0.

Отсюда легко получается искомое оптимальное решение:

В формуле (9.21) прослеживается близкое родство с формулой (9.8), описывающей оптимальный в статистическом смысле фильтр Винера. Действительно, если

учесть, что (XXT )−1 K дает оценку корреляционной матрицы сигнала, полученную по одной реализации эргодического случайного процесса путем временного усреднения, а X d / K является аналогичной оценкой взаимных корреляций между образцовым сигналом и содержимым линии задержки фильтра, то формулы (9.8) и (9.21) совпадут.

Адаптивный алгоритм RLS

В принципе, в процессе приема сигнала можно на каждом очередном шаге пересчитывать коэффициенты фильтра непосредственно по формуле (9.21), однако это связано с неоправданно большими вычислительными затратами. Действительно, на каждом шаге увеличивается размер матрицы X; кроме того, необходимо каждый

раз заново вычислять обратную матрицу (X XT )−1 .

Сократить вычислительные затраты можно, если заметить, что на каждом шаге к матрице X добавляется лишь один новый столбец, а к вектору d — один новый элемент. Это дает возможность организовать вычисления рекурсивно. Соответствующий алгоритм получил название рекурсивного метода наименьших квадратов (Recursive Least Squares, RLS).

На k-м шаге оптимальный вектор коэффициентов фильтра, согласно (9.21), равен

w(k) = (X(k)XT (k))−1 X(k)d(k) . (9.22)

Посмотрим, какие новые данные добавляются к X(k) и d(k) на следующем шаге и как мы можем использовать их для обновления рассчитанного ранее решения — ведь именно в таком постоянном обновлении состоит суть работы адаптивного фильтра.

К матрице X на каждом шаге добавляется новый столбец x(k + 1), а к вектору d — новый элемент d(k + 1):

В алгоритме RLS производится рекурсивное обновление оценки обратной корреляционной матрицы

P(k) = (X(k)XT (k))−1 .

При переходе к следующему шагу мы имеем

которое является справедливым для произвольных квадратных невырожденных матриц A и C и произвольных матриц B и D совместимых размеров.

Установим между (9.24) и (9.25) следующее соответствие:

A = X(k)XT (k) = P−1(k) , A−1 = P(k) (квадратная матрица);

B = x(k + 1) (вектор-столбец);

C = 1 (скаляр);

D = xT (k +1) (вектор-строка).

В результате равенство (9.24) можно записать следующим образом:

P(k +1) = P(k) − P(k)x(k +1) (1 + xT (k +1)P(k)x(k +1))−1 xT (k +1)P(k) . (9.26)

Здесь следует заметить, что фрагмент выражения в круглых скобках, возводимый в минус первую степень, представляет собой скаляр. Далее используем (9.23) и (9.26) в выражении (9.22) для коэффициентов фильтра:

= P(k +1) ( X(k)d(k) + x(k +1)d(k +1)) =

Раскроем скобки в этом выражении:

w(k +1) = P(k)X(k)d(k) − P(k)x(k +1)xT (k +1)P(k) X(k)d(k) + 1 + xT (k +1)P(k)x(k +1)

+ P(k)x(k +1)d (k +1) − P(k)x(k +1)xT (k +1)P(k) x(k +1)d(k +1). 1 + xT (k +1)P(k)x(k +1)

Первое слагаемое в полученной формуле, согласно (9.22), представляет собой коэффициенты оптимального фильтра для k-го шага — w(k). Этот же вектор можно выделить в качестве множителя во втором слагаемом. Что касается третьего и четвертого слагаемых, в них можно выделить общий множитель P(k) x(k + 1) d(k + 1):

Выполнив вычитание в круглых скобках, получаем

Теперь видно, что второе и третье слагаемые имеют общий множитель, который можно вынести за скобки:

Далее заметим, что произведение xT (k +1) w(k) есть результат обработки нового (поступившего на (k + 1)-м шаге) входного сигнала фильтром со старыми (т. е. имеющимися на данный момент) коэффициентами w(k). Иными словами, это произведение представляет собой выходной сигнал адаптивного фильтра y(k + 1). Соответственно, получается, что разность, стоящая в скобках, — не что иное, как ошибка фильтрации e(k + 1):

Наконец, в качестве последнего штриха введем обозначение для вынесенного за скобки векторного множителя:

С учетом этого формула для коэффициентов фильтра примет вид

Вектор K(k + 1) называют вектором коэффициентов усиления.

ЗАМЕЧАНИЕ

Структура формулы (9.27) совпадает со структурой формулы (9.14), описывающей обновление коэффициентов фильтра для LMS-алгоритма — к текущему вектору коэффициентов фильтра добавляется слагаемое, рассчитываемое как произведение сигнала ошибки на некий вектор коэффициентов усиления. Различие состоит лишь в способе получения этого вектора — в алгоритме LMS он пропорционален содержимому линии задержки фильтра, а в алгоритме RLS рассчитывается более сложным образом.

Итак, при использовании адаптивного алгоритма RLS необходимо на каждом временном такте выполнить следующие шаги:

1.При поступлении новых входных данных x(k) производится фильтрация сигнала с использованием текущих коэффициентов фильтра w(k – 1) и вычисление величины ошибки воспроизведения образцового сигнала:

y(k) = x(k)T w(k −1) , e(k) = d (k) − y(k) .

2.Рассчитывается вектор-столбец коэффициентов усиления (обратите внимание на то, что при вычислениях не используются предыдущие значения — вектор K всякий раз рассчитывается заново, т. е. вычисления не являются рекурсивными; кроме того, еще раз отметим, что знаменатель дроби в следующих двух формулах является скаляром, а не матрицей):

3. Выполняется обновление оценки обратной корреляционной матрицы сигнала:

4. Наконец, производится обновление коэффициентов фильтра: w(k) = w(k −1) + K(k)e(k) .

Осталось разобраться с начальными значениями рекурсивно обновляемых матрицы P и вектора w. Вектор коэффициентов фильтра w перед началом работы алгоритма обычно заполняется нулями. Что касается матрицы P, то строгий анализ показыва-

606 Глава 9

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

На практике диагональ матрицы заполняется большими положительными значениями. Например, в [13] рекомендуется величина 100σ2x .

По сравнению с LMS-алгоритмом алгоритм RLS требует значительно большего числа вычислительных операций. При оптимальной организации вычислений для

обновления коэффициентов фильтра необходимо 2,5N2 + 4N пар операций «ум- ножение—сложение» [13]. Под оптимальной организацией здесь, в частности, подразумевается учет симметрии матрицы P. Таким образом, число операций в алгоритме RLS квадратично возрастает с увеличением порядка фильтра.

Зато алгоритм RLS сходится значительно быстрее, чем алгоритм LMS. Строго говоря, само понятие сходимости применимо здесь лишь условно, поскольку данный алгоритм не является алгоритмом последовательного приближения. Алгоритм RLS на каждом шаге дает оптимальные коэффициенты фильтра, соответствующие формуле (9.21), и переходный процесс в начале работы связан лишь с рекурсивным расчетом оценки матрицы P и постепенным заполнением линии задержки фильтра отсчетами входного сигнала.

Экспоненциальное забывание

В формулах (9.18) и (9.20) значениям ошибки на всех временных тактах придается одинаковый вес. В результате, если статистические свойства входного сигнала со временем изменяются, это приведет к ухудшению качества фильтрации. Чтобы дать фильтру возможность отслеживать нестационарный входной сигнал, можно применить в (9.18) экспоненциальное забывание (exponential forgetting), при котором вес прошлых значений сигнала ошибки экспоненциально уменьшается:

K −1

J (w) = ∑λK −1−k | e(k) |2 ,

k=0

где λ — коэффициент забывания (forgetting factor), 0 <λ ≤ 1.

При использовании экспоненциального забывания формулы (9.28) и (9.29) принимают следующий вид:

Применение адаптивных фильтров

Адаптивные фильтры в настоящее время нашли применение во многих системах обработки сигналов. В данном разделе мы рассмотрим лишь некоторые из возможных областей их использования.

Идентификация систем

Все способы применения адаптивных фильтров так или иначе сводятся к решению задачи идентификации, т. е. определения характеристик, некоторой системы. Возможны два варианта идентификации — прямая и обратная. В первом случае адаптивный фильтр включается параллельно с исследуемой системой (рис. 9.4, а).

Рис. 9.4. Идентификация систем с помощью адаптивного фильтра: а — прямая, б — обратная

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

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

Теперь перейдем от этих обобщенных схем к рассмотрению более конкретных вариантов.

Линейное предсказание

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

Рис. 9.5. Линейное предсказание с помощью адаптивного фильтра

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

Подавление шума

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

Рис. 9.6. Подавление шума с помощью адаптивного фильтра

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

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

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

Система речевой связи не является единственно возможным объектом применения рассматриваемой системы шумоподавления. Возможны и иные области использования данной идеи.

Выравнивание частотной характеристики канала связи

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

При использовании адаптивного фильтра в качестве эквалайзера возникает проблема получения образцового сигнала. Эта проблема решается путем передачи специального настроечного сигнала перед началом передачи данных. В качестве такого настроечного сигнала обычно используется псевдослучайная последовательность символов. Алгоритм формирования этого сигнала известен приемной стороне, поэтому образцовый сигнал может быть сгенерирован автономно и использован для обучения адаптивного фильтра. Этот режим работы называется режимом обучения (training mode) (рис. 9.7).

Рис. 9.7. Выравнивание частотной характеристики канала связи с помощью адаптивного фильтра

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

Оно используется в качестве образцового сигнала, а разность между этим значением и принятым сигналом дает сигнал ошибки, используемый для адаптации. На рис. 9.8 это иллюстрируется применительно к 16-позиционной квадратурной модуляции.

Im

Принятый сигнал

Сигнал ошибки

Re

Ближайшее допустимое значение выбирается в качестве образцового сигнала

Рис. 9.8. Формирование образцового сигнала и сигнала ошибки в режиме оценивания

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

ЗАМЕЧАНИЕ

В некоторых, в частности многопользовательских, системах связи передавать настроечный сигнал не представляется возможным. В этом случае может использоваться слепая (blind) адаптация, рассмотрение которой выходит за рамки данной книги. Некоторые алгоритмы слепой адаптации приведены в [35].

Эхоподавление

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

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

Для правильной работы системы эхоподавления необходимо, чтобы передаваемый и принимаемый сигналы были некоррелированы. Для обеспечения некоррелиро-

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

Эхоподавление, осуществляемое согласно схеме рис. 9.10, используется во всех современных модемах.

Рис. 9.10. Эхоподавление, осуществляемое с помощью адаптивного фильтра

Объекты адаптивной фильтрации пакета Filter Design

В версии 4.7 пакета Filter Design (R2010a) адаптивные фильтры реализованы только в виде объектов MATLAB, а существовавшие ранее функции адаптивной фильтрации не только объявлены устаревшими, но и удалены.

Общие сведения об использовании объектов в MATLAB были приведены в главе 4. Для создания объекта адаптивного фильтра служит функция adaptfilt. При ее вызове необходимо указать требуемый алгоритм адаптивной фильтрации — для этого после имени функции ставится точка и указывается соответствующий идентификатор метода (конструктора):

ha = adaptfilt.algorithm(…)

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

которые имеются в справочной системе MATLAB, а также литературу, упомянутую в начале этой главы.

Полный список реализованных алгоритмов адаптивной фильтрации и имена соответствующих конструкторов объектов adaptfilt приведены в табл. 9.1.

Набор входных параметров конструктора зависит от реализуемого алгоритма. Рассмотрим более подробно использование объектов адаптивных фильтров для тех алгоритмов, которые были ранее описаны в этой главе. Вызов конструкторов для них имеет следующий вид:

ha = adaptfilt.lms(L, step, leakage, coeffs, states)

ha = adaptfilt.nlms(L, step, leakage, offset, coeffs, states) ha = adaptfilt.sd(L, step, leakage, coeffs, states)

ha = adaptfilt.se(L, step, leakage, coeffs, states) ha = adaptfilt.ss(L, step, leakage, coeffs, states) ha = adaptfilt.rls(L, lambda, invcov, coeffs, states)

Назначение входных параметров расшифровано в табл. 9.2. Там же приведены и значения по умолчанию, которые имеются для всех параметров конструкторов.

Таблица 9.2. Параметры конструкторов рассматриваемых адаптивных фильтров

ЗАМЕЧАНИЕ

В табл. 9.2 и далее в табл. 9.3 упоминается не рассматривавшийся в теоретической части главы параметр LMS-алгоритма — коэффициент утечки (leakage). Это коэффициент, лежащий в диапазоне от 0 до 1, на который может умножаться старый вектор коэффициентов w(k) в формуле (9.14) для улучшения стабильности алгоритма.

Список свойств объектов адаптивных фильтров для рассматриваемых алгоритмов приведен в табл. 9.3. Считывать и задавать значения этих свойств можно с помощью функций get и set соответственно (см. разд. «Общая информация об объектах MATLAB» главы 4).

Обработка сигнала адаптивным фильтром осуществляется с помощью функции

filter:

[y, e] = filter(ha, x, d)

Здесь ha — объект адаптивного фильтра, x — вектор отсчетов входного сигнала, d — вектор отсчетов образцового сигнала. Результатами работы функции являются векторы отсчетов выходного сигнала y и сигнала ошибки фильтрации e.

Как показывает синтаксис вызова функции фильтрации, образцовый сигнал должен быть известен заранее, что делает невозможным применение данных функций в режиме адаптации по оценке сигнала (decision-directed mode). Кроме того, реализация данной функции подразумевает на каждом временном такте сдвиг данных

влинии задержки фильтра на один отсчет. Это затрудняет применение объектов

втех ситуациях, когда данные обновляются иным образом (это имеет место, например, в так называемых дробных (fractionally-spaced equalizer) эквалайзерах).

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

ЗАМЕЧАНИЕ

Возможность адаптации по оценке сигнала реализована в пакете расширения Communications Toolbox, где для этого имеются специальные объекты эквалайзеров.

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

Аналогичным образом объекты адаптивных фильтров поддерживаются и визуализатором фильтров FVTool:

fvtool(ha)

К методам, специфическим именно для адаптивных фильтров, относятся три функции, перечисленные далее:

[mumax, mumaxmse] = maxstep(ha, x) — данная функция производит оценку максимально допустимой величины размера шага µ для объекта адаптивного фильтра ha (функция применима для следующих вариантов алгоритма LMS: lms, nlms, blms, blmsfft, se). Входной параметр x — пример возможного входного сигнала фильтра (для алгоритма nlms задавать этот параметр не нужно). Этот параметр может быть матрицей, тогда ее столбцы рассматриваются как различные реализации входного случайного сигнала. Выходной параметр mumax — предельная величина µ, при которой обеспечивается сходимость алгоритма в среднем (см. формулу (9.12)), выходной параметр mumaxmse — предельная величина µ, при которой обеспечивается сходимость алгоритма в среднем квадрате (см. формулу (9.15));

[mmse, emse, meanw, mse, tracek] = msepred(ha, x, d, m) — данная функция

производит оценку среднего квадрата ошибки фильтрации для объекта адаптивного фильтра ha (функция применима для следующих вариантов алгоритма LMS: lms, nlms, blms, blmsfft, se). Входные параметры x и d — векторы отсчетов входного и образцового сигнала соответственно. Выходные параметры: mmse — средний квадрат ошибки фильтрации в установившемся режиме, emse — избыточная ошибка фильтрации (см. формулу (9.17)), meanw — матрица предсказанных математических ожиданий коэффициентов фильтра для всех моментов времени, mse — вектор зависимости среднего квадрата ошибки от времени, tracek — вектор зависимости от времени для суммарного квадратичного отклонения коэффициентов фильтра от оптимальных. При расчете трех последних выходных параметров может производиться прореживание по времени с коэффициентом m (по умолчанию m = 1);

[mse, meanw, w, tracek] = msesim(ha, x, d, m) — данная функция по назначе-

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

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

Сказанному соответствуют и итоговые импульсные характеристики фильтров (см. рис. 9.11, справа) — характеристика LMS-фильтра хорошо совпадает с экспоненциально затухающими колебаниями импульсной характеристики анализируемой системы, а характеристика RLS-фильтра выглядит абсолютно случайным набором значений.

ЗАМЕЧАНИЕ

В данном примере не использовалось экспоненциальное забывание, дающее RLSалгоритму возможность отслеживать изменения статистических свойств обрабатываемого сигнала. Чтобы использовать забывание, полю ForgettingFactor объекта RLS-фильтра в данном примере необходимо присвоить значение, меньшее единицы.

Линейное предсказание

Реализуем с помощью RLS-фильтра линейное предсказание сигнала, сформированного с помощью авторегрессионной модели. В качестве коэффициентов модели используем знаменатель фильтра Баттерворта 8-го порядка с частотой среза, равной 0,3 от частоты Найквиста:

>>ha = adaptfilt.rls(16); % создание RLS-фильтра

>>% линейное предсказание с помощью адаптивного фильтра

>>[y1, e] = filter(ha, y(1:end-1), y(2:end));

>>subplot(2, 2, 1)

Результаты работы программы показаны на рис. 9.12. Видно, что коэффициенты адаптивного фильтра совпадают с коэффициентами использованной авторегрессионной модели (при сравнении верхних графиков на рис. 9.12 следует учесть, что первый элемент вектора a, равный единице, не входит в число коэффициентов модели; кроме того, коэффициенты модели и коэффициенты предсказывающего фильтра имеют разные знаки (см. рис. 5.13 в главе 5)). На нижнем графике рис. 9.12, демонстрирующем поведение ошибки фильтрации во времени, хорошо видно, что после начального переходного процесса устанавливается стационарный режим работы. Согласно идее линейного предсказания (см. разд. «Авторегрессионная модель» главы 5), сигнал ошибки фильтрации при этом является оценкой белого шума, возбуждающего авторегрессионную модель.

Шум noise1 мы добавим к речевому сигналу. На графике (см. рис. 9.13, в центре) видно, что полезный сигнал визуально не выделяется над шумом:

Завершая рассмотрение данного примера, следует отметить, что сигнал был более чем в четыре раза (более чем на 6 дБ) слабее по мощности, чем шум:

На рис. 9.14 показаны результаты работы программы. В левом столбце — зависимость сигнала ошибки от номера шага, в правом — график выходного сигнала

Рис. 9.14. Компенсация искажений, вносимых каналом связи, с помощью различных вариантов LMS-алгоритма: слева — ошибка фильтрации, справа — выходной сигнал

Медиана для N значений лёгкая

// облегчённый вариант медианы для N значений
// предложен Виталием Емельяновым, доработан AlexGyver
// фильтр лёгкий с точки зрения кода и памяти, но выполняется долго
// возвращает медиану по последним NUM_READ вызовам

#define NUM_READ 10

int findMedianN(int newVal) {
  static int buffer[NUM_READ];  // статический буфер
  static byte count = 0;    // счётчик
  buffer[count] = newVal;
  if (++count >= NUM_READ) count = 0;  // перемотка буфера
  
  int buf[NUM_READ];    // локальный буфер для медианы
  for (byte i = 0; i < NUM_READ; i++) buf[i] = buffer[i];  
  for (int i = 0; i <= (int) ((NUM_READ / 2) + 1); i++) { // пузырьковая сортировка массива (можно использовать любую)
    for (int m = 0; m < NUM_READ - i - 1; m++) {
      if (buf[m] > buf[m + 1]) {
        int buff = buf[m];
        buf[m] = buf[m + 1];
        buf[m + 1] = buff;
      }
    }
  }
  int ans = 0;
  if (NUM_READ % 2 == 0) {             // кол-во элементов в массиве четное (NUM_READ - последний индекс массива)
    ans = buf[(int) (NUM_READ / 2)];   // берем центральное
  } else {
    ans = (buf[(int) (NUM_READ / 2)] + buf[((int) (NUM_READ / 2)) + 1]) / 2; // берем среднее от двух центральных
  }
  return ans;
}

Медиана для N значений лёгкая оптимальная

// облегчённый вариант медианы для N значений
// предложен Виталием Емельяновым, доработан AlexGyver
// возвращает медиану по последним NUM_READ вызовам
// НАВЕРНОЕ ЛУЧШИЙ ВАРИАНТ!

#define NUM_READ 10  // порядок медианы

// медиана на N значений со своим буфером, ускоренный вариант
float findMedianN_optim(float newVal) {
  static float buffer[NUM_READ];  // статический буфер
  static byte count = 0;
  buffer[count] = newVal;

  if ((count < NUM_READ - 1) and (buffer[count] > buffer[count + 1])) {
    for (int i = count; i < NUM_READ - 1; i++) {
      if (buffer[i] > buffer[i + 1]) {
        float buff = buffer[i];
        buffer[i] = buffer[i + 1];
        buffer[i + 1] = buff;
      }
    }
  } else {
    if ((count > 0) and (buffer[count - 1] > buffer[count])) {
      for (int i = count; i > 0; i--) {
        if (buffer[i] < buffer[i - 1]) {
          float buff = buffer[i];
          buffer[i] = buffer[i - 1];
          buffer[i - 1] = buff;
        }
      }
    }
  }
  if (++count >= NUM_READ) count = 0;
  return buffer[(int)NUM_READ / 2];
}

float _err_measure = 0.8;  // примерный шум измерений
float _q = 0.1;   // скорость изменения значений 0.001-1, варьировать самому

float simpleKalman(float newVal) {
  float _kalman_gain, _current_estimate;
  static float _err_estimate = _err_measure;
  static float _last_estimate;

  _kalman_gain = (float)_err_estimate / (_err_estimate + _err_measure);
  _current_estimate = _last_estimate + (float)_kalman_gain * (newVal - _last_estimate);
  _err_estimate =  (1.0 - _kalman_gain) * _err_estimate + fabs(_last_estimate - _current_estimate) * _q;
  _last_estimate = _current_estimate;
  return _current_estimate;
}

Особенности использования

• Хорошо фильтрует и постоянный шум, и резкие выбросы
• Делает только одно измерение за раз, не блокирует код на длительный период
• Слегка запаздывает, как бегущее среднее
• Подстраивается в процессе работы
• Чем чаще измерения, тем лучше работает
• Алгоритм весьма тяжёлый, вычисление длится ~90 мкс при системной частоте 16 МГц

1. ^ Thakor, N.V.; Zhu, Yi-Sheng (1991-08-01). «Applications of adaptive filtering to ECG analysis: noise cancellation and arrhythmia detection». IEEE Transactions on Biomedical Engineering. 38 (8): 785–794. doi:10.1109/10.83591. ISSN 0018-9294. PMID 1937512. S2CID 11271450.
2. ^ Widrow, Bernard; Stearns, Samuel D. (1985). Adaptive Signal Processing (1st ed.). Prentice-Hall. p. 329. ISBN 978-0130040299.
3. ^ Widrow p 304
4. ^ Widrow p 212
5. ^ Widrow p 313
6. ^ Widrow p 100
7. ^ ^{a b} Widrow p 103
8. ^ Danilo Comminiello; José C. Príncipe (2018). Adaptive Learning Methods for Nonlinear System Modeling. Elsevier Inc. ISBN 978-0-12-812976-0.
9. ^ M.Poluektov and A.Polar. Urysohn Adaptive Filter. 2019.
10. ^ «Nonlinear Adaptive Filtering». ezcodesample.com.
11. ^ Weifeng Liu; José C. Principe; Simon Haykin (March 2010). Kernel Adaptive Filtering: A Comprehensive Introduction (PDF). Wiley. pp. 12–20. ISBN 978-0-470-44753-6.