如何在R中对时间序列进行频率滤波并返回时域结果

基础疑问解答 #

• 原生fft的逆变换不需要额外函数，只需调用fft(fft_result, inverse = TRUE)，输出后取实部并除以序列长度归一化，即可得到时域结果
• 原signal包维护更新较慢，推荐使用功能更稳定的替代包gsignal，语法几乎兼容，运行结果更可靠

需求匹配说明 #

针对AirPassengers月度数据的季节成分提取需求，带通滤波的参数设置逻辑：

• 月度数据的季节周期为12个月，对应频率为1/12 ≈ 0.083 周期/月
• 带通范围设置为0.07~0.09 周期/月，刚好保留12个月周期的季节波动，同时剔除低于该范围的趋势项、高于该范围的高频噪声

方案1：基于gsignal包的封装实现（推荐入门使用） #

步骤1：安装加载依赖包 #

# 安装包（首次运行执行）
install.packages("gsignal")
# 加载包
library(gsignal)

步骤2：数据预处理与滤波参数设置 #

# 取对数后的目标序列
y <- log(AirPassengers)
# 序列采样频率：月度数据为1次/月
fs <- 1
# 带通上下限（单位：周期/月）
low <- 0.07
high <- 0.09
# 设计4阶巴特沃斯带通滤波器
bf <- butter(4, c(low, high), type = "pass", fs = fs)
# 零相位滤波，避免相位偏移
seasonal_component <- filtfilt(bf, y)

步骤3：结果可视化验证 #

# 绘制原始对数序列
plot(log(AirPassengers), main = "对数AirPassengers序列与提取的季节成分")
# 叠加提取的季节成分，叠加均值偏移方便观测对比
lines(time(AirPassengers), seasonal_component + mean(log(AirPassengers)), col = "red", lwd = 2)
legend("topleft", legend = c("原始序列", "提取的季节成分"), col = c("black", "red"), lwd = c(1,2))

方案2：手动实现FFT→滤波→逆FFT流程 #

如果需要完全控制傅里叶变换的每一步，可以用原生fft函数实现：

y <- log(AirPassengers)
n <- length(y)
low <- 0.07
high <- 0.09
# 1. 执行傅里叶变换
y_fft <- fft(y)
# 2. 生成频率轴
freq <- seq(0, (n-1)/n, length.out = n)
# 3. 构建滤波掩码：实序列FFT结果共轭对称，需同时保留对称的高低频段
mask <- (freq >= low & freq <= high) | (freq >= (1-high) & freq <= (1-low))
# 4. 对频谱做滤波
y_fft_filtered <- y_fft * mask
# 5. 逆傅里叶变换，取实部并归一化得到时域结果
seasonal_component_manual <- Re(fft(y_fft_filtered, inverse = TRUE)) / n

该方法输出的结果与封装包结果基本一致，适合需要自定义滤波规则的场景。

内容的提问来源于stack exchange，提问作者Sebastian