Фильтрация 1D-массивов в Python

Question

Объяснение:

У меня есть два массива: dataX и dataY, и я пытаюсь отфильтровать каждый массив, чтобы уменьшить шум.Изображение, показанное ниже, показывает фактические входные данные (синие точки) и пример того, как я хочу, чтобы это было (красные точки).Мне не нужно, чтобы отфильтрованные данные были идеальными, как в примере, но я хочу, чтобы они были максимально прямыми.Я предоставил пример данных в коде.

Что я пробовал:

Во-первых, вы можете видеть, что данные не являются «непрерывными», поэтому я сначаларазделил их на отдельные «сегменты» (4 из них в этом примере), а затем применил фильтр к каждому «сегменту».Кто-то предложил мне использовать фильтр Савицкого-Голея .Полный, исполняемый код приведен ниже:

import scipy as sc
import scipy.signal
import numpy as np
import matplotlib.pyplot as plt

# Sample Data
ydata = np.array([1,0,1,2,1,2,1,0,1,1,2,2,0,0,1,0,1,0,1,2,7,6,8,6,8,6,6,8,6,6,8,6,6,7,6,5,5,6,6, 10,11,12,13,12,11,10,10,11,10,12,11,10,10,10,10,12,12,10,10,17,16,15,17,16, 17,16,18,19,18,17,16,16,16,16,16,15,16])
xdata = np.array([1,2,3,1,5,4,7,8,6,10,11,12,13,10,12,13,17,16,19,18,21,19,23,21,25,20,26,27,28,26,26,26,29,30,30,29,30,32,33, 1,2,3,1,5,4,7,8,6,10,11,12,13,10,12,13,17,16,19,18,21,19,23,21,25,20,26,27,28,26,26,26,29,30,30,29,30,32])

# Used a diff array to find where there is a big change in Y. 
# If there's a big change in Y, then there must be a change of 'segment'.
diffy = np.diff(ydata)

# Create empty numpy arrays to append values into
filteredX = np.array([])
filteredY = np.array([])

# Chose 3 to be the value indicating the change in Y
index = np.where(diffy >3)

# Loop through the array
start = 0
for i in range (0, (index[0].size +1) ):
# Check if last segment is reached
    if i == index[0].size:
        print xdata[start:]
        partSize = xdata[start:].size
# Window length must be an odd integer
        if partSize % 2 == 0:
            partSize = partSize - 1

        filteredDataX = sc.signal.savgol_filter(xdata[start:], partSize, 3)
        filteredDataY = sc.signal.savgol_filter(ydata[start:], partSize, 3)
        filteredX = np.append(filteredX, filteredDataX)
        filteredY = np.append(filteredY, filteredDataY)

    else:
        print xdata[start:index[0][i]]
        partSize = xdata[start:index[0][i]].size
        if partSize % 2 == 0:
            partSize = partSize - 1
        filteredDataX = sc.signal.savgol_filter(xdata[start:index[0][i]], partSize, 3)
        filteredDataY = sc.signal.savgol_filter(ydata[start:index[0][i]], partSize, 3)
        start = index[0][i] 
        filteredX = np.append(filteredX, filteredDataX)
        filteredY = np.append(filteredY, filteredDataY)

# Plots
plt.plot(xdata,ydata, 'bo', label = 'Input Data')
plt.plot(filteredX, filteredY, 'ro', label = 'Filtered Data')

plt.xlabel('X')
plt.ylabel('Y')
plt.title('Result')
plt.legend()
plt.show()

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

Я также попробовал эти фильтры, среди нескольких других:

• scipy.signal.medfilt
• scipy.ndimage.filters.uniform_filter1d

Но пока ни один изфильтры, которые я пробовал, были близки к тому, что я действительно хотел.Каков наилучший способ фильтрации данных, таких как этот?Ждем вашей помощи.

Answer 1

Один из способов получить что-то похожее на ваш идеал - кластеризация + линейная регрессия.

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

import numpy as np
from scipy import cluster, stats

ydata = np.array([1,0,1,2,1,2,1,0,1,1,2,2,0,0,1,0,1,0,1,2,7,6,8,6,8,6,6,8,6,6,8,6,6,7,6,5,5,6,6, 10,11,12,13,12,11,10,10,11,10,12,11,10,10,10,10,12,12,10,10,17,16,15,17,16, 17,16,18,19,18,17,16,16,16,16,16,15,16])
xdata = np.array([1,2,3,1,5,4,7,8,6,10,11,12,13,10,12,13,17,16,19,18,21,19,23,21,25,20,26,27,28,26,26,26,29,30,30,29,30,32,33, 1,2,3,1,5,4,7,8,6,10,11,12,13,10,12,13,17,16,19,18,21,19,23,21,25,20,26,27,28,26,26,26,29,30,30,29,30,32])

def split_to_lines(x, y, k):
    yo = np.empty_like(y, dtype=float)
    # get the cluster centers and the labels for each point
    centers, map_ = cluster.vq.kmeans2(np.array((x, y * 2)).T.astype(float), k)
    # for each cluster, use the labels to select the points belonging to
    # the cluster and do a linear regression
    for i in range(k):
        slope, interc, *_ = stats.linregress(x[map_==i], y[map_==i])
        # use the regression parameters to construct y values on the
        # best fit line
        yo[map_==i] = x[map_==i] * slope + interc
    return yo

import pylab
pylab.plot(xdata, ydata, 'or')
pylab.plot(xdata, split_to_lines(xdata, ydata, 4), 'ob')
pylab.show()