Поскольку пассивный радиолокатор (Passive Coherent Locator) и радиопеленгатор — родственные души, посмотрим как может пригодится адаптация в кольцевой антенной решетке. Что есть главное зло для PCL? Оригинальный сигнал базовой станции. Что есть главное зло для РП? Переотражения от местных предметов, или как любят говорить связисты, многолучевое распространение. И в том и другом случае, антенная решетка, обладая направленными свойствами, создает максимально благоприятные условия для полезного сигнала и по мере возможности — как там лягут боковые лепестки диаграммы направленности (ДН) — давит помеху.
Вот в том то и дело, что по мере возможности. Помеха может встать по азимуту на хороший боковой лепесток ДН, и прощай динамический диапазон и точность. Возникает задача: можно ли сформировать диаграмму таким образом, чтобы и пропустить полезный сигнал без потерь, и подавить помеховый? Другими словами: можно ли так «деформировать» классическую ДН, чтобы в направлении полезного сигнала смотрел главный лепесток, а в направлении помехи — лепестка не было вообще, то есть полный провал?
Структура тракта обработки для адаптивного алгоритма
Данную задачу решает определенный класс адаптивных алгоритмов. Структурно все выглядит достаточно просто, как на этой схеме. На рисунке — кольцевая антенная решетка радиусом R, состоящая из N элементов. В общем, форма решетки может быть произвольной; кольцевую мы выбрали из соображений практичности.
Кольцевая антенная решетка с весовыми коэффициентами
Как видно из рисунка, мы обходимся минималистким набором оборудования. Сделаем решетке один единственный выход, в который сведем все элементы, пропущенные через блоки W.
W — это элементы управления амплитудно-фазовым раскрывом с коэффициентами передачи в виде комплексных чисел w. Модуль комплексного числа — коэффициент усиления, аргумент — доворот фазы. В общем все как обычно — набору элементов антенной решетки соответствует набор весовых коэффициентов. На выходе сумматора будем наблюдать сигнал x, который зависит от направления прихода и параметров волны и значений весовых коэффициентов w. Замечу что нас совершенно не волнует геометрия решетки — пусть об этом заботятся сами коэффициенты w, которые при желании могут сотворить из кольцевой решетки линейную или еще какую нибудь (для определенного угла прихода, конечно).
Тренировка решетки — берем количеством
Итак, начнем. «Поскольку диссертации никто не читает, примем w=0«. Здесь эта крылатая фраза абсолютно к месту, поскольку в общем то все равно, какие значения весовые коэффициенты будут принимать поначалу. Наша цель — найти их значения в процессе адаптации. Адаптация сродни магии: процесс очень простой, но почему получается именно такой результат — это загадка. Кстати, какой результат мы хотим получить? А вот какой: положим, полезный сигнал приходит с направления 0°, и хотим максимум ДН в этом направлении. Или, другими словами, на выходе сумматора решетки мы рассчитываем получить желаемый сигнал d=1.
Итак, поехали. Включаем источник в направлении 0°, принимаем на антенную решетку, обрабатываем весами которые на первом этапе у нас w=0 и получаем некий сигнал x. Есть все основания полагать, что мы получили совсем не то что нам нужно, но это нас не беспокоит. Отличие полученного x от желаемого d называется сигналом ошибки e = d — x. Задача адаптивного алгоритма — непрерывно минимизировать эту ошибку. Не углубляясь в теорию и доверяя первоисточникам, поверим тому, что нужно сформировать поправку к весовым коэффициентам w и величина поправки равна
ΔW = μ ⋅ S ⋅ e*
Не забываем о том, что имеем дело с комплексными переменными, и значок * рядом с переменной ошибки соответствует комплексно сопряженной величине. S — вектор сигналов элементов решетки, W — вектор весовых коэффициентов. μ — некоторый параметр, который будет определять качество алгоритма. Чем он меньше, тем лучше решение будет сходиться к ожидаемому результату. Если этот параметр больше, то решение будет находиться быстрее, но есть вероятность что задача «пойдет вразнос».
Есть максимум
Вот собственно и все! Что делаем дальше? Много-много раз повторяем эту процедуру: вносим поправки в весовые коэффициенты, находим ошибку, вычисляем поправку и так далее. Такой процесс называется тренировкой антенной решетки. В результате получим вектор весовых коэффициентов, который будет задавать вот эту синюю диаграмму (на красный цвет пока не смотрим):
Диаграмма направленности адаптивной кольцевой антенной решетки радиопеленгатора и пассивного радиолокатора в линейных координатах
Синим цветом показана ДН решетки, повернутая в направлении полезного сигнала по азимуту 0°. Первая часть задачи решена! Заметим, что в постановке задачи не использовались данные по геометрии антенной решетки — адаптивный алгоритм все сделал сам.
Давим помеху сохраняя полезный сигнал
Теперь усложним задачу. Поставим помеху в направлении к примеру 35° и попробуем ее подавить, исходя из того что желаемый сигнал примет значение d = 0. Если мы запустим адаптивный алгоритм в такой формулировке, то он свою задачу выполнит, но довольно своеобразным способом: все весовые коэффициенты получат значение ноль, на выходе для направления помехи тоже ноль (а вы разве не этого хотели?), и как вы догадались со всех остальных направлений тоже ноль (а вот этого мы никак не ожидали). Что же делать? В этом случае мы используем другую стратегию тренировки решетки, комбинируя в определенной последовательности измерения как полезного сигнала, так и помехи с разных углов. Несмотря на справедливое ощущение, что в этом случае мы сваливаем все в одну непонятную кучу, адаптация проходит замечательно: на определенном цикле веса выстраиваются таким образом, чтобы и максимум сохранить — в направлении полезного сигнала 0°, и подавить помеху в направлении 35°.
Модифицированная ДН показана на рисунке красным цветом. Видно, что она деформирована таким образом, что в направлении помехи возникает существенный провал. Привожу здесь значения этих волшебных весовых коэффициентов, которые делают это:
Заметим на разброс коэффициентов передачи: для решения такой задачи простого доворота фазы для каждого элемента решетки уже недостаточно, и в игру вступают модули комплексных чисел.
В заключение приведем диаграмму в логарифмических декартовых координатах, где минимум в направлении помехи 35° выражен особенно наглядно:
Диаграмма направленности адаптивной кольцевой антенной решетки радиопеленгатора и пассивного радиолокатора в логарифмических координатах
Подавление помех и многолучевого распространения позволяет улучшить точность пеленгования. В настоящее время в отечественных РП эти методы не используются. Разработка осталась на уровне начала 2000 годов. В свое время, например «Платан» DF-2000 был хорошим радиопеленгатором, но время неумолимо идет вперед. С учетом новых методов обработки и появления новой компонентной базы он отстал уже лет на 15. Заказчику — Госкорпорации по ОрВД до сих пор поставляются РП пятнадцатилетней давности.
Вот зачем мне спрашивается карточка графического ускорителя в компе? Игры не нужны, а на вычислительный процесс это не влияет. В последнем ноуте оказалась графика Radeon, и так я и воспринимал свой комп — как вместилище неиспользованного ресурса.
Оказалось это не совсем так. Когда понадобилась визуализация сигналов в реальном времени при моделировании на Питоне, выяснилось что пакет matplotlib страшно медленный. Чтобы отрисовать двумерную функцию неопределенности, ждать приходилось минуты. И тут в зоне видимости возник пакет PyQtGraph — замечательная штука для отображения всяческих осциллограмм в нужных точках.
Для тех кто освоился в логике фреймворка Qt, все достаточно просто. Работает эта штука на OpenGL очень шустро. Вот, можете посмотреть сами:
https://youtu.be/2nfbk2CYLPE
Так что рано или поздно все неиспользуемое пригождается. Так же хорошо этот пакет работает и на другом компе с Nvidia. Будьте осторожны с клипом: завораживает и затягивает, можно смотреть беспрерывно, как на горящий огонь. Например на этот:
Обращаю ваше внимание: тандыр авторской работы, изготовлен лично мной, true hand made, любое копирование и воспроизведение запрещено без согласия автора )
Наконец подошло время и для альфа — бета фильтра. Этот фильтр представляет собой более простую реализацию фильтра Калмана и использует два параметра — альфа и бета.
Фактически, эти параметры — весовые коэффициенты к ошибкам прогноза позиционирования и скорости объекта. Альфа и бета остаются постоянными во время работы фильтра и задаются изначально как константы.
Работа 2D фильтра моделировалась в предположении, что движение объекта, или цели в двумерной системе координат является независимым, то есть цель может двигаться совершенно произвольно. Если сделать другое предположение, что движение цели обусловлено, то подозреваю можно улучшить работу фильтра. Некоторые соображения на эту тему я выскажу в конце статьи.
Как работает Альфа — Бета фильтр
Поскольку фильтрация по каждой из осей X, Y происходит независимо, проследим работу фильтра по одной из осей. Для нас важен физический смысл того, что происходит, поскольку математическая реализация фильтра довольно проста и занимает вместе с графикой не больше 100 строк на Питоне.
Итак, цель движется в горизонтальной плоскости, и мы наблюдаем ее положение скажем по оси X. Наблюдение сводится к измерению x координаты цели в каждый дискретный момент времени с равномерными промежутками. Возьмем отдельно взятый отсчет по координате x, допустим пусть это будет 1523 метра. В этот момент альфа — бета фильтр делает предположение, или прогноз, или предсказание о том, где цель очутится в следующий момент времени. Как он это делает?
Предположим, что в тот момент когда мы достигли отсчета 1523м, у нас уже было предсказание, сделанное в предыдущей точке о том, где мы будем находится. Скажем, до текущего отсчета альфа — бета фильтр предсказал нам что наша текущая позиция (а для него тогда это была наша будущая позиция) составит 1535м. Этого для нас, точнее для фильтра достаточно, и теперь внимательно наблюдаем за тем что он будет делать с настоящим и прошлым. Он из этого будет делать будущее.
Во-первых, фильтр тут же узнает о том, насколько он ошибся в своем прошлом прогнозе. Ошибка составит 12м. Фильтр постоянно старается минимизировать эту ошибку и улучшить точность предсказания. Для этого он берет долю альфа от этой ошибки и добавляет эту долю к своему прошлому прогнозу. После добавления это значение становится предсказанием будущего расстояния.
То же самое фильтр проделывает со скоростью движения, только для корректировки скорости используется параметр бета. Прогноз скорости также используется для определения будущего расстояния, поскольку цель проходит его за промежуток между двумя отсчетами.
Предсказание положения цели и траекторная обработка
Обратимся к диаграмме в начале статьи. На ней показаны результаты моделирования движения цели в горизонтальной плоскости. Альфа — бета фильтр использовался в двух каналах X, Y независимо. Измеренные, или реальные положения цели обозначены зелеными треугольниками, а предсказанные фильтром положения цели — красными кружочками. Цель движется с левого верхнего угла к правому нижнему углу диаграммы.
Поскольку в начале работы альфа — бета фильтра предыдущее положение цели неизвестно, ее координаты задаются нулевыми. Соответственно, с самого начала работы фильтра возникает максимальная ошибка, которую обозначает самый большой эллипс на рисунке. Эта ошибка — разница между текущим положением цели и нулевой координатой. На втором шаге фильтр уже существенно улучшает прогноз и ставит его в левую верхнюю часть диаграммы, но ошибка обозначенная эллипсом поменьше по прежнему высока. Видно, как за первые четыре итерации значения ошибок неуклонно снижаются и к началу выхода траектории цели на движение вдоль оси X фильтр предсказывает будущее значение красным кружком практически там же, где будет находиться цель. По мере дальнейшего движения цели и роста крутизны маневра, начинает наблюдаться отставание предсказанного значения от реального. Для того чтобы сократить это отставание, фильтру не хватает работы по второй производной движения — по ускорению. Ускорение учитывается в более сложном альфа — бета — гамма фильтре, где параметр гамма определяет вклад ускорения цели в предсказание положения.
На основе предсказания строится траекторная обработка. Цель такой обработки — логически связать вместе все отсчеты (plots) и назначить их принадлежащему одному треку цели (track). Если каждый плот попадает в предсказанную область, то его с высокой степенью вероятности можно считать принадлежащему этому треку.
Зависимое по координатам движение
Как я заметил выше, X и Y каналы фильтра работают независимо друг от друга. Это значит, что погрешность предсказания, полученная в одном канале, никак не учитывается в другом. Однако, если мы сделаем дополнительные предположения о характере движения цели, мы можем сделать каналы зависимыми и улучшить точность предсказания. Одно из таких предположений может заключаться в следующем. Допустим, цель совершает маневр — движение по кругу. Сразу в двух каналах возникает ошибка предсказания, которая будет полностью скомпенсирована в альфа — бета фильтре только после выхода на прямолинейную траекторию движения. Однако, если мы введем зависимый параметр — длину вектора, ортогональные составляющие которого образованы ошибкой в каждом из каналов, и убедимся в том что длина вектора неизменна (цель сохраняет скорость и только поворачивает), то мы можем обеспечить дополнительную компенсацию ошибки. Такая межканальная компенсация может заключаться в следующем: очевидно, что при повороте цели ошибка в каналах будет расти, но с разными знаками. Вводя дополнительные коэффициенты, можно создать упреждающее предсказание, основанное на предположении о маневре цели.
В общем надо бы попробовать это сделать и сравнить с обычным альфа — бета фильтром )
В верхней части — SDR радиостанции серии 4200 Rohde&Schwarz, в нижней — радиостанции серии 2000 АО Азимут
В эпоху дискретных компонентов — ламп и транзисторов — техника радиоприема была необычайно изощренной. Рефлексные схемы приема, регенеративные и сверхрегенеративные приемники, и наконец классические супергетеродины бились за место под солнцем по показателям чувствительности и избирательности. В любительских конструкциях во главу угла было «поймать то что нужно», и только умудренные опытом профессионалы знали, что этого мало — нужно еще не поймать то, что не нужно. Другими словами, речь идет об ЭМС, электромагнитной совместимости (немного неуклюжий термин, но что поделаешь ).
Если забыть про ЭМС, то стать профи в радиоприеме может любой программист: современные АЦП (аналого-цифровые преобразователи) позволяют оцифровать радио напрямую на скорости до 500 млн. выборок в секунду, где данные перемещаются в домен программного обеспечения и где с ними можно делать все что пожелаешь. Так обычно описывают Software Defined Radio, SDR. Не заставили себя ждать и другие доступные компоненты: SDR радиоприемники размером с флешку и собственно программируемый конструктор радиоблоков GnuRadio.
Снижение порога вхождения на рынок средств радиосвязи стимулировало большой поток желающих работать в этой области. И это здорово, только если не забывать о том, что SDR — это еще не весь радиоприемник. В чем же отличие? В той же самой ЭМС. Поясним на простом примере.
Чего не хватает SDR приемнику
Положим, приемник настроен на рабочую частоту 124.8 МГц, а на частоте 125 МГц присутствует сильная помеха. С точки зрения программиста, проблемы нет: программно определяемый цифровой фильтр формирует полосу пропускания 25 кГц с центром на рабочей частоте и с высокой крутизной характеристики, и помеха остается за бортом. Все просто и очевидно. С точки зрения инженера, проблема есть: при величине помехи +80 дБ замечательный 16-разрядный АЦП по отношению к сигналу будет работать всего на 2 — 3 разряда, а при более мощной помехе полезного сигнала мы не увидим вообще. И никакой программный цифровой фильтр этому не поможет. Есть и другие негативные эффекты, связанные с оцифровкой, такие как алиасинг.
В сухом остатке: профессиональные тракты радиоприема содержат в своем составе, помимо собственно SDR, еще входные узлы (фильтры), которые ограничивают полосу входного сигнала (и не только, но здесь мы фокусируемся только на этом). Если говорить не только о приемниках, но также и о радиопередатчиках, то аналогичные узлы используются в выходном тракте для улучшения спектра излучаемого сигнала.
Профессиональная реализация Software Defined Radio
Теперь смотрим, как SDR радиостанция выглядит в профессиональной реализации. На коллаже сверху — сравнительные фотографии SDR станций двух производителей: компании Rohde&Schwarz и АО Азимут. Те блоки которые поуже это радиоприемники; которые пошире — радиостанции, т.е. содержат в своем составе еще и передатчик. Радиостанции обоих поставщиков предназначены для радиосвязи диспетчер — борт в системе управления воздушным движением и работают практически в одинаковом частотном диапазоне (как минимум в диапазоне 118 — 136 МГц). Они могут использоваться как автономно, так и в самом распространенном варианте — в составе многоканальных автоматизированных приемо — передающих центров, АППЦ.
Бросающаяся в глаза схожесть дизайнов навевает мысли о клонировании немецких устройств в стиле China Copy ) Однако, радиостанция АО Азимут — вполне самостоятельное изделие, полученное в результате заказной разработки. Приемный тракт разрабатывала компания НЭЛС, усилитель мощности, компоненты ВЧ тракта и антенны — АО Сапфир.
Теперь посмотрим, какие фильтры используют поставщики радиостанции для служб управления воздушным движением.
Реализация от Rohde&Schwarz
Автоматический фильтр Rohde&Schwarz, используемый в АППЦ
Rohde&Schwarz комплектует свои АППЦ автоматическими фильтрами FU221. В фирменном описании этих устройств с немецкой точностью и педантичностью перечислено, для чего они нужны. Этот список нужно каждый раз показывать для иллюстрации того, что радиосредства должны разрабатывать инженеры, а не программисты ) Чтобы не возиться с переводом, я даю текст на русском языке, в привычной отечественной терминологии.
Назначение фильтров на входе:
подавление внеполосных продуктов интермодуляции внесением дополнительного затухания для тех сигналов и их гармоник, которые находятся за пределами диапазона приема и которые, вследствие нелинейных характеристик входных цепей могут вызвать продукты интермодуляции, которые создают помехи приему;
подавление продуктов кросс — модуляции третьего порядка уменьшением помехового сигнала, который может воздействовать своей модуляцией на полезный сигнал;
подавление зеркального канала;
предотвращение смешивания (на самом деле — перемножения) помехового сигнала с шумовыми составляющими опорного генератора, что приведет появлению в полосе приема составляющих, уменьшающих чувствительность радиоприемника.
Назначение фильтров на выходе:
подавление шумов опорного генератора и усилителя мощности;
уменьшение гармоник, которые возникают в выходных цепях передатчика;
подавление продуктов интермодуляции третьего и более высоких порядков, вызванных одновременной работой нескольких передатчиков на одну и ту же антенну.
Замечу, что с перечисленным при радиоприеме после АЦП разбираться бесполезно: по спектральным характеристикам помеха становится подобной сигналу и фильтровать ее бессмысленно. Это к с слову о том, что «АЦП оцифрует сигнал с антенны, а там в программе все порешаем» ) При этом неважно, содержат ли входные цепи SDR смеситель с переносом частоты или оцифровка выполняется на несущей — входные компоненты, в том числе и входные цепи АЦП, имеют изначально нелинейную природу, и если вы подали на вход f1, то ждите также других гармоник, а если пришло только f1 и f2, то взаимные продукты этих частот и их гармоник не заставят себя ждать к вам в гости. Хотя как программист вы имеете право думать что имеете дело с линейной моделью )
Вот если бы как в старые добрые времена, ламповые триоды с их линейной характеристикой…
Внимательный читатель заметит, что есть паразитные эффекты которые возникают на системном уровне (уровень многоканального АППЦ). Помимо работы нескольких передатчиков на одну нагрузку, я бы добавил сюда одновременную работу нескольких приемников с одной приемной антенны. Ни за что не догадаетесь, но и в этом случае вас ждет засада: паразитные излучения опорных генераторов и синтезаторов каждого из приемников будут пролезать на входы друг друга и понятное дело качеству радиоприема это совсем не способствует )
Автоматический фильтр: найти мотор в радиоэлектронике Rohde&Schwarz
На примере автоматического фильтра FU221 Rohde&Schwarz я дал картинку того, как в профессиональной области формируются тракты радиоприема и передачи. И здесь дополнительно я хотел упомянуть один важный момент, который проявит себя в сравнении этого принципа построения с построением АППЦ АО Азимут. Этот малозаметный момент кроется в слове «автоматический» в названии фильтра.
Слово «автоматический» здесь совсем не случайно, это совсем не дань моде, как в сокращении АППЦ. Чтобы оценить всю степень последовательности и дотошности инженеров Rohde&Schwarz, вспомним как диспетчеры работают с АППЦ. В самом простом случае, за каждым диспетчером закрепляется своя частота из набора «Старт», «Подход», «Круг» и так далее. Каждой частоте соответствует своя радиостанция АППЦ. Логично предположить, что современная радиостанция позволяет выбрать произвольную частоту из заданной сетки частот в диапазоне; и так оно и есть, это можно сделать как с клавиатуры (как на фото), так и через систему дистанционного управления.
Если возникнет задача поменять частоты у диспетчеров, перестройка радиостанций произойдет достаточно просто по одной команде. Если больше ничего не трогать, то нас ждет неприятное открытие: для входных фильтров новые частоты будут восприняты как помеховые и будут безжалостно подавлены. Кроме этого, помехи на частотах, соответствующих настройке фильтра, будут спокойно проходить на вход радиоприемника. Следовательно, фильтры также нужно перестраивать! Не удержусь от того, чтобы продемонстрировать как это было сделано в ламповых приемниках )
Двухсекционный конденсатор переменной емкости обеспечивает одновременную перестройку входного фильтра и гетеродина
На фото — конденсатор переменной емкости, КПЕ, состоящий из двух секций. Одна секция работает в колебательном контуре гетеродина, читай — синтезатора сигнала. Другая секция работает во входном фильтре приемника. Параметры контуров и фильтра, в частности индуктивности, подобраны таким образом, чтобы вне зависимости от частоты настройки, которая определяется поворотом оси КПЕ, всегда сохранялась постоянная разница между частотой входного фильтра и частотой гетеродина. Собственно эта разница и есть промежуточная частота, которая должна быть постоянной.
Конечно в наше время забавно наблюдать такой девайс, в котором непонятно чего больше — механического или электронного. Однако, если мы изучим FU221 поподробнее, то обнаружим, что этот фильтр содержит аналогичное устройство — механически перестраиваемый конденсатор. Только вместо ручного управления, движение обеспечивает электрический мотор. Очевидно, что механическая перестройка фильтра производится каждый раз, когда меняется частота радиостанции, причем прецезионный механизм обеспечивает точное следование частоты фильтра частоте радиостанции.
Зачем же нужно было создавать устройство таких габаритов, да которое еще требует механической перестройки? Ответ дает еще одна функция, которую я не упомянул: это обеспечение развязки между работающими приемниками и передатчиками. Объемные резонаторы, которые для этого используются, хорошо проявляют себя только в линейных размерах, соизмеримых с длиной волны.
Шкаф для размещения фильтров в реализации АО Азимут
Реализация от АО Азимут:
АППЦ TRS 2000
Отечественный поставщик также комплектует свои АППЦ TRS 2000 аналогичными фильтрами. На сайте и в сети не удалось найти описание этих устройств, но зато точно известно что они есть )
Под эти устройства выделен отдельный шкаф (кто в домике живет?). Точно также как в АППЦ Rohde&Schwarz, эти фильтры перестраиваются механически. На заводе. В процессе регулировки изделия, на заранее известные частоты. В эксплуатации перестройка этих фильтров не предусмотрена. Если вдруг случится необходимость поменять частоту канала АППЦ TRS-2000, эксплуатирующая организация вызывает регулировщика и он по прибытии на место производит перестройку фильтра. Вот такая своеобразная реализация моторчика в отечественном исполнении.
С таким темпом перестройки фильтров (недели для удаленных позиций) возникает сомнение в целесообразности оперативной перестройки частот радиостанции. Известно, что ценой за такую перестройку, с учетом поддержания частоты синтезатора с требуемой точностью во всем диапазоне, является в конечном счете стоимость. С такими фильтрами было бы гораздо проще и органичнее использовать радиостанции на фиксированные частоты и менять частоту в момент перестройки фильтров также механической заменой частотно — зависимых узлов. Отсутствие оперативной перестройки фильтров в АППЦ Азимут ведет к забавной ситуации: сменный инженер через систему дистанционного управления или непосредственно на передней панели радиостанции конечно же может перестроить частоту, но это не означает что АППЦ TRS 2000 на этой частоте будет работать ) И если фильтры не перестраиваются, то зачем в радиостанциях функционал оперативной перестройки частоты? Избыточная стоимость радиостанции, помноженная на количество каналов.
На самом деле, необходимость оперативной перестройки частоты на позициях возникает достаточно редко. Единственная ситуация которая приходит мне в голову — это замена отказавших каналов. Но и в этом случае скорость перестройки является критичной.
Извлекаем уроки
Но не будем сильно придираться к этой реализации: в конце концов, если заказчика — Госкорпорацию по ОрВД все устраивает, значит так тому и быть. Для нас из этих двух коротких примеров, важнее сделать определенные выводы. Для себя я например сделал такие:
SDR приемник — это еще не радио;
отдельные радиокомпоненты, собранные вместе — еще не система;
заказная разработка модулей так и остается набором модулей, если нет комплексного решения / системной интеграции;
даже если заказчик и молчит, это не значит что ему все нравится )
Классический, моностатический радар использует антенную систему с совмещенными функциями передачи зондирующего сигнала и приема отраженного. Бистатический радар содержит разнесенные антенны, одна из которых в передающем тракте Tx передает зондирующий сигнал, другая — в приемном тракте Rx принимает отраженный сигнал.
Частным случаем бистатического радара является пассивный радар (Passive Coherent Locator — PCL), который не имеет собственного передающего тракта, как и передающей антенны, и использует излучающие источники радиоокружения. Выигрышные стороны PCL — это отсутствие собственного излучения, что затрудняет его подавление, простая антенная система, а также малая потребляемая мощность (до 100 Вт) и объем аппаратуры.
По аналогии с PCL, в данном случае рассматривается реализация бистатического радара, использующего непрерывный излучающий сигнал. Данный режим имеет принципиальное отличие от традиционного импульсного режима классических радиолокаторов, когда последние принимают отраженный импульс в моменты времени, когда излучение мощного зондирующего импульса уже завершено. В нашем случае, бистатический радар ведет прием непрерывного отраженного сигнала на фоне мощной помехи, которая представляет собой зондирующий сигнал.
Если для PCL отсутствует возможность влиять на радиоокружение для получения наиболее выигрышного соотношения отраженного сигнала и помехи, для бистатического радара такая возможность имеется. Поэтому бистатический радар, по сравнению с PCL, позволяет выполнить следующую дополнительную оптимизацию:
достичь максимального подавления в прямом направлении на собственный сигнал излучения путем управления диаграммой направленности как в точке излучения, так и в точке приема;
выбрать форму излучающего сигнала для получения лучшей функции неопределенности по доплеровскому сдвигу и времени задержки, приближенную к дельта — функции.
Бистатический радар воздушного базирования
Рассмотрим концепт построения бистатического радара воздушного базирования. Радар предназначен для обнаружения маневрирующих воздушных целей с ожидаемым направлением подлета и определения их координат (Рисунок 1). Оборудование бистатического радара базируется на двух БПЛА (беспилотных летательных аппаратах). На нижем по рисунке БПЛА размещена передающая часть радара Tx, на верхнем БПЛА — приемная часть Rx.
Рисунок 1. Бистатический радар UAV базирования, состоящий из передающей части Tx и приемной части Rx
Для частотного диапазона от 100 до 500 МГц размах крыла БПЛА может составлять от 1 до 3 м.
Как передающая, так и приемная часть радара содержат идентичные направленные антенны, которые имеют нелинейную, в частности логопериодическую структуру, которая может формироваться как физически, так и электронно из антенных элементов. Антенна может быть частью конструкции радиопрозрачного крыла БПЛА.
Максимум диаграммы направленности антенн формируется в сторону ожидания поступления цели со стороны одного борта БПЛА, причем БПЛА двигаются в направлении примерно перпендикулярном направлению движения цели. Расстояние между БПЛА может составлять от 1 до100 км.
При типовом значении подавления в антенной системе в направлении соседнего БПЛА 40 дБ, развязка по помехе между Tx и Rx составит 80 дБ, что ставит бистатический радар в неизмеримо более выигрышную по соотношению сигнал/помеха ситуацию по сравнению с PCL.
Характеристики радара могут быть существенно улучшены применением синтеза апертуры антенны.
Для работы приемной части Rx радара, помимо отраженного сигнала, также используется прямой (помеховый) опорный сигнал, поступающий напрямую от Tx. Принцип приемного тракта основан на построении функции неопределенности принимаемого отраженного сигнала на основе опорного в двумерной плоскости координат: по доплеровской частоте и по временной задержке (Рисунок 2).
Рисунок 2. Функция неопределенности сигнала бистатического радара
На рисунке показана функция неопределенности сигнала бистатического радара. Функция содержит главный пик на нулевой отметке (Delay = 0 us, Doppler = 0 Hz), который соответствует прямому распространению сигнала от Tx к Rx (помеха). Другие пики соответствуют движущимся целям (эхо — сигнал), в частности отметка цели с доплеровским сдвигом Doppler = 100 Hz. Сопоставимые уровни помехи и эхо-сигналов обясняются подавлением помехового сигнала в тракте обработки.
Для корректного определения координат целей требуется определение координат каждого БПЛА в течение времени обнаружения и сопровождения цели. Данное время составляет около 1 с, что позволяет достичь очень узкой полосы доплеровского фильтра в 1 Гц. Данные значения ограничены стабильностью несущей частоты передатчика и взаимной флуктуацией расстояния между БПЛА в течение 1 секунды.
Только для обнаружения цели, без определения координат знания координат БПЛА и расстояния между ними не требуется. В этом случае информационным параметром, который будет передавать бистатический радар, будет оповещение типа «приближение объекта с радиальной скоростью 900 км/ч на приблизительном удалении 100 км».
Энергетический бюджет радиолинии
Для оценки уровней принимаемого полезного отраженного эхо — сигнала и помех воспользуемся уравнением бистатического радара, приведенного в статье «Bistatic and Multistatic Radar Systems»:
Для типовых, практически значимых значений параметров, в частности: излучаемая мощность передатчика Tx 10 кВт, расстояние до цели от приемника Rx и передатчика Tx 50 км и 80 км соответственно, ЭПР 10 кв.м, частота 100 МГц искомые значения Pr для сигнала и помех сведены в таблицу (таблица приведена в цитируемом источнике).
Тип сигнала
Уровень на входе Rx
Относительный уровень
Тепловой шум
-121,0 dBm
+8,8 dB
Прямой сигнал Tx->Rx
-37,7 dBm
+92,1 dB
Электромагнитный шум
-90,0 dBm
+30 – 40 dB
Переотражения
— 99,8 dBm
+39,8 dB
Эхо — сигнал цели
-129,8 dBm
0,0 dB
Как видно из таблицы, при отражении от цели мощность эхо — сигнала в точке приема составит -129,8 dBm, что для типовой нагрузки 50 Ом составит 0,1 мкВ. Электромагнитный и тепловой шум не являются проблемой, поскольку в процессе обработки в доплеровском фильтре происходит сжатие сигнала до полосы 1 Гц, что дает увеличение отношения сигнал/шум на величину около 50 дБ. Вследствие этого отношение сигнал/шум для данных видов помех составит около +10 дБ.
Гораздо более серьезным и интенсивным источником помехи является прямой сигнал, идущий напрямую от передатчика Tx к приемнику Rx минуя цель. Точно также как и для переотражений, которые являются аналогичной помехой, единственными различающими параметрами для эхо — сигнала являются приобретенные доплеровский частотный сдвиг и запаздывание сигнала по времени.
Как следует из таблицы, а также из практических приложений, превышение помехи над сигналом может составлять 90 — 100 дБ. Задача подавления прямого сигнала для PCL решается комбинацией следующих методов:
формирование минимума диаграммы антенной системы в направлении помехи: 20 — 30 дБ;
использование в тракте обработки адаптивного фильтра, который настроен на подавление сигнала базовой станции, а также более сложных реализаций лестичных адаптивнных фильтров с ортогонализацией Грамма — Шмидта: 30 — 40 дБ;
различение на уровне боковых лепестков функции неопределенности: 30 — 40 дБ.
Бистатический радар с собственным передатчиком Tx, как отмечалось выше, обладает еще большими возможностями по подавлению собственного прямого сигнала. За счет управления диаграммами приемника и передатчика достигается дополнительное подавление +40 дБ; использование собственного шумоподобного сигнала существенно улучшает работу адаптивного фильтра и уменьшает уровень боковых лепестков функции неопределенности, что дает дополнительный выигрыш в 10 — 20 дБ.
Таким образом, бистатический радар является реализуемым по отношению сигнал/шум для приведенных исходных данных, и существуют возможности по улучшению этого параметра для более сложных условий окружения, что требует дальнейших исследований.
Ключевые компоненты и программное обеспечение
Реализация бистатического радара требует разработки ключевых компонентов с высокой плотностью интеллектуальной составляющей. В частности, это касается радиотракта с широким динамическим диапазоном и высокой когерентностью в прямом и отраженном канале, а также быстродействующих средств обработки сигнала на базе DSP.
Такие компоненты представлены на рисунке 3.
Рисунок 3. Радиомодуль, содержащий прямой и отраженный канал, и DSP процессоры
Тракт обработки, представленный на фото, работает с различными сигналами. В частности, на рисунке 4 представлено сечение функции неопределенности сигнала по частоте для широкобазового сигнала с линейной частотной модуляцией (ЛЧМ):
Рисунок 4. Сечение функции неопределенности ЛЧМ сигнала по частоте
Появились первые результаты моего проекта по разработке аэродромного обзорного радиолокатора (АОРЛ) S- диапазона. Главный партнер в проекте — небольшая итальянская компания Argos — состоит из сотрудников, которые долгое время работали в Selex по радиолокационной тематике. На первом этапе проекта партнерами разработана конструкторская документация на антенну и функционально проработана структура изделия.
Пьедестал антенной системы
За то время, когда прорабатывалась структура локатора и велась разработка, я посетил Рим раз наверное десять. Визиты были такими насыщенными, что мне предложили выделить кабинет в римском офисе для работы ) Работали быстро и эффективно. В команде проекта с итальянской стороны — три специалиста в предметной области, два программиста. Для антенны — один радиоинженер и один конструктор, которые кстати после Selex’а работают в другой маленькой компании — разрабатывают, изготавливают и продают спутниковые тарелочки.
Чувствуется, что у старой гвардии, которая ушла в малый бизнес, по прежнему есть порох в пороховницах. В который раз убеждаешься в том, что небольшие профессиональные команды на голову эффективнее традиционных НИИ — монстров, состоящих из кучи отделов и секторов. Small is beautiful, как говорит Гринпис.
Справа вместо фото представлен рисунок другого компонента антенной системы — пьедестала, в котором располагается вращпереход, редуктор с двигателями и устройство юстировки антенны. Предполагалось, что этот чисто механический узел (за исключением вращперехода) будет изготавливаться на нашей территории. Однако, все идет к тому, что пьедестал будет делать маленький симпатичный заводик в Неаполе, который заточен под радиолокационные СВЧ компоненты.
Элементы волноводного тракта, разрабатываемые и изготавливаемые на экспериментальном заводе в Неаполе
Гибкий волновод требует высокой точности изготовления и специальной технологии
Антенну решили делать со смещенным фокусом, чтобы не затенять рефлектор и как следствие получить дополнительный выигрыш в энергетической эффективности. Как результат — более жесткие требования к изготовлению рефлектора, который опирается на горизонтальную несущую раму и вертикальные ребра жесткости. Антенна формирует два луча: основной (Main Beam Tx/Rx) и вспомогательный (Aux Beam Rx) и содержит поляризатор для работы в условиях дождя. Рабочий диапазон антенны: 2700 — 2900 МГц.
Argos разрабатывал конструкторскую документацию на антенну в САПР Pro Engineer. По условиям контракта, документация должна быть в формате Solid Works. В процессе экспорта выяснилось, что при преобразовании форматов теряется информация о преобразованиях которым надо подвергнуть деталь (например, гибка). В Solid Works согнутый лист выглядел как отлитый 🙁 Поэтому фактически всю КД пришлось перерабатывать уже в Солиде. Вообще, надо отметить, что Solid Works — это какое-то специфическое отечественное применение; западные разработчики его практически не используют.
Изготовление рефлектора, главной деталью которого является металлическая сетка, требует высокой точности и соответственно специальной технологической остнастки. Для проверки характеристик антенны необходимо разворачивать целый полигон и смотреть диаграмму в дальнем поле. Когда мы выбирали площадку для полигона, итальянцы очень придирчиво рассматривали особенности местности. Вплоть до того, какие кустики растут на подстилающей поверхности 🙂
Специфичные компоненты антенны, такие как волноводы и облучатель изготавливает маленький симпатичный экспериментальный заводик в Неаполе. Неаполь запомнился хаотичным автомобильным движением и автомобильными гудками, совсем как у нас на кавказском юге. Видимо, южане везде одинаковы ) На фото которое я сделал с горки видно, как Неаполь огибает море. Картинка сильно напоминает амфитеатр. К слову, VIP — места этого амфитеатра, а именно квартиры окна которых обращены к морю, стоят существенно дороже 😎
В процессе экспериментов с малышкой RPi на просторах инета мне попался скрипт motion.py, который ведет непрерывную видеозапись с камеры и сохраняет в выходном видео только движущиеся изображения. Это очень удобно как с точки зрения экономии места на флешке, так и для просмотра: очень быстро надоест проматывать запись с автомобильной стоянки за всю ночь, когда нужна буквально пара — другая информационных кадров. Если быть точным, то это не совсем видео в современном смысле, а скорее видео в классическом: RPi беспрерывно делает фотографии, которые склеиваются в формат motion jpeg.
Детектор движения у нас есть, сделаем полноценную систему видеонаблюдения того, что происходит за окном. И не только за окном: ночью в магазине, офисе, на складе, дома когда вы в отъезде — применений просто огромное количество. Устроено все очень просто: Raspberry Pi с видеокамерой и WiFi донглом для доступа в интернет, и пара скриптов: motion.py и conv.sh.
Скрипт motion.py: детектор движения
Питоновский скрипт motion.py я немного доработал под систему:
Python
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
#!/usr/bin/python
# /home/pi/bin/motion.py
# original script by brainflakes, improved by pageauc, peewee2 and Kesthal
# type "sudo apt-get install python-imaging-tk" in an terminal window to do this
importStringIO
importsubprocess
importos
importtime
fromdatetimeimportdatetime
fromPIL importImage
# Motion detection settings:
# Threshold - how much a pixel has to change by to be marked as "changed"
# Sensitivity - how many changed pixels before capturing an image, needs to be higher if noisy view
# ForceCapture - whether to force an image to be captured every forceCaptureTime seconds, values True or False
# filepath - location of folder to save photos
# filenamePrefix - string that prefixes the file name for easier identification of files.
# diskSpaceToReserve - Delete oldest images to avoid filling disk. How much byte to keep free on disk.
# cameraSettings - "" = no extra settings; "-hf" = Set horizontal flip of image; "-vf" = Set vertical flip; "-hf -vf" = both horizontal and vertical flip
threshold=10
sensitivity=20
forceCapture=True
forceCaptureTime=60*60# Once an hour
# filepath = "/home/pi/picam"
filepath="/home/pi/cache"
filenamePrefix="img"
diskSpaceToReserve=40*1024*1024# Keep 40 mb free on disk
cameraSettings=""
# settings of the photos to save
saveWidth=1296
saveHeight=972
saveQuality=15# Set jpeg quality (0 to 100)
# Test-Image settings
testWidth=100
testHeight=75
# this is the default setting, if the whole image should be scanned for changed pixel
testAreaCount=1
testBorders=[[[1,testWidth],[1,testHeight]]]# [ [[start pixel on left side,end pixel on right side],[start pixel on top side,stop pixel on bottom side]] ]
# testBorders are NOT zero-based, the first pixel is 1 and the last pixel is testWith or testHeight
# with "testBorders", you can define areas, where the script should scan for changed pixel
# for example, if your picture looks like this:
#
# ....XXXX
# ........
# ........
#
# "." is a street or a house, "X" are trees which move arround like crazy when the wind is blowing
# because of the wind in the trees, there will be taken photos all the time. to prevent this, your setting might look like this:
# testAreaCount = 2
# testBorders = [ [[1,50],[1,75]], [[51,100],[26,75]] ] # area y=1 to 25 not scanned in x=51 to 100
forzinxrange(0,testAreaCount):# = xrange(0,1) with default-values = z will only have the value of 0 = only one scan-area = whole picture
forxinxrange(testBorders[z][0][0]-1,testBorders[z][0][1]):# = xrange(0,100) with default-values
foryinxrange(testBorders[z][1][0]-1,testBorders[z][1][1]):# = xrange(0,75) with default-values; testBorders are NOT zero-based, buffer1[x,y] are zero-based (0,0 is top left of image, testWidth-1,testHeight-1 is botton right)
Скрипт работает так. В цикле через утилиту raspistill, которая работает с камерой, непрерывно делаются фотоснимки, после чего сравнивается каждый текущий и предыдущий кадр. Сравнение выполняется в модуле Imagemagic. Если кадры различаются начиная с некоторого уровня порога, значит в зоне наблюдения присутствует движение и снимки начинают записываться в директорию cache SD — карты RPi в формате jpeg, причем название файла образовано от отметки текущего времени.
Скрипт motion.py стартует с загрузкой RPi и работает непрерывно. Частота появления файлов в директории cache будет зависеть от движения в кадре: если картинка сильно меняется то снимки будут появляться постоянно; если движения нет — то новых файлов не будет.
Далее, с определенной периодичностью, нужно смонтировать полученные кадры в видеоролик и отправить его в сеть. Эту задачу выполняет скрипт conv.sh.
Скрипт conv.sh: формирование видеороликов и выкладка их на файловый сервер
Задача второго скрипта — conv.sh — периодически запускать на обработку полученные снимки. Для этого он выполняет следующую последовательность действий:
собирает снимки, накопленные к данному времени скриптом motion.py;
вставляет в снимки титры с указанием даты и времени снимка;
переименовывает снимки для того, чтобы к ним можно было применить групповую операцию;
преобразует последовательность снимков в видеоролик;
стирает обработанные снимки;
выкладывает видеоролик на файловый сервер типа box.net или Яндекс.
Скрипт conv.sh содержит следующую последовательность команд:
Первым делом, составляется список всех снимков накопленных в директории cache и каждый снимок из списка пропускается через утилиту convert, которая вставляет в снимок отметку времени. Сразу замечу, что это накладная процедура, поэтому лучше обходиться без нее. В этом случае надо просто закомментировать вызов convert и раскомментировать вызов команды mv. И в том и другом случае файл снимка получает новое имя — возрастающее число. Это необходимо сделать для того, чтобы конвертор ffmpeg мог выполнить групповую операцию над файлами с маской %06d.jpg.
Поскольку ffmpeg из пакетов RPi имеет ограниченную функциональность, я компилировал утилиту из исходников и разместил полную версию в каталоге /home/pi/bin. После монтажа снимков в видеоролик, который будет временно располагаться в директории /tmp/*.avi, другая команда — curl отправляет ролик на файловый сервер. В нашем случае это сервер dav.box.com.
После этого подчищаются обработанные файлы и сам отправленный ролик.
Скрипт conv.sh запускается каждые 15 минут по cron. За это время не должно быть переполнения файловой системы; ну и сам размер клипа должен быть вменяемым для передачи через сеть.
В результате, мы можем залогиниться на box.net и посмотреть нарезку роликов, которые будут появляться каждые 15 минут. Само собой, это время можно изменить в настройках cron. Косвенно, размер каждого ролика определяет активность в кадре: чем больше движений тем длиннее ролик, а если ничего не происходит и картинка статичная, длина ролика будет вплоть до нулевой.
Во всей этой истории есть слабое место, а именно сборка клипа с помощью ffmpeg. Поскольку он не использует аппаратное ускорение, то загрузка ARM идет под 100%, и это весьма печально; хуже всего что при интенсивном движении в кадре и как следствие большом количестве файлов сборка может не завершиться к тому моменту когда conv.sh должен быть запущен в очередной раз. С этим приходилось мириться, пока я не нашел это решение, заточенное под GPU.
RPi style: как малышка Raspberry Pi уделает Intel при монтаже клипа
Как мы знаем, несмотря на свою миниатюрность, RPi обладает мощным ресурсом: это GPU, Graphic Processig Unit, в котором реализована аппаратная поддержка такого известного кодека как H.264. На программном уровне с GPU работает программная прослойка OpenMax, а в свою очередь реализация мультимедийной утилиты gstreamer на RPi использует OpenMax и соответственно ресурсы GPU.
Не буду описывать установку дополнительных пакетов, как это описано по ссылке: я поставил их без проблем. Что нам нужно сделать — это теперь выкинуть из скрипта conv.sh утилиту ffmpeg и заменить ее следующей строкой:
Теперь все стало любо — дорого: выходной ролик вместо mjpeg получил навороченный формат h.264, а время монтажа сократилось с минут до нескольких секунд! Вот она, мощь рациональной архитектуры по сравнению с много-много CPU GHz.
В сети можно найти сравнительные тесты обработки видео в Raspberry Pi с использованием аппаратного ускорения и с использованием платформы Intel. Посмотрите, не поленитесь, результаты вас удивят.
Просматривая записи в архиве интернета, который кстати почему-то заблокирован, нашел старые снимки своего сайта (смотреть через прокси — напрямую ссылка работать не будет). Мои старые записи касаются разработки АРП «Платан», которая была начата в конце 90-х. Было интересно сравнить то что задумано и что появилось в результате, после чего я решил написать этот пост. Свои неоконченные записи я даю здесь без изменений, как они были сделаны в 2000 — 2006гг, под названием «Разработка драйвера реального времени FreeBSD для радиопеленгатора Платан». Под ссылкой достаточно подробные тексты про железо, которое предполагалось использовать, и исходный текст драйвера, с которого начался весь радиопеленгатор.
Копировать и использовать текст драйвера, а также приведенный здесь софт без разрешения автора, то есть меня, не допускается 😎
Запуск проекта: прелюдия
В бурные 90-е я покинул свою контору, которая окончательно перестала платить зарплату сотрудникам, и переехал в Москву. Позади была моя радиопеленгационная история, которая началась аж с 1983 года: запуск в серийное производство двухканального аналогового АРП-80 «Тополь-К», участие в разработке и испытаниях ряда цифровых унифицированных радиопеленгаторов АРП-85 «Пихта», дальше уже свои собственные разработки в ипостаси главного конструктора: судовой радиопеленгатор «Пихта-2С» и встраиваемый в мобильный радиолокационный комплекс «Комар-2» радиопеленгатор. С судовым РП была забавная история: для того чтобы продвинуть свою разработку, пришлось буквально «захватить» судно на котором шли испытания старого изделия. РП «Комар-2» был предназначен для Tesla Pardubice, Чехия, которая делала радиолокатор для наших ВВС, и это тоже была длинная интересная история, которую я как-нибудь напишу.
Поскольку в новом царстве-государстве ex-USSR инженерная подготовка никому не была нужна, я прошел полную профессиональную перезагрузку и начал работать в области телекоммуникаций и разработки ПО для этого рынка. За это время в конторе сменилось руководство, она как-то оправилась, появились деньги. Но не было самого главного — продукта, который надо было выдвинуть на рынок. Контору уже начинали душить приподнявшиеся конкуренты: челябинск выкатил новый радиопеленгатор АРП-90; сделан он был страшненько, на базе стандарной PC, но в аэродромной эксплуатации и того не было.
Чтобы было ясно, какова цена вопроса, забегая вперед, скажу: к настоящему времени в российские аэропорты поставлено более 150 РП «Платан», который был разработан в этом проекте. А это ни много ни мало 8 — 10 млн. рублей за одну позицию.
Новое руководство отыскало меня в Москве в 2000 году, когда я работал в компании «ИнкомА«. Тогда она располагалась на пятачке на Ленинградском проспекте, по иронии судьбы связанным с авиацией: бывший аэровокзал, Госкорпорация. Компенсация в Инкоме была отличная: пресловутая «штука» в месяц. Это не нынешние легковесные 1000 долларов — в те времена, «штука» в месяц означала покупку московской квартиры в перспективе полутора — двух лет. Отмечаю эти подробности для того, чтобы упомянуть Инкому добрым словом: несколько месяцев спустя после ухода из нее, когда я уже работал над новым пеленгатором, мне позвонили с Инкомы. Звонок озадачил: с работы я уволился, дела сдал, человека на свое место подготовил. Что же вы думаете? Меня приглашают подъехать и забрать бонус по результатам прошедшего года.
Честно говоря, я был обескуражен. Выплатить существенную сумму своему бывшему сотруднику когда он ни о чем не просит и не требует — для этого нужна сила высшего порядка )
Я иногда рассказываю эту историю как пример того, что не все так плохо в бизнесе по отношению к персоналу. Есть и приличные фирмы. ИнкомА, шлю тебе респект, вспоминая далекий 2000 год!
Риски проекта
На проект я согласился не сразу. Хотя предложенные условия и оплата были вменяемыми, мне надо было оценить вероятность его успешного завершения. Сроки поджимали, права на ошибку не было. Самая главная проблема была это технологический разрыв, который произошел за это время в конторе. Последние разработанные РП работали на отечественной PDP-подобной платформе 588/1806 — когда еще получалось резать импортные кристаллы и их копировать. Операционной системы не было, разработка велась на ассемблере. Несмотря на то, что возможности предыдущих технологий были ограниченными (попробуйте сейчас написать что нибудь для 64К оперативной памяти?), тем не менее это были системы жесткого реального времени. Предстояло решить вопрос с аппаратной платформой и операционной системой.
Для того чтобы снять риски я взял маленький таймаут для того, чтобы поиграться с ОС. Windows Embedded отпала сразу, по понятным причинам. Был серьезный кандидат на разработку — QNX. Изучив возможности этой системы, в частности возможности по realtime, я понял что накладные расходы будут слишком велики. Сюда входят компетенция программистов и возможности по обновлению системы. Кроме этого, меня не интересовал программисткий реалтайм: в отличие от программистов, инженеры работают в двух доменах — аппаратном и программном, и всегда можно сделать жизнь программы проще, задействуя возможности в аппаратном домене. Много времени спустя, к командировке у коллег, которые решили подсесть на QNX, поинтересовался как идут дела. Поделились печалью: купили диск емкостью в два раза больше, пытаемся подсоединить. Драйвер в QNX старый, система обновляется медленно… Портировать бы весь софт на Linux, но затраты по времени кусаются: год — полтора. Такова цена своевременно сделанных решений или отсутствия оных.
С таким подходом, мне подошла бы unix- подобная ОС. Linux тогда не был так хорошо распространен как сейчас, и выбор пал на FreeBSD. Был еще один аргумент, который в то время никто особо не принимал во внимание, но который сейчас очень дорого стоит: это тип лицензии. Во FreeBSD очень либеральная лицензия, она не обязывает раскрывать код проприетарной разработки. В Linux работает GPL, и грамотные юристы рано или поздно заставят вас вытащить все исходные тексты изделия на свет божий. А это уже — потеря продукта и технологий.
Чтобы подстраховаться, я написал разработчикам FreeBSD и рассказал про свой проект и спросил каковы риски раскрытия кода. Ответ был: делай что хочешь, хоть закрывай, хоть открывай, условий лицензии BSD это не нарушает )
Что и говорить, с драйвером пришлось повозиться. Загружаемые модули тогда были не в моде; драйвер жестко компилился в ядро, и не счесть многочисленных часов, потраченных на перезагрузки после очередного kernel panic. Одновременно пришлось хорошо продрать ядро, вплоть до функций инициализации контроллера прерываний. Но в результате фря не подвела. Первый опытный образец АРП «Платан», который был поставлен в Пулково в 2001г., как работал так и работает до сих пор. Не знаю, какой там сейчас uptime, но сообщений и претензий по сбоям не было. Так что будете пролетать над Питером и заприметите внизу, рядом с ВПП радиопеленгатор с характерной кольцевой решеткой — знайте, что там крутится FreeBSD 4.3.
С аппаратной платформой все было проще. Появились одноплатные встраиваемые компьютеры (Single Board Computers — SBC) на основе малопотребляющего процессора Geode с тактовой частотой порядка 200 МГц. По той же ссылке есть обзор аппаратных платформ, который я делал для того чтобы определить какую платформу использовать в пеленгационном проекте. SBC который был применен в изделии обладал всеми характерными признаками встраиваемых систем: отсутствие движущихся частей, таких как вентилятор охлаждения и жесткий диск, присутствие flash памяти на борту: микросхемы Disk-On-Chip емкостью целых 256 МБ (непозволительная роскошь, не правда ли?), шина совместимая с периферией, сторожевой таймер, который обнаруживает потенциальное «зависание» SBC и перезагружает его.
Single Board Computer (SBC) в конструктиве Mini-PC. Применяется в АРП «Платан» для обработки сигнала и управления
Все было замечательно, осталось решить последний принципиальный вопрос: где будет граница между программной и аппаратной частью в тракте обработки сигнала и как следствие каким будет драйвер этой аппаратной части. Очевидно, что это должно быть результатом компромисса: чем проще железо тем жестче требования к realtime драйвера, и наоборот: ставя FreeBSD в комфортные условия, мы увеличиваем объем предварительной обработки которая связана с аппаратными затратами.
Тракт обработки АРП «Платан»: программно-аппаратная реализация
В отечественной фазовой радиопеленгации с использованием квазидоплеровского эффекта сложилась следующая схема построения антенно-приемного тракта. Она идет от «бедности», а точнее от того, что для того чтобы принять сигнал как с коммутируемого кольцевого, так и с центрального вибратора одновременно — нужно два одновременных работающих радиоприемника. Поскольку в то время разработка такого приемника — это отдельная сложная ОКР, использовались покупные. Покупные приемники были как правило бортовыми: это приемники «Полет» в АРП-80 «Тополь-К», Р-872 в гражданском варианте АРП-85 «Пихта» и радиостанция Р-863 в военном варианте АРП-85, наконец «Юрок» который и использовался в Платане. Раз бортовые — значит дорогие. Поэтому, для снижения объема аппаратуры и стоимости зачинателями отечественной радиопеленгации был придуман способ, как запихнуть оба сигнала в один приемный тракт.
Суть этого трюка заключается в следующем. В антенной системе сигнал кольцевого вибратора сносится по спектру на частоту скажем 5550 Гц с помощью однополосного модулятора. После этого он суммируется с сигналом кольцевого вибратора и дальше идет на вход приемника. В приемном тракте возникают биения частоты модуляции 5550 Гц, которые выделяются амплитудным детектором. В детектированном сигнале сохраняются все фазовые соотношения между кольцевыми и центральным вибратором, что является принципиальным для фазовой угломерной системы.
Как видно, частота несущей не так велика: сейчас никого не удивишь на скорости обработки 500 млн. выборок в секунду. Однако использование подобного рода компонентов в то время было избыточным. Оказалось достаточным определить амплитудно — фазовые соотношения детектированного сигнала для каждого вибратора, и темп поступления данных с частотой коммутации вибраторов оказался подъемным для разрабатываемого драйвера.
Секция обработки АРП «Платан». Справа — плата демодулятора и опорных сигналов формата 6U
Так и сложился тракт обработки: антенна с коммутатором и однополосным модулятором, стандартный приемник, плата демодулятора и SBC на котором крутится FreeBSD и драйвер. Временная диаграмма показала, что в среднем каждую половину электронного оборота антенны система проводит в драйвере, и только вторая половина отдается в userspace. Однако этого оказалось достаточным для работы приложений.
Тестирование системы с драйвером привело к отказу от первоначальной идеи запуска потока частотного канала в момент появления сигнала: выяснилось, что накладные расходы на создание процесса настолько велики, что драйверу в этот момент не хватает времени на обработку сигнала. Поэтому была принята схема загрузки приложений каналов в момент старта системы по количеству частотных каналов.
Люди и технологии
Если с концептуальными подходами к архитектуре нового радиопеленгатора стало более-менее понятно, следующей проблемой были люди. На период проекта, чтобы средние и не очень начальники не сильно мешали разработке, меня сделали главным конструктором предприятия и дали отдельный кабинет. Это было забавно, но так и есть: для многих отношения в первую очередь определяются твоим кабинетом и должностью, а не тем что ты знаешь и умеешь )
В конторе появилось нормальное измерительное оборудование, новая молодежь, даже стали платить командировочные за выезд на полигон и вовремя предоставлять транспорт: это вообще меня удивило, вспоминая как в стародавние времена можно было часами ждать на жаре на оградой предприятия, когда водитель соизволит подать свой ПАЗик. Но не было главного: носителей опыта, знаний; специалистов в предметной области, или инженеров старой школы. Не осталось ни одного человека, которые хорошо ориентировались в предметной области радиопеленгации: продвинутые спецы в 90-х разъехались кто по другим городам и странам, а те кто остались, стали начальниками. Поэтому я главным образом работал с молодой порослью. Первым делом надо было привить то, что я называю unix- культурой, которая напрочь отсутствовала. Это подход и стиль мышления основанный на восприятии компьютера как лаборатории с многочисленным инструментарием разработчика, а не как машины конечного пользователя. И само собой, использование FreeBSD как среды разработки. Молодежь с удовольствием открыла для себя новую поляну, в которой вся FreeBSD со своими исходными текстами как на ладони: готовый конструктор, настольный toolbox инженера.
Следующее нововведение в конторе было сделать сложнее, поскольку оно затрагивало ее организационную структуру. В НИИ и КБ советских времен сложилась устойчивая ассоциация отдела/сектора разработки с блоком. Смотришь на опытный образец и сразу видишь что блоки разрабатывали разные подразделения: светодиоды разные, кнопки тоже, гравировка отличается. Блоки были напичканы железом под завязку: так, в АРП-85 «Пихта» блоки были выполнены в УБНК, блок обработки содержал 7 плат, аппаратура передачи информации по 10 км линии на КДП, фактически — модем, состояла из двух блоков, один из которых состоял из 15 плат, генератор опорных и управляющих напряжений — тоже блок, в котором было около 10 плат. И соответственно секторы, отделы, комнаты: отдел электропитания, отдел передачи данных и еще много-много отделов.
Эта рамка — что изделие делают отделы, которые делают блоки — прочно сидела в мозгах и выносить ее пришлось силовым способом, ломая оргструктуру конторы. Новый радиопеленгатор разрабатывался в соответствии с COTS подходом, и не было никакого смысла заниматься разработкой тех узлов, которые в массовом порядке присутствуют на рынке. Решения, которые были приняты в самом начале, определили облик радиопеленгатора и быструю разработку: опытный образец был разработан и изготовлен в течение буквально одного года. Вот какие были эти решения:
источники электропитания НЕ разрабатываем а покупаем. Все покупаем — и те которые стоят на платах и сетевые. До этого разработкой занимался целый отдел, и сколько я себя помню в конторе — они постоянно отказывали;
используем SBC с операционной системой для обработки, управления и контроля: есть многозадачность, IP стек для передачи данных, резко упрощается сопровождение;
портирование FreeBSD до системы реального времени, разработка в дружественной среде на языке Си;
вместо двух блоков передачи данных — одна плата 100% резервированного модема с покупным QPSK чипом, на четырехпроводной 10-км линии АРП — КДП поднимаем IP/PPP, никаких проблем с маршрутизацией и созданием дополнительных точек управления и контроля;
индикаторы пеленга, которых может быть подключено множество — до 20, работают по одной линии RS-485, а не подключаются отдельными кабелями к блоку сопряжения как это было раньше;
сами индикаторы пеленга приобрели приличный вид, появилось новое отображение уровня и точности сигнала;
существенна упрощена аппаратная часть тракта: вместо одного блока одна плата демодулятора, вместо другого блока опорных напряжений — плата формирователя опорных сигналов;
из-за резкого сокращения количества блоков и соединений ушли габаритные разъемы и кабели, монтаж выполняется плоскими кабелями;
скользящий резерв радиоприемных устройств;
создана автоматизированная система контроля, в том числе антенной системы.
Гладко выбритый главный конструктор АРП «Платан» со своим детищем. Тестирование опытного образца на полигоне в аэропорту Махачкалы
На фото я собственной персоной в дутом плаще, который остался у меня с 90-х годов. В эти времена они были дешевы и очень популярны; теперь старой одежде нашлось применение. Видно что единственный шкаф пеленгатора практически пустой: в верхней части 6U секция обработки, место которое занимает дисплей отдано под радиоприемные тракты в зависимости от канальности РП. Внизу секция 4-х радиоприемников с управлением, еще ниже — усилитель — распределитель.
Индикатор пеленга. Справа внизу – характерная бабочка, которая осталась только у первых партий Платанов: отображает уровень сигнала и его точность
Единственное, что руководство конторы помешало мне сделать — это поставить в аппаратной топчан. В предыдущих изделиях он также использовался для хранения ЗИПа. Топчан очень удобен — на нем можно сидеть, лежать, его можно использовать как подручную полку. Топчану были рады все — эксплуатация всех пеленгаторов, и особенно военные, не говоря уже о нас, которые проводили в аппаратных кучу времени. И вот теперь топчана нет: красиво, эстетично, но нефункционально. Печалька.
Параллельно шла разработка программного обеспечения. Весь софт, а также функциональную проработку критичных узлов, я взял на себя. Второстепенные функции были переданы помощнице, которую я взял сразу же как занял свой кабинет. Для конторы это все было необычно: компьютеры на столах считались признаком хайтека и престижа, а не печатными машинками, за которыми самому работать совсем не обязательно. К слову, многие организационные проблемы таких контор есть следствие низкой оплаты квалифицированного персонала. Если бы зарплата инженера была от 5000 Евро, то начальники десять раз подумали бы, стоит ли тратить дорогое инженерное время на написание разных бумаг, служебных записок и отчетности.
Вибратор антенной системы АРП с измерительным стендом. Жесткие требования к фазовой идентичности вибраторов удалось снизить за счет использования корректирующего ПО в тракте обработки
Благодаря хорошей архитектурной проработке изделия, а также тому, что в Платане было задействовано минимальное количество собственных разрабатываемых компонентов, удалось существенно снизить риски. В опытном образце практически не было ни одной серьезной ошибки, которая потребовала доработки. Первоначально, как и планировалось, мы изготавливали один образец. Видимо, руководство так поверило в то, что результат скоро будет, что по ходу разработки заказало сразу три образца. Потом их стало пять, и к концу проекта, с колес мы настраивали уже 10 новых Платанов, под которые были контракты. Все эти образцы пошли к заказчикам практически без доработок. В результате получилось, что критичной по срокам была не собственно разработка и даже изготовление, а приобретение комплектации. Со многими поставщиками мы начинали работать впервые, не были отработаны механизмы входной приемки, но потом все наладилось.
Предварительные испытания проводились в аэропорту Махачкалы, а опытный образец был поставлен в Пулково, куда мы выехали целой бригадой на проведение испытаний. Новый пеленгатор не подвел — все испытания завершились успешно. Моя пеленгаторная команда набралась опыта, заматерела, и можно было спокойно оставлять на нее поддержку и сопровождение Платана.
Развертывание антенной системы опытного образца АРП Платан в аэропорту Махачкалы, 2003 год. Помогает лично руководство базы ЭРТОС )
Как я уже говорил, сейчас Платанов по России поставлено более 150 комплектов. Не так давно Пулково купило второй радиопеленгатор. Как не грезили теоретики от гражданской авиации, что РП скоро сойдут на нет и им нет места в управлении воздушным движением, а пеленгатор особенно любим диспетчерами и никуда уходить не собирается. Любит диспетчер знать совершенно точно, с каким бортом говорит: а ведь метка вторичного локатора на экране может подвести. А спокойствие диспетчера дорого стоит.
Опытный образец АРП Платан, развернутый в Пулково в 2003г. Высокая посадка аппаратной — защита от снега. Аппаратная располагается непосредственно под антенной
После успешного окончания проекта, я покинул контору и вернулся в Москву. Много времени спустя я наконец осознал, чем мне был интересен этот проект: в нем я сделал пеленгатор именно таким, каким хотел его видеть. Естественно, на то время. А дальше, просто пребывать в конторе в статусной должности было уже неинтересно — и я окунулся в другой поток событий.
В основу программного обеспечения АРП «Платан» легли исследования, которые я обобщил в своей диссертации. Точнее будет сказать что в диссертации просто закреплен теоретический и практический материал, который я наработал занимаясь пеленгационной тематикой. Приведенный здесь софт тесно связана с математической задачей принятия решения о наличии сигнала в шуме и оценке угломестных параметров.
Первый и самый критичный компонент в цепочке обработки сигнала — это обнаружитель.
/* printf("(%d)ph0=%f ph1=%f w=%f", i, ph0[i], ph1[i], work);*/
ph1[i]=work;
}
(*counter)++;
/*printf("\ncounter=%d\n", *counter);*/
}
else{
/* no signal */
for(i=0;i<DIPOLS;i++){
ph1[i]=ph0[i];
}
*counter=0;
}
/* printf("d=%f factor=%f\n", d, factor);*/
returnd;
}
Обнаружитель определяет параметры ложной тревоги и пропуска сигнала и в конечном счете дальность действия радиопеленгатора. В предыдущих поколениях РП использовались амплитудные обнаружители. Приложение dfdetect.c работает по фазовому принципу и использует для оценки параметра обнаружения данные о фазовом раскрыве за два последовательных оборота антенной системы. Поскольку за межпериодный интервал электронного оборота корреляционные связи между фазами сигналов вибраторов, в том числе вызванные многолучевым распространением практически отсутствуют, используется метрика определяющая степень схожести двух фазовых раскрывов. Степень схожести задается как dot функция.
Весьма специфический алгоритм, характерный для радиопеленгатора — разрешение фазовой неоднозначности и восстановление доплеровской огибающей. Для современных РП характерна работа в частотных диапазонах, когда длина волны принимаемого сигнала существенно меньше базы антенной системы. Например, для АРП «Платан» для нижней точки диапазона 118 МГц λ =2.5 м, что с учетом радиуса антенны R = 1.6 м вызывает фазовую девиацию Ψ = 2π (R / λ) = 0.6 π. Если мы возьмем верхний диапазон 400 МГц, то Ψ = 2π (R / λ) = 2.1 π.
Картинка доплеровской фазовой огибающей с размахом ± 2π хороша только в книжке. Реально вы ее никогда не увидите, потому что фазовый детектор не в состоянии определить набег фазы больше 2π и выдает значения только в диапазоне 0..2π. Это значит, что огибающая будет нещадно разломана по высоте, и разглядеть в этих обломках отрезок синуса, да еще с учетом количества пространственных выборок = 16 (количество вибраторов антенной системы) будет весьма проблематично.
Есть и другая проблема — это вероятность аномального разрешения неоднозначности, т.е. построение огибающей сигнала с несуществующим азимутом, который может отличаться от истинного на любую величину. В общем, решать все эти проблемы призвано приложение dfres.c
dfres.c: обнаружитель пеленгационного сигнала
C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
/*
* dfres.c
* Copyright (c) 2002
* Author: Nazim Aliev, alnaz@mail.ru
*
* resolve phase ambiguity and restore Doppler modulation
dev, dev_x, q, q_x, s, min, dev_xstore, q_xstore);
*/
}/* for q */
}/*for dev */
for(i=0;i<DIPOLS;i++){
work=0;
fX=bufmin[i];
fY=phi[i];
/*fprintf(stderr, "phi=%f ", 180*fY/PI);*/
while(1){
if(fX>0)
fAbs=fX;
else
fAbs=-fX;
if(fAbs<PI/2)
break;
else
if(fX>0){
fX-=PI;
fY-=PI;
work++;
}
else{
fX+=PI;
fY+=PI;
work--;
}
}/* while */
pho[i]=fY;
/*fprintf(stderr, "work=%d\n", work);*/
}/* for */
returnmin;
}
Приложение dfres.c делает предположение об азимуте и девиации сигнала, создавая матрицу ошибок — отличия полученной огибающей от гипотетической. Для минимальной ошибки полученная огибающая восстанавливается (дополняется значениями π где это необходимо) по гипотетической огибающей. В результате выходной сигнал представляет собой восстановленную фазовую огибающую, по которой уже можно определить азимут как ее сдвиг по времени.
И наконец собственно оценка азимута.
dfdetect.c: обнаружитель пеленгационного сигнала
C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
/*
* dfangle.c
* Copyright(c) 2002
* Author: Nazim Aliev, alnaz@mail.ru
*
* estimate azimuth angle by imput phases
* f - input phases
* return value - Q
*
*/
#include "df.h"
#include "dfchan.h"
#include <math.h>
#define REFS 4
floatdfangle(float*f)
{
inti;
floats,c;
floatrs[REFS];
floatrc[REFS];
for(i=0;i<REFS;i++){
rs[i]=sin((DPI*i)/DIPOLS);
rc[i]=cos((DPI*i)/DIPOLS);
}
s=
(f[1]+f[7]-f[9]-f[15])*rs[1]+
(f[2]+f[6]-f[10]-f[14])*rs[2]+
(f[3]+f[5]-f[11]-f[13])*rs[3]+
f[4]-f[12];
c=
f[0]-f[8]+
(f[1]-f[7]-f[9]+f[15])*rc[1]+
(f[2]-f[6]-f[10]+f[14])*rc[2]+
(f[3]-f[5]-f[11]+f[13])*rc[3];
returnatan2(s,c);
/* e = sqrt(s*s + c*c); */
}
Приложение определяет первую гармонику входной фазовой последовательности методом FFT разлагая ее на комплексные составляющие. В результате получается вектор с амплитудой, соответствующей фазовой девиации, или углу места на источник сигнала, и углом поворота который соответствует пеленгу на источник сигнала.
Когда я смотрю на этот арктангенс в конце программы все время вспоминаю, как писали софт для 1806ВМ3 на ассемблере. Ведь эта тригонометрическая функция была единственным местом, где позарез требовалась математическая библиотека с плавающей точкой — а ее не было. Приходилось добывать целочисленный арктангенс целым алгоритмом, докручивая вектор в область линейного приближения. Знания добывали самостоятельно, потому что Гугла, как и интернета вообще — не было.
Данные приложения определяют ключевые параметры РП — дальность, точность и быстродействие; кроме этого есть неявные факторы: устойчивость обработки к перемодуляции и искажениям амплитудно-фазового раскрыва, вызванным многолучевым распространением.
Есть и другие программы, которые содержат вкусные вещи: подавление «вертолетного эффекта» — паразитной частотной модуляции при работе вертолетной радиостанции, защита от амплитудной перемодуляции, адаптивная обработка и другие. Я также планирую выложить эти алгоритмы.
Лабораторная сеть, в которой живут разные ARM/DSP/FPGA имеет фиксированные IP адреса и в общем случае не предназначена для выхода в инет. Однако часто возникает такая необходимость, например чтобы поставить новые пакеты. Это делается в несколько строк.
Конфигурация сети следующая:
лабораторная сетка 192.168.77.0;
в ней плата Zed Board с адресом 192.168.77.5 которой нужен выход в инет;
Wi-Fi сеть 192.168.1.0 которая имеет выход в инет через маршрутизатор провайдера;
комп который смотрит в лабораторную сеть интерфейсом eth1 с адресом 192.168.77.7 и в сеть Wi-Fi интерфейсом wlan1 с адресом 192.168.1.10 (эти адреса для настройки нам не понадобятся).
Делаем следующее, все команды — от root.
На плате Zed Board:
Shell
1
# route add default gw 192.168.77.7
На компе:
Shell
1
2
3
# echo 1 > /proc/sys/net/ipv4/ip_forward
# iptables -t nat -A POSTROUTING -o wlan1 -j MASQUERADE
# iptables -A FORWARD -i eth1 -j ACCEPT
Проверяем на Zed Board:
Shell
1
2
3
4
# ping www.ru
PING www.ru(217.112.35.75)56(84)bytes of data.
64bytes from v76-u.valuehost.ru(217.112.35.75):icmp_seq=1ttl=56time=49.2ms
64bytes from v76-u.valuehost.ru(217.112.35.75):icmp_seq=2ttl=56time=49.1ms
Для тех, кто пользует Arch Linux, обновляемся и ставим то что недостает на Zed Board:
Разбирая свойства 2D альфа — бета фильтра для траекторной обработки радарных данных, мне постоянно приходилось уходить в сторону ортогонализации входных данных. Поначалу не очень хотелось вникать в этот Gram–Schmidt process, чтобы не вдаваясь в детали попользоваться его результатами. Однако со временем стало ясно что без погружения не получится. Поэтому я отложил в сторону свой 2D Alpha-Beta Filter и решил сделать простое описание процесса Грама — Шмидта как памятку на будущее.
В интернете и литературе математическая сторона дела описана досконально; только ускользает физический смысл: зачем собственно эта ортогонализация нужна. Поэтому я сделал изложение с упором на рассмотрение конкретного примера, причем само изложение идет на трех уровнях: алгебраическом (кому нравится восприятие на абстрактном уровне), на уровне собственно примера с конкретными числовыми данными (для тех кто не любит теорию не подкрепленную практическими задачами) и также на уровне программной реализации. Для последнего я как обычно выбрал Python, который отлично подходит для инженерных задач; при этом можно было ограничиться нотацией питона для записи матричных и векторных уравнений. Тем не менее, формульная запись была сохранена для того чтобы оставить общность рассмотрения.
Сосед с перфоратором
Сформулируем исходную задачу на следующем примере. Мы решили заняться вокалом и записать свое собственное музыкальное произведение. После первых композиций активизируется ваш сосед и включает свой перфоратор, атакуя нашу общую стенку. Будущий сингл века под угрозой!
У вас есть организационно-административное и инженерное решение. Сразу опускаем первый вариант как самый скучный и непосредственно переходим к самому интересному, второму варианту. Очевидно, что можно расположить второй микрофон (а первый вы используете для записи своего голоса) непосредственно у стены и сделать запись шума перфоратора в надежде, что удастся скомпенсировать его в тракте первого, солирующего микрофона.
И эти надежды не беспочвенны. Для того чтобы двигаться дальше, начнем давать обозначения нашей конфигурации. Наша запись, сделанная с солирующего микрофона, подпорченная перфорирующими звуками будет называться X. Это будет вход фильтра, который мы собираемся использовать для подавления шума. Наша конечная цель — убрать из входного сигнала X этот шум подчистую. Другими словами, нам надо знать эту шумовую составляющую D, которая содержится в X, которая в адаптивной фильтрации именуется desired signal.
Здесь я должен вас предостеречь вот от чего. Признайтесь, в вашей голове уже мелькнула шальная мысль взять сигнал с «перфораторного» микрофона и просто вычесть его из записи X. Ну и что, если интенсивность помехи в полезном сигнале и у стенки разные: подкрутим уровень потенциометром, и всего делов. Сразу скажу, что потенциометров понадобится много, в том числе для компенсации сдвига фазы, учета расстояния, перемещения соседа с дрелью за стеной, изменения свойств среды распространения звука по причине допустим заливания соседей сверху… И крутить эти потенциометры придется непрерывно. Собственно, этим и занимается адаптивный фильтр: он непрерывно ищет измеренный сигнал помехи во входном сигнале и старается минимизировать значение ошибки. И если это так, то зачем нам делать за него всю эту работу?
Теперь мы понимаем, что у стены мы записываем шум как некоторую оценку D^ помеховой составляющей D. Эти сигналы связаны, похожи друг на друга, обладают схожими параметрами, но они не идентичны. При этом D^ нам известен в записи как помеха перфоратора у стены, а нашей конечной целью является D. А еще точнее — на самом деле, зная D, мы получим отфильтрованный сигнал X^ = X — D, очищенный от помехи насколько это возможно (в смысле минимума среднеквадратической ошибки, конечно!)
Микрофон включен… поехали
В нашем распоряжении сейчас следующие записи:
будем считать что это нота «ля» нашего вокала, искаженная перфоратором;
и скажем следующая «подпорченная» нота «си — бемоль».
Алгоритму все равно, каким образом представлены входные сигналы: здесь мы демонстрируем пятимерные векторы. Значок транспонирования, который был бы уместным чтобы не записывать векторы вертикально, для удобства опущен.
Соответственно, сигнал помехи D^ который мы «подслушали» у соседской стенки:
Теперь подготовим наш главный инструмент — язык Python. Работать будем прямо из командной строки Питона, набирая операторы по мере надобности. Запускаем интерпретатор (у меня в Ubuntu 14.04 стоит python3, а у вас?):
Команда импортирует модуль numpy, который мы будем интенсивно использовать для работы с векторами и матрицами.
Зададим три наших исходных вектора:
Python
1
2
3
4
>>>x1=np.array([-15,35,50,50,25])
>>>x2=np.array([-20,30,40,30,-5])
>>>de=np.array([-10,10,0,10,-10])
>>>
Векторы в какой — то части заданы совершенно произвольно, я только старался чтобы X1 не слишком сильно отличался от X2 и это различие как-то ассоциировалось с D^.
Замечу, что здесь я не придумал, как обозначить d с крышечкой сверху, и поэтому добавил к d букву e: estimated.
Все, для работы нам больше ничего не нужно. Следующий вопрос, на который надо ответить, будет таким: существует ли строгое теоретическое решение задачи нахождения desired signal D?
Классическое решение…
Строгое теоретическое решение существует, и оно представляет собой одну из форм выражения Винера — Хопфа (это известная тема и ее легко найти в интернете). Вот это соотношение:
(1)
Искомым значением является вектор оптимальных весовых коэффициентов W, с помощью которого по входному сигналу будет найдено значение помехи D:
Rxx представляет собой автокорреляционную функцию сигнала X и для нашего примера будет выглядеть следующим образом:
где выражение <X1 X2> представляет собой inner product, или dot product, определяемый вот так:
Попросим Питон посчитать автокоррелационную матрицу для сигналов, которые мы уже ввели в командной строке выше:
Python
1
2
3
4
5
6
7
8
9
>>>rxx=np.zeros((2,2))
>>>rxx[0,0]=np.dot(x1,x1)
>>>rxx[0,1]=np.dot(x1,x2)
>>>rxx[1,0]=np.dot(x2,x1)
>>>rxx[1,1]=np.dot(x2,x2)
>>>rxx
array([[7075.,4725.],
[4725.,3825.]])
>>>
Запишем автокорреляционную матрицу нашего сольного выступления, подпорченного помехой, выглядит следующим образом:
О чем говорит автокорреляция? Значения по главной диагонали 7075 и 3825 фактически представляют собой мощности сигналов X1 и X2 соответственно; по составу этих векторов мы и так видим, что X1 будет помощнее (нота «ля» особенно удалась). Гораздо интереснее значения вне главной диагонали, в нашем случае оно всего одно 4725, и характеризует степень взаимозависимости сигналов X1 и X2. Отметим это в уме и пойдем дальше.
Вектор Pdx показывает взаимную корреляцию между сигналом X и оценочным значением помехи D^:
Снова обращаемся к Питону:
Python
1
2
3
4
5
6
7
>>>pdx=np.zeros((2,1))
>>>pdx[0]=np.dot(x1,de)
>>>pdx[1]=np.dot(x2,de)
>>>pdx
array([[750.],
[850.]])
>>>
Теперь у нас есть все для нахождения вектора весовых коэффициентов из выражения (1).
В недрах модуля numpy Питона есть целый раздел функций для работы с линейной алгеброй, в том числе функция для решения подобных матричных уравнений. Воспользуемся ею и получим W:
Python
1
2
3
4
5
>>>w=np.linalg.solve(rxx,pdx)
>>>w
array([[-0.24228029],
[0.52150963]])
>>>
В результате, находим помеху D, а заодно и наш сигнал, очищенный от помехи:
Запишем красиво результат нашей работы, слегка округлив результаты:
Вот собственно и все, поздравляю с успешным решением задачи — наша песня снова звучит, очищенная от назойливого соседского перфоратора. И только один вопрос не дает покоя: при чем здесь собственно ортогонализация Грама — Шмидта?
… которое так и останется на страницах учебников
Если бы этот мир был идеальным, то все было бы просто. Одна строчка np.linalg.solve(a, b) на Питоне, и дело в шляпе. Но поскольку мы еще и инженеры, нас начинает интересовать вопрос производительности. Насколько быстро все это будет работать?
Для оценки быстродействия Питон любезно предлагает нам воспользоваться модулем datetime. Чтобы не создавать исходный сигнал вручную, эксперименты будем проводить над тестовой матрицей размером M, которую мы будем заполнять случайными данными. Делая отметки времени до и после выполнения функции np.linalg.solve(a, b), мы определим время, которое тратится на решение матричного уравнения.
Вот кусок кода, который выводит значение времени для матрицы размерностью в 1000 строк и столбцов:
Python
1
2
3
4
5
6
7
M=1000
a=5*np.random.random_sample((M,M))-5
b=5*np.random.random_sample((M,1))-5
n1=dt.datetime.now()
z=np.linalg.solve(a,b)
n2=dt.datetime.now()
print(M,n2-n1)
Результат теста:
Python
1
10000:00:00.119958
Результат в 120 миллисекунд наверное неплох. Последовательно наращиваем размер тестовой матрицы с 2000 до 8000 получаем:
Python
1
2
3
4
20000:00:00.816943
40000:00:05.747642
60000:00:18.522508
80000:00:43.691010
Представляете, сколько мне пришлось ждать решения уравнения с 8000 размером тестовой матрицы? Целых 44 секунды. Результат несколько обескураживает, но нас всегда предупреждали, что обращение матриц — это дорогая вычислительная операция.
Не стоит грешить на то, что Питон является скриптовым языком: основная работа идет в Си библиотеках модуля numpy (и есть еще фортрановские библиотеки — не ожидали?), поэтому Питон делает все максимально быстро. У меня обычный комп с тактовой частотой в несколько GHz. Конечно, результат будет лучше, если numpy backend Питона будет работать на числодробилке DSP или Neon (есть и такие реализации). Однако, каким бы не было быстродействие платформы, оно просто убивается увеличением размера выборки.
Хорошо заметно, что вычислительная сложность растет в геометрической прогрессии и это ставит крест на прямом использовании классических теоретических решений.
Ортогонализация как избавление от зависимости
Как существенным образом сократить трудоемкость вычислительного процесса? Все портит автокорреляционная матрица Rxx, в которой учтены все возможные зависимости элементов входного сигнала. Позволим себе немного пофантазировать и представить, что мы работаем с сигналом, для которого автокорреляционная матрица — диагональная, то есть все остальные элементы, не лежащие на главной диагонали, равны нулю. Это говорит о том, что корреляция между элементами сигнала отсутствует, что ведет нас к важному выводу: каждый элемент вектора D может быть оценен только по одному соответствующему элементу вектора X. Это резко сокращает вычислительный объем, который будет линейно зависеть от входной выборки. При этом векторы, описывающий сигнал, будут взаимно ортогональными.
Если опустить промежуточные выкладки, сделанные в предположении что входные сигналы взаимно некоррелированы и сответственно матрица Rxx является диагональной, оптимальные весовые коэффициенты могут быть определены независимо в следующей форме:
Соответственно, найденное решение D будет выглядеть следующим образом:
(2)
Здесь мы плавно подходим к процессу ортогонализации, результат которого примечателен тем, что в нем разрушаются зависимости, в общем случае присутствующие во входном сигнале и автокорреляционная матрица будет представлена как раз в диагональной форме. Однако, этот деструктивный процесс исказит собственно наш полезный сигнал, а это нам надо? С другой стороны, ортогонализация могла бы пригодится только для нахождения весовых коэффициентов, а сигнал уже можно оценить по исходной выборке, которая еще не подверглась ортогонализации.
Все это означает то, что нам нужен ответ на следующий вопрос: возможно ли привести задачу линейного оценивания набора произвольных векторов к задаче оценивания набора взаимно ортогональных векторов и получить идентичные результаты? К счастью, ответ на этот вопрос существует и он — утвердительный. Поэтому мы переходим к манипуляциям с нашими сигнальными векторами X1, X2. Для кого теория это кошмарный сон наяву, следующий параграф можно смело пропустить и продолжать рассмотрение нашего музыкального примера дальше.
Разложение как решение
Здесь наш анализ переходит в геометрическую плоскость, потому что на самом деле векторам, как произвольным так и ортогональным, положено комфортно чувствовать себя в геометрии. Чтобы не ломать голову над пятимерным пространством, рассмотрим пример на обычной плоскости с двумя сигнальными векторами X1, X2. Сейчас мы соединим алгебраическую задачу и ее понятную геометрическую интерпретацию (смотрите рисунок):
Ортогональное разложение сигнального вектора
Предположим, мы последовательно получили два сигнальных вектора X1, X2. После получения X1, мы выполнили оценку и теперь хотим улучшить ее, пользуясь тем что у нас есть еще и X2. Используя вектор X1 как основу новой, повернутой системы координат, разложим вектор X2 на ортогональные составляющие в этой системе координат. Проще говоря, опустим от конца X2 перпендикуляр на линию на которой лежит X1, и получим проекцию X2 на X1 которую обозначим projx1X2. Это будет одна ортогональная составляющая. Вторая составляющая, само собой, будет находится под прямым углом, и ее мы обозначим как V2.
Что касается X1, то будем считать что этот вектор разложен сам на себя и имеет единственную проекцию V1 = X1.
Данный процесс и носит название ортогонализации Грама — Шмидта, когда в результате разложения вектора X2 в системе координат вектора X1 получаем ортогональную пару векторов V1, V2. В общем случае процесс Грама — Шмидта конечно не ограничивается двумя векторами: если мы возьмем третий вектор X3, то он будет раскладываться на векторы X1, X2 и так далее. Если помимо алгебраического и геометрического слегка затронем схемотехнический аспект, то такое последовательное разложение по ортогональным составляющим приведет нас к схеме лестничного фильтра (Lattice Filter), но это уже совсем другая история, которой мы займемся позже.
Очевидно, что если бы не было помех, мы всегда получали одинаковый сигнал и в нашем случае выполнялось бы равенство X2 = X1. Поскольку присутствует помеха, векторы X1 и X2 отличаются по длине и по направлению. При этом геометрически и интуитивно понятно, что V2 является мерой отклонения X2 от X1 и есть не что иное, как величина ошибки.
Не полагаясь на интуицию, попробуем выразить высказанную мысль в математической форме. Из курса линейной алгебры известно, что проекцию одного вектора на другой можно найти следующим способом:
Сравниваем значение проекции projx1X2с выражением (2) которое было приведено ранее, и приходим к выводу, что это не что иное как оценка вектора X2 (который выполняет в этом выражении роль вектора D^) с использованием X1, поэтому мы запишем данное выражение со следующим представлением левой части:
Далее, из геометрических соображений очевидно что
и соответственно если выразить V2 из правой части то получим
Как следует из этого выражения, V2 есть не что иное, как разница между вектором и его оценкой, или другими словами — значение ошибки оценивания.
Жизнь в новом базисе
То, что мы сделали в предыдущем параграфе — заменили набор сигнальных векторов X1, X2 на их ортогональных наследников V1, V2. С вектором V1 все ясно:
Посчитаем что из себя представляет вектор V2:
где по исходным условиям нашей задачи (нота си-бемоль, не забыли?)
Вновь призываем Питон, который у нас в консольной строке:
Пока мы еще не отошли от Питона, проверим ортогональность векторов V1, V2:
Python
1
2
3
>>>np.dot(v1,v2)
5.6843418860808015e-13
>>>
Не будем слишком придирчивыми и посчитаем число в минус тринадцатой степени за ноль:
Это было так, к слову, что-то вроде проверочного упражнения.
Теперь наступает момент истины: до этого, для выражения (2) мы сделали утверждение, что во-первых в новом базисе V1, V2 оценка D будет считаться независимо, а во вторых — ее значение для этого базиса будет таким же, как и для теоретического метода. Проверим это: нам нужен результат выражений
Действительно, первое утверждение справедливо: для нахождения D(V1) используется только V1 и не используется V2, и наоборот (это следствие диагональности новой автокорреляционной матрицы!) Напомню, что исходный вектор D^ выглядит следующим образом:
Для нахождения значений оценок призываем наш безотказный Питон:
И это точно такое же значение ошибки, полученное теоретическим методом! Теперь у нас есть такой же отфильтрованный сигнал
только сформированный итерационным способом буквально из нескольких арифметических действий.
В результате: оргононализация разрушает корреляционные связи в сигнале, что позволяет представить сложную матричную задачу, которая учитывает все эти связи, в виде отдельных независимых уравнений, решение которых не представляет вычислительной сложности.
Теперь с спокойной душой можно вернуться к 2D alpha — beta фильтру.
(1)
(2)
Last comments