时间/位置列表上的Python低通滤波器

Question

让我先解释一下我的目标和问题。 我有一个物体在bocal（准确地说是一只老鼠）内循环移动，我需要知道它在一分钟内旋转了多少圈。这个运动不是一个完美的圆，而且可以是混乱的。它可以停止几秒钟，然后重新开始移动。在

但我知道的是，我每分钟旋转3到20次（0.05赫兹到0.33赫兹）。在

我用一个30 fps的相机在30分钟内拍摄这个物体。 使用Python和OpenCV，我设法提取（X，Y）坐标。由于30fps对于我所寻找的频率来说有点高，所以我选择了15帧，以便采样频率为2hz。在

第一个问题是，有时候，我会漏掉一点，因为OpenCV找不到对象。这是一个我无法解决的问题，因为我用HSV色罩跟踪的点有时是隐藏的。在

事实上，我要过滤的数据是：

• 帧，帧编号的列表。在符号[i]处，我有对应于X[i]和Y[i]的帧编号
• 十、 X坐标的列表
• Y、 Y坐标的列表

例如，我有这样的列表：

• 框架=[15，90，165，180，195，210，225，300，315，375，405]
• X=[395、487、389、389、358、382、373、389、397、403、446]
• Y=[167、211、166、173、180、185、191、223、241、231、238]

This is the king of movement I get, as you can see, I need to filter the high frequency.

现在我的问题是：即使我有一些过滤和编码的基础知识，我也不知道怎么做，也不知道该用哪个库。 我读过scipy应该有所有必需的功能，但我不知道该用哪一个，以及如何使用。在

我不确定，但我想我应该做这样的事情：

• 富里埃 转化
• 低通滤波
• 反变换

你能告诉我我是不是对的，并给我指出编码部分的正确方向吗？在

谢谢你的帮助

马修

Answer 1

下面是我用你告诉我的代码： 我创建了两个函数：

• 插值，如果帧 间隔不超过6个间隔。否则它将插值 按段划分。在
• 过滤也在对段进行过滤。在

目前，我还没有尝试改变a的值来找到正确的值，但它确实给了我很好的结果。在

Preliminary results

def Interpolation(Frame, X, Y, Frequence = 2):
    """Fonction réalisant l'interpolation cubique des points manquants si possible, avec le découpage en segment si nécessaire."""

    # Détermination de la séparation entre les Frames pour la fréquence d'échantillonage choisie
    FPS = 30
    DT = FPS / Frequence

    # Initialisation
    i = 0

    # Détermination du début de l'interpolation
    try:
        while Frame[i+2] - Frame[i] != 2 * DT:
            i += 1
            continue
    except:
        print ("Erreur : Il est impossible de commencer l'interpolation.")

    # Suppresion des points antérieurs
    Frame = Frame[i:]
    X = X[i:]
    Y = Y[i:]

    k = 0
    Frame_id = list()
    # Recherche des segments
    for i in range(len(Frame) - 1):
        if Frame[i+1] - Frame[i] > 6 * DT:
            Frame_id.append((k, i))
            k = i + 1

    # Ajout du dernier segment
    if k != len(Frame) - 1:
        Frame_id.append((k, len(Frame)-1))

    # Suppresion des segments de moins de 3 points
    for k in reversed(range(len(Frame_id))):
        i = Frame_id[k][0]
        j = Frame_id[k][1]

        if j - i < 2:
            del Frame_id[k]
            continue

    Abs_inter, X_inter, Y_inter = [], [], []
    for index in Frame_id:
        i = index[0]
        j = index[1]

        # Génération des abscisses sur chaque segment, ainsi que des fonctions interpolés
        Abscisse = np.arange(Frame[i], Frame[j] + DT, DT)

        fX = interp1d(Frame[i:j+1], X[i:j+1], kind = 'cubic')
        fY = interp1d(Frame[i:j+1], Y[i:j+1], kind = 'cubic')

        # Génération des nouvelles coordonnées sur le segment
        Xnew = fX(Abscisse)
        Ynew = fY(Abscisse)

        Abs_inter += Abscisse.tolist()
        X_inter += Xnew.tolist()
        Y_inter += Ynew.tolist()

    # Création des listes résultats finaux
    Frame_final, X_final, Y_final = [], [], []

    Time = 0
    while Time <= Abs_inter[len(Abs_inter)-1]:
        if Time in Frame:
            Frame_final.append(Time)
            X_final.append(X[Frame.index(Time)])
            Y_final.append(Y[Frame.index(Time)])
            Time += DT
            continue
        elif Time in Abs_inter:
            Frame_final.append(Time)
            X_final.append(X_inter[Abs_inter.index(Time)])
            Y_final.append(Y_inter[Abs_inter.index(Time)])
            Time += DT
            continue
        else:
            Time += DT

    return (Frame_final, X_final, Y_final)

def Filtrage2(Frame, X, Y, DT = 15, a = 0.1):
    """Fonction réalisant un filtrage passe-bas d'ordre 1 du signal."""
    # Initialisation
    X_temp, Y_temp = [], []
    X_filter, Y_filter = [], []
    i = 1
    X_temp.append(X[0])
    Y_temp.append(Y[0])

    # Filtrage par morceau
    while i < len(Frame):

        while Frame[i] - Frame[i-1] == DT:
            Xnew = a * X[i-1] + (1 - a) * X_temp[len(X_temp)-1]
            Ynew = a * Y[i-1] + (1 - a) * Y_temp[len(Y_temp)-1]
            X_temp.append(Xnew)
            Y_temp.append(Ynew)
            if i < len(Frame)-1:
                i += 1
            else:
                break

        X_filter += X_temp
        Y_filter += Y_temp
        X_temp, Y_temp = [], []
        X_temp.append(X[i])
        Y_temp.append(Y[i])
        i += 1

    return (X_filter, Y_filter)

我也尝试过在scipy库中使用butterworth过滤器，但是我无法使其工作。以下是我无法理解的代码：

^{pr2}$

如果你知道如何让这个巴特沃斯工作，我很乐意听到。在

Answer 2

在深入研究傅立叶变换之前，您可以只应用一阶或二阶低通滤波器。你可以先线性插值数据，这样你就可以有一个恒定的2Hz频率。然后可以对数据点应用一阶低通滤波器。在

y_k = a * x_k + (1-a) * y_km1, a in [0,1]

x是你的观察值，y是你的过滤后的估计值。在

然后，如果第一个原型产生一些有用的结果，您就可以使用一些物理模型来获得更好的估计值。Kalman滤波器可以更好地利用你潜在的物理现实。但要做到这一点，你首先需要知道如何模拟物理现实。在

https://en.wikipedia.org/wiki/Kalman_filter

在这里，你也许还可以找到一个更接近计算机视觉跟踪运动的例子： http://www.diss.fu-berlin.de/docs/servlets/MCRFileNodeServlet/FUDOCS_derivate_000000000473/2005_12.pdf