43

我是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应该有什么值的提示，但最初的 octave 文档至少提到了让我前进的弧度。我的原始图表中的值没有选择任何特定的频率，所以我最终得到了以下不那么简单的频率：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。

4

8 回答 8

31

我最近遇到了类似的问题，并没有发现这里的答案特别有用。这是另一种方法。

让我们从定义问题中的示例数据开始：

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 值使其更复杂？对我来说，该应用程序是可视化的，其中黄土/平滑样条曲线是合适的方法。这两种方法都假设多项式关系，不同之处在于 loess 也使用更高次多项式更灵活，而三次样条始终是三次 (x^2)。使用哪一个取决于数据集中的趋势。也就是说，下一步是使用loess()or对数据集应用低通滤波器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)

在此处输入图像描述

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

于 2014-04-08T12:31:39.957 回答
17

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

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

红线显示 filtfilt 的结果

于 2014-10-20T06:43:08.863 回答
8

一种方法是使用fast fourier transformR 中实现的 as 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')

y 噪声和 fft 频谱

### 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')

在此处输入图像描述

于 2016-06-24T04:06:30.027 回答
7

根据 OP 的要求：

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

于 2011-08-19T08:26:11.387 回答
3

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

biostattmat.com

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

于 2011-08-18T10:35:33.270 回答
3

我还努力弄清楚黄油函数中的 W 参数如何映射到过滤器截止值，部分原因是过滤器和 filtfilt 的文档在发布时不正确（它表明 W = .1 将导致 10信号采样率 Fs = 100 时与 filtfilt 组合使用的 Hz lp 滤波器，但实际上，它只是一个 5 Hz 的 lp 滤波器——使用 filtfilt 时半幅度截止为 5 Hz，但半功率截止为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
于 2018-04-18T17:08:12.497 回答
2

FastICACRAN上有一个包名为只有 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 = "")
于 2014-11-08T21:55:25.473 回答
1

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

于 2013-03-01T12:43:30.063 回答