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

數列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)。

vector<float> x, a, y; void evaluate() { for (int i=0; i<x.size(); ++i) { float sum = 0; for (int j=0; j<a.size(); ++j) if (i-j >= 0) sum += x[i-j] * a[j]; y[i] = sum; } }

求解（反向通過濾波器）

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時刻，兩串不同數列的點積）。

vector<float> correlation(vector<float> x, vector<float> y) { vector<float> xy; xy.resize(x.size()); for (int i=0; i<x.size(); ++i) // gap { float sum = 0; for (int j=0; i+j<y.size(); ++j) sum += x[j] * y[i+j]; xy[i] = sum; } return xy; }

反函數（求反濾波器）

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

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

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

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

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

延伸閱讀：linear time-invariant的定義

linear time-invariant（連續函數）與linear shift-invariant（離散數列）是指同時滿足三種特性：加性、倍性、移位性。

   function of a sequence
   F(x(n)) = y(n)

1. additive function
   F(x₁(n) + x₂(n)) = y₁(n) + y₂(n)

2. homogeneous function of order 1
   F(k x(n)) = k y(n)

3. shiftable function / shift-invariant function
   F(x(n + m)) = y(n + m)

       def
linear === 1 and 2
                       def
linear shift-invariant === 1 and 2 and 3

移位性：輸入數列移位，等同於輸出數列移位。

曾訊號學者提出shiftable這個詞彙。然而大家不採用。

訊號學者採用time-invariant（連續函數）與shift-invariant（離散數列）這兩個詞彙。然而X-invariant有兩種定義。這兩個詞彙屬於第二種定義。

1. F(X(x(n))) = F(x(n)) = y(n)
2. F(X(x(n))) = X(y(n))

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(0)   ⎤   ⎡ xx(1)   ⎤
⎢ xx(1)   xx(0)   ... xx(k-2) ⎥ ⎢ a(1)   ⎥   ⎢ 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)   ⎦

vector<float> x; // 遞迴數列 const int N = 10; // 線性遞迴函數的次方數 float a[N]; // 線性遞迴函數的係數 float xx[N + 1]; // 間距為t的自相關數的總和 void autocorrelation() { for (int i=0; i<=N; i++) { xx[i] = 0; for (int j=0; i+j<x.size(); j++) xx[i] += x[j] * x[i+j]; } } float autoregression() { // 預先計算自相關數的總和 autocorrelation(); float error = xx[0]; for (int i=0; i<N; i++) { // 計算倍率（reflection coefficient） float scale = xx[i + 1]; for (int j=0; j<i; j++) scale -= a[j] * xx[i - j]; scale /= error; // 計算這個階段的係數。 a[i] = scale; for (int j=0; j<i/2; j++) { float d1 = a[j] - scale * a[i-1 - j]; float d2 = a[i-1 - j] - scale * a[j]; a[j] = d1; a[i-1 - j] = d2; } if (i % 2) a[i / 2] -= scale * a[i / 2]; error *= (1.0 - scale * scale); } return error; }

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ᵢ⁻¹ᵀ