如何对 R 中的数据点运行高通或低通滤波器?
我是 R 的初学者,我尝试查找有关以下内容的信息,但没有找到任何内容。
图中的绿色图形是由红色和黄色图形组成的。但假设我只有绿色图之类的数据点。如何使用 低通< 提取低频/高频(即大约红色/黄色图表) /a>/高通滤波器?
更新:该图是通过
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
应该具有什么值的提示,但是 原始八度文档至少提到弧度 这让我去。我的原始图表中的值没有考虑到任何特定的频率,因此我最终得到了以下不那么简单的频率: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)。
I am a beginner in R and I have tried to find information about the following without finding anything.
The green graph in the picture is composed by the red and yellow graphs. But let's say that I only have the data points of something like the green graph. How do I extract the low/high frequencies (i.e. approximately the red/yellow graphs) using a low pass/high pass filter?
Update: The graph was generated with
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)
Update 2: Using the Butterworth filter in the signal
package suggested I get the following:
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)
The calculations was a bit work, signal.pdf gave next to no hints about what values W
should have, but the original octave documentation at least mentioned radians which got me going. The values in my original graph was not chosen with any specific frequency in mind, so I ended up with the following not so simple frequencies: f_low = 1/500 * 2 = 1/250
, f_high = 1/500 * 2*10 = 1/25
and the sampling frequency f_s = 500/500 = 1
. Then I chose a f_c somewhere inbetween the low and high frequencies for the low/high pass filters (1/100 and 1/50 respectively).
如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱,如果其他用户或者作者回复了您的评论,将不能在第一时间通知您!
发布评论
评论(8)
我不确定是否有任何过滤器最适合您。实现这一目标更有用的工具是快速傅里叶变换。
I am not sure if any filter is the best way for You. More useful instrument for that aim is the fast Fourier transformation.
我最近遇到了类似的问题,但没有发现这里的答案特别有帮助。这是一种替代方法。
让我们首先定义问题中的示例数据:
所以绿线是我们想要的数据集到低通和高通滤波器。
旁注:这种情况下的直线可以使用三次样条函数 (
spline(x,y, n = length(x))
) 表示为函数,但对于现实世界的数据,这很少会出现这种情况,所以我们假设不可能将数据集表示为函数。我遇到的平滑此类数据的最简单方法是使用
loess
或smooth.spline
以及适当的span
/spar
代码>.根据统计学家的说法,loess/smooth.spline 可能不是正确的方法,因为它并没有真正呈现这个意义上的数据的定义模型。另一种方法是使用广义加法模型(mgcv 包中的gam()
函数)。我在这里使用黄土或平滑样条的理由是,它更容易,并且没有什么区别,因为我们对可见的结果模式感兴趣。现实世界的数据集比本例中的数据集更复杂,并且找到用于过滤多个相似数据集的定义函数可能很困难。如果可见拟合良好,为什么要使用 R2 和 p 值使其变得更复杂?对我来说,该应用程序是可视化的,黄土/平滑样条是合适的方法。这两种方法都假设多项式关系,不同之处在于黄土也使用更高次数的多项式更灵活,而三次样条始终是三次 (x^2)。使用哪一种取决于数据集中的趋势。也就是说,下一步是使用loess()
或smooth.spline()
对数据集应用低通滤波器:红线是平滑样条滤波器,蓝色是黄土滤波器。如您所见,结果略有不同。我想使用 GAM 的一个论点是找到最佳拟合,如果数据集之间的趋势确实如此清晰且一致,但对于此应用程序,这两种拟合对我来说都足够好了。
找到合适的低通滤波器后,高通滤波就像从
y
中减去低通滤波值一样简单:这个答案来晚了,但我希望它可以帮助其他遇到类似问题的人。
I bumped into similar problem recently and did not find the answers here particularly helpful. Here is an alternative approach.
Let´s start by defining the example data from the question:
So the green line is the dataset we want to low-pass and high-pass filter.
Side note: The line in this case could be expressed as a function by using cubic spline (
spline(x,y, n = length(x))
), but with real world data this would rarely be the case, so let's assume that it is not possible to express the dataset as a function.The easiest way to smooth such data I have came across is to use
loess
orsmooth.spline
with appropriatespan
/spar
. According to statisticians loess/smooth.spline is probably not the right approach here, as it does not really present a defined model of the data in that sense. An alternative is to use Generalized Additive Models (gam()
function from package mgcv). My argument for using loess or smoothed spline here is that it is easier and does not make a difference as we are interested in the visible resulting pattern. Real world datasets are more complicated than in this example and finding a defined function for filtering several similar datasets might be difficult. If the visible fit is good, why to make it more complicated with R2 and p values? To me the application is visual for which loess/smoothed splines are appropriate methods. Both of the methods assume polynomial relationships with the difference that loess is more flexible also using higher degree polynomials, while cubic spline is always cubic (x^2). Which one to use depends on trends in a dataset. That said, the next step is to apply a low-pass filter on the dataset by usingloess()
orsmooth.spline()
:Red line is the smoothed spline filter and blue the loess filter. As you see results differ slightly. I guess one argument of using GAM would be to find the best fit, if the trends really were this clear and consistent among datasets, but for this application both of these fits are good enough for me.
After finding a fitting low-pass filter, the high-pass filtering is as simple as subtracting the low-pass filtered values from
y
:This answer comes late, but I hope it helps someone else struggling with similar problem.
使用 filtfilt 函数代替滤波器(封装信号)来消除信号偏移。
Use filtfilt function instead of filter (package signal) to get rid of signal shift.
一种方法是使用 R 中实现的
快速傅立叶变换
作为fft
。这是高通滤波器的示例。从上图中,本示例实现的想法是从绿色系列(您的真实数据)开始获取黄色系列。One method is using the
fast fourier transform
implemented in R asfft
. Here is an example of a high pass filter. From the plots above, the idea implemented in this example is to get the serie in yellow starting from the serie in green (your real data).根据OP的请求:
信号包包含各种过滤器用于信号处理。其中大部分与Matlab/Octave中的信号处理功能相当/兼容。
Per request of OP:
The signal package contains all kinds of filters for signal processing. Most of it is comparable to / compatible with the signal processing functions in Matlab/Octave.
查看此链接,其中有用于过滤(医疗信号)的 R 代码。作者是 Matt Shotwell,该网站充满了有趣的 R/stats 信息以及医学倾向:
biostattmat.com
fftfilt 包包含许多过滤算法,应该也有帮助。
Check out this link where there's R code for filtering (medical signals). It's by Matt Shotwell and the site is full of interesting R/stats info with a medical bent:
biostattmat.com
The package fftfilt contains lots of filtering algorithms that should help too.
我还努力弄清楚黄油函数中的 W 参数如何映射到过滤器截止,部分原因是过滤器和 filtfilt 的文档在发布时不正确(它表明 W = .1 将导致 10当信号采样率 Fs = 100 时,与 filtfilt 结合使用 Hz lp 滤波器,但实际上,它只是一个 5 Hz lp 滤波器——使用 filtfilt 时半幅截止为 5 Hz,但当您仅应用一次滤波器(使用滤波器功能)时,半功率截止为 5 Hz。我发布了一些我在下面编写的演示代码,这些代码帮助我确认这一切是如何工作的,并且您可以使用它来检查过滤器是否正在执行您想要的操作。
I also struggled to figure out how the W parameter in the butter function maps on to the filter cut-off, in part because the documentation for filter and filtfilt is incorrect as of posting (it suggests that W = .1 would result in a 10 Hz lp filter when combined with filtfilt when signal sampling rate Fs = 100, but actually, it's only a 5 Hz lp filter -- the half-amplitude cut-off is 5 Hz when use filtfilt, but the half-power cut-off is 5 Hz when you only apply the filter once, using the filter function). I'm posting some demo code I wrote below that helped me confirm how this is all working, and that you could use to check a filter is doing what you want.
CRAN 上有一个名为 FastICA 的包,它计算独立源信号的近似值,但是为了计算这两个信号,您需要一个至少包含 2xn 混合观测值的矩阵(对于本例),这算法无法仅用 1xn 向量确定两个独立信号。请参阅下面的示例。希望这可以帮助你。
there is a package on CRAN named
FastICA
, this computes the approximation of the independent source signals, however in order to compute both signals you need a matrix of at least 2xn mixed observations (for this example), this algorithm can't determine the two indpendent signals with just 1xn vector. See the example below. hope this can help you.