一文学透Crane DSP预测算法
云
孟凡杰，余宇飞约 12 分钟
孟凡杰，余宇飞约 12 分钟孟凡杰，腾讯云容器技术专家，FinOps产品研发负责人。
余宇飞，腾讯云专家工程师，专注云原生可观测性、成本优化等领域，Crane 核心开发者，现负责 Crane 资源预测、推荐落地、运营平台建设等相关工作。
（一）真实世界的时间序列
监控系统 采集CPU利用率等系统指标时，通常按固定采样频率，如每30秒或60秒采集一次。 这样采集形成的时间序列是一个个离散的点，在Prometheus中存储格式为Timestamp: Value键值对。 Crane代码（ 代码仓库地址： https://github.com/gocrane/crane ）中保存了我们在真实生产环境拉取的典型业务监控指标数据，以pkg/prediction/dsp/test_data/input0.csv 为例，其数据片段如下： Times tamp,Value
1628640000,39.722330 1628640060,39.821460 1628640120,40.021330 ...你可以运行signal_test.go中的TestSignal_Plot()将input0.csv画出如下的图，需要注意的是，图像看起来是一条连续的曲线，实际只是我们的绘图工具将上述离散的点连起来造成的假象。
大部分规律性的在线业务，时间序列都如下图一样有规律的潮汐变动。 那么，如何对这种复杂的时间序列变动进行预测呢？ 200年前的数学家傅里叶已经给出了答案。
（二）时序分解
下面的动图有助于直观理解，当指标随时间规律波动时，多个不同频率和振幅的信号如何组合起来生成不同的时域图。可以看到，不同频率、相位、振幅的信号组合，最终能生成不同形状的时序曲线。
下图更直观的展现了从时域和频域两个不同角度观测方波的结果，从时域角度看过去我们能看到上述代码绘制出的方波图像；从频域角度看过去，我们能看到在不同频率下面的振幅。无论如何复杂的曲线，在变换到频域以后，都是频率（波动的快慢）、振幅（波峰的高低）、相位（起始的偏移量）的组合。观察频域，是否一种世界突然静止了的感觉？
import numpy as npimport matplotlib.pylab as pltx = np.linspace(-np.pi, np.pi, 100)# 振幅amp=4/np.pi# 基波频率为 1/2πbase = amp*np.sin(x)# 协波分量频率为3, 5, 7harmonic1 = (amp/3)*np.sin(3*x)harmonic2 = (amp/5)*np.sin(5*x)harmonic3 = (amp/7)*np.sin(7*x)merged=base+harmonic1+harmonic2+harmonic3plt.plot(x,base)plt.plot(x,harmonic1)plt.plot(x,harmonic2)plt.plot(x,harmonic3)plt.plot(x,merged)plt.show()
运行上述python代码，可以得到与动图完全一致的图像，我们用简单的四个波的叠加画出了下图中紫色曲线。叠加的协波越多，最终的曲线越接近一个理想的方波。通过不同频率、振幅、相位的波的组合，我们可以拟合任何周期曲线。
前文提到，在现实世界中，我们获得的信号基本都是经过采样后的离散信号，所以处理采样信号时应用更广泛算法的是离散傅里叶变换 （Discrete Fourier Transform，DFT），它是一种将离散时间信号从时域转换到频域的方法。
（三）离散傅里叶变换DFT
以最简单的二维空间为例，任意一个向量都可以被表示为一对数，如下图的（1,2) 代表该向量向x轴的投影为1个单位长度，向y轴的投影为2个单位长度。这就是一种最简单的的变换。
傅里叶变换与二维空间中的变换类似，本质上是将时序空间中的信号投影到不同频率空间上去。那么要如何确定投影的频率，以及如何计算每个频率的振幅和相位呢？
一个复杂的周期信号可能会含有许多不同频率的分量，图4的方波就是四个不同频率的正弦波的叠加。其中，周期最长（也就是频率最小）的那个正弦波（蓝色曲线）被称之为基波（它的周期和叠加信号的周期相等）。而那些频率是基波频率整数倍的分量则被称为谐波，频率是几倍就被称为几次谐波，因此图4中有3次，5次和7次三个谐波分量。在频域坐标系中，基波频率和协波频率就构成了水平坐标轴上的频率单位。离散傅里叶变换就是将时域信号投射到这些频率上去。
在时域中，x轴代表时间，周期信号沿x轴以基波周期不断循环往复。生活中有一个非常常见的工具来记录这种时间的循环，对的，就是钟表。
可以把基波想象成在一个周期内，围绕这个表盘走一圈。而3次谐波、5次谐波和7次谐波的频率是基波的3倍、5倍和7倍，因此在一个周期时间内，会分别绕表盘走3圈、5圈和7圈。
在数学的世界里也有一个类似的表盘，那就是复平面。这里快速回顾下复数的概念。复数，作为实数的延伸，它使任一多项式方程都有解。复数中的虚数单位i，定义为-1的平方根。任一复数都可表示为a + bi，其中a及b皆为实数，分别称为复数的实部和虚部，对应复平面的实轴Re和虚轴Im上的投影。
复平面中半径为1的圆形叫做单位圆，单位圆上的任意一点到圆心的向量可以用自然常数e的指数表示，并且可以转换成为正弦和余弦函数的表示法。φ是该向量的夹角，e的iφ次方代表向量逆时针旋转，e的−iφ次方代 代表向量顺时针旋转。
复平面中半径为1的圆形叫做单位圆，单位圆上的任意一点到圆心的向量可以用自然常数e的指数表示，并且可以转换成为正弦和余弦函数的表示法。φ是该向量的夹角，e的iφ次方代表向量逆时针旋转，e的−iφ次方代 代表向量顺时针旋转。
既然傅里叶变换能将任何时域信号转换成为正弦波的组合，正弦函数又可以转换为自然常数的对数形式（参见图7中的欧拉公式），那么自然的傅里叶变换公式可以以自然对数的形式表示。是时候祭出离散傅里叶变换的定义了：
我们先找出最熟悉的部分，自然常数底数的部分。我们令
，在复平面中代表的意义是什么呢？其实就是将单位圆N等分，那么每一个等分的角度就是ω。
假设某个时序数据在一个周期内有8个采样点，也就是基频信号绕复平面转一圈的过程中会有8个采样点，每次采样时间间隔完全一样。那么在复平面上，我们是不是就是把单位圆切分成了8等分，每一等分的夹角为 2 π / 8 那么要提取信号频率特征，我们是否就可以通过振幅和复指数的乘积就可以计算出每个数据数据在ω0-ω7的投影。
基频的投影逻辑也可以一样套用到协波频率上去。以2次谐波为例，2次谐波的频率是基波的两倍，也就是在相同时间内，二次谐波绕着单位圆走了2圈。每一个采样间隔，走了2倍的夹角，因此，2次谐波落在了ω2，ω4，... ω14这8个点上。我们同样可以通过振幅和复指数的乘积就计算出每个数据数据在ω2-ω14的投影。
最终，我们将这投影值按照采样顺序求和，就得到了每个频率的特征信息，这就是对离散傅里叶变换的朴素表述。我们可以用 F=W*x 矩阵计算来表示傅里叶变换。x是原始输入信号，也就是监控数据，W是如下的矩阵：
从上面的矩阵乘法看出，DFT计算的时间复杂度为O(n2)，在数据样本较大的时候几乎无法得出实时结果，因此应用并不广，直到出现了快速傅里叶变换。
（四）快速傅里叶变换FFT
快速傅里叶变换是通过巧妙的数学技巧加速傅里叶变换计算的方法。 我们先取W矩阵中的第二行，与x（假设x=[1,5,7,1,2,3,2,3])相乘，则可做如下推导：
你可能
注意到
了，通过这样的分解，我们将一个高阶的多项式可以拆分成两个低阶多项式，多项式中有很多重复的部分，而且。
而被拆解出来的两个多项式，我们可以继续分解。
从程序实现的角度思考，是不是有了递归分解以及动态规划减少重复运算的感觉？
对的，这就是FFT的本质。
请注意是复平面上面的单位圆上被N等分的点，这些点有如下一些特性：
| 公式 | 解释 |
|---|---|
|
|
一个点的平方等于将该点绕复平面原点旋转两倍夹角 |
|
|
对称性，一个点绕复平面原点转半圈得到的点与原始点相反 |
|
|
共轭，即实部相等，虚部相反 |
|
|
一个点绕一圈以后与原点重合 |
这些特性使得我们采用分治的方法快速计算傅里叶变换，因为基于递归降低了复杂度，基于复数的特性，使得无论计算的多少次方，事实上都是在单位圆上的被N等分的点上反复计算和取值，利用这些特性，就能实现减少大量重复计算的目的，而递归的算法也是时间复杂度从O(n2)变成了O(n*log2(n))。
（一）位反转
记得我们如何通过多次多项式分解将高阶多项式分解成奇偶两个部分吗？该步骤的目的是将多项式不断分解为Fe、Fo，直到每组中只剩下一对样本，此步骤的本质上是交换样本位置。
方法是将奇数项和偶数项归类，然后再对分解后的低阶多项式做重复动作。针对8个元素组成的多项式，经过3次分解，原有的多项式中的0-7的顺序就被交换成为了0,4,2,6,1,5,3,7。我们观察二进制表示，可以看出，位反转的本质是元素序号二进制表示的位反转，比如100反转后为110，011反转后为011。
| 原始位（二进制） | 第一次反转 | 第二次反转 |
|---|---|---|
| 0 (000) | 0 (000) | 0 (000) |
| 1 (001) | 2 (010) | 4 (100) |
| 2 (010) | 4 (100) | 2 (010) |
| 3 (011) | 6 (110) | 6 (110) |
| 4 (100) | 1 (001) | 1 (001) |
| 5 (101) | 3 (011) | 5 (101) |
| 6 (110) | 5 (101) | 3 (011) |
| 7 (111) | 7 (111) | 7 (111) |
（二）蝶形变换
对每组样本执行DFT计算，这个计算过程可被简化为加减法和乘法的基本运算组合；使用刚刚执行的计算结果，组合成样本对，作为下一阶段的输入，直到得出最终答案。
下图是8个样本的FFT的计算过程。
（三）Crane对FFT的应用
Crane调用了go-dsp完成FFT的计算，仔细分析会发现具体过程与上述流程完全一致。下图展示了该过程的详细过程，并附上了8点数据的debug过程。
-
将输入采样数据转为复数
-
判断是否采样数据为2的指数，如果不是则在采样数组后面补0
-
调用getRadix2Factors，依据样本长度计算
- 调用reorderData进行位反转，重新排列样本
-
起动多个worker线程并发执行蝶形变换的计算逻辑，并汇总结果
回到开篇的场景，针对任何运行在Kubernetes集群上的作业我们没有计算函数，只有检测系统采集上来的时域数据，如何从这些观测数据中提取特征呢。
input0.csv作为在线应用的系统监控指标，数据有如下主要特征：
-
长度有限：无论是一天、一个月还是一年的数据，都是有限长度。
-
数据离散：指标采集通常是按固定时间周期采集，比如Prometheus scrape_interval默认配置为一分钟，也就是每分钟采集一次数据，每天1440个离散的数据点。
-
周期波动：如图一所示，业务量呈现周期性，并且在一个周期内的不同时段有明显的峰谷特征。
-
数据预处理：
- 主周期判断
-
分为两个阶段：
（2）计算序列的自相关函数（Auto Correlation Function, ACF）。ACF是一个信号于其自身在不同时间点的互相关。通俗的讲，它就是两次观察之间的相似度对它们之间的时间差的函数。相差k个采样间隔的离散时间序列x的自相关函数定义如下：
可以想象，如果M是一个序列的周期，它的自相关函数在M点的取值一定是一个局部的高点。根据这个特性，我们对第一阶段得到的候选周期在ACF图上进行确认，最终选出位于「最高峰」的点作为序列的主周期（也就是基波周期）。
（
3）预测
调用FFT函数将时域指标数据转换为频域数据，过滤掉噪音，并调用逆快速傅里叶变换（IFFT），将频域信号转换成时域信号，作为下一个周期的预测数据。
将预测后的指标数据的下标，按采样频率转换成timestamp，与对应的预测值组合成TimeSeries，并写回TimeSeriesPrediction.Status
DSP 参数说明与配置
下面的代码片段展示了DSP算法的核心参数，参数说明如下，在使用时可依据您的业务特性酌情配置：
dsp:sampleInterval: "60s" //采样时间间隔，应与监控数据保持一致，Prometheus默认为60historyLength: "14d" // predictor 从监控系统拉取指标的历史长度，4d代表拉取过去十四天的指标，historyLength 至少需要预测周期的两倍及以上。estimators:fft:- marginFraction: "0.2" // 预留资源余量，0.2代表给出的预测值是同时间段历史指标的1.2倍// 高于最高频率阈值、切振幅低于此阈值的信号会被忽略掉，lowAmplitudeThreshold: "1.0" // 最低振幅阈值highFrequencyThreshold: "0.05" //最高频率阈值// 最大最小频率数，时域信号转为频域信号时，频率数会介于二者之间，如下面的配置代表频率数不低于10，不高于20minNumOfSpectrumItems: 10maxNumOfSpectrumItems: 20
Crane支持多种预测算法，且算法可灵活扩展。本文从DSP算法实现层面展开，希望能从数学原理层面将预测原理讲解清楚。
参考文献
1.傅里叶级数- 维基百科，自由的百科全书 2.离散傅里叶变换- 维基百科，自由的百科全书 3.复平面- 维基百科，自由的百科全书 4.复数(数学) - 维基百科，自由的百科全书 5.The Fast Fourier Transform (FFT): Most Ingenious Algorithm Ever? 6.DIT FFT algorithm | Butterfly diagram | Digital signal processing点击下方空白 ▼ 查看明日开发者黄历