时序顶会基础创新知识点-傅立叶变换篇

前言

很长一段时间以来，我在读论文的时候经常看到时序研究中，会运用傅立叶变换做初步的处理，然后基于处理结果，进行后续的建模研究。又因为“傅立叶变换”这个名字太熟悉了，但是对于它的推理过程、应用步骤，每次似乎都没有彻底吃透，于是占据时间序列分析建模重要地位的傅立叶变换，似乎既熟悉又陌生。

我回顾了自己学习傅立叶变换的过程，觉得没能学好的原因大概在于两方面：一是自己数学基础确实比较弱，一到数学推理部分就比较吃力；二是大量的科普内容要么侧重于可视化展示，要么侧重于严谨的数学推导；三是应用场景并不集中在时间序列领域。但是对于我而言，我更希望的理解过程是从代码的角度，给出数据流向、算法处理过程、以及对结果的解释。打通这个过程，基本上我就知道算法该怎么用、得到什么结果，此外我希望从时序研究的角度理解和应用。

带着这样的思路，本篇文章试图总结以下几件事：

• 什么是傅立叶变换，为什么要（可以）做傅立叶变换？
• 如何基于python做傅立叶变换，得到的结果如何解读？
• 时序顶会论文是如何运用傅立叶变换增加创新点的？

什么是傅立叶变换？

傅里叶变换（Fourier Transform）能够将一个函数（在时间序列问题中，通常是离散的点）从时域转换到频域。它可以把一个复杂的随时间变化的信号分解成许多不同频率的正弦波和余弦波的组合。

其实这里就隐含了两个要点：“把随时间变化的函数（或离散点）转换为不同三角函数组合”，即从时域转换到频域。这里面值得思考的问题很多，比如：变换的本质是什么？为什么要使用三角函数？如果你要深入理解背后的数学原理，我推荐你看“喵星考拉”、“我们都来自星尘”和博客主“Yacan Man”的讲解（本文代码参考了以上作者）。但简而言之，可以理解为傅立叶变换就是一种坐标系变换，和我们在三维空间对其中的点，做坐标系变换并无本质区别。既然是坐标系变换，自然引出了基向量概念，首先看三维空间中的(x,y,z)基向量，一个最重要的性质就是正交，基坐标内积为0，而三角函数也恰恰有这样的性质，我个人理解：这种正交性是傅立叶变换使用正余弦函数作为基向量的数学基础。

此外，根据傅里叶级数，任何周期函数都可以用一系列正弦函数和余弦函数的无穷级数来表示。那么，对于非周期函数，怎么办的呢？可以将其看作是周期趋于无穷大的周期函数。当周期趋于无穷时，傅里叶级数就演变成了傅里叶变换。

为什么要做傅立叶变换？

我们把应用场景放到时间序列研究领域，我认为做傅立叶变换至少有三个好处：频率成分分析、滤波噪声、数据压缩，下面我们逐一来看。

频率成分分析

在时间序列研究中，许多时序数据包含复杂的周期性和非周期性成分。傅里叶变换可以将时间序列从时域转换到频域，从而清晰地揭示出数据中隐藏的频率成分。例如，在分析电力消耗的时间序列时，通过傅里叶变换可以发现日周期（由于人们的日常活动模式导致白天和夜晚用电量不同）和周周期（工作日和休息日用电模式差异）以及季度周期等周期性成分。通过把混合数据分解为季度+周+日不同的频率维度，有助于理解数据背后的规律。

下图右下角就是三个频率叠加后的结果，假如我们拿到了这个时序数据，通过傅立叶变换我们可以反向得到不同频率分量，这些不同频率分量就相当于我们上面说的日周期、周周期、季周期，这样可解释性明显提升了。下图是一个例子，前三个曲线合成为第四个曲线，反之可以分解出来。

从另一个方面看，原本在时域上被掩盖的特征，从频域的角度就能看的很清晰，下图是nips24的一篇文章，从图中我们看到了一条非平稳时间序列，但是它的统计特征，比如均值和方差却都是不变的，但是把这条序列做傅立叶变换后，我们发现它们的频域分量是不同的。

时序数据滤波

由于测量误差或短暂干扰，时间序列中可能存在高频噪声，可以在频域中去除这些高频成分，通过设置频率阈值来实现滤波，然后通过反傅里叶变换将数据恢复到时域，得到滤波后的时间序列，下图就是时序研究中常用的电力数据集，从红线中去除高频杂波，然后还原到时域得到“纯净“序列的过程。

这里我特别想说明一点，尽管上面的滤波例子去掉了高频分量，但这并不意味着高频分量不重要！实际上高频分量描述了时序数据的细节，现在就有一些研究反其道而行之，选择过滤低频分量保留高频分量。总而言之，不同的频率成分可能代表了时间序列中的不同特征，例如在股票市场时间序列分析中，低频成分可能与市场的长期趋势有关，而高频成分可能反映了短期的波动和噪声，我们可以有选择的过滤高频或低频分量。

数据压缩和降维

在频域中，通常只有少数频率成分包含了大部分的信息。可以通过保留这些主要频率成分，忽略次要频率成分，实现数据的压缩。对于长周期的时间序列数据，傅里叶变换后的频域表示可以用于数据压缩。

如何基于python做傅立叶变换

数据导入

这一部分我们使用时间序列建模分析领域的经典电力数据集，从数据导入、可视化分析、傅立叶变换、逆变换等几个方面，介绍如何基于 python做傅立叶变换。首先导入所需的数据集，由于数据集本身很大，我们截取其中的500个点，同时数据包含多列特征，我们这里以第‘0’列为例进行后续处理，具体代码如下：

import pandas as pd
import torch
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, ifft

df_electricity = pd.read_csv('electricity.csv')
df_electricity = df_electricity[0:500]

df0 = df_electricity['0']
electricity_time = pd.to_datetime(df_electricity['date'])

plt.figure(figsize=(16,8))
plt.plot(electricity_time, df0)

结果如下图所示，我们发现数据还是有一定周期规律的，但是看起来又没那么规整，我们的目的就是建模这种规律，剔除不规整的地方。

傅立叶变换

下面的代码就是如何进行傅立叶变换，有一点需要注意，scipy库在实现离散傅里叶变换时，没有在内部进行除以N这一步操作。这意味着scipy返回的结果与标准公式有一个归一化的差异。所以，当我们使用scipy得到 DFT 结果后，如果要和理论上的傅里叶变换结果保持一致，就需要手动除以N，乘以 2 是因为由于复数的引入，同一个振幅被分配至两个共轭复数上。

df_fft = fft(df0.to_list())
df_fft_norm = df_fft/len(df_fft)*2

fig, ax = plt.subplots(1, 2, figsize=(10, 4), dpi=150)

ax[0].stem([i for i in range(len(df_fft))], np.abs(df_fft), 'b', markerfmt=' ', basefmt='-b')
ax[1].stem([i for i in range(len(df_fft_norm))], np.abs(df_fft_norm), 'b', markerfmt=' ', basefmt='-b')

ax[1].set_xlim(-1, 50)

下图就分别是原始图和经过归一化后的傅立叶变换结果图：

过滤高频分量

freq_clean = pd.Series(df_fft).apply(lambda x: x if np.abs(x) > 1000 else 0).to_numpy()

fig, ax = plt.subplots(figsize=(10, 4))
ax.stem([i for i in range(len(df_fft))], np.abs(freq_clean), 'b', markerfmt=' ', basefmt='-b')
ax.set_xlim(-1, 100)
plt.show()

这里我们简单的选择过滤超过1000的高频分量，并展示过滤后的傅立叶变换结果图。

傅立叶逆变换

# 逆傅里叶变换
ix = ifft(freq_clean)

说起来很复杂，调用函数包就是一句命令的事情～～这时候，逆变换回去的结果实际已经把高频分量，也就是噪声过滤掉了，我们把滤波前后的结果同时展示出来。

# 绘制信号
fig, ax = plt.subplots(figsize=(12, 3))

ax.plot([i for i in range(len(df_fft))], ix, linewidth=.5, color='c', label='IFFT')
ax.plot([i for i in range(len(df_fft))], df0, linewidth=.5, color='r', linestyle='-.', label='Raw signal', alpha=0.7)

ax.set_xlabel('Sampling Time')
ax.set_ylabel('Amplitude')
ax.legend()

plt.show()

怎么样，是不是感觉傅立叶变换虽然原理挺难理解，但是从代码的角度，实现起来其实一点都不难！python封装的实在太好了，只要知道输入和输出，用起来就是几行代码。那么上文我们就算是把傅立叶变换的基本用法学会了。下面来看几篇经典顶会论文是如何使用傅立叶变换的。

顶会论文如何使用傅立叶变换

FITS: Modeling Time Series with 10k parameters (ICLR 2024 Spotlight)

def forward(self, x):
        # RIN，做归一化
        x_mean = torch.mean(x, dim=1, keepdim=True)
        x = x - x_mean
        x_var=torch.var(x, dim=1, keepdim=True)+ 1e-5
        x = x / torch.sqrt(x_var)

        low_specx = torch.fft.rfft(x, dim=1)
        low_specx[:,self.dominance_freq:]=0 # LPF

        low_specx = low_specx[:,0:self.dominance_freq,:] # LPF
        if self.individual:
            low_specxy_ = torch.zeros([low_specx.size(0),int(self.dominance_freq*self.length_ratio),low_specx.size(2)],dtype=low_specx.dtype).to(low_specx.device)
            for i in range(self.channels):
                low_specxy_[:,:,i]=self.freq_upsampler[i](low_specx[:,:,i].permute(0,1)).permute(0,1)
        else:
            low_specxy_ = self.freq_upsampler(low_specx.permute(0,2,1)).permute(0,2,1)
        low_specxy = torch.zeros([low_specxy_.size(0),int((self.seq_len+self.pred_len)/2+1),low_specxy_.size(2)],dtype=low_specxy_.dtype).to(low_specxy_.device)
        low_specxy[:,0:low_specxy_.size(1),:]=low_specxy_ # zero padding
        low_xy=torch.fft.irfft(low_specxy, dim=1)

        low_xy=low_xy * self.length_ratio # energy compemsation for the length change
        xy=(low_xy) * torch.sqrt(x_var) +x_mean
        return xy, low_xy* torch.sqrt(x_var)

这篇文章发表于ICLR24，是之前Dlinear模型的延续。论文核心思路是结合傅立叶变换与低通滤波，在复频域进行插值操作。模型先对时间序列数据预处理，经频域转换、低通滤波、线性变换与上采样后再逆变换回时域得到预测结果。它参数少（约1万）、效率高，适合资源受限的边缘设备，且在参数大幅减少情况下性能优越。

看上面的代码，x是输入数据，其维度是[batch_size, seq_len, fea_size]，经过归一化后，通过torch.fft进行傅立叶变换，变换过程和我们上面给出的python代码并无本质区别。然后，把超过dominance_freq的频域分量置为零，这里起到了过滤作用。再然后经过线性层变换，既上面的upsampler()，上采样和0填充，最后转回时域，这篇文章傅立叶变换起到滤波作用。

Frequency Adaptive Normalization For Non-stationary Time Series Forecasting（NIPS2024）

def main_freq_part(x, k, rfft=True):
    # freq normalization
    # start = time.time()
    if rfft:
        xf = torch.fft.rfft(x, dim=1)
    else:
        xf = torch.fft.fft(x, dim=1)

    k_values = torch.topk(xf.abs(), k, dim = 1)  
    indices = k_values.indices

    mask = torch.zeros_like(xf)
    mask.scatter_(1, indices, 1)
    xf_filtered = xf * mask

    if rfft:
        x_filtered = torch.fft.irfft(xf_filtered, dim=1).real.float()
    else:
        x_filtered = torch.fft.ifft(xf_filtered, dim=1).real.float()


    norm_input = x - x_filtered
    # print(f"decompose take:{ time.time() - start} s")
    return norm_input, x_filtered

这篇文章发表于NIPS24，基本思路是：考虑到时间序列傅立叶变换分解中，低频分量一般对应趋势项，高频分量一般对应季节项。通过傅立叶变换之后，可以分别对趋势项和季节项进行建模，对应代码中的norm_input和x_filtered，然后合并建模结果。

看代码，main_freq_part函数的主要目的是对输入的张量x进行基于频率的分解和滤波操作。它选择x在频域中的k个最大幅值对应的频率成分，将其余频率成分置零，然后通过逆傅里叶变换得到滤波后的信号x_filtered，并计算原始信号x与滤波后信号的差值norm_input。

mask = torch.zeros_like(xf)
mask.scatter_(1, indices, 1)
xf_filtered = xf * mask

创建一个与xf形状相同的全零张量mask，然后使用scatter_函数根据indices将对应位置的值设为1，从而创建一个掩码。这个掩码用于选择k个最大幅值对应的频率成分，将xf与mask相乘实现滤波，得到xf_filtered。

x_filtered = torch.fft.irfft(xf_filtered, dim=1).real.float()
norm_input = x - x_filtered

对滤波后的频域信号xf_filtered使用irfft（逆实值快速傅里叶变换）函数在维度1上进行逆变换，得到时域的滤波后信号x_filtered。最后，计算原始信号x与滤波后信号x_filtered的差值norm_input。返回norm_input和x_filtered进行下一步建模。