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

【问题标题】：How do I run a high pass or low pass filter on data points in R?如何对 R 中的数据点运行高通或低通滤波器？
【发布时间】：2011-10-29 16:13:34
【问题描述】：

我是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 包中的 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 应该具有什么值，但original octave documentation 至少提到了radians，这让我继续前进。我的原始图表中的值没有选择任何特定的频率，所以我最终得到了以下不那么简单的频率：f_low = 1/500 * 2 = 1/250、f_high = 1/500 * 2*10 = 1/25 和采样频率f_s = 500/500 = 1。然后我为低通/高通滤波器（分别为 1/100 和 1/50）选择了介于低频和高频之间的 f_c。

【问题讨论】：

如果您给我们reproducible example，例如您用于图表的数据/代码，人们将能够更轻松地为您提供帮助。这将有助于向我们展示您到目前为止所做的尝试。
添加：信号包包含所有类型的过滤器：cran.r-project.org/web/packages/signal/signal.pdf
无论如何，这都是一个过于宽泛的编程问题。您至少应该指定要使用的过滤器。有很多选项可能对您的真实数据有意义，也可能没有意义。
@Joris 请把你对信号的评论变成答案，我会接受的。这正是我一直在寻找的东西（尽管我发现我必须对多年前所学的过滤器进行大量的重新学习......）。
我认为你已经错过了傅里叶分析的整个领域。正确应用的分析应该能够提取出只有两个正弦信号的事实。

标签： r signal-processing frequency-analysis

【解决方案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 is probably not the right approach here 的说法，因为它并没有真正呈现这种意义上的数据定义模型。另一种方法是使用广义加法模型（mgcv 包中的gam() 函数）。我在这里使用黄土或平滑样条曲线的论点是它更容易并且不会产生影响，因为我们对可见的结果模式感兴趣。现实世界的数据集比这个例子中的更复杂，找到一个定义的函数来过滤几个相似的数据集可能很困难。如果可见拟合良好，为什么要使用 R2 和 p 值使其更复杂？对我来说，该应用程序是可视化的，其中黄土/平滑样条曲线是合适的方法。这两种方法都假设多项式关系，不同之处在于 loess 也使用更高次多项式更灵活，而三次样条始终是三次 (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)

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

【讨论】：

谢谢，很好的回答。下次遇到这种问题我会记住的。

【解决方案2】：

使用 filtfilt 函数代替过滤器（封装信号）来消除信号偏移。

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

【讨论】：

小心使用这个函数，因为在它的文档中有“......所以这个函数还需要一些工作 - 并且处于 2000 年版本的 Octave 代码的状态。”跨度>

【解决方案3】：

一种方法是使用在 R 中实现的fast fourier transform 作为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')

【讨论】：

【解决方案4】：

根据 OP 的请求：

signal package 包含用于信号处理的各种过滤器。大部分与Matlab/Octave中的信号处理功能相当/兼容。

【讨论】：

【解决方案5】：

查看此链接，其中有用于过滤的 R 代码（医疗信号）。由 Matt Shotwell 撰写，该网站充满了有趣的 R/stats 信息，具有医学倾向：

biostattmat.com

fftfilt 包包含许多过滤算法，它们也应该有所帮助。

【讨论】：

有一个包。复制一个非常基本的过滤器的手动实现，您不知道它是否会实际执行，这不是一个好主意。

【解决方案6】：

我还努力弄清楚黄油函数中的 W 参数如何映射到过滤器截止值，部分原因是过滤器和 filtfilt 的文档在发布时不正确（它表明 W = .1 会导致在信号采样率 Fs = 100 时与 filtfilt 组合使用 10 Hz lp 滤波器，但实际上，它只是一个 5 Hz lp 滤波器 - 使用 filtfilt 时半幅度截止为 5 Hz，但半功率截止-off 是 5 Hz，当您只应用一次过滤器时，使用过滤器功能）。我将发布一些我在下面编写的演示代码，帮助我确认这一切是如何工作的，并且您可以使用它来检查过滤器是否正在执行您想要的操作。

#Example usage of butter, filter, and filtfilt functions
#adapted from https://rdrr.io/cran/signal/man/filtfilt.html

library(signal)

Fs <- 100; #sampling rate

bf <- butter(3, 0.1);       
#when apply twice with filtfilt, 
#results in a 0 phase shift 
#5 Hz half-amplitude cut-off LP filter
#
#W * (Fs/2) == half-amplitude cut-off when combined with filtfilt
#
#when apply only one time, using the filter function (non-zero phase shift),
#W * (Fs/2) == half-power cut-off


t <- seq(0, .99, len = 100)   # 1 second sample

#generate a 5 Hz sine wave
x <- sin(2*pi*t*5)

#filter it with filtfilt
y <- filtfilt(bf, x)

#filter it with filter
z <- filter(bf, x)

#plot original and filtered signals
plot(t, x, type='l')
lines(t, y, col="red")
lines(t,z,col="blue")

#estimate signal attenuation (proportional reduction in signal amplitude)
1 - mean(abs(range(y[t > .2 & t < .8]))) #~50% attenuation at 5 Hz using filtfilt

1 - mean(abs(range(z[t > .2 & t < .8]))) #~30% attenuation at 5 Hz using filter

#demonstration that half-amplitude cut-off is 6 Hz when apply filter only once
x6hz <- sin(2*pi*t*6)

z6hz <- filter(bf, x6hz)

1 - mean(abs(range(z6hz[t > .2 & t < .8]))) #~50% attenuation at 6 Hz using filter


#plot the filter attenuation profile (for when apply one time, as with "filter" function):

hf <- freqz(bf, Fs = Fs);

plot(c(0, 20, 20, 0, 0), c(0, 0, 1, 1, 0), type = "l", 
 xlab = "Frequency (Hz)", ylab = "Attenuation (abs)")

lines(hf$f[hf$f<=20], abs(hf$h)[hf$f<=20])

plot(c(0, 20, 20, 0, 0), c(0, 0, -50, -50, 0),
 type = "l", xlab = "Frequency (Hz)", ylab = "Attenuation (dB)")

lines(hf$f[hf$f<=20], 20*log10(abs(hf$h))[hf$f<=20])

hf$f[which(abs(hf$h) - .5 < .001)[1]] #half-amplitude cutoff, around 6 Hz

hf$f[which(20*log10(abs(hf$h))+6 < .2)[1]] #half-amplitude cutoff, around 6 Hz

hf$f[which(20*log10(abs(hf$h))+3 < .2)[1]] #half-power cutoff, around 5 Hz

【讨论】：

【解决方案7】：

CRAN 上有一个名为 FastICA 的包，它计算独立源信号的近似值，但是为了计算这两个信号，您需要一个至少包含 2xn 混合观察的矩阵（对于本示例），该算法可以't 仅用 1xn 向量确定两个独立信号。请参见下面的示例。希望这可以帮到你。

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)
######################################################
library(fastICA)
S <- cbind(y,noise1)#Assuming that "y" source1 and "noise1" is source2
A <- matrix(c(0.291, 0.6557, -0.5439, 0.5572), 2, 2) #This is a mixing matrix
X <- S %*% A 

a <- fastICA(X, 2, alg.typ = "parallel", fun = "logcosh", alpha = 1,
method = "R", row.norm = FALSE, maxit = 200,
tol = 0.0001, verbose = TRUE)

par(mfcol = c(2, 3))
plot(S[,1 ], type = "l", main = "Original Signals",
xlab = "", ylab = "")
plot(S[,2 ], type = "l", xlab = "", ylab = "")
plot(X[,1 ], type = "l", main = "Mixed Signals",
xlab = "", ylab = "")
plot(X[,2 ], type = "l", xlab = "", ylab = "")
plot(a$S[,1 ], type = "l", main = "ICA source estimates",
xlab = "", ylab = "")
plot(a$S[, 2], type = "l", xlab = "", ylab = "")

【讨论】：

【解决方案8】：

我不确定是否有任何过滤器是最适合您的方式。更有用的工具是快速傅里叶变换。

【讨论】：