平滑循环数据

Smooth circular data

我有一个数据数组 Y 使得 Y 是自变量 X(另一个数组)的函数。

X 中的值从 0 到 360 不等,带有回绕。

Y 中的值从 -180 到 180 不等,也有回绕。

(也就是说,这些值是围绕一个圆的角度。)

有谁知道 Python 中的任何函数(在 numpyscipy 等中)能够将我的 Y 值作为函数进行低通滤波X?

如果这有点令人困惑,这里有一个示例数据图:

这是一个使用 pandas 做移动平均线的解决方案。首先 unwrap 数据(需要转换为弧度并返回),因此没有不连续性(例如,从 180 跳到 -179)。然后计算移动平均值,如果需要,最后转换回包装数据。另外,使用 np.convolve().

检查这个 numpy cookbook recipe
import numpy as np
import pandas as pd

# generate random data
X = pd.Series([(x  + 5*np.random.random())%360       for x in range(-100, 600, 15)])
Y = pd.Series([(y  + 5*np.random.random())%360 - 180 for y in range(-200, 500, 15)])

# 'unwrap' the angles so there is no wrap around
X1 = pd.Series(np.rad2deg(np.unwrap(np.deg2rad(Y))))
Y1 = pd.Series(np.rad2deg(np.unwrap(np.deg2rad(Y))))

# smooth the data with a moving average
# note: this is pandas 17.1, the api changed for version 18
X2 = pd.rolling_mean(X1, window=3)
Y2 = pd.rolling_mean(Y1, window=3)

# convert back to wrapped data if desired
X3 = X2 % 360
Y3 = (Y2 + 180)%360 - 180

假设你从

开始
import numpy as np

x = np.linspace(0, 360, 360)
y = 5 * np.sin(x / 90. * 3.14) + np.random.randn(360)

plot(x, y, '+');

要执行循环卷积,您可以执行以下操作:

yy = np.concatenate((y, y))
smoothed = np.convolve(np.array([1] * 5), yy)[5: len(x) + 5]

这在每个点都使用前 5 个点(含)的循环平均值。当然,还有其他方法。

>>> plot(x, smoothed)

您可以使用 scipy.signal 中的 convolve2D。这是一个函数,它将平滑应用于 numpy 数组 a。如果 a 有不止一个维度,平滑应用于最里面(最快)的维度。

import numpy as np
from scipy import signal

def cyclic_moving_av( a, n= 3, win_type= 'boxcar' ):
  window= signal.get_window( win_type, n, fftbins=False ).reshape( (1,n) )
  shp_a= a.shape
  b= signal.convolve2d( a.reshape( ( np.prod( shp_a[:-1], dtype=int ), shp_a[-1] ) ), 
                        window, boundary='wrap', mode='same' )
  return ( b / np.sum( window ) ).reshape( shp_a )

例如可以这样使用

import matplotlib.pyplot as plt

x = np.linspace(0, 360, 360)
y1 = 5 * np.sin(x / 90. * 3.14) + 0.5 * np.random.randn(360)
y2 = 5 * np.cos(0.8 * x / 90. * 3.14) + 0.5 * np.random.randn(360)

y_av=  cyclic_moving_av( np.stack((y1,y2)), n=10 )  #1

plt.plot(x, y1, '+')
plt.plot(x, y2, '+')
plt.plot(x, y_av[0])
plt.plot(x, y_av[1])
plt.show()

这导致

#1等同于

y_av[0]=  cyclic_moving_av( y1, n=10 )
y_av[1]=  cyclic_moving_av( y2, n=10 )

win_type= 'boxcar' 对具有相同权重的邻居进行平均。有关其他选项,请参阅 signal.get_window