如何在R中的数据点上运行高通或低通滤波器？

Question

我是R的初学者,我试图在没有找到任何内容的情况下找到有关以下内容的信息.

图中的绿色图由红色和黄色图组成.但是,假设我只有绿色图形的数据点.如何使用低通/高通滤波器提取低/高频率(即大约红/黄图)？

更新:图表是使用生成的

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包中的Butterworth过滤器建议我得到以下内容:

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应该有什么值的提示,但原始的八度文档至少提到了让我走的弧度.在我原来的图中的值是没有考虑任何特定的频率选择,所以我结束了以下并非如此简单频率: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/ spar.根据统计学家的说法,loess/smooth.spline可能不是正确的方法,因为在这个意义上它并没有真正呈现数据的定义模型.另一种方法是使用广义加法模型(gam()函数来自包mgcv).我在这里使用黄土或平滑样条曲线的论点是,它更容易并且没有区别,因为我们对可见的结果模式感兴趣.真实世界数据集比此示例中更复杂,并且找到用于过滤若干类似数据集的已定义函数可能是困难的.如果可见拟合良好,为什么要使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)

这个答案来得很晚,但我希望它可以帮助其他人解决类似的问题.

Answer 2

使用filtfilt函数代替滤波器(封装信号)来摆脱信号转换.

library(signal)
bf <- butter(2, 1/50, type="low")
b1 <- filtfilt(bf, y+noise1)
points(x, b1, col="red", pch=20)

Answer 3

根据OP的要求:

该信号包包含一个有信号处理各种滤波器.其中大部分与Matlab/Octave中的信号处理功能相当/兼容.

Answer 4

一种方法是使用fast fourier transformR中实现的fft.以下是高通滤波器的示例.从上面的图中可以看出,本例中实​​现的想法是从绿色系列(您的真实数据)开始将系列设为黄色.

# I've changed the data a bit so it's easier to see in the plots
par(mfrow = c(1, 1))
number_of_cycles = 2
max_y = 40
N <- 256

x = 0:(N-1)
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)

### Apply the fft to the noisy data
y_noise = y + noise1
fft.y_noise = fft(y_noise)


# Plot the series and spectrum
par(mfrow = c(1, 2))
plot(x, y_noise, type='l', main='original serie', col='green4')
plot(Mod(fft.y_noise), type='l', main='Raw serie - fft spectrum')

### The following code removes the first spike in the spectrum
### This would be the high pass filter
inx_filter = 15
FDfilter = rep(1, N)
FDfilter[1:inx_filter] = 0
FDfilter[(N-inx_filter):N] = 0
fft.y_noise_filtered = FDfilter * fft.y_noise

par(mfrow = c(2, 1))
plot(x, noise1, type='l', main='original noise')
plot(x, y=Re( fft( fft.y_noise_filtered, inverse=TRUE) / N ) , type='l', 
     main = 'filtered noise')