Filter (Sequence)

Filter

「濾波器」。數列到數列的函數。輸入數列、輸出數列。

每個輸出變數分開來看,濾波器其實是由許多函數組成。

簡單的例子是每項加1的濾波器、每項延遲1時刻的濾波器。

當全部函數都相同,僅索引值(時刻)不同,可以簡化成一個函數。

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

Linear Time-Invariant Filter可以寫成卷積形式

數列末端補零K個。

( x(0) ~ x(n)     ) · ( a(0)                           ) = y(0)
( x(0) ~ x(n)     ) · ( a(1) a(0)                      ) = y(1)
        :                                :                   :
( x(0) ~ x(n)     ) · (      a(k) ~~~~~~~~~~ a(0)      ) = y(n-1)
( x(0) ~ x(n)     ) · (           a(k) ~~~~~~~~~~ a(0) ) = y(n)
—————————————————————————————————————————————————————————————————
( x(0) ~ x(n)     ) · (                a(k) ~~~~~ a(1) ) = 0
        :                                          :      :
( x(0) ~ x(n)     ) · (                           a(k) ) = 0

    (x append 0s)   ∗                 a                  = (y append 0s)

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(N²)。

滑動視窗,取加權平均數。時間複雜度O(NK)。

求解(恢復原數列)

resolution (deconvolution): f(?) = y , find x

寫成矩陣形式,化作常對角矩陣求解。時間複雜度O(N²)。

滑動視窗,解加權平均數方程式。時間複雜度O(NK)。

迴歸(求濾波器)

regression (Wiener filter): ?(x) = y , find f

寫成另一種矩陣形式。等式太多,通常無解。改求平方誤差最小的解。套用「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

迴歸(求濾波器):旁門左道解法

數列末端補零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時刻,兩串不同數列的點積)。

反函數(求反濾波器)

inversion: f(x) = y , x = ?(y) , find f⁻¹

濾波器通常不是一對一函數,通常不存在反函數。什麼時候存在反函數呢?我查不到相關資料!

我也沒找到反函數演算法!目前有兩種做法:

一、隨便找一組輸出輸入,對調一下,然後迴歸。

二、濾波器,轉換到頻域,取倒數,轉換回時域。

Filter (Recurrence)

Linear Recurrence

遞迴數列可以視作數列暨濾波器:輸入與輸出是同一數列,但是輸出延遲1時刻。

Linear Recurrence and Linear Time-Invariant Filter

遞迴數列可以做為濾波器的輸入。

Linear Recurrence and Linear Time-Invariant Filter運算

evaluation:                  f(?) = ? >> 1   , find x
regression (autoregression): ?(x) = x >> 1   , find f
resolution (Kalman filter):  { f(?) = ? >> 1 , find x
                             { h(?) = y

三種運算:求值(求遞迴數列)、迴歸(求遞迴函數)、求解(恢復原遞迴數列)。

三種運算應該是首次歸類在一起。

求值(求遞迴數列)

evaluation: f(?) = ? >> 1 , find x

前面章節介紹過了。線性遞迴函數K項,求數列第N項:

一、動態規劃,第0項到第N項。O(NK)。

二、相伴矩陣,矩陣的N次方。O(KωlogN)。

三、多項式乘法、多項式餘數。O(K²logN)或O(KlogKlogN)。

四、生成函數,手工推導公式解。

當K是常數,一變成O(N),二三變成O(logN)。

迴歸(求遞迴函數)

regression (autoregression): ?(x) = x >> 1 , find f

每一個數值,皆是先前緊鄰的K個數值的加權總和,權重是定值。反覆套用線性遞迴函數、代入數列最後N個數值,就能預測下一個即將出現的數值。此行為稱作「線性預測Linear Prediction」。

一串長長的數列,壓縮成一個線性遞迴函數,只需儲存K個函數係數、K個初始數值。反覆套用函數、代入數列最後K個數值,解壓縮得到一串長長的數列。此行為稱作「線性預測編碼Linear Predictive Coding」。

旁門左道的解法,得到Yule-Walker Equation。

視作普通卷積、常對角矩陣,時間複雜度簡單寫成O(N²)。

視作循環卷積、循環矩陣,時間複雜度簡單寫成O(NlogN)。

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)   ]

求解(恢復原遞迴數列)

resolution (Kalman filter): { f(?) = ? >> 1 , find x
                            { h(?) = y

使用測量儀器h,時時收集數據y,時時推敲現況x。推定現況是自迴歸濾波器f,合乎自然規律。自創測量儀器是濾波器h,合乎電路原理。根據結果,推敲原因,此行為稱作「線性二次估計Linear Quadratic Estimation」。

輸出誤差,通過反濾波器,得到輸入誤差,用以修正輸入。

為什麼要大費周章,將觀察值的誤差通過反濾波器,而不是直截了當,將觀察值通過反濾波器?我也不清楚。

h(x) = y
狀態 x₀ x₁ x₂ ...   未知(一串數列,線性遞迴關係)
函數 h₀ h₁ h₂ ...   已知(數列每一項皆有一個函數)
觀察 y₀ y₁ y₂ ...   已知(數列每一項皆有一個輸出)

函數可以一樣               h₀ = h₁ = h₂ = ...
                          
狀態是自迴歸的LTI filter   xₙ = f(xₙ₋₁, xₙ₋₂, ..., xₙ₋ₖ)
                           xₙ = a₁xₙ₋₁ + a₂xₙ₋₂ + ... + aₚxₙ₋ₖ
                
函數可以推廣成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ᵢ⁻¹ᵀ