Filter
Filter
「濾波器」。輸入一串數列、輸出一串數列的函數。
每個輸出變數分開來看,濾波器其實是由許多函數組成。
簡單的例子是每項加1的濾波器、每項延遲1時刻的濾波器。
當全部函數都相同,僅索引值(時刻)不同,可以簡化成一個函數。
Linear Time-Invariant Filter
Linear Time-Invariant Filter(LTI Filter)
「線性非時變濾波器」。濾波器是相同的線性函數。
數學家和計算學家擅長線性函數,容易計算。
數列x,通過濾波器a,得到數列y。
Linear Time-Invariant Filter可以寫成線性組合形式
濾波器視作權重。
y(n) = a(0)x(n) + a(1)x(n-1) + a(2)x(n-2) + ... + a(k)x(n-k)
Linear Time-Invariant Filter可以寫成矩陣形式
濾波器視作矩陣。
[ a(0) 0 0 ... 0 0 0 ] [ x(0) ] [ y(0) ] [ a(1) a(0) 0 ... 0 0 0 ] [ x(1) ] [ y(1) ] [ a(2) a(1) a(0) ... 0 0 0 ] [ : ] [ : ] [ : : : : : : ] [ : ] = [ : ] [ 0 0 0 ... a(0) 0 0 ] [ : ] [ : ] [ 0 0 0 ... a(1) a(0) 0 ] [ : ] [ : ] [ 0 0 0 ... a(2) a(1) a(0) ] [ x(n) ] [ y(n) ] A x y
數列x視作矩陣。
[ x(0) 0 ... 0 0 ] [ y(0) ] [ x(1) x(0) ... 0 0 ] [ a(0) ] [ y(1) ] [ : : : : ] [ a(1) ] [ : ] [ : : : : ] [ : ] = [ : ] [ : : : : ] [ : ] [ : ] [ x(n-1) x(n-2) ... x(n-k) x(n-k-1) ] [ a(k) ] [ : ] [ x(n) x(n-1) ... x(n-k+1) x(n-k) ] [ y(n) ] X a y
Linear Time-Invariant Filter可以寫成卷積形式
數列末端補零K個。
( x(0) .. x(n) 0 .. 0 ) · ( a(0) ) = y(0) ( x(0) .. x(n) 0 .. 0 ) · ( a(1) a(0) ) = y(1) : : : ( x(0) .. x(n) 0 .. 0 ) · ( a(k) ....... a(0) ) = y(n-1) ( x(0) .. x(n) 0 .. 0 ) · ( a(k) ....... a(0) ) = y(n) ————————————————————————————————————————————————————————————————— ( x(0) .. x(n) 0 .. 0 ) · ( a(k) .. a(1) ) = 0 : : : ( x(0) .. x(n) 0 .. 0 ) · ( a(k) ) = 0 x' ∗ a = y'
Linear Time-Invariant Filter運算
evaluation (moving average): f(x) = ? , find y resolution (deconvolution): f(?) = y , find x regression (Wiener filter): ?(x) = y , find f inversion: f(x) = y , x = ?(y) , find f⁻¹
四種運算:求值(通過濾波器)、求解(反向通過濾波器)、迴歸(求濾波器)、反函數(求反濾波器)。
四種運算都有時域和頻域兩種解法。頻域解法的時間複雜度較低,但是沒人用!一、濾波器長度通常很短(函數項數很少),傅立葉轉換反而浪費時間。二、濾波器無法完美轉換到頻域。
求值(通過濾波器)
evaluation (moving average): f(x) = ? , find y
滑動視窗,取加權平均數。時間複雜度O(NK)。
求解(反向通過濾波器)
resolution (deconvolution): f(?) = y , find x
滑動視窗,解加權平均數方程式。時間複雜度O(NK)。
迴歸(求濾波器)
regression (Wiener filter): ?(x) = y , find f
數列x視作矩陣。等式太多,通常無解。改為找到平方誤差最小的解,套用「Normal Equation」。時間複雜度O(N³)。
[ x(0) 0 ... 0 0 ] [ y(0) ] [ x(1) x(0) ... 0 0 ] [ a(0) ] [ y(1) ] [ : : : : ] [ a(1) ] [ : ] [ : : : : ] [ : ] = [ : ] [ : : : : ] [ : ] [ : ] [ x(n-1) x(n-2) ... x(n-k) x(n-k-1) ] [ a(k) ] [ : ] [ x(n) x(n-1) ... x(n-k+1) x(n-k) ] [ y(n) ] X a y 對應項平方誤差總和最小,即是滿足 Xᵀ X a = Xᵀ y。 a = (Xᵀ X)⁻¹ Xᵀ y
旁門左道的解法是Wiener–Hopf Equation。
數列末端補零K個,讓矩陣變漂亮。答案不對,但是算得快。
一、建立矩陣:卷積。O(NK)甚至O(NlogN)。
二、矩陣求解:循環矩陣求解。O(K²)甚至O(KlogK)。
[ x(0) ] [ y(0) ] [ : x(0) ] [ : ] [ : : ] [ a(0) ] [ : ] [ : : x(0) ] [ : ] [ : ] [ : : ..... : x(0) ] [ : ] = [ : ] [ x(n) : : : ] [ : ] [ y(n) ] [ x(n) : : ] [ a(k) ] [ 0 ] [ : : ] [ : ] [ x(n) : ] [ : ] [ x(n) ] [ 0 ] X a y [ xx(0) xx(1) ... xx(k) ] [ a(0) ] [ xy(0) ] [ xx(1) xx(0) ... xx(k-1) ] [ a(1) ] [ xy(1) ] [ : : : ] [ : ] = [ : ] [ xx(k-1) xx(k-2) ... xx(1) ] [ : ] [ : ] [ xx(k) xx(k-1) ... xx(0) ] [ a(k) ] [ xy(k) ] Xᵀ X a Xᵀ y xx(t) = ∑ x(i-t)⋅x(i) autocorrelation of x(i-t) and x(i). xy(t) = ∑ x(i-t)⋅y(i) cross-correlation of x(i-t) and y(i). 對應項平方誤差總和最小,即是滿足 Xᵀ X a = Xᵀ y。 Xᵀ X 和 Xᵀ y 乘起來,稱作Wiener–Hopf Equation。 xx(t):數列x延遲t時刻、數列x,兩者的自相關數(相差t時刻,兩串相同數列的點積)。 xy(t):數列x延遲t時刻、數列y,兩者的互相關數(相差t時刻,兩串不同數列的點積)。
Linear Recurrence
Linear Recurrence
遞迴數列可以視作數列暨濾波器:輸入與輸出是同一數列,但是輸出延遲1時刻。
Linear Recurrence運算
evaluation: f(?) = ? >> 1 , find x autoregression: ?(x) = x >> 1 , find f
兩種運算:求值(求遞迴數列)、自迴歸(求遞迴函數)。
求值(求遞迴數列)
evaluation: f(?) = ? >> 1 , find x
線性遞迴函數K項,求數列第N項,總共三種方法:
一、動態規劃,第0項算到第N項。O(NK)。
二、同伴矩陣的N次方。O(K³logN)。假設矩陣相乘O(K³)。
三、xᴺ模特徵多項式。O(K²logN)甚至O(KlogKlogN)。
當K是常數,一變成O(N),二三變成O(logN)。
自迴歸(求遞迴函數)
autoregression: ?(x) = x >> 1 , find f
每一個數值,皆是先前緊鄰的K個數值的加權總和,權重是定值。反覆套用線性遞迴函數、代入數列最後N個數值,就能預測下一個即將出現的數值。此行為稱作「線性預測Linear Prediction」。
一串長長的數列,壓縮成一個線性遞迴函數,只需儲存K個函數係數、K個初始數值。反覆套用函數、代入數列最後K個數值,解壓縮得到一串長長的數列。此行為稱作「線性預測編碼Linear Predictive Coding」。
演算法是Berlekamp–Massey Algorithm。時間複雜度O(NK)。
旁門左道的解法是Yule–Walker Equation。
數列末端補零K個,讓矩陣變漂亮。答案不對,但是算得快。
一、建立矩陣:卷積。O(NK)甚至O(NlogN)。
二、矩陣求解:循環矩陣求解。O(K²)甚至O(KlogK)。
x(n) = a(1)x(n-1) + a(2)x(n-2) + ... + a(k)x(n-k) [ xx(0) xx(1) ... xx(k-1) ] [ a(1) ] [ xx(1) ] [ xx(1) xx(0) ... xx(k-2) ] [ a(2) ] [ xx(2) ] [ : : : ] [ : ] = [ : ] [ xx(k-2) xx(k-1) ... xx(1) ] [ a(k-1) ] [ xx(k-1) ] [ xx(k-1) xx(k-2) ... xx(0) ] [ a(k) ] [ xx(k) ]
Kalman Filter
Linear Recurrence and Linear Time-Invariant Filter
遞迴數列可以作為濾波器的輸入。
Linear Recurrence and Linear Time-Invariant Filter運算
resolution (Kalman filter): { f(?) = ? >> 1 , find x { h(?) = y
求解。求得原本遞迴數列。
Kalman Filter
使用測量儀器h,時時收集數據y,時時推敲現況x。推定現況是自迴歸濾波器f,合乎自然規律。自創測量儀器是濾波器h,合乎電路原理。根據結果,推敲原因,此行為稱作「線性二次估計Linear Quadratic Estimation」。
輸出誤差,通過反濾波器,得到輸入誤差,用以修正輸入。
為什麼要大費周章,將觀察值的誤差通過反濾波器,而不是直截了當,將觀察值通過反濾波器?我也不清楚。
h(x) = y 狀態 x₀ x₁ x₂ ... 未知(線性遞迴數列) 函數 h₀ h₁ h₂ ... 已知(數列每一項皆有一個函數) 觀察 y₀ y₁ y₂ ... 已知(數列每一項皆有一個輸出) 狀態是自迴歸的LTI filter xₙ = f(xₙ₋₁, xₙ₋₂, ..., xₙ₋ₖ) xₙ = a₁xₙ₋₁ + a₂xₙ₋₂ + ... + aₚxₙ₋ₖ 函數可以一樣 h₀ = h₁ = h₂ = ... 函數可以推廣成LTI filter yₙ = hₙ(xₙ, xₙ₋₁, xₙ₋₂, ..., xₙ₋ₚ) yₙ = bₙ₀xₙ + bₙ₁xₙ₋₁ + ... + bₙₚxₙ₋ₚ 從 x₁ 開始一個一個算到 xᵢ,現在要算 xᵢ 是多少: 1. x̂ᵢ = f(xᵢ₋₁) x 數列估計下一項。 (x₀ = 0) 此處以一階函數為例。 2. ŷᵢ = hᵢ(x̂ᵢ) y 數列估計下一項。 盜版 ŷᵢ 通常不等於正版 yᵢ。 這表示 x̂ᵢ 不夠準,需要修正。 3. kᵢ⁻¹ 求得誤差的反濾波器 kᵢ⁻¹。稱作Kalman gain。 kᵢ(x̂ᵢ - xᵢ) = (ŷᵢ - yᵢ) (1) pᵢ = (x̂ᵢ - xᵢ)² 為了計算反濾波器,需要 xᵢ 的誤差的平方。 = cov(x̂ᵢ - xᵢ) 即是誤差的covariance。 sᵢ = (ŷᵢ - yᵢ)² 不過由於我們還沒求出 xᵢ,所以不能這樣算。 = cov(ŷᵢ - yᵢ) (2) p̂ᵢ = f²(pᵢ₋₁) p 數列估計下一項。 (p₀ = 0) = f pᵢ₋₁ fᵀ 原理是 Var[k⋅X] = k² ⋅ Var[X]。 (3) ŝᵢ = hᵢ²(p̂ᵢ) s 數列估計下一項。 = hᵢ p̂ᵢ hᵢᵀ 如果 f h 各有誤差項:常態分布 q r,平均數0。 那麼 p̂ᵢ ŝᵢ 各加上 var(q) var(r)。 (4) kᵢ⁻¹ = p̂ᵢ hᵢᵀ ŝᵢ⁻¹ 反濾波器是變異數的比值。這是什麼鬼啦? 4. xᵢ = x̂ᵢ - kᵢ⁻¹(ŷᵢ - yᵢ) 修正 x。 yᵢ 的誤差,通過反濾波器,還原成 xᵢ 的誤差。 x̂ᵢ 補上誤差,答案就對了。 5. pᵢ = p̂ᵢ - kᵢ⁻²(ŝᵢ) 修正 p。 = p̂ᵢ - kᵢ⁻¹ ŝᵢ kᵢ⁻¹ᵀ