Мой ввод представляет собой временной ряд 2d (x, y) точки, движущейся по экрану для программного обеспечения трекера. У него есть шум, который я хочу удалить с помощью фильтра Калмана. Кто-нибудь может указать мне код Python для фильтра Kalman 2d? В кулинарной книге scipy я нашел только 1 пример: http://www.scipy.org/Cookbook/KalmanFiltering я видел там это реализация фильтра Калмана в OpenCV, но не удалось найти примеры кода. Спасибо!
Kalman 2d фильтр в питоне
Ответы (2)
Вот моя реализация фильтра Калмана на основе уравнений, приведенных в Википедии. Имейте в виду, что мое понимание фильтров Калмана очень рудиментарно, поэтому, скорее всего, есть способы улучшить этот код. (Например, он страдает от численной нестабильности, обсуждаемой здесь. Насколько я понимаю, это влияет на численную стабильность только тогда, когда Q
, шум движения, очень мал. В реальной жизни шум обычно не мал, поэтому, к счастью (по крайней мере, для моей реализации) на практике числовая нестабильность не проявляется.)
В приведенном ниже примере kalman_xy
предполагает, что вектор состояния представляет собой 4-кортеж: 2 числа для местоположения и 2 числа для скорости. Матрицы F
и H
были определены специально для этого вектора состояния: если x
является состоянием из 4 кортежей, то
new_x = F * x
position = H * x
Затем он вызывает kalman
, который является обобщенным фильтром Калмана. Он является общим в том смысле, что он по-прежнему полезен, если вы хотите определить другой вектор состояния — возможно, 6-кортеж, представляющий местоположение, скорость и ускорение. Вам просто нужно определить уравнения движения, указав соответствующие F
и H
.
import numpy as np
import matplotlib.pyplot as plt
def kalman_xy(x, P, measurement, R,
motion = np.matrix('0. 0. 0. 0.').T,
Q = np.matrix(np.eye(4))):
"""
Parameters:
x: initial state 4-tuple of location and velocity: (x0, x1, x0_dot, x1_dot)
P: initial uncertainty convariance matrix
measurement: observed position
R: measurement noise
motion: external motion added to state vector x
Q: motion noise (same shape as P)
"""
return kalman(x, P, measurement, R, motion, Q,
F = np.matrix('''
1. 0. 1. 0.;
0. 1. 0. 1.;
0. 0. 1. 0.;
0. 0. 0. 1.
'''),
H = np.matrix('''
1. 0. 0. 0.;
0. 1. 0. 0.'''))
def kalman(x, P, measurement, R, motion, Q, F, H):
'''
Parameters:
x: initial state
P: initial uncertainty convariance matrix
measurement: observed position (same shape as H*x)
R: measurement noise (same shape as H)
motion: external motion added to state vector x
Q: motion noise (same shape as P)
F: next state function: x_prime = F*x
H: measurement function: position = H*x
Return: the updated and predicted new values for (x, P)
See also http://en.wikipedia.org/wiki/Kalman_filter
This version of kalman can be applied to many different situations by
appropriately defining F and H
'''
# UPDATE x, P based on measurement m
# distance between measured and current position-belief
y = np.matrix(measurement).T - H * x
S = H * P * H.T + R # residual convariance
K = P * H.T * S.I # Kalman gain
x = x + K*y
I = np.matrix(np.eye(F.shape[0])) # identity matrix
P = (I - K*H)*P
# PREDICT x, P based on motion
x = F*x + motion
P = F*P*F.T + Q
return x, P
def demo_kalman_xy():
x = np.matrix('0. 0. 0. 0.').T
P = np.matrix(np.eye(4))*1000 # initial uncertainty
N = 20
true_x = np.linspace(0.0, 10.0, N)
true_y = true_x**2
observed_x = true_x + 0.05*np.random.random(N)*true_x
observed_y = true_y + 0.05*np.random.random(N)*true_y
plt.plot(observed_x, observed_y, 'ro')
result = []
R = 0.01**2
for meas in zip(observed_x, observed_y):
x, P = kalman_xy(x, P, meas, R)
result.append((x[:2]).tolist())
kalman_x, kalman_y = zip(*result)
plt.plot(kalman_x, kalman_y, 'g-')
plt.show()
demo_kalman_xy()
Красные точки показывают зашумленные измерения положения, зеленая линия показывает положения, предсказанные Калманом.
P = np.matrix(np.eye(4))*1000
. Зачем умножать на 1000?
- person Sophia; 16.03.2019
demo_kalman_xy()
R действительно должна быть матрицей 2x2 (ковариация шума измерения), например R = np.matrix([[r,0], [0,r]])
. Значение Q по умолчанию в kalman_xy()
также, вероятно, слишком велико, чтобы можно было легко увидеть эффекты регулировки R.
- person hyperspasm; 03.01.2021
Для моего проекта мне нужно было создать интервалы для моделирования временных рядов, и, чтобы сделать процедуру более эффективной, я создал tsmoothie: библиотека Python для сглаживания временных рядов и обнаружения выбросов векторизованным способом.
Он предоставляет различные алгоритмы сглаживания вместе с возможностью вычисления интервалов.
В случае KalmanSmoother
вы можете управлять сглаживанием кривой, объединяя различные компоненты: уровень, тренд, сезонность, длинная сезонность.
import numpy as np
import matplotlib.pyplot as plt
from tsmoothie.smoother import *
from tsmoothie.utils_func import sim_randomwalk
# generate 3 randomwalks timeseries of lenght 100
np.random.seed(123)
data = sim_randomwalk(n_series=3, timesteps=100,
process_noise=10, measure_noise=30)
# operate smoothing
smoother = KalmanSmoother(component='level_trend',
component_noise={'level':0.1, 'trend':0.1})
smoother.smooth(data)
# generate intervals
low, up = smoother.get_intervals('kalman_interval', confidence=0.05)
# plot the first smoothed timeseries with intervals
plt.figure(figsize=(11,6))
plt.plot(smoother.smooth_data[0], linewidth=3, color='blue')
plt.plot(smoother.data[0], '.k')
plt.fill_between(range(len(smoother.data[0])), low[0], up[0], alpha=0.3)
Я также указываю, что tsmoothie может выполнять сглаживание нескольких временных рядов векторизованным способом.