一、簡介

在實際的工程應用中，經常會遇到初始結果噪聲太多的問題，比如信號強度抖動的太厲害，比如視頻流中的bbox抖動的太厲害，比如光譜信號抖動的太厲害等等，這時候就需要一些簡單的滑動平均算法。滑動平均其實是一個很樸素的方法，但是要與實際結合，構造出合適的平滑方式，是需要一些思考的。下面我將分別介紹滑動平均法(Moving Average)、指數滑動平均法(Exponential Mean Average)、SG濾波法(Savitzky Golay Filter)。

二、滑動平均法

簡單來說，滑動平均法把前后時刻的一共2n+1個觀測值做平均，得到當前時刻的濾波結果。這是一個比較符合直覺的平滑方法，在生活中、工作中很經常會用到，但是很少去思考這么做的依據是什么，下面我就來仔細分析一下其中的原理。

對于一個觀測序列，我們有這樣的假設：每一次的觀測值是帶有噪聲的，而我們期望噪聲的均值為0，方差為

，觀測值和真實值之間的關系如下：

(1)

其中，

為觀測值，

為真實值，

為噪聲。為了降低噪聲的影響，我們把相鄰時刻的觀測值相加后平均，公式如下：

(2)

表示

時刻的濾波結果，

表示

時刻的觀測值，

代表滑動窗口半徑。將公式(1)代入公式(2)，可以得到

(3)

前面說到了，我們假設噪聲的均值為0，所以

為0，那么我們得到的結果就是:

(4)

當觀測數據的真實值變化較小時，或者變化為線性時，可以近似認為：

(5)

從上面的分析過程我們可以看到，當滑動窗口內的真實數據變化不大的時候，我們可以抑制掉很大一部分噪聲，濾波結果近似真實值；當滑動窗口內的真實值變化較大時，這種濾波方式就會損失一部分精確度，濾波結果接近真實值的平均期望。所以，窗口的大小會對濾波結果有很大影響。窗口越大，濾波結果越平滑，但會一定程度上偏離真實值；窗口越小，濾波結果越接近觀測值，但噪聲偏大。

滑動平均法還有一個升級版本，也就是加權滑動平均法。實際場景中，每個觀測值的重要程度不同，忽略每個觀測值的置信度直接平均不能得到精確的結果，所以就需要給觀測值加權。加權滑動平均法的公式如下：

(6)

為

時刻的權重。(6)式表示的是把每個觀測值乘以權重后再平均。這種方法適用于觀測值本身帶有置信度的情況。注意，這里有一個小問題，如果置信度的取值范圍是0到1之間，那么加權之后計算得到的觀測值往往小于真實值，我來解釋一下為什么。首先，我們假設觀測值和真實值的均值是相等的，也就是

(7)

當我們把觀測值乘以權重了之后，觀測值和真實值的均值就不相等了，因為真實值的權重均值為1，而觀測值的權重均值為

，是小于等于1的，最終的預測值也是小于等于真實值的，而且大概率是小于。所以我們需要對(6)式增加一個修正：

(8)

這樣，得到的預測值就會更加合理了。

小結：滑動平均法使用的前提是，噪聲的均值為0，真實值變化不大或線性變化的場景。如果真實值有較高頻率的非線性突變的話，滑動平均法的效果就不夠好了。同時，滑動平均法的窗口選取很重要，需要根據具體數據來選擇。如果需要使用在線版本的滑動平均，那么就要把窗口前移，也就是把當前時刻的前n個觀測值進行平均，但這樣得到的結果會明顯滯后于當前觀測值，窗口越大，滯后的現象越嚴重。

class MovAvg(object):

def __init__(self, window_size=7):

self.window_size = window_size

self.data_queue = []

def update(self, data):

if len(self.data_queue) == self.window_size:

del self.data_queue[0]

self.data_queue.append(data)

return sum(self.data_queue)/len(self.data_queue)鼠標軌跡的滑動平均效果一維數據的滑動平均效果

三、指數滑動平均法

指數滑動平均法相當于加權滑動平均法的變體，主要區別在于，指數滑動平均法的權重是固定的、隨時間推移呈指數衰減。指數滑動平均法的公式如下：

(9)

表示預測值，

表示衰減權重，通常我們設為固定值0.9，

表示觀測值，這是一個遞推公式。前面說了，指數滑動平均法的權重是隨時間推移呈指數衰減的，那么上面的這個遞推公式的指數體現在哪里呢？我們把(9)式進行延伸：

(10)

將(9)和(10)兩式子聯立，可得

(11)

發現沒有，在(11)式中

與

的關系是

倍，而在(9)式中

與

的關系是

倍，呈指數衰減關系。同時，在初始時刻有如下關系：

(12)

根據這一關系和上述的遞推公式，我們就能夠得到整個算法的公式了。

由于這種指數衰減的特性，指數滑動平均法會比滑動平均法的實時性更強，更加接近當前時刻的觀測值。在實際場景下，如果目標的波動較大時，指數滑動平均法會比滑動平均更加接近當前的真實值。那么是不是就說明，指數滑動平均法在任意場景下都比滑動平均法更好呢？不一定。我們來分析一下指數衰減法的誤差項，這里為了簡便表示，設定

，同時，將(1)式和(12)式代入公式(11)，可得到誤差項：

(13)

所以誤差項也是呈指數衰減的，越接近當前時刻的誤差項權重越大。假如在當前的工程場景中，誤差是固定的分布，不受目標的觀測值大小影響的話，那么指數滑動平均法會更接近真實值；假如誤差會受目標觀測值影響，比如我們觀測的是一個連續運動的目標，中間突然出現了一個偏離很遠的觀測點，那么這個點為誤檢的概率相當大，也就是該觀測值的誤差比之前其他點的誤差要大得多，此時指數加權平均法的結果就會波動較大，結果就不如滑動平均了。

小結：當誤差不受觀測值大小影響的話，指數滑動平均比滑動平均好；當誤差隨觀測值大小變化時，滑動平均比指數滑動平均更好。

class ExpMovAvg(object):

def __init__(self, decay=0.9):

self.shadow = 0

self.decay = decay

self.first_time = True

def update(self, data):

if self.first_time:

self.shadow = data

self.first_time = False

return data

else:

self.shadow = self.decay*self.shadow+(1-self.decay)*data

return self.shadow鼠標軌跡的指數滑動平均效果一維數據的指數滑動平均效果

四、SG濾波法

SG濾波法(Savitzky Golay Filter)的核心思想也是對窗口內的數據進行加權濾波，但是它的加權權重是對給定的高階多項式進行最小二乘擬合得到。它的優點在于，在濾波平滑的同時，能夠更有效地保留信號的變化信息，下面我來介紹一下其原理。

我們同樣對當前時刻的前后一共2n+1個觀測值進行濾波，用k-1階多項式對其進行擬合。對于當前時刻的觀測值，我們用下面的公式進行擬合：

(14)

同樣，對于前后時刻(如t-1、t+1、t-2、t+2等時刻)的預測值，我們同樣可以用(14)式來計算，這樣一共得到2n+1個式子，構成一個矩陣(似乎發不了矩陣，我放個圖片吧)：

要使得整個矩陣有解，必須滿足 2n+1>k，這樣我們才能夠通過最小二乘法確定參數

、

、

...

。我們把上面的矩陣簡化表示為下面公式：

(15)

各個參數下標表示它們各自的維度，如

表示有k行1列的參數。通過最小二乘法，我們可以求得

的解為：

(16)

上標trans表示轉置。那么，模型的濾波值為：

(17)

最終可以得到濾波值和觀測值之間的關系矩陣：

(18)

算出了B矩陣，我們就能夠快速的將觀測值轉換為濾波值了。

小結：SG濾波法對于數據的觀測信息保持的更好，在一些注重數據變化的場合會比較適用。

class SavGol(object):

def __init__(self, window_size=11, rank=2):

assert window_size % 2 == 1

self.window_size = window_size

self.rank = rank

self.size = int((self.window_size - 1) / 2)

self.mm = self.create_matrix(self.size)

self.data_seq = []

def create_matrix(self, size):

line_seq = np.linspace(-size, size, 2*size+1)

rank_seqs = [line_seq**j for j in range(self.rank)]

rank_seqs = np.mat(rank_seqs)

kernel = (rank_seqs.T * (rank_seqs * rank_seqs.T).I) * rank_seqs

mm = kernel[self.size].T

return mm

def update(self, data):

self.data_seq.append(data)

if len(self.data_seq) > self.window_size:

del self.data_seq[0]

padded_data = self.data_seq.copy()

if len(padded_data) < self.window_size:

left = int((self.window_size-len(padded_data))/2)

right = self.window_size-len(padded_data)-left

for i in range(left):

padded_data.insert(0, padded_data[0])

for i in range(right):

padded_data.insert(

len(padded_data), padded_data[len(padded_data)-1])

return (np.mat(padded_data)*self.mm).item()鼠標軌跡的SG濾波效果一維數據的SG濾波效果

附錄

本文圖片制作的相關代碼。

import cv2

import numpy as np

import matplotlib.pyplot as plt

import imageio

# 一維數據濾波

ma, ema, sg = MovAvg(), ExpMovAvg(), SavGol()

data_list, data_ma, data_ema, data_sg = [], [], [], []

for i in range(200):

data = i+np.random.randint(-50, 50)

data_list.append(data)

data_ma.append(ma.update(data))

data_ema.append(ema.update(data))

data_sg.append(sg.update(data))

plt.plot(data_list, label='raw')

plt.plot(data_ma, label='ma')

plt.plot(data_ema, label='ema')

plt.plot(data_sg, label='sg')

plt.legend()

plt.show()

# 鼠標軌跡濾波

ma_x, ma_y = MovAvg(), MovAvg()

ema_x, ema_y = ExpMovAvg(), ExpMovAvg()

sg_x, sg_y = SavGol(), SavGol()

def draw_circle(event, x, y, flags, param):

if event == cv2.EVENT_MOUSEMOVE:

sx = np.random.randint(-50, 51)

sy = np.random.randint(-50, 51)

cv2.circle(show, (x+sx, y+sy), 5, (255, 255, 255), -1)

x, y = ma_x.update(x+sx), ma_y.update(y+sy)

cv2.circle(show, (int(x), int(y)), 5, (0, 0, 255), -1)

x, y = ema_x.update(x+sx), ema_y.update(y+sy)

cv2.circle(show, (int(x), int(y)), 5, (0, 255, 0), -1)

x, y = sg_x.update(x+sx), sg_y.update(y+sy)

cv2.circle(show, (int(x), int(y)), 5, (255, 0, 0), -1)

show = np.zeros((1024, 1024, 3), np.uint8)

cv2.namedWindow('image')

buff = []

while True:

cv2.setMouseCallback('image', draw_circle)

cv2.imshow('image', show)

save = cv2.resize(show, (512, 512))

save = cv2.cvtColor(save, cv2.COLOR_BGR2RGB)

buff.append(save)

if cv2.waitKey(100) == ord('q'):

break

cv2.destroyAllWindows()

imageio.mimwrite('test.gif', buff, 'GIF', duration=0.1)

總結

以上是生活随笔為你收集整理的python 数据平滑_数据平滑方法的原理和应用的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。