平滑的圆形数据答案

【问题标题】：Smooth circular data平滑的圆形数据
【发布时间】：2020-05-29 17:47:08
【问题描述】：

我有一个数据数组Y，这样Y 是自变量X（另一个数组）的函数。

X 中的值从 0 到 360 不等，带有环绕。

Y 中的值从 -180 到 180 不等，也有环绕。

（也就是说，这些值是以度为单位的圆的角度。）

有谁知道 Python 中的任何函数（在 numpy、scipy 等中）能够低通过滤我的 Y 值作为 X 的函数？

如果这完全令人困惑，这里有一个示例数据图：

【问题讨论】：

你想如何处理环绕？如果 X 从 359 度环绕。到001度，Y值是否应该与001度附近的其他数据一起平滑？还是应该将 Y 值视为 361 度？同样，你想如何处理 Y 的环绕？
问题是有两个数组 X 和 Y，其中 Y[i] 是 X[i] 的函数。我假设我代表时间。我的问题是如果 Y 有一个看起来像 [..., 179, 180, -179, -178, 177, ...] 的序列，从字面上看是否更有意义，或者像它包含 [... , 179, 180, 181, 182, 183, ...]。对于 X 也是如此。在下面的回答中，我在平滑数据之前“解包”数据。

标签： python numpy scipy filtering smooth

【解决方案1】：

假设你开始

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)

【讨论】：

请注意，这不处理 y 值中的环绕 - 但在这里没关系，因为它们永远不会获得接近环绕点 (+/- 180) 的值。

【解决方案2】：

这是一个使用 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

【讨论】：

我看不出平滑 X1 的意义。那是我的自变量。那里没有不确定性。上面，您将X 和Y 视为时间的函数，但这不是正确的思考方式。就像我说的，Y 是X 的函数。 X 是这里时间的替身。
附注感谢您指出np.unwrap。我没有意识到这一点。

【解决方案3】：

您可以从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。

【讨论】：