SG濾波原理

 

作者：藏雲閣主
鏈接：https://zhuanlan.zhihu.com/p/78848809

 

 

一、簡介

在實際的工程應用中，經常會遇到初始結果噪聲太多的問題，比如信號強度抖動的太厲害，比如視頻流中的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)