1. 说明

 用傅立叶变换预测时序数据,原理是把时域数据转换到频域,再转换回来.python的numpy和scipy里面都有现成的转换工具fft()和ifft(),但使用时会遇到一个问题:比如25天的数据转到频域再转回时域,还是25天,虽然拟合了数据,但没法直接预测未来,本篇介绍用它实现预测的方法.

2. 傅立叶变换

(1) 相关知识

 之前写过关于傅立叶变换原理的文档,这次就不再重复了,具体请见:https://www.jianshu.com/p/9e786be6dccb 本篇只从程序的角度看如何使用它. 经过FFT转换的数据和转换前长度一致,每个数据分为实部和虚部两部分,假设时序时数长度为N(N最好是2的整数次幂,这样算起来更快),用fft()转换后:下标为0和 N /2的两个复数的虚数部分为0,下标为i和 N - i 的两个复数共辄,也就是其虚部数值相同、符号相反。再用ifft()从频域转回时域之后,出现了由误差引起的很小的虚部,用np.real()取其实部即可.
 由于一半是另一半的共轭,因此只需要关心一半数据.fft转换后下标为0的实数表示时域信号中的直流成分(不随时间变化),下标为i的复数 a + bj,其中a表示余弦成分,b表示其正弦成分.

(2) 示例功能

 数据是航空乘客数据"AirPassengers.csv",可以从CSDN下载,其中包括从1949-1960年,每月旅客的数量,程序预测未来几年的旅客数据.

 如图所示,数据为非平稳数据,其趋势向上,且波动加俱,为将其变为平稳数据, 先对其做了对数和差分处理.

(3) 示例代码

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# 函数功能:将频域数据转换成时序数据
# bins为频域数据,n设置使用前多少个频域数据,loop设置生成数据的长度

def fft_combine(bins, n, loops=1):
 length = int(len(bins) * loops)
 data = np.zeros(length)
 index = loops * np.arange(0, length, 1.0) / length * (2 * np.pi)
 for k, p in enumerate(bins[:n]):
   if k != 0 : p *= 2 # 除去直流成分之外, 其余的系数都 * 2
   data += np.real(p) * np.cos(k*index) # 余弦成分的系数为实数部分
   data -= np.imag(p) * np.sin(k*index) # 正弦成分的系数为负的虚数部分
 return index, data

if __name__ == '__main__':
 data = pd.read_csv('AirPassengers.csv')
 ts = data['Passengers']

# 平稳化
 ts_log = np.log(ts)
 ts_diff = ts_log.diff(1)
 ts_diff = ts_diff.dropna()
 print(fy[:10]) # 显示前10个频域数据
 fy = np.fft.fft(ts_diff)
 conv1 = np.real(np.fft.ifft(fy)) # 逆变换
 index, conv2 = fft_combine(fy / len(ts_diff), int(len(fy)/2-1), 1.3) # 只关心一半数据
 plt.plot(ts_diff)
 plt.plot(conv1 - 0.5) # 为看清楚,将显示区域下拉0.5
 plt.plot(conv2 - 1)
plt.show()

(4) 运行结果

[ 1.34992672+0. j -0.09526905-0.14569535j -0.03664114-0.12007802j
 -0.2670005 +0.24512406j -0.10075074+0.0314084 j -0.26409417+0.04197159j
 0.14411338+0.18703009j 0.07467991+0.05367644j -0.26663142+0.15324939j
 0.03248223+0.14130114j]

(5). 示例分析

 输出的是fft转换后的数据,只显示了前十个,形式为复数.复数模(绝对值)的两倍为对应频率的余弦波的振幅;复数的辐角表示对应频率的余弦波的相位。第0个元素表示直流分量,虚部为0.在数据中的位置标记了频率大小,值标记了振幅大小.
 图中显示的三条曲线分别为原始数据,做了fft以及ifft逆变换后的数据,以及fft后自己实现算法还原并预测了未来的数据,从图中可见,基本拟合了原始曲线,预测曲线看起来也比较合理.
 上述方法可实现用傅里叶变换预测时序数据.与ARMA算法相比,它没有明显衰减,更适合长时间的预测.
 对于随时间变化的波形,比如语音数据,一般使用加窗后做傅立叶变量的方法拟合数据.

3. 小波变换

(1) 相关知识

 有了傅立叶变换,为什么还用小波呢?上面提到,如果波型随时间变化,就需要对波型加窗分段后再处理,而且有时候需要大窗口,有时候需要小窗口,处理起来就更加麻烦.于是引入了更灵活的小波.
 傅立叶变换的基是正余弦函数,而小波的基是各种形状小波,也就是说它把整个波形看成是多个位置和宽度不同的小波的叠加.小波有两个变量:尺度a和平移量t,尺度控制小伸的伸缩,平移量控制小波的平移,它不需将数据切分成段,就可以处理时变数据.尤其对突变信号处理得更好.
 下图是几种常见的小波.

 离散小波变换,Discrete Wavelet Transformatio (dwt),可以说是小波变换中最简单的一种。这里使用Python调用pywt库实现最简单的功能.
 经过变换之后的返回值:cA:Approximation(近似), cD:Detail(细节),其中近似cA是周期性有规律的部分,可以被模拟和预测,而cD可看做是噪声。换言之,用此方法可以拆分周期性数据,和其上的扰动数据。

(2) 示例功能

 示例使用的仍然是乘客数据,下面代码是将细节D设为0,然后还原。

(3) 示例代码

import pywt
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

if __name__ == '__main__':
 data = pd.read_csv('AirPassengers.csv')
 ts = data['Passengers']
 # 平稳化
 ts_log = np.log(ts)
 ts_diff = ts_log.diff(1)
 ts_diff = ts_diff.dropna()

 cA,cD = pywt.dwt(ts_diff, 'db2')
 cD = np.zeros(len(cD))
 new_data = pywt.idwt(cA, cD, 'db2')

 plt.plot(ts_diff)
 plt.plot(new_data - 0.5)

 plt.show()

(4) 运行结果

(5) 示例分析

 可以看到,用小波拟合的效果也还可以,一般可以使用小波拟合cA,使用ARMA拟合cD部分,两种方法配合使用.

posted on 2018-03-28 19:35  xieyan0811  阅读(114)  评论(0编辑  收藏  举报