基于R语言的低通滤波算法

参考链接：CRAN Task View: Time Series Analysis，bwfilter: Butterworth filter of a time series in mFilter: Miscellaneous Time Series Filters

版本相关信息：

Maintainer:	Rob J Hyndman
Contact:	Rob.Hyndman at monash.edu
Version:	2021-12-20
URL:	https://CRAN.R-project.org/view=TimeSeries

属于R语言的时间序列分析

原文：

• Filters and smoothing : filter() in stats provides autoregressive and moving average linear filtering of multiple univariate time series. The robfilter package provides several robust time series filters. smooth() from the stats package computes Tukey's running median smoothers, 3RS3R, 3RSS, 3R, etc. sleekts computes the 4253H twice smoothing method. mFilter implements several filters for smoothing and extracting trend and cyclical components including Hodrick-Prescott and Butterworth filters. smoots provides nonparametric estimation of the time trend and its derivatives.

滤波器和平滑：统计中的filter()函数提供了多个单变量时间序列的自回归和移动平均线性滤波。robfilter包提供了几个强有力的时间序列滤波器。来自stats package的smooth()函数计算Tukey的运行中位数平滑，3RS3R，3RSS，3R等。sleekts(笔者按：这里不确定是平滑函数还是平滑包)计算4253H两次平滑方法。mFilter实现了几个过滤器来平滑和提取趋势和循环组件，包括霍德里克-普雷斯科特和巴特沃斯滤波器。smoots(笔者按：这里不确定是平滑函数还是平滑包)提供了时间趋势及其导数的非参数估计。

由此可知，第一优选包为mFilter，点击链接中mFilter发现包链接CRAN - Package mFilter，证明此包可下载

install.packages("mFilter")并下载其参考说明书，发现其中的BWFILTER（巴特沃斯滤波器）命令如下所示：

bwfilter(x,freq=NULL,nfix=NULL,drift=FALSE)

这里我们用一个案例来说明bwfilter的应用方案

rm(list = ls())#清空工作区，可在括号内加装变量等定向清除某个变量等
library(mFilter)
library(TSA)#larain数据集所在包
data(larain)#注意，larain是一个时间序列，其输出图见封面图1

opar <- par(no.readonly=TRUE)
#函数par()在进行参数设置之前会自动的打开一个新绘图设备。
#par(no.readonly=TRUE)都可以获取当前的各个绘图参数。opar存储了当前的绘图参数
larain.bw <- bwfilter(larain)
plot(larain.bw)#在不设置截止频率和参数n的情况下，可以看到仅能滤除很少的噪声

#如封面图2所示

larain.bw1 <- bwfilter(larain, drift=TRUE)
larain.bw2 <- bwfilter(larain, freq=8,drift=TRUE)
larain.bw3 <- bwfilter(larain, freq=10, nfix=3, drift=TRUE)
larain.bw4 <- bwfilter(larain, freq=10, nfix=4, drift=TRUE)

par(mfrow=c(2,1),mar=c(3,3,2,1),cex=.8)
plot(larain.bw1$x,
     main="Butterworth filter of larainloyment: Trend,
     drift=TRUE",col=1, ylab="")
lines(larain.bw1$trend,col=2)
lines(larain.bw2$trend,col=3)
lines(larain.bw3$trend,col=4)
lines(larain.bw4$trend,col=5)
legend("topleft",legend=c("series", "freq=10, nfix=2",
                          "freq=8, nfix=2", "freq=10, nfix=3", "freq=10, nfix=4"),
       col=1:5, lty=rep(1,5), ncol=1)

#如封面图3所示，不同截止频率和阶数n的趋势分量，注意观察，n=2时，初始的数据滤波后失真##较小，n=3和4失真较大。这部分涉及到时间序列的分解，具体的说：（参考##时间序列分解算法——简述 - 知乎）

#对于一个时间序列{y(t)}，假设它是加性模型（an additive decomposition），则可以写成　

　　　　　　　　　　　　　　　#y(t)=S(t)+T(t)+R(t),

#其中S(t)、T(t)、R(t)分别是周期成分（seasonal component）、趋势成分（trend-cycle #component）、残差成分（remainder component）。这里巴特沃斯滤波器输出的就是趋势成分##和循环成分，是将T（t）进行进一步分解的结果。

#类似地，一个乘性模型可以写成

　　　　　　　　　　　　　　　#y(t)=S(t)×T(t)×R(t)，

#对于乘性模型，可以取对数（当然是有意义的前提下），将其转化为加性模型。

plot(larain.bw1$cycle,
     main="Butterworth filter of larainloyment: Cycle,drift=TRUE",
     col=2, ylab="", ylim=range(larain.bw3$cycle,na.rm=TRUE))
lines(larain.bw2$cycle,col=3)
lines(larain.bw3$cycle,col=4)
lines(larain.bw4$cycle,col=5)
## legend("topleft",legend=c("series", "freq=10, nfix=2", "freq=8,
## nfix=2", "freq## =10, nfix=3", "freq=10, nfix=4"), col=1:5,
## lty=rep(1,5), ncol=1)

#如封面图3所示，不同截止频率和阶数n的循环分量，注意观察，直接滤波bw1线为0，证明没有循环分量，n=2时，初始的数据滤波后失真##较小，n=3和4失真较大。由此可见，想要获得理想的滤波信号，需要对滤波参数反复修正才可以。

par(opar)