在R中对时间序列执行FFT：拟合谐波、绘图及数据点预测

嘿，我看你想用FFT来拟合时间序列的波形，还要绘制不同谐波的图像，再预测后续n个数据点——刚好你用的是@catastrophic-failure的代码，我帮你把这个思路补全、解释明白，让你能直接用起来：

完整实现函数 #

先把补全后的完整代码放出来，这个函数包含了拟合、预测和绘图的功能：

nff = function(y = NULL, n = NULL, up = 10L, plot = TRUE, add = FALSE, main = NULL, ...){
  # 对原始时间序列执行FFT变换
  # 注意：FFT结果第一个元素是直流分量，其余频率分量是对称重复的
  dff = fft(y)
  
  # 原始时间序列的时间点索引
  t = seq(from = 1, to = length(y))
  
  # 构建包含预测点的升采样时间轴：up控制平滑度，n是要预测的点数
  nt = seq(from = 1, to = length(y) + n, length.out = length(y)*up + n*up)
  
  # 处理FFT系数：只保留前半部分（后半部分是共轭对称的）
  len_y = length(y)
  half_len = ceiling(len_y/2)
  
  # 初始化升采样后的FFT系数数组
  dff_up = numeric(length(nt))
  # 复制前半部分的正频率系数
  dff_up[1:half_len] = dff[1:half_len]
  
  # 偶数长度的序列需要单独处理中间的频率分量
  if(len_y %% 2 == 0){
    dff_up[length(nt) - half_len + 2] = dff[half_len + 1]
  }
  
  # 填充后半部分的共轭对称系数
  dff_up[(length(nt) - half_len + 2):length(nt)] = Conj(rev(dff[2:half_len]))
  
  # 逆FFT变换得到拟合/预测序列，取实部并归一化
  y_fit = Re(fft(dff_up, inverse = TRUE)) / length(dff_up)
  
  # 绘制原始数据和拟合曲线
  if(plot){
    if(!add){
      plot(t, y, type = "p", main = ifelse(is.null(main), "FFT拟合与预测结果", main), ...)
    }
    lines(nt, y_fit, col = "red", lwd = 2)
  }
  
  # 返回所有结果供后续分析
  list(original = y, fitted = y_fit, time_original = t, time_fitted = nt, fft_coeff = dff)
}

关键步骤拆解 #

• FFT核心变换：fft(y)把时域的时间序列转换成频域的系数，第一个系数是直流分量（信号的平均水平），后面的系数对应不同频率的谐波成分，且正负频率是对称的。
• 升采样时间轴：nt不仅包含原始数据的时间点，还扩展了n个预测点的时间范围，up参数越大，拟合出来的曲线越平滑，适合观察波形细节。
• 系数对称处理：因为FFT的结果是对称的，我们只需要保留前半部分的有效系数，再通过共轭对称填充后半部分，这样逆变换后就能得到平滑的拟合波形。
• 逆FFT还原：逆变换后取实部（因为原始信号是实数），再除以长度做归一化，就能得到和原始信号同尺度的拟合/预测序列。

绘制单个谐波的波形 #

如果想单独看每个谐波对原始信号的贡献，可以用这个小工具函数：

# 绘制指定次数的谐波波形
plot_single_harmonic = function(y, harmonic_idx){
  len_y = length(y)
  dff = fft(y)
  
  # 初始化系数数组，只保留目标谐波的分量
  dff_single = numeric(len_y)
  dff_single[harmonic_idx + 1] = dff[harmonic_idx + 1]
  
  # 非直流分量需要保留对应的负频率共轭系数
  if(harmonic_idx != 0){
    dff_single[len_y - harmonic_idx + 1] = dff[len_y - harmonic_idx + 1]
  }
  
  # 逆变换得到该谐波的时域波形
  harmonic_wave = Re(fft(dff_single, inverse = TRUE)) / len_y
  t = seq(1, len_y)
  
  # 绘图
  plot(t, harmonic_wave, type = "l", main = paste("第", harmonic_idx, "次谐波"), 
       col = "steelblue", lwd = 2, xlab = "时间点", ylab = "幅值")
}

# 示例：绘制直流分量（0次）、基波（1次）、2次谐波
plot_single_harmonic(y, 0)
plot_single_harmonic(y, 1)
plot_single_harmonic(y, 2)

• 0次谐波就是信号的平均值，是一条水平线
• 1次谐波是信号的基波，对应原始信号的基本周期
• 更高次谐波是叠加在基波上的高频“波纹”，这些谐波组合起来就构成了原始信号的复杂波形

用函数预测n个数据点 #

调用nff函数时，只需要指定n参数就能得到后续n个点的预测结果：

# 先生成一个测试用的时间序列（带噪声的正弦波）
set.seed(123)
y = sin(seq(1, 20, by = 0.5)) + rnorm(40, 0, 0.1)

# 预测后续10个点，设置up=10保证曲线平滑
result = nff(y, n = 10, up = 10, main = "FFT拟合与10点预测", xlab = "时间点", ylab = "幅值")

# 查看最后10个预测值
tail(result$fitted, 10)

运行后，你会看到黑色散点是原始数据，红色曲线是拟合和预测的结果——红色曲线的后半段就是你要的n个预测点。

内容的提问来源于stack exchange，提问作者m.inam