时间序列分析方法：离散傅里叶逼近

前言

在时间序列的分析当中，人眼直接观察到的往往是数据的趋势、形状，这些属于数据中变化缓慢的成分，即数据的低频成分，离散傅里叶近似的思想即利用傅里叶系数前一部分系数逼近原始数据的大体形状和趋势。

---

一、理论部分

1.傅里叶变换

讲傅里叶变换的文章有很多，我就不在这里献丑啦，把我自己学习过程中的一些内容总结总结~
傅里叶变换公式：
$$X(f)=\int x(t)e^{-j2\pi ft}dt$$
离散傅里叶变换公式：
$$X(k)=\sum _{n=0}^{N-1}x(n)e^{-j\frac{2\pi }{N}nk}$$
其中$N$为原始序列长度
对比一下连续和离散公式的关系：
$$2\pi ft\Rightarrow 2\pi nk/N_{fft}\Rightarrow \{\begin{array}{r}t\to n/fs\\ f\to k\cdot df\end{array}$$
其中$fs$为信号的采样频率，$f$为对应的频率，$N_{fft}$为傅里叶变换点数，$df$为频谱分辨率，$df=fs/N_{fft}$
之前在学习的过程中对于傅里叶变换点数、数据点数这些弄得很迷糊，做毕设的时候又好好捯饬了一下，把它们写成矩阵形式可能更好理解：
$$\left[\begin{array}{c}X(0)\\ X(1)\\ ⋮\\ X(N_{fft}-1)\end{array}\right]=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ 1& W^{1\cdot 1}& \cdots & W^{1\cdot (N-1)}\\ ⋮& ⋮& \ddots & ⋮\\ 1& W^{(N_{fft}-1)\cdot 1}& \cdots & W^{(N_{fft}-1)(N-1)}\end{array}\right]\left[\begin{array}{c}x(0)\\ x(1)\\ ⋮\\ x(N-1)\end{array}\right]$$
其中$W=e^{-j2\pi /N}$
这样看起来就更直观一些了，傅里叶变换点数即矩阵的行数，对应k的选择范围，在matlab或者其他编程语言的函数中，不刻意设定的话基本都是$N_{fft}=N$。
根据这个矩阵表达形式，也可以调整一些参数来实现限制傅里叶变换得到的频率范围。
有点跑题了，总而言之，傅里叶变换利用三角基函数的组合对原始函数、数据进行了表示。

2.离散傅里叶逼近

对原始数据进行离散傅里叶变换，该算法选择离散傅里叶系数中的前$\omega$个，之后还原至时域表示，可以得到的原始数据大体形状趋势。但是需要注意，对于文章提出的离散傅里叶逼近方法还有一些细节：

• 索引为0的频率成分可以选择是否删除，即对应时域的均值；
• 离散傅里叶变换得到的数据是复数序列，并且是共轭对称的，当我们选择了前$\omega$系数后直接傅里叶反变换是肯定不能得到我们期望的结果的，需要将这些系数补零至原始长度，再将这些系数取共轭后反向排列放置在末端，即：
  $$[c_1,c_2,...,c_{\omega }]\Rightarrow [c_1,c_2,...,c_{\omega },0,...,0]\Rightarrow [c_1,c_2,...,c_{\omega },0,...,0,{\bar{c}}_{\omega },...,{\bar{c}}_2,{\bar{c}}_1]$$
  这样再进行傅里叶反变换得到的结果就使我们希望的结果了，但是也要注意因为精度等问题会使得结果中仍然存在很小很小的虚数部分，可以取实部来去除掉。

二、实战

1.Matlab自编代码

function [dftcoeff,out] = DFT(x,coeff)
    n = length(x);
    realindex = 2:2:coeff;
    imagindex = 3:2:coeff;
    %提取正半轴频率
    tempxf = fft(x);
    %将实部虚部重新组合
    if mod(n/2,2) == 0
        realPart = real(tempxf(1:round(n/2)));
        imagPart = imag(tempxf(1:round(n/2)));
        newsub = reshape([realPart;imagPart],[1,n]);  
    else
        realPart = real(tempxf(1:round(n/2)+1));
        imagPart = imag(tempxf(1:round(n/2)+1));
        newsub = reshape([realPart;imagPart],1,length(realPart)+length(imagPart));
    end
    newsub = [newsub(1),newsub(3:end-1)];
    newsub = newsub(1:coeff);
    if coeff/2 == round(coeff/2)
        newsub = [newsub(1),newsub(realindex)+1i*[newsub(imagindex),zeros(1,1)]];
    else
        newsub = [newsub(1),newsub(realindex)+1i*[newsub(imagindex)]];
    end
    newsub(1) = 0;%根据需求去均值化
    dftcoeff = real(newsub);
    newsub = [newsub,zeros(1,length(x)-2*length(newsub)+2),conj(newsub(end-1:-1:2))];%取共轭对称
    out = real(ifft(newsub));
end

2.测试结果

如图为选择不同数量傅里叶系数还原的时域图像以及误差（误差采用欧氏距离），数据为UCR数据库中GunPoint数据集的其中一个。

---

总结

离散傅里叶逼近算法比较简单，文章主要就是捋了一下自己学习过程中的一些思考和总结，如果有什么错误烦请指出，也希望能和大家多多学习多多交流。

参考文献

Schäfer P, Högqvist M. SFA: a symbolic fourier approximation and index for similarity search in high dimensional datasets[C]//Proceedings of the 15th international conference on extending database technology. 2012: 516-527.