r - 如何对 R 中的数据点运行高通或低通滤波器？

Question

我是 R 的初学者我试图查找有关以下内容的信息，但没有找到任何内容。

图中的绿色图形是由红色和黄色图形组成的。但假设我只有绿色图之类的数据点。如何使用 low pass 提取低频/高频(即大约红色/黄色图表)/high pass filter ？

更新:图表是用

生成的

number_of_cycles = 2
max_y = 40

x = 1:500
a = number_of_cycles * 2*pi/length(x)

y = max_y * sin(x*a)
noise1 = max_y * 1/10 * sin(x*a*10)

plot(x, y, type="l", col="red", ylim=range(-1.5*max_y,1.5*max_y,5))
points(x, y + noise1, col="green", pch=20)
points(x, noise1, col="yellow", pch=20)

更新 2:使用 signal 包中的巴特沃斯滤波器建议我得到以下结果:

library(signal)

bf <- butter(2, 1/50, type="low")
b <- filter(bf, y+noise1)
points(x, b, col="black", pch=20)

bf <- butter(2, 1/25, type="high")
b <- filter(bf, y+noise1)
points(x, b, col="black", pch=20)

计算有点麻烦，signal.pdf 几乎没有提示 W 应该具有什么值，但 original octave documentation至少提到radians这让我继续前进。我的原始图表中的值没有考虑到任何特定的频率，因此我最终得到了以下不那么简单的频率:f_low = 1/500 * 2 = 1/250， f_high = 1/500 * 2*10 = 1/25 和采样频率 f_s = 500/500 = 1。然后我为低/高通滤波器选择了低频和高频之间的 f_c(分别为 1/100 和 1/50)。

Answer 1

我最近遇到了类似的问题，但没有发现这里的答案特别有帮助。这是一种替代方法。

让我们首先定义问题中的示例数据:

number_of_cycles = 2
max_y = 40

x = 1:500
a = number_of_cycles * 2*pi/length(x)

y = max_y * sin(x*a)
noise1 = max_y * 1/10 * sin(x*a*10)
y <- y + noise1

plot(x, y, type="l", ylim=range(-1.5*max_y,1.5*max_y,5), lwd = 5, col = "green")

所以绿线是我们想要低通和高通滤波的数据集。

旁注:本例中的直线可以使用三次样条函数 (spline(x,y, n = length(x))) 表示为函数，但对于现实世界的数据，此这种情况很少发生，所以我们假设不可能将数据集表示为函数。

我遇到的平滑此类数据的最简单方法是使用 loess 或 smooth.spline 以及适当的 span/晶石。据统计学家loess/smooth.spline is probably not the right approach here ，因为它并没有真正呈现这个意义上的数据的定义模型。另一种方法是使用广义加法模型(mgcv 包中的gam() 函数)。我在这里使用黄土或平滑样条的理由是，它更容易，并且不会产生任何影响，因为我们对可见的结果模式感兴趣。现实世界的数据集比本例中的数据集更复杂，并且找到用于过滤多个相似数据集的定义函数可能很困难。如果可见拟合良好，为什么要使用 R2 和 p 值使其变得更复杂？对我来说，该应用程序是可视化的，黄土/平滑样条是合适的方法。这两种方法都假设多项式关系，不同之处在于黄土使用更高次数的多项式更灵活，而三次样条始终是三次 (x^2)。使用哪一种取决于数据集中的趋势。也就是说，下一步是使用 loess() 或 smooth.spline() 对数据集应用低通滤波器:

lowpass.spline <- smooth.spline(x,y, spar = 0.6) ## Control spar for amount of smoothing
lowpass.loess <- loess(y ~ x, data = data.frame(x = x, y = y), span = 0.3) ## control span to define the amount of smoothing

lines(predict(lowpass.spline, x), col = "red", lwd = 2)
lines(predict(lowpass.loess, x), col = "blue", lwd = 2)

红线是平滑样条滤波器，蓝色是黄土滤波器。如您所见，结果略有不同。我猜想使用 GAM 的一个论点是找到最佳拟合(如果数据集之间的趋势确实如此清晰且一致)，但对于此应用程序，这两种拟合对我来说都足够好了。

找到合适的低通滤波器后，高通滤波就像从y中减去低通滤波值一样简单:

highpass <- y - predict(lowpass.loess, x)
lines(x, highpass, lwd =  2)

这个答案来晚了，但我希望它可以帮助其他遇到类似问题的人。