Хабрахабр

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

А именно основами адаптивной фильтрации. Развивая тему конспектов по магистерской специальности "Communication and Signal Processing" (TU Ilmenau), продолжить хотелось бы одной из основных тем курса "Adaptive and Array Signal Processing".

Для кого в первую очередь была написана эта статья:

1) для студенческой братии родной специальности;
2) для преподавателей, которые готовят практические семинары, но ещё не определились с инструментарием — ниже будут примеры на python и Matlab/Octave;
3) для всех, кто интересуется темой фильтрации.

Что можно найти под катом:

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

В общем, добро пожаловать и давайте разбирать всё по пунктам.

Фильтр его имени мы, по большей части, и будем изучать. Задумчивый человек на фото — знакомый многим, я думаю, Норберт Винер. Однако, нельзя не упомянуть и о нашем соотечественнике — Андрее Николаевиче Колмогорове, чья статья 1941 года также внесла значительный вклад в развитие теории оптимальной фильтрации, которая даже в англоязычных источниках так и называется Kolmogorov-Wiener filtering theory.

Что рассматриваем?

1). Сегодня мы рассматриваем классический фильтр с конечной импульсной характеристикой (КИХ, FIR — finite impulse response), описать который можно следующей простой схемой (рис.

1. Рис. c. Схема КИХ-фильтра для изучения фильтра Винера.[1. 117]

Запишем в матричном виде, что именно будет на выходе данного стенда:

e(n) = d(n) - \hat(n|\mathcal{U}_n) = d(n) - \mathbf{w}^H\mathbf{u} \qquad(1)

Расшифруем обозначения:

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

Что будем оптимизировать?

Оптимизировать, а точнее минимизировать мы будем не просто ошибку, среднюю квадратичную ошибку (MSE — Mean Sqared Error):

MSE: J(\mathbf{w}) = E\{e(n)^2\} \qquad (2)

ожидание. где J(\mathbf{w}) обозначает функцию издержек (cost function) от вектора коэффициентов фильтра, а E\{*\} обозначает мат.

Квадрат в данном случае очень приятен, так как он означает, что перед нами задача выпуклого программирования (я нагуглил только такой аналог английскому convex optimization), что, в свою очередь, подразумевает единственный экстремум (в нашем случае минимум).

2. Рис. Поверхность средней квадратичной ошибки.

121]: У нашей функции ошибки есть каноническая форма [1, c.

J(\mathbf{w}) = \sigma^2_d - \mathbf{w}^H\mathbf{p} - \mathbf{p}^H\mathbf{w} + \mathbf{w}^H\mathbf{R}\mathbf{w} \qquad (3)

где \sigma^2_d — это дисперсия ожидаемого сигнала, \mathbf{p} = E\{\mathbf{u}(n)d^*(n)\} — это вектор кросс-корреляции между входным вектором и ожидаемым сигналом, а \mathbf{R} = E\{\mathbf{u}(n)\mathbf{u}^H(n)\} — это вектор автокорреляции входного сигнала.

Вывод данной формулы здесь (старался понагляднее).

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

\frac{\delta J(\mathbf{w})}{\delta w^*} = - \mathbf{p} + \mathbf{R}\mathbf{w} \qquad (4)

В оптимальном случае (\mathbf{w} = \mathbf{w}_{opt}), ошибка должна быть, конечно же, минимальной, а значит приравниваем производную к нулю:

\mathbf{R}\mathbf{w}_{opt} = \mathbf{p} \qquad (5)

Собственно, вот она, наша печка, от которой мы будем плясать дальше: перед нами система линейных уравнений.

Как будем решать?

Однако, разбор таких вот упрощенных примеров — это лучшее, что можно придумать для понимания основных подходов. Нужно отметить сразу, что оба решения, которые мы рассмотрим ниже, в данном случае являются теоретическими и учебными, так как \mathbf{R} и \mathbf{p} заранее известны (то есть у нас была предполагаемая возможность собрать достаточную статистику для вычисления оных).

Аналитическое решение

Решать эту задачу можно, так сказать, в лоб — с помощью обратных матриц:

\mathbf{w}_{opt} = \mathbf{R}^{-1}\mathbf{p} \qquad (6)

Такое выражение называется уравнением Винера-Хопфа (Wiener–Hopf equation) — оно нам ещё пригодится в качестве некого эталона.

не с ^{-}, а с ^+ (псевдо-инвертирование):
\mathbf{R}^+ = \mathbf{R}^H(\mathbf{R}\mathbf{R}^H)^{-1} Конечно, если быть совсем дотошным, то, наверное, правильнее было бы записать это дело в общем виде, т.е.

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

в нашем случае MMSE — minimum mean square error): Из данного уравнения аналитически можно вывести, чему будет равняться минимальное значение функции издержек (т.е.

J_{min} = J(\mathbf{w}_{opt}) = \sigma^2_d - \mathbf{p}^H\mathbf{R}^{-1}\mathbf{p} \qquad (7)

Вывод формулы здесь (тоже старался покрасочнее).

Хорошо, одно решение есть.

Решение итеративным методом

Для этой цели рассмотрим родной и понятный метод градиентного спуска (method of steepest/gradient descent). Однако, да, решать систему линейных уравнений можно и без инвертирования автокорреляционной матрицы — итеративно (to save computations).

Суть алгоритма можно свести к следующему:

  1. Выставляем искомую переменную в какое-то значение по умолчанию (например, \mathbf{w}(0) = \mathbf{0})
  2. Выбираем некоторый шаг \mu (как именно выбираем, мы поговорим ниже).
  3. И далее, как бы, спускаемся вниз по нашей исходной поверхности (в нашем случае, это поверхность MSE) с заданным шагом \mu и определенной скоростью, определяемой величиной градиента.

Отсюда и название: gradient — градиентный или steepest — пошаговый descent — спуск.

220]). Градиент в нашем случае уже известен: по сути, мы нашли его, когда дифференцировали функцию издержек (поверхность же вогнутая, сравните с [1, с. 220]: Запишем как будет выглядеть формула для итеративного обновления искомой переменной (коэффициентов фильтра) [1, с.

\mathbf{w}(n+1) = \mathbf{w}(n) - \mu [-\mathbf{p} + \mathbf{R}\mathbf{w}(n)] \qquad (8)

где n — это номер итерации.

Теперь давайте поговорим о выборе величины шага.

Перечислим очевидные предпосылки:

  • шаг не может быть отрицательным или нулевым
  • шаг не должен быть слишком большим, иначе алгоритм не сойдется (будет, как бы, перескакивать от края до края, не попадая в экстремум)
  • шаг, конечно, может быть очень маленьким, но и это тоже не совсем желательно — алгоритм будет тратить большее количество времени

222-226]: Относительно фильтра Винера ограничения, конечно, уже были давно найдены [1, с.

0 < \mu < \frac{2}{\lambda_{max}} \qquad (9)

где \lambda_{max} — это наибольшее собственное число автокорреляционной матрицы \mathbf{R}.

Под это дело есть даже целый Eigen filter (см. Кстати говоря, собственные числа и вектора — это отдельная интересная тема в контексте линейной фильтрации. Приложение 1).

Есть и оптимальное, адаптивное решение: Но и это, к счастью, не все.

\mu(n) = \frac{\mathbf{\gamma}(n)^H\mathbf{\gamma}(n)}{\mathbf{\gamma}(n)^H\mathbf{R}\mathbf{\gamma}(n)} \qquad(10)

Как видно из формулы, шаг пересчитывается в каждую итерацию, то есть адаптируется. где \mathbf{\gamma}(n) = \mathbf{p} - \mathbf{R}\mathbf{w}(n) — это отрицательный градиент.

Вывод формулы здесь (много математики - смотреть только таким же отъявленным ботанам, как я).

Окей, под второе решение мы тоже подготовили почву.

А нельзя ли на примерах?

Использовать будем Python 3. Наглядности ради проведем небольшое моделирование. 4. 6.

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

import numpy as np import matplotlib.pyplot as plt
from scipy.linalg import toeplitz def convmtx(h,n): return toeplitz(np.hstack([h, np.zeros(n-1)]),\ np.hstack([h[0], np.zeros(n-1)])) def MSE_calc(sigmaS, R, p, w): w = w.reshape(w.shape[0], 1) wH = np.conj(w).reshape(1, w.shape[0]) p = p.reshape(p.shape[0], 1) pH = np.conj(p).reshape(1, p.shape[0]) MSE = sigmaS - np.dot(wH, p) - np.dot(pH, w) + np.dot(np.dot(wH, R), w) return MSE[0, 0] def mu_opt_calc(gamma, R): gamma = gamma.reshape(gamma.shape[0], 1) gammaH = np.conj(gamma).reshape(1, gamma.shape[0]) mu_opt = np.dot(gammaH, gamma) / np.dot(np.dot(gammaH, R), gamma) return mu_opt[0, 0]

Наш линейный фильтр мы применим для задачи выравнивания канала (channel equalization), основной целью которого является нивелировать различные воздействия этого самого канала на полезный сигнал.

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

Модель системы

Допустим, есть антенная решетка (её мы уже рассматривали в статье про MUSIC).

3. Рис. 32]. Ненаправленная линейная антенная решетка (ULAA — uniform linear antenna array) [2, с.

Определим исходные параметры решетки:

M = 5 # количество элементов решетки (number of sensors)

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

4. Рис. 29]. Модель широкополосного канала при n фиксированных задержках.[3, c. Как вы понимаете, конкретные обозначения роли не играют — далее мы будем использовать немного другие.

Модель принимаемого сигнала для одного сенсора выразим следующим образом:

x(n) = \sum_{l=0}^Lh(l)s(n-l) + \nu(n)

В данном случае n обозначает номер отсчета, h(l) — это отклик канала по l-ому лучу, L — количество регистров задержки, s — передаваемый (полезный) сигнал, \nu(n) — аддитивный шум.

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

\mathbf{x}(n) = \mathbf{H}\mathbf{s}(n) + \mathbf{\nu}(n)

где \mathbf{x}(n) и \mathbf{\nu}(n) — имеют размерность M \times 1, размерность \mathbf{H} равна M \times (M-L), а размерность \mathbf{s}(n) равняется (M-L) \times 1.

Матрица \mathbf{H} в нашем случае будет сверточной матрицей для вектора откликов по каждому лучу. Предположим, что каждый сенсор принимает сигнал тоже с какой-то задержкой, в силу падения волны под каким-то углом. Думаю, в коде будет более понятно:

h = np.array([0.722-1j*0.779, -0.257-1j*0.722, -0.789-1j*1.862])
L = len(h)-1 # number of signal sources
H = convmtx(h,M-L) print(H.shape)
print(H)

Выводом будет:

>>> (5, 3)
>>> array([[ 0.722-0.779j, 0. +0.j , 0. +0.j ], [-0.257-0.722j, 0.722-0.779j, 0. +0.j ], [-0.789-1.862j, -0.257-0.722j, 0.722-0.779j], [ 0. +0.j , -0.789-1.862j, -0.257-0.722j], [ 0. +0.j , 0. +0.j , -0.789-1.862j]])

Далее зададим исходные данные для полезного сигнала и шума:

sigmaS = 1 # мощность полезного сигнала (the signal's s(n) power)
sigmaN = 0.01 # мощность шума (the noise's n(n) power)

Теперь переходим к корреляциям.

Rxx = (sigmaS)*(np.dot(H,np.matrix(H).H))+(sigmaN)*np.identity(M) p = (sigmaS)*H[:,0]
p = p.reshape((len(p), 1))

Вывод формул здесь (тоже простыня для самых отчаянных).

Найдем решение по Винеру:

# Solution of the Wiener-Hopf equation:
wopt = np.dot(np.linalg.inv(Rxx), p)
MSEopt = MSE_calc(sigmaS, Rxx, p, wopt)

Теперь перейдем к методу градиентного спуска.

формулу (9)): Найдем наибольшее собственное число, чтобы из неё потом вывести верхнюю границу шага (см.

lamda_max = max(np.linalg.eigvals(Rxx))

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

coeff = np.array([1, 0.9, 0.5, 0.2, 0.1]) mus = 2/lamda_max*coeff # different step sizes

Определим максимальное количество итераций:

N_steps = 100

Запускаем алгоритм:

MSE = np.empty((len(mus), N_steps), dtype=complex)
for mu_idx, mu in enumerate(mus): w = np.zeros((M,1), dtype=complex) for N_i in range(N_steps): w = w - mu*(np.dot(Rxx, w) - p) MSE[mu_idx, N_i] = MSE_calc(sigmaS, Rxx, p, w)

Теперь сделаем то же самое, но уже для адаптивного шага (формула (10)):

MSEoptmu = np.empty((1, N_steps), dtype=complex)
w = np.zeros((M,1), dtype=complex)
for N_i in range(N_steps): gamma = p - np.dot(Rxx,w) mu_opt = mu_opt_calc(gamma, Rxx) w = w - mu_opt*(np.dot(Rxx,w) - p) MSEoptmu[:, N_i] = MSE_calc(sigmaS, Rxx, p, w)

Получить должны нечто такое:

Отрисовка

x = [i for i in range(1, N_steps+1)] plt.figure(figsize=(5, 4), dpi=300) for idx, item in enumerate(coeff): if item == 1: item = '' plt.loglog(x, np.abs(MSE[idx, :]),\ label='$\mu = '+str(item)+'\mu_{max}$') plt.loglog(x, np.abs(MSEoptmu[0, :]),\ label='$\mu = \mu_{opt}$') plt.loglog(x, np.abs(MSEopt*np.ones((len(x), 1), dtype=complex)),\ label = 'Wiener solution')
plt.grid(True)
plt.xlabel('Number of steps')
plt.ylabel('Mean-Square Error')
plt.title('Steepest descent')
plt.legend(loc='best')
plt.minorticks_on()
plt.grid(which='major')
plt.grid(which='minor', linestyle=':')
plt.show()

5. Рис. Кривые обучения (learning curves) для шагов разных размеров.

Закрепления ради проговорим основные моменты по градиентному спуску:

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

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

А есть что-то более близкое к реальности?

Мы уже несколько раз проговорили, что сбор статистики (т.е. Конечно! Однако, человечество приспособилось и к этим сложностям: вместо подхода детерминированного на практике используются подходы адаптивные. вычисление корреляционных матриц и векторов) в real-time системах — это далеко не всегда доступная роскошь. 246]: Разделить их можно на две большие группы [1, c.

  • вероятностные (stochastic) (например, SG — Stochastic Gradient)
  • и на основе метода наименьших квадратов (например, LMS — Least Mean Squares или RLS — Recursive Least Squares)

Тема адаптивных фильтров неплохо представлена в рамках open-source сообщества (примеры для python):

Однако, будьте внимательны! Во втором примере мне особенно нравится документация. Возможно, такие же проблемы могут возникнуть и при работе с какими-то другими реализациями. Когда я тестировал пакет padasip, я столкнулся со сложностями в обработке комплексных чисел (по дефолту там подразумеваются float64).

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

Коротко взглянем на примеры: рассматривать будем три уже упомянутых нами алгоритма SG, LMS и RLS (моделировать будем на языке MATLAB — каюсь, заготовки уже были, а переписывать всё на питон единообразия ради… ну...).

Описание алгоритмов LMS и RLS можно подсмотреть, например, в доке по padasip.

Описание SG можно посмотреть здесь.

Основное отличие от градиентного спуска состоит в изменяемом градиенте:

\mathbf{w}[n] = \mathbf{w}[n-1] + \mu \left(\mathbf{\hat{p}}[n] - \mathbf{\hat{R}}_{xx}[n]\mathbf{w}[n-1]\right)

при

\mathbf{\hat{R}}_{xx}[n] = \frac{1}{n}\left( (n-1) \mathbf{\hat{R}}_{xx}[n-1] + \mathbf{x}[n]\mathbf{x}[n]^H\right)

\mathbf{\hat{p}}[n] = \frac{1}{n}\left( (n-1) \mathbf{\hat{p}}[n-1] + \mathbf{x}[n]d[n]^*\right)

1) Случай, аналогичный рассмотренному выше.

Исходники (MatLab/Octave).

Исходники можно скачать здесь.

6. Рис. Кривые обучения (learning curves) для LMS, RLS и SG.

Самый быстрый результат даёт RLS, но и он может засбоить при относительно небольшом факторе забывания (forgetting factor). Сразу можно отметить, что при своей относительной простоте алгоритм LMS может в принципе и не прийти к оптимальному решению при относительно большом шаге. Пока молодцом держится SG, но давайте посмотрим на другой пример.

2) Случай, когда канал меняется во времени.

Исходники (MatLab/Octave).

Исходники можно скачать здесь.

7. Рис. Кривые обучения (learning curves) для LMS, RLS и SG (канал изменяется во времени).

Кто бы мог подумать. А вот тут картина уже гораздо интересней: при резком изменении отклика канала самым надежным решением уже кажется LMS. Хотя и RLS при правильном forgetting factor тоже обеспечивает приемлемый результат.

Пара слов о производительности.

То есть разница, в общем-то, сравнимая. Да, конечно каждый алгоритм имеет свою определенную вычислительную сложность, однако по моим замерам моя старенькая машина справляется с одним ансамблем примерно за 120 мкс в итерацию в случае с LMS и SG и примерно за 250 мкс в итерацию в случае RLS.

Спасибо всем, кто заглянул! А на этом у меня сегодня всё.

Литература

  1. S. Haykin S. – Pearson Education India, 2005. Adaptive filter theory.

  2. Handbook on array processing and sensor networks. Haykin, Simon, and KJ Ray Liu. 63. Vol. pp. John Wiley & Sons, 2010. 102-107

  3. (2015). Arndt, D. On Channel Modelling for Land Mobile Satellite Reception (Doctoral dissertation).

Приложение 1

Eigen filter

Основная цель такого фильтра — это максимизировать отношение сигнал-шум (SNR — Signal-to-Noise Ratio).

Но судя по наличию в расчетах корреляций, это тоже больше теоретический конструкт, нежели практическое решение.

Теги
Показать больше

Похожие статьи

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»
Закрыть