Моя реализация фильтра Калмана, кроме сглаженных координат, дает еще вектор скорости, который мне удалось перевести из угловых величин в метры за секунду. Теперь, если даже передвигаться по Wi-Fi точкам, от которых локейшен-менеджер не может получить скорость, в моем компоненте скорость будет доступна и довольно достоверна.
Скоро напишу подробнее о ФК и его реализации на Objective-C.
Я тут нашел, как сымитировать движение по собственному треку на симуляторе и хочу показать вам Фильтр Калмана в действии.
Я подписал там характерные моменты и добавил музыку, чтобы не скучно было смотреть. Видео немного поможет вам почувствовать это волшебство, когда приемник передает позицию то там, то здесь, а булавочка все равно идет вместе с тобой.
Я переписал фильтр Калмана для матриц размером 4x4, чтобы его можно было использовать для уточнения долготы и широты. Кроме того, в таком виде он сможет вычислять и скорость движения более точно. После отладки фильтра на Питоне, я переписал его на Objective-C++, выбрав матричную библиотеку GLM.
Матричный Калман показал себя очень хорошо. Вечером с Кеней ходили вокруг квартала 85% по вайфаю — все равно показывает фактически, как идешь, но с отставанием метров на 30-40. При этом измеряемый локейшен прыгает вокруг тебя чудовищно, от здания к зданию, а трек все равно выглядит, как по точному GPS-сигналу.
Я шел и не верил своим глазам. Такое впечатление, что это какая-то магия. Когда приемник переходит на спутниковый точный сигнал, отставание фильтрованной позиции от измеренной уходит и идешь, как обычно.
Вот и настал момент, когда пришлось с ним познакомиться поближе. Вкратце, это математический аппарат, который позволяет сглаживать данные на лету, не накапливая их для анализа. Используется для обработки показаний датчиков и любых приборов. Как правило, такие показания подвержены шумам и отклонениям, их нужно отсекать. Фильтра Калмана (ФК) позволяет отбрасывать пики и видеть усредненную, более вероятную картину процесса.
Рудольф Калман. Что бы мы без него делали. А ведь когда он придумал свой фильтр в 1960 году, то встретился с таким скептицизмом, что был вынужден опубликовать первые работы о нем в журнале, связанном с механикой, хотя сам занимался электротехникой.
Применений множество. Представьте датчик топлива в баке автомобиля. Он же болтается там внутри, если просто записать его показания, а потом попытаться проанализировать, получится, что в какой-то момент в баке было 30 л, через секунду 50, а потом 27. А это просто волны внутри гуляют. Сколько же в этот момент было в баке? ФК может дать более вероятный ответ.
Я сейчас делаю навигационные компоненты для iOS и столкнулся с проблемой биения сигнала и резкими скачками позиции при включении приемника. Переключение со спутникового сигнала на вайфай тоже дает неприятные перепады. Кроме того, моментальную и среднюю скорость тоже неплохо бы фильтровать. И начал разбираться.
Материалов по ФК много, мне показались лучшими эти документы: раз и два. Но ими лучше пользоваться потом, когда суть работы фильтра уже станет понятной.
Если очень кратко, то суть в том, что большинство явлений подчиняются Гауссовскому распределению вероятностей.
Например, мы думаем, что находимся в точке x. Но мы ведь можем ошибаться, потому, например, что отмеряли расстояние до x каким-нибудь прибором, а он, конечно, врет. Поэтому мы считаем, что находимся в точке x с какой-то вероятностью p. А может, мы сейчас в точке x', но вероятность этого меньше — p'. Ну такое вот распределение вероятностей.
Наиболее вероятный x посредине — это μ, математическое ожидание. Ширина этой "шляпы" зависит от σ2 — это дисперсия, квадрат стандартного отклонения. Чем больше отклонение, тем шляпа шире и неопределенность больше.
У этой кривой такое уравнение:
(1)
Но его бояться не надо, т.к. нам надо искать середину, самый вероятный x. А в середине он равен μ, в числителе степени разница между x и μ дает 0, это многое потом упрощает.
Что же делает фильтр? Он помнит только одно значение x1 и его дисперсию σ2. Когда поступает новое значение x2 со своим отклонением r2, он умножает их (по правилам распределения, см. формулу 2) друг на друга и в результате будет новое значение x3 с погрешностью меньшей, чем у двух предыдущих!
Например, текущее значение фильтра — синее распределение. Поступает новое значение, со своим отклонением (красное). Фильтр перемножает их по правилам распределения (см. формулу 2) и получает новое распределение (зеленое). Оно гораздо точнее двух изначальных, т.к. его дисперсия меньше (см. формулу 3) и становится текущим значением фильтра.
Умножение распределений делается так:
(2)
Дисперсия тоже уточняется:
(3)
И тоже становится текущей. Поступает еще значение, состояние еще больше уточняется и так далее.
Я разобрался с принципом работы фильтра, пройдя этот онлайн-курс. Там дядечка рассказывает все довольно понятно, на каком-то этапе видео останавливается и ты должен ответить на вопросы или дописать кусок функции (на питоне), чтобы подтвердить, что ты усвоил материал.
В процессе прохождения этого курса я написал скалярную и матричную версию фильтра, выкладываю их здесь.
# Multi-dimensional Kalman Filter from online course # https://www.udacity.com/course/viewer#!/c-cs373/l-48723604/m-48670988 # written by Pavel Malinnikov, May 2013
def__init__(self, value): self.value= value self.dimx=len(value) self.dimy=len(value[0]) if value ==[[]]: self.dimx=0
def zero(self, dimx, dimy): # check if valid dimensions if dimx <1or dimy <1: raiseValueError,"Invalid size of matrix" else: self.dimx= dimx self.dimy= dimy self.value=[[0for row inrange(dimy)]for col inrange(dimx)]
def identity(self, dim): # check if valid dimension if dim <1: raiseValueError,"Invalid size of matrix" else: self.dimx= dim self.dimy= dim self.value=[[0for row inrange(dim)]for col inrange(dim)] for i inrange(dim): self.value[i][i]=1
def show(self): for i inrange(self.dimx): printself.value[i] print' '
def__add__(self, other): # check if correct dimensions ifself.dimx!= other.dimxorself.dimy!= other.dimy: raiseValueError,"Matrices must be of equal dimensions to add" else: # add if correct dimensions
res = matrix([[]])
res.zero(self.dimx,self.dimy) for i inrange(self.dimx): for j inrange(self.dimy):
res.value[i][j]=self.value[i][j] + other.value[i][j] return res
def__sub__(self, other): # check if correct dimensions ifself.dimx!= other.dimxorself.dimy!= other.dimy: raiseValueError,"Matrices must be of equal dimensions to subtract" else: # subtract if correct dimensions
res = matrix([[]])
res.zero(self.dimx,self.dimy) for i inrange(self.dimx): for j inrange(self.dimy):
res.value[i][j]=self.value[i][j] - other.value[i][j] return res
def__mul__(self, other): # check if correct dimensions ifself.dimy!= other.dimx: raiseValueError,"Matrices must be m*n and n*p to multiply" else: # subtract if correct dimensions
res = matrix([[]])
res.zero(self.dimx, other.dimy) for i inrange(self.dimx): for j inrange(other.dimy): for k inrange(self.dimy):
res.value[i][j] +=self.value[i][k] * other.value[k][j] return res
def transpose(self): # compute transpose
res = matrix([[]])
res.zero(self.dimy,self.dimx) for i inrange(self.dimx): for j inrange(self.dimy):
res.value[j][i]=self.value[i][j] return res
# Thanks to Ernesto P. Adorio for use of Cholesky and CholeskyInverse functions
def Cholesky(self, ztol =1.0e-5): # Computes the upper triangular Cholesky factorization of # a positive definite matrix.
res = matrix([[]])
res.zero(self.dimx,self.dimx)
for i inrange(self.dimx):
S =sum([(res.value[k][i]) ** 2for k inrange(i)])
d =self.value[i][i] - S ifabs(d)< ztol:
res.value[i][i]=0.0 else: if d <0.0: raiseValueError,"Matrix not positive-definite"
res.value[i][i]= sqrt(d) for j inrange(i + 1,self.dimx):
S =sum([res.value[k][i] * res.value[k][j]for k inrange(self.dimx)]) ifabs(S)< ztol:
S =0.0
res.value[i][j]=(self.value[i][j] - S) / res.value[i][i] return res
def CholeskyInverse(self): # Computes inverse of matrix given its Cholesky upper Triangular # decomposition of matrix.
res = matrix([[]])
res.zero(self.dimx,self.dimx)
# Backward step for inverse. for j inreversed(range(self.dimx)):
tjj =self.value[j][j]
S =sum([self.value[j][k] * res.value[j][k]for k inrange(j + 1,self.dimx)])
res.value[j][j]=1.0 / tjj ** 2 - S / tjj for i inreversed(range(j)):
res.value[j][i]= res.value[i][j]= -sum([self.value[i][k] * res.value[k][j]for k inrange(i + 1,self.dimx)]) / self.value[i][i] return res
def inverse(self):
aux =self.Cholesky()
res = aux.CholeskyInverse() return res
def__repr__(self): returnrepr(self.value)
########################################
result =[]
def kalman(x, P): for n inrange(len(positions)):
# measurement update
z = matrix([[positions[n]]])
y = z - H * x
S = H * P * H.transpose() + R
K = P * H.transpose() * S.inverse()
x = x + (K * y)
P =(I - K * H) * P
# prediction
x = F * x + u
P = F * P * F.transpose()
positions =[x + random() * randint(-1,1) * 5for x in data]
x = matrix([[0.],[0.]])# initial state (location and velocity)
P = matrix([[1000.,0.],[0.,1000.]])# initial uncertainty
u = matrix([[0.],[0.]])# external motion
F = matrix([[1.,1.],[0,1.]])# next state function
H = matrix([[1.,0.]])# measurement function
R = matrix([[4.]])# measurement uncertainty
I = matrix([[1.,0.],[0.,1.]])# identity matrix
А это скалярная, упрощенная реализация. Я сразу попробовал ее на раздельном сглаживании широты и долготы. Тренировался на треках, которые люди загружают на OpenStreetMap. Вот она:
# One-dimensional Kalman Filter testing for geo-track data smoothing # Tracks were taken from http://www.openstreetmap.org/traces # Pavel Malinnikov, May 2013
import matplotlib.pyplotas plt importxml.etree.ElementTreeas ET
tree = ET.parse('/Users/pavelmalinnikov/Downloads/1464857.gpx')
root = tree.getroot()
Яндекс придумал хорошую вещь — возможность для сайтов вставлять виджеты с частью функциональности в результат поиска:
Мне нравится эта идея, в первую очередь тем, что она должна заставить сайты при разработке подобных виджетов привести свою основную функциональность в порядок.
Магазинам актуализировать данные об ассортименте и остатках, таксопаркам о количестве свободных машин и т.д.
Боюсь только, что виджеты будут делать только те сайты, у которых и так все в порядке. А остальное большинство охламонов никакими островами не проймешь.
Меня тут часто стали спрашивать, чем я занимаюсь сейчас и вообще. Я вам сделаю небольшой обзорчик некоторых последних проектов.
Вот, например, одна фирма, разрабатывающая софт для Европы, заказывает у нас компоненты для iOS-приложений. Проекты небольшие, на 2 — 3 недели, но заковыристые. Нужно делать не совсем тривиальные вещи, требующие много исследований.
Чтобы дать представление об этих компонентах, покажу пару пояснительных записочек к некоторым из них.
Вот компоненты для создания навигационных программ:
А вот листалка страниц с предварительным кешированием:
Видео демки:
Пусть люди делают клевые приложения, мне будет приятно знать, что в них есть мои компоненты. Они тоже клевые, за них мне не стыдно, вот тебе и источник хорошего настроения.
У меня тоже есть программа в эппловском Аппсторе. "The Tracker" называется. Показывает друзей на карте.
Их можно не только в телефоне видеть, кстати, но и на сайте GPSov.net. Пользователей около 5000 человек уже, но активно пользуются около 10%:
Как-то для английской фирмы Calibr8 Ltd я делал софтину для работы с их оборудованием через ком-порт. Несмотря на то, что программа под Виндовс, делал я ее в макоси. Англичанам была нужна поддержка от Виндовс 98 до Виндовс 7. Сделал на С++ (MFC) — красота. Англичане довольны: "We are very pleased with the software", говорят. Ну а чего же не быть pleased, если там все четко, асинхронно, безопасно, как я люблю. Тоже есть видео разработки, кстати.
Недавно обнаружил, что они на Ebay ее продают по 55 фунтов. Надо было с них больше денег брать :-)
Ну в общем, все в таком духе. Но это проекты, которые мы делаем для внешних заказчиков.
А вообще основная работа происходит по решению задач предприятия. Автоматизация, документооборот, в торговой компании много всего. Оракл, Flex, яваскрипт, php, C++, всего не перечислишь. Вот сейчас синхронизацию с ToodleDo делаю, завтра еще что-то будет.
Сделал сайт для нашей школы. Скромно заявляю, что теперь это самый лучший сайт школы в Луганске. Страницы адаптируются под разные разрешения, можно смотреть с телефона.
Материалы, которая предоставила школа, размещены там самым человечным способом. Вроде бы просто, надо взять шаблон поприличнее, допилить его и аккуратно опубликовать материалы. Почему большинство сайтов школ сделаны так коряво — загадка.
Вот сайт нашей соседней школы, №60. Даже текст читается с трудом, не говоря уже о том, чтобы телефон школы вообще узнать:
Дорогие веб-дизайнеры! Пора задуматься о том, что нужно побеждать свои проявления кустарности и местечковости. Пользуйтесь готовыми шаблонами, но развивайте вкус и аккуратность!
P.S. Адрес сайта для 59-й школы я сделал school59.lg.ua (школа 59 - Луганск - Украина). А сайт 60-й школы почему-то расположен в коммерческой зоне (school60.com.ua) — я понимаю, что это не намеренный намек, но все равно забавно.
После перехода на Lion обнаружилось, что оракловый клиент под макосью больше нормально не работает. Вызов sqlpus вызывает «Segmentation fault: 11». Установка 32-битной версии клиента позволяет пользоваться плюсом (по слухам) но не решает проблему доступа к ораклу из php через oci8. Расширение oci8 ведь должно компилироваться для php, как 64-битная библиотека и с сдк от 32-битного инстанс-клиента ее сбилдить не получится.
Глобальное решение давно назрело и обкатано на работе — виртуальные машины.
Вот и дома я решил настроить какой-нибудь линукс в виртуалке, в качестве вебсервера. Кроме того, что его всегда можно будет откатить к любому сохраненному состоянию, единожды настроенную машину можно будет использовать на других тачках, да и если придется переустановить всю систему, не придется настраивать веб-сервер заново.