============================‼️注意‼️============================
由于 jupyter nbconvert 的后端限制，无法渲染使用有控件的图片，所以文章存在很多"空白“的图像。读者不妨自己写一下。
============================‼️注意‼️============================

滑动滤波

1
2
3
4
5
6
7
8
9
10
11
%matplotlib widget

import matplotlib.pyplot as plt
import numpy as np
import myfilter.MovingAverageFilter as MAF
import myfilter.FilterPerformance as FP

import ipywidgets as widgets
from ipywidgets import FloatSlider

figsize=(24, 8)

在第一节中，我们讨论了对于静态目标的测量噪声的抑制方法。我们认为测量过程中的所有信息都是有效且平权的，所以采取的最直接的方法————累计平均（Cumulative Average）。但实际应用中，我们往往面对着动态过程的测量，比如测量一辆匀加速运动中小车的实时速度。这种情况下，使用累计平均来进行滤波显然是不可靠的。Why？

  1. 从数学模型上来说，CA滤波器在对系统建模时认为其是完全不变的 ————> 同一滤波器对不同模型的所表现出的性能是不同的。
  2. 从性能分析上来说，滤波器的性能可以分为静态性能(对噪声的抑制能力)与动态性能(对真值变化的跟踪能力) ————> 选择滤波器时应综合考虑。

如果在我自己浅薄的知识范围内对所有滤波器的工作原理进行分析的话应该是————运用合理的数学方法对信息进行提取与融合。
每次测量的结果都可以视为真值+噪声真值+噪声对于单一结果的话我们无法区分每种成分所占的比例(我们所拥有的信息太少了)。但是，每一次新的测量都会带来新的信息(新息序列)，利用合理的数学方法将其中的信息提取出来的过程便是滤波。从过去的数据中提取出信息，再与新息进行融合，从而得到一个更可靠的数值。

本教程大多讨论滤波器的时域性能而非频域性能，一般来说滤波器的时域性能与频域性能是互斥的。而我写这个教程的初衷是希望新生能够对传感器的数据(加速度、角速度、磁场和电压等等物理量)进行有效提取，而非对频域上的信号进行分离，故FIR与IIR滤波器不做重点讨论。

运动状态的建模

递推模型的建立

设某单位质量的质点做匀加速直线运动，加速度为 aa，速度为 vv，考虑常见的运动学方程 v=at+v0v=at+v_0 。
我们的目标是得到从 n 时刻到 n+1 时刻, 每个时刻间隔为 Δt\Delta t 的递推表达式，记 t 时刻的速度与加速度分别为 vt,atv_t, a_t 显然有

\begin{cases}
v_{n+1} = v_{n} + \Delta t \cdot a_{n} \
a_{n+1} = a_{n}
\end{cases}

记 xn=[vnan]\mathbf{x_n} = \begin{bmatrix} v_{n}\\a_{n}\end{bmatrix} (称为递推变量), 于是矩阵形式为

xn+1=[1Δt01]xn\mathbf{x_{n+1}} = \begin{bmatrix} 1 & Δt \\ 0 & 1 \end{bmatrix} \mathbf{x_{n}}

另记 A=[1Δt01]\mathcal{A} = \begin{bmatrix} 1 & Δt \\ 0 & 1 \end{bmatrix} 。

噪声的引入

在现实情况中，一个运动物体往往会因周围环境的改变，而受到各种干扰力(地面的摩擦力、空气阻力或自身输出动力的不平稳等等)。由上节内容可知，我们可以认为这些力的合力满足正态分布，即

f=ma∼N(μ,σ2)f=ma\sim N(\mu,\sigma ^2)

设 w∼N(μ,σ2)w \sim N(\mu,\sigma ^2), 于是递推方程应该被修改为

xn+1=Axn+[0w]\mathbf{x_{n+1}} = \mathcal{A} \mathbf{x_{n}} + \begin{bmatrix} 0\\w \end{bmatrix}

下述代码展示了仅使用模型进行的纯推断过程。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
### 模型参数 ###
# 采样参数
v0 = 2.
a0 = 5.
dt = 0.01
TotalTime = 3.
xSequence = np.arange(0., TotalTime, dt)

A = np.array([[1., dt], 
			  [0., 1.]])
B = np.array([[0.],
			  [0.]])
xn = np.array([[v0],
			   [a0]])

### 仿真 ###
ModelSequence = np.empty((2, xSequence.size))
RealSequence  = np.empty((2, xSequence.size))
ModelSequence[:, 0] = xn.T
RealSequence[:, 0]  = xn.T

for i in range(1, xSequence.size):
	ModelSequence[:, i] = (A @ ModelSequence[:, i-1])				
	B[1, 0] = np.random.normal(0, 1)
	RealSequence[:, i] = (A @ RealSequence[:, i-1]) + B[:, 0]
	if i % int(xSequence.size/5) == 0:			
		RealSequence[0, i] += np.random.normal(0, 1)		# 模拟碰撞或其他某种原因导致的突变

### 绘图展示 ###
plt.figure(figsize=figsize)
plt.plot(xSequence, ModelSequence[0,:], '.-', label="Model", markersize=5)
plt.plot(xSequence, RealSequence[0,:], 	'*-', label="Actual", markersize=5)
plt.legend(loc='lower right') #绘制图例
<matplotlib.legend.Legend at 0x77cff0f66a20>
Figure

至此，模型估计值在本章就杀青了。将模型推导过程放在这里只是为了展示逻辑过程与数学过程(仿真要用到)，本章教程内容是以观测为核心的。
等到模型估计值回到我们视野的时候，可以用一句诗来形容：“此去泉台招旧部，旌旗十万斩阎罗”。

观测过程

不难发现，即使噪声标准差为 1, 纯理论上的估计也很可能会产生极大的误差。在长时间的工作下，(模型估计值-真实值)必然是发散的。
因此我们需要进行观测，即引入新息。
一般的，设在 n 时刻的被观测量为 zn\mathbf{z_{n}}, 与上述过程类似，建立被观测量与递推变量 xn\mathbf{x_{n}} 的联系

zn=H∗xn+v\mathbf{z_{n}} = \mathcal{H} * \mathbf{x_{n}} + \mathbf{v}

大多数观测过程都是近似线性的，所以一个线性方程就足够了。
举个栗子，我们想知道某串联电路流过的电流，于是测量电路上一个 2Ω 电阻的电压，那么 U = 2I 。

对于本例，我们直接认为被测量为速度，并且测量误差 vv 为标准正态噪声 N(0,1)N(0,1) ，即

H=[1000]v=[v0]\mathcal{H}=\begin{bmatrix}1&0\\0&0\end{bmatrix} \qquad \mathbf{v}=\begin{bmatrix}v\\0\end{bmatrix}\\

与上文类似，下述代码展示了仅使用传感器进行的纯观测过程。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
def linear_regression(xSequence, ySequence):    
    N = len(xSequence)
    sumx = sum(xSequence)
    sumy = sum(ySequence)
    sumx2 = sum(xSequence**2)
    sumxy = sum(xSequence*ySequence)

    A = np.array([[N, 	 sumx],
                  [sumx, sumx2]])
    b = np.array([sumy, sumxy])

    return np.linalg.solve(A, b)

MeasSequence = np.empty((2, xSequence.size))
MeasSequence[0, :] = RealSequence[0, :] + np.random.normal(0, 2, xSequence.size)
MeasSequence[1, :] = 0

plt.figure(figsize=figsize)

plt.plot(xSequence, MeasSequence[0,:], '*', label="Measure", markersize=6)
plt.plot(xSequence, RealSequence[0,:], '*', label="Actual", markersize=6)

b, k = linear_regression(np.arange(0, TotalTime, dt), MeasSequence[0, :])

y = k * xSequence + b
plt.plot(xSequence, y, 'ro-', markersize=3,label="Regression") 	# 线性回归拟合
plt.title("y = {}x + {} ".format(k.round(2), b.round(2)))

plt.legend(loc='lower right') #绘制图例
<matplotlib.legend.Legend at 0x77cff11e7860>
Figure

不同的滑动滤波策略

在动态系统中，较早时刻的测量数据很可能已无法反映当前系统的状态，若对所有数据一视同仁，易引入累积偏差，导致滤波结果滞后或失真；相比之下，近期数据通常包含更多与当前状态相关的新息，具有更高的信息价值。
因此，以平均滤波为基础，提升滤波器的动态响应性能的一种直接的思路是：引入遗忘机制，逐步舍弃历史数据中时效性较弱的部分，而更重视近期观测数据的作用。
从数学形式看，这种逐步遗忘历史数据的机制可表述为某种加权策略。（人话：通过赋予不同权重，使新数据的信息更重要)

实际上就是一个卷积过程嘛， t 时刻的响应等于全部时刻激励对该时刻的激发之和。而我们所说的某种加权策略对应着卷积核的构造。

简单滑动滤波(SMA, Simple Moving Average)

介绍

一种最简单的方法是只对最近 k 个时刻的简单平均。一个卷积核长度为 n 的 SMA 滤波器，其数学形式为

x^n=1k⋅∑i=0k−1xn−i\hat{x}_n = \frac{1}{k} \cdot \sum_{i=0}^{k-1}x_{n-i}

python 实现

SMA class 在 python 中的实现如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
class SimpleMovingAverageFilter:
    def __init__(self, window_size):
        self.window_size = window_size

    def update(self, new_value):
        self.window.append(new_value)
        if len(self.window) > self.window_size:
            self.window.pop(0)
        return sum(self.window) / len(self.window)

    def estimate(self, MeasSequence):
        EstimateSequence = []
        for i in range(len(MeasSequence)):
            EstimateSequence.append(self.update(MeasSequence[i]))
        return EstimateSequence

仅供参考，从没学过 python ，就是凭着感觉随手写写+面向谷哥编程。把面向对象语言当面向过程写，丑陋至极😓。

性能分析

仿个真先

1
2
3
4
5
6
7
8
9
moveFilter = MAF.SimpleMovingAverageFilter(5)
EstimateSequence = moveFilter.estimate(MeasSequence[0,:])

### 绘图展示 ###
plt.figure(figsize=figsize)
plt.plot(xSequence, RealSequence[0,:],	'*-', label="Actual", markersize=5)
plt.plot(xSequence, MeasSequence[0,:],	'+-',  label="Measure", markersize=5)
plt.plot(xSequence, EstimateSequence[:],	'o-', label="Estimate", markersize=5)
plt.legend(loc='lower right') #绘制图例
<matplotlib.legend.Legend at 0x77cff0fa2bd0>
Figure

很容易有一个感性认识：滤波效果是不错的，对噪声有明显的抑制。但问题是：有多好，如何给出一个客观比较？
我们需要把感性认识上升为理性认识，即通过某种数学方法来描述出滤波前后两组数据的误差程度。

首先，对于任何一个滤波器来说，在不考虑代价的时候，我们希望其一定是将噪声尽可能的去除，同时对于信号有一个很好的跟踪能力。
也就是说衡量一个滤波器可以分为两个部分，静态降噪能力(平滑)与动态跟踪性能(保真)。但是如果进行一些简单的分析，就会发现这两个目标是互斥的，如果增强平滑性能，那么跟踪能力就会下降。

傅里叶变换、或者上文的分析都可以简单得到这个结论。工程师要做的就是在诸多目标中找到一个平衡点，以实现原定需求。

在我们的需求下，滤波器应该是工作在稳态，所以我们舍弃前 30 个点。回顾整个建模过程，我们一共引入了两次噪声。一次是影响物体运动的过程噪声，一次是影响观测的测量噪声。这两种噪声同为正态分布，则总体噪声应该满足分布 N(Δt⋅μ1+μ2,Δt2σ12+σ22)N(\Delta t\cdot \mu_1+\mu_2, \Delta t^2 \sigma^2_1+\sigma^2_2) 。
显然，我们所希望得到的噪声应该服从 N(0,σ‘2)σ‘2<Δt2σ12+σ22N(0, \sigma`^2)\quad\sigma`^2<\Delta t^2 \sigma^2_1+\sigma^2_2 。所以我们对其残差进行分析，问题：

  1. 残差是否服从正态分布？ \qquad --> \qquad 采用 Q-Q 图与 P-P 图进行数据展示，并使用 Jarque-Bera 检验量化分析是否符合正态分布。
  2. 残差是否具有无偏性？ \qquad \quad --> \qquad 单样本 t 检验量化分析是否具有无偏性。
  3. 方差是否明显减小？ \qquad --> \qquad 使用信噪比量化降噪能力。(肯定减小啊，不然在干嘛?偷懒😁，就不进行 F 检验了)

一般认为信噪比至少要达到 10dB 才可用。

代码如下所示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
class FilterPerformance:
	def __init__(self, figsize=(10,5)):
		self.figsize = figsize
	
	def CalculateSNRdB(self, RealSequence, SignalSequence):
		signal  = np.array(RealSequence) 
		noise0  = np.array(SignalSequence  - RealSequence)
		snr_dB = 10 * np.log10(np.var(signal) / np.var(noise0))
		return snr_dB


	def CalculateGain(self, RealSequence, MeasSequence, EstimateSequence):
		snr_dB0 = self.CalculateSNRdB(RealSequence, MeasSequence)
		snr_dB1 = self.CalculateSNRdB(RealSequence, EstimateSequence)
		return snr_dB1 - snr_dB0

	def Analysis(self, RealSequence, MeasSequence, EstimateSequence):
		residual = EstimateSequence - RealSequence
		# 正态性检验
		fig, axs = plt.subplots(nrows=1, ncols=2, figsize=self.figsize)
		axs[0].title.set_text('Q-Q Plot')
		axs[1].title.set_text('P-P Plot')
		sm.qqplot(residual, line='45', ax=axs[0])
		sm.ProbPlot(residual, loc=0, scale=1, fit=False).ppplot(line='45', ax=axs[1])

		jb_statistic, p_value = stats.jarque_bera(residual)
		print(f"JB Statistic: {jb_statistic}， P-value: {p_value}")

		if p_value < 0.05:
			print("拒绝原假设，数据显著不服从正态分布")
		else:
			print("接受原假设，数据服从正态分布")

		# 单样本 t 检验
		t_statistic, p_value = stats.ttest_1samp(residual, 0)
		print(f"t-statistic: {t_statistic}, p-value: {p_value}")
		if p_value < 0.05:
			print("拒绝原假设，均值显著不等于0")
		else:
			print("接受原假设，均值等于0")

		# 降噪增益
		gain = self.CalculateGain(RealSequence, MeasSequence, EstimateSequence)
		fig.suptitle(f"gain: {gain:.2f} dB")
		
		plt.tight_layout()
		plt.show()

下面展示一下分析效果。

1
2
3
performance = FP.FilterPerformance(figsize=figsize)
performance.Analysis(xSequence, RealSequence[0,:], MeasSequence[0,:], EstimateSequence[:])
plt.show()
Figure 
JB Statistic: 54.5257580270603， P-value: 1.4450473690815274e-12
拒绝原假设，数据显著不服从正态分布
t-statistic: -1.3454773078886768, p-value: 0.17949059236742446
接受原假设，均值等于0

然后展示一下，滤波核大小变化对数据的影响。(拖动滑动条可以改变滤波核，并重新计算性能表现)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
performance = FP.FilterPerformance(figsize=figsize)
performance.PreAnalysis()
def PerformanceDemo(SMA_KernalSize):
	moveFilter = MAF.SimpleMovingAverageFilter(int(SMA_KernalSize))
	EstimateSequence = moveFilter.estimate(MeasSequence[0,:])
	performance.ReAnalysis(xSequence ,RealSequence[0,:], MeasSequence[0,:], EstimateSequence)
	plt.show()

p_max = 0
max = 0
for i in range(1,51,1):
	moveFilter = MAF.SimpleMovingAverageFilter(i)
	EstimateSequence = moveFilter.estimate(MeasSequence[0,:])
	t = performance.CalculateGain(RealSequence[0,30:], MeasSequence[0,30:], EstimateSequence[30:])
	if t < max:
		max = t
		p_max = i
widgets.interact(PerformanceDemo, 
                 SMA_KernalSize = FloatSlider(p_max, min=1, max=50 , step=1, description='SMA_KernalSize'))
print(f"最大降噪增益: {max} dB 出现在卷积核大小: {p_max}")
Figure 
interactive(children=(FloatSlider(value=10.0, description='SMA_KernalSize', max=50.0, min=1.0, step=1.0), Outp…


最大降噪增益: -8.761633964492876 dB 出现在卷积核大小: 10

加权滑动滤波(WMA, Weithed Moving Average)

介绍

在简单滑动滤波中，我们认为滤波核内的信息是同等重要的。而加权滑动滤波则给予滤波核内数据不同的权值。一般常见的权值构造为三角形。

1
2
3
4
5
6
7
8
9
10
11
window_size = 20
kernel = []
normal_Coefficient = window_size * (window_size + 1) / 2
for i in range(window_size):
    kernel.append(i / normal_Coefficient)

fig = plt.figure(figsize=figsize)
ax = fig.add_subplot(1, 1, 1)
ax.plot(kernel, '-o', label="kernel", markersize=10)
ax.legend(loc='lower right')
plt.suptitle('WMA kernel')
Text(0.5, 0.98, 'WMA kernel')
Figure

python 实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
class WeightedMovingAverageFilter:
    def __init__(self, window_size):
        self.window_size = window_size
        self.kernal = []
        self.window = []
        normal_Coefficient = window_size * (window_size + 1) / 2
        for i in range(window_size):
            self.kernal.append((i + 1) / normal_Coefficient)       

    def update(self, new_value):
        self.window.append(new_value)
        if len(self.window) > self.window_size:
            self.window.pop(0)
        return sum([self.window[i] * self.kernal[i] for i in range(len(self.window))])

性能分析

1
2
3
4
5
6
7
8
9
moveFilter = MAF.WeightedMovingAverageFilter(5)
EstimateSequence = moveFilter.estimate(MeasSequence[0,:])

### 绘图展示 ###
plt.figure(figsize=figsize)
plt.plot(xSequence, RealSequence[0,:],	'*-', label="Actual", markersize=5)
plt.plot(xSequence, MeasSequence[0,:],	'+-',  label="Measure", markersize=5)
plt.plot(xSequence, EstimateSequence[:],	'o-', label="Estimate", markersize=5)
plt.legend(loc='lower right') #绘制图例
<matplotlib.legend.Legend at 0x77cff0d57350>
Figure

如果初值没有设为 0 的话，应该可以明显注意到前几个点存在迅速下降的巨大误差。这是因为滤波器初始化时窗口内部全零填充导致的，这段时间被称为滤波器的启动时间/warm-up时间。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
performance = FP.FilterPerformance(figsize=figsize)
performance.PreAnalysis()
def PerformanceDemo(SMA_KernalSize):
	moveFilter = MAF.WeightedMovingAverageFilter(int(SMA_KernalSize))
	EstimateSequence = moveFilter.estimate(MeasSequence[0,:])
	performance.ReAnalysis(xSequence ,RealSequence[0,:], MeasSequence[0,:], EstimateSequence)
	plt.show()

p_max = 0
max = 0
for i in range(1,201,1):
	moveFilter = MAF.WeightedMovingAverageFilter(i)
	EstimateSequence = moveFilter.estimate(MeasSequence[0,:])
	t = performance.CalculateGain(RealSequence[0,30:], MeasSequence[0,30:], EstimateSequence[30:])
	if t < max:
		max = t
		p_max = i
widgets.interact(PerformanceDemo, 
                 SMA_KernalSize = FloatSlider(p_max, min=1, max=50 , step=1, description='SMA_KernalSize'))
print(f"最大降噪增益: {max} dB 出现在卷积核大小: {p_max}")
Figure 
interactive(children=(FloatSlider(value=15.0, description='SMA_KernalSize', max=50.0, min=1.0, step=1.0), Outp…


最大降噪增益: -9.105796862143208 dB 出现在卷积核大小: 15

指数滑动滤波(EMA, Exponential Moving Average)

介绍

表达式的推导

指数滑动滤波是以指数形式构造其滤波核的滤波器。各数值的加权影响力随时间而指数式递减，越近期的数据加权影响力越重，但较旧的数据也给予一定的加权值。根据定义直接写出其权重序列 {wi}\{w_i\} 及估值表达式 x^\hat{x}

wi=(1−α)i∑i=0n(1−α)iw_i = \frac{(1-\alpha)^i}{\sum_{i=0}^n{(1-\alpha)^i}}

x^n+1=(1−α)0xn+1+(1−α)1xn+(1−α)2xn−1+⋯+(1−α)nx1∑i=0n(1−α)i\hat{x}_{n+1} = \frac{(1-\alpha)^0x_{n+1} + (1-\alpha)^1x_{n} + (1-\alpha)^2x_{n-1} + \cdots + (1-\alpha)^nx_{1}}{\sum_{i=0}^n{(1-\alpha)^i}}

注意到

lim⁡n−>∞∑i=0n(1−α)i=lim⁡n−>∞1⋅(1−(1−α)n)1−(1−α)=1α\lim_{n->\infty}\sum_{i=0}^n{(1-\alpha)^i} = \lim_{n->\infty} \frac{1 \cdot (1 - (1-\alpha)^n)}{1 - (1-\alpha)} = \frac {1}{\alpha}

所以在极限情况下，估值表达式变为

lim⁡n→∞x^n+1=α⋅(xn+1+(1−α)1xn+⋯+(1−α)nx1)lim⁡n→∞x^n+1=α⋅xn+1+(1−α)⋅[α⋅(xn+(1−α)1xn−1+⋯+(1−α)n−1x1)]lim⁡n→∞x^n+1=α⋅xn+1+(1−α)⋅x^nlim⁡n→∞x^n+1=x^n+α⋅(xn+1−x^n)\begin{align} \nonumber \lim_{n\to \infty} \hat{x}_{n+1} &= \alpha \cdot (x_{n+1} + (1-\alpha)^1x_{n} + \cdots + (1-\alpha)^nx_{1})\\ \nonumber \lim_{n\to \infty} \hat{x}_{n+1} &= \alpha \cdot x_{n+1} + (1-\alpha) \cdot [ \alpha \cdot (x_{n} + (1-\alpha)^1x_{n-1} + \cdots + (1-\alpha)^{n-1}x_{1})]\\ \nonumber \lim_{n\to \infty} \hat{x}_{n+1} &= \alpha \cdot x_{n+1} + (1-\alpha) \cdot \hat{x}_{n}\\ \nonumber \lim_{n\to \infty} \hat{x}_{n+1} &= \hat{x}_{n} + \alpha \cdot (x_{n+1} - \hat{x}_{n}) \end{align}

至此我们已经得到了具有应用价值的指数滤波的递推表达式。

半衰期

1−(1−α)k=12(1−α)k=12k=−ln⁡2/ln⁡(1−α)\begin{align} \nonumber 1 - (1-\alpha)^k & = \frac{1}{2}\\ \nonumber (1-\alpha)^k & = \frac{1}{2}\\ \nonumber k & = -\ln {2} / \ln{(1-\alpha)}\\ \end{align} 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
k = []
com = []
for i in range(1, 100, 1):
    k.append(- np.log(2) / np.log(1 - i/100))
    com.append((1 - i/100) / i/100)

fig = plt.figure(figsize=figsize)
ax = fig.add_subplot(111)
ax.plot([i/100 for i in range(1, 100, 1)], k, '.-g', label="k", markersize=5)
ax.legend(loc='lower left')

ax2 = ax.twinx()
ax2.plot([i/100 for i in range(1, 100, 1)], com, '*-', label="com", markersize=5)
ax2.legend(loc='lower right')

plt.grid(True)
plt.suptitle("Half-life of EMA & Center of Mass")
Text(0.5, 0.98, 'Half-life of EMA & Center of Mass')
Figure
质心 (COM, Center Of Mass)

由质心公式有

∑i=1nximi∑i=1nmi=∑i=0∞i⋅(1−α)i=xCOM\frac{\sum_{i=1}^{n} x_{i} m_{i}}{\sum_{i=1}^{n} m_{i}} = \sum_{i=0}^{\infty} i \cdot (1-\alpha)^i = x_{COM}

由错位相减有

xCOM−(1−α)xCOM=∑i=0∞i⋅(1−α)i−∑i=1∞(i−1)⋅(1−α)i+1α⋅xCOM=(1−α)−lim⁡i→∞(i−1)⋅(1−α)ixCOM=(1−α)α\begin{align} \nonumber x_{COM} - (1-\alpha) x_{COM} & = \sum_{i= 0}^{\infty} i \cdot (1-\alpha)^i - \sum_{i = 1}^{\infty} (i-1) \cdot (1-\alpha)^{i+1}\\ \nonumber \alpha \cdot x_{COM} & = (1-\alpha) - \lim_{i\to\infty}(i-1) \cdot (1-\alpha)^i\\ \nonumber x_{COM} & = \frac{(1-\alpha)}{ \alpha} \end{align}

质心表征了“数据的重心”, 显然质心越小说明估计值越接近测量值，半衰期越小。但问题是质心不够直观，为了让其具有更浅显的物理意义，我们将其映射为 SMA 滤波核长度 L 。
假定 SMA 与 EMA 滤波核具有相同质心，则有

xCOM,EMA=xCOM,SMA1−αα=L2α=2L+2\begin{align} \nonumber x_{COM, EMA} &= x_{COM, SMA}\\ \nonumber \frac{1-\alpha}{\alpha} &= \frac{L}{2}\\ \nonumber \alpha &= \frac{2}{L+2} \end{align} 

1
2
3
4
5
6
alpha = []
for i in range(1, 21):
    alpha.append(2 / (i + 2))

fig = plt.figure(figsize=figsize)
plt.plot(range(1, 21), alpha, '-o', label="alpha", markersize=5)
[<matplotlib.lines.Line2D at 0x77cff0981fa0>]
Figure

python 实现

1
2
3
4
5
6
7
8
9
10
11
class ExponentialMovingAverageFilter:
    def __init__(self, alpha):
        self.alpha = alpha
        self.previous_value = None

    def update(self, new_value):
        if self.previous_value is None:
            self.previous_value = new_value
        else:
            self.previous_value = self.alpha * new_value + (1 - self.alpha) * self.previous_value
        return self.previous_value

性能分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
performance = FP.FilterPerformance(figsize=figsize)
performance.PreAnalysis()
def PerformanceDemo(alpha):
	moveFilter = MAF.ExponentialMovingAverageFilter(alpha)
	EstimateSequence = moveFilter.estimate(MeasSequence[0,:])
	performance.ReAnalysis(xSequence ,RealSequence[0,:], MeasSequence[0,:], EstimateSequence)

p_max = 0.
max = 0.
for i in range(1, 100, 1):
	moveFilter = MAF.ExponentialMovingAverageFilter(float(i)/100.0)
	EstimateSequence = moveFilter.estimate(MeasSequence[0,:])
	t = performance.CalculateGain(RealSequence[0,30:], MeasSequence[0,30:], EstimateSequence[30:])
	if t < max:
		max = t
		p_max = i
widgets.interact(PerformanceDemo, 
                 alpha = FloatSlider(p_max/100, min=0.01, max=1 , step=0.01, description='EMA_alpha'))
print(f"最大降噪增益: {max} dB 出现在alpha: {p_max/100.}")
Figure 
interactive(children=(FloatSlider(value=0.15, description='EMA_alpha', max=1.0, min=0.01, step=0.01), Output()…


最大降噪增益: -8.924687523432539 dB 出现在alpha: 0.15

二阶指数滑动滤波(DEMA, Double Exponential Moving Average)

就是两层指数滑动滤波，不再赘述。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
performance = FP.FilterPerformance(figsize=figsize)
performance.PreAnalysis()
def PerformanceDemo(alpha):
	moveFilter = MAF.DoubleExponentialMovingAverageFilter(alpha)
	EstimateSequence = moveFilter.estimate(MeasSequence[0,:])
	performance.ReAnalysis(xSequence ,RealSequence[0,:], MeasSequence[0,:], EstimateSequence)
	plt.show()

p_max = 0.
max = 0.
for i in range(1, 100, 1):
	moveFilter = MAF.DoubleExponentialMovingAverageFilter(float(i)/100.0)
	EstimateSequence = moveFilter.estimate(MeasSequence[0,:])
	t = performance.CalculateGain(RealSequence[0,30:], MeasSequence[0,30:], EstimateSequence[30:])
	if t < max:
		max = t
		p_max = i/100
widgets.interact(PerformanceDemo, 
                 alpha = FloatSlider(p_max, min=0.01, max=1 , step=0.01, description='DEMA_alpha'))
print(f"最大降噪增益: {max} dB 出现在alpha: {p_max}")
Figure 
interactive(children=(FloatSlider(value=0.06, description='DEMA_alpha', max=1.0, min=0.01, step=0.01), Output(…


最大降噪增益: -10.279182019578673 dB 出现在alpha: 0.06

总结

不难发现，当获得最大降噪增益时，残差往往无法通过正态性检验，且均值也不为 0 。这是由滞后性导致的，接下来我们将着手解决这件事情。

参考文献

Han-Jia Ye. Time Series Analysis. 2022
Diederik P. Kingma and Jimmy Ba, Adam: A Method for Stochastic Optimization[C], 2017, 1412(6980)
Brown. Robert G. et al, The Fundamental Theorem of Exponential Smoothing[J], 1961, 9(5): 673–687