引言
訊號學家自創一個詞彙system identification。
標題本來應該是system parameter estimation。
硬要區分的話嘛：
identification是找到系統模型。就是建模！
estimation是找到系統參數。就是迴歸！
system model
FIR system / IIR system
開迴路系統、閉迴路系統
訊號學家自創兩個同義詞彙，就是這樣而已。
脈衝響應分成兩種：
1. finite impulse response (FIR)
有限脈衝響應。輸入脈衝函數，輸出很快歸零。時間長度有限。
2. infinite impulse response (IIR)
無限脈衝響應。輸入脈衝函數，輸出永不歸零。時間長度無限。
系統分為兩種款式：
1. FIR system = open-loop system
開迴路->輸入只取幾項->有限脈衝響應
2. IIR system = closed-loop system
閉迴路->輸入包含輸出->輸出強行展開->輸入取所有項->無限脈衝響應
LTI system
系統分為兩種款式：
1. LTI FIR system = MA model
開迴路->輸入的加權總和->輸入只取幾項->有限脈衝響應
2. LTI IIR system = ARMA model
閉迴路->輸入與輸出的加權總和->輸出強行展開->輸入取所有項->無限脈衝響應
stochastic system
系統分為兩種款式：
1. LTI system
沒有雜訊／干擾->整體視作一個LTI system->system operation
2. stochastic LTI system
追加雜訊／干擾->考慮各種system diagram->system identification
system identification
stochastic LTI system
stochastic LTI system
├ stochastic LTI FIR system
│ └ MAX model Y = AX + E skip (z)
└ stochastic LTI IIR system
├ output error model Y = (A/B)X + E
├ ARX model Y = (A/B)X + (1/B)E
├ ARMAX model Y = (A/B)X + (C/B)E
└ Box–Jenkins model Y = (A/B)X + (C/D)E
system model identification
系統模型識別的步驟如下：
一、判斷系統是線性非時變系統／不是線性非時變系統：
依序檢查因果性、時間不變性、加性、倍性。
令輸入訊號是特定函數，
觀察輸出訊號是否不受控制、隨之延遲、相加、翻倍。
二、判斷系統是有限脈衝響應／無限脈衝響應：
令輸入訊號是脈衝函數，觀察輸出訊號。
甲、輸出訊號迅速歸零：FIR system。
乙、輸出訊號永不歸零：IIR system。
一般使用脈衝響應。再不濟，矩形響應、三角形響應。
三、決定系統模型：
stochastic LTI FIR system只有一種基礎模型。
stochastic LTI IIR system擁有四種基礎模型。
四種基礎模型通通嘗試一遍，看看哪種誤差較少。
你也可以自己發明新模型。
四、決定系統參數：
各種系統參數數量通通嘗試一遍，看看哪種誤差最少。
另外還要檢查系統延遲時間、穩定性。
詳情請見講義：
http://mocha-java.uccs.edu/ECE5560/ECE5560-Notes04.pdf
system parameter estimation
(1) convolution kernel: sequence deconvolution
y = f * x
Xᵀ X f = Xᵀ y
f = (Xᵀ X)⁻¹ Xᵀ y
(2) transfer function: polynomial division & interpolation
F(exp(𝑖ω)) = Y(exp(𝑖ω)) / X(exp(𝑖ω))
系統參數估計的演算法，原理只有兩種，時域和頻域。
一、時域卷積核：
已知輸入訊號、輸出訊號，求得系統參數。
如果已知系統參數數量，那麼系統參數數值有唯一解。
因為對象是LTI system，所以形成linear equation。
高斯消去法可以求解。
現實世界的訊號數值，無法完美精確地測量，總是有雜訊／干擾。
大家習慣改用least squares method，找到平方誤差最小的解。
因為對象是LTI system，所以形成linear least squares。
normal equation可以求解。
二、頻域轉移函數／系統頻譜：
已知輸入訊號頻譜、輸出訊號頻譜，求得系統頻譜。
如果已知系統參數數量，那麼系統頻譜有唯一解。
多項式除法與多項式內插可以求解。
訊號頻譜有雜訊／干擾，那麼改用least squares method。
最佳化演算法可以求解。
experimental data
input signal / system response:
輸入特殊訊號，直接量頻譜。
(1) chirp (swept sine): 弦波頻率漸增，依序得到輸出訊號頻譜每個bin。
(2) white noise: 訊號頻譜是常數函數，方便計算系統的gain。
(3) pseudorandom binary sequence：針對數位訊號。功能類似white noise。
steady state / initial state:
一、進行實驗之時，確保輸出訊號已經抵達穩態，才做測量。
二、進行實驗之前，確保系統內部狀態已經恢復初始值，才做測量。
連續進行實驗的情況下，
輸入訊號需要插入足夠多個零。
甚至切斷電源重開機。
尤其是stochastic LTI IIR system。
model selection / model validation
模型選擇：找到最符合的系統模型與系統參數。
指標有AIC、BIC。方法有cross-validation。此處省略。
模型驗證：承上，接著檢查該系統模型與系統參數。
嘗試各種輸入訊號，檢查實際系統與估計系統的輸出訊號是否足夠相符。
指標有平均數、變異數。方法有cross-validation。此處省略。
stochastic LTI FIR system
system model
MAX model:
e (zero-mean white noise)
╷
┌─────┐ ↓
x ───→│ A │──→⊕──→ y
└─────┘
y[n] = sum { aₖ x[n-k] + e[n] }
k=0⋯∞
MAX model = Moving Average model with eXogenous inputs
移動平均模型MA附帶外生輸入X
（此處的外生輸入是指zero-mean white noise）
system parameter estimation
system parameter estimation:
(1) correlation: least squares estimation
input signal has time-invariant autocorrelation.
e.g. weakly stationary process
system spectrum estimation:
(2) correlation spectrum: H₁ estimate and H₂ estimate
input signal has time-invariant autocorrelation.
e.g. weakly stationary process
(3) frequency response: QAM estimate
input signal is sinusoid.
e.g. cosine wave
(4) transfer function: empirical transfer function estimate
input signal is purpose-built.
e.g. white noise
針對LTI FIR model，
可以直接估計系統參數，
也可以間接估計系統頻譜。
LTI FIR model = MA model。
卷積核恰是系統參數。
系統頻譜做逆向傅立葉轉換得到卷積核。
least squares estimation
correlation
autocorrelation function:
rₓₓ[k] = sum { x[n+k] x[n] }
n=0⋯∞
cross-correlation function:
rₓ[k] = sum { x[n+k] y[n] }
n=0⋯∞
實務上訊號長度有限。訊號長度是N，加總運算範圍是n=0⋯N-1。
property:
(1) rₓ[k] ≠ rₓ[k] not commute
(2) rₓ[k] = rₓ[-k] however negative index is not defined
(3) rₓₓ ∗ a = rₓ a is LTI system that y = a ∗ x
(4) -⃡a ∗ rₓₓ ∗ a = r a is LTI system that y = a ∗ x
proof of property (3):
y[n] = sum { aₖ x[n-k] }
k=0⋯∞
rₓ[k] = sum { y[n+k] x[n] }
n
= sum { sum { aᵢ x[n+k-i] } x[n] }
n i
= sum { sum { aᵢ x[n+k-i] x[n] } }
n i
= sum { sum { aᵢ x[n+k-i] x[n] } }
i n
= sum { aᵢ sum { x[n+k-i] x[n] } }
i n
= sum { aᵢ rₓₓ[k-i] }
i
least squares estimation
MAX model:
e (zero-mean white noise)
╷
┌─────┐ ↓+
x ───→│ A │──→⊕──→ y
└─────┘
y[n] = sum { aₖ x[n-k] + e[n] }
k=0⋯∞
assumption:
assume e and x are independent.
theorem:
independent => uncorrelated.
rₑₓ = 0
rₑₓ[k] = sum { e[n+k] x[n] } = 0
correlation:
假設e與x獨立、不相關，就可以用correlation消除雜訊影響。
rₓ = rₓₓ ∗ a + rₑₓ
= rₓₓ ∗ a
rₓ[k] = sum { y[n+k] x[n] }
= sum { sum { aᵢ x[n+k-i] + e[n+k] } x[n] }
= sum { sum { aᵢ x[n+k-i] x[n] + e[n+k] x[n] } }
= sum { sum { aᵢ x[n+k-i] x[n] } + sum { e[n+k] x[n] } }
= sum { aᵢ rₓₓ[k-i] } ^^^^^^^^^^^^^^^^^^^
= 0
linear regression:
解一次方程組得到系統參數。
即是先前章節介紹的方法。
given rₓ and rₓₓ, solve a.
assume aₜ = 0 when t ≥ N
⎡ rₓ[0] ⎤ ⎡ rₓₓ[0] ... rₓₓ[N-1] ⎤ ⎡ a₀ ⎤
⎢ : ⎥ = ⎢ : : ⎥ ⎢ : ⎥
⎣ rₓ[N-1] ⎦ ⎣ rₓₓ[-(N-1)] ... rₓₓ[0] ⎦ ⎣ aɴ₋₁ ⎦
rₓ Toeplitz(rₓₓ) a
solution:
a = Toeplitz(rₓₓ)⁻¹ rₓ
H₁ estimate and H₂ estimate
correlation spectrum
符號：z轉換或者拉普拉斯轉換。
數值：傅立葉轉換。
(1) symbol: z-transform / Laplace transforms
autocorrelation spectrum:
Rₓₓ(z) = rₓₓ[0] z⁰ + rₓₓ[1] z⁻¹ + rₓₓ[2] z⁻² + ...
= sum { rₓₓ[k] z⁻ᵏ }
k=0⋯∞
cross-correlation spectrum:
Rₓ(z) = rₓ[0] z⁰ + rₓ[1] z⁻¹ + rₓ[2] z⁻² + ...
= sum { rₓ[k] z⁻ᵏ }
k=0⋯∞
(2) numeral: discrete-time Fourier transform
autocorrelation spectrum:
Rₓₓ(ω) = sum { rₓₓ[k] exp(-𝑖ωk) }
k=-∞⋯+∞
cross-correlation spectrum:
Rₓ(ω) = sum { rₓ[k] exp(-𝑖ωk) }
k=-∞⋯+∞
property (derived from convolution theorem):
(1) Rₓ(z) ≠ Rₓ(z) not commute
(2) Rₓ(z) = Rₓ(1/z) however not being calculated
(3) Rₓₓ(z) A(z) = Rₓ(z) A is LTI system that Y = AX
(4) A(1/z) Rₓₓ(z) A(z) = R(z) A is LTI system that Y = AX
H₁ estimate and H₂ estimate
e
╷
┌─────┐ ↓+
x ───→│ A │──→⊕──→ y
└─────┘
rₑₑ[k] = σ² δ[k] since e is white
Rₑₑ(ω) = sum { σ² δ[k] exp(-𝑖ωk) } = σ² since e is white
r[k] = aₖ ∗ rₓ[k] + σ² δ[k]
R(ω) = A(exp(𝑖ω)) Rₓ(ω) + Rₑₑ(ω) where Rₑ(ω) = σ²
d e
╷ ╷
↓+ ┌─────┐ ↓+
x ──→⊕──→│ A │──→⊕──→ y
└─────┘
H₁ estimate: Ĥ₁ = Rₓ / Rₓₓ = (A Rₓₓ) / Rₓₓ + Rdd
H₂ estimate: Ĥ₂ = R / Rₓ = (A Rₓ + Rₑₑ) / Rₓ
inequality: |Ĥ₁| ≤ |A| ≤ |Ĥ₂|
https://dsp.stackexchange.com/questions/71811/
QAM estimate【查無正式學術名稱】
frequency response
理論上是輸入複數exp波，實務上是輸入實數cos波。
輸入餘弦波，頻率ω、振幅α、相位0。
調整輸入振幅α，避免輸出訊號太弱太強而測量不到。
e
╷
α cos(ωn) ┌─────┐ ↓+
x ──────────→│ A │──→⊕──→ y
└─────┘
x[n] = α cos(ωn)
y[n] = α |A(exp(𝑖ω))| cos(ωn + ∠A(exp(𝑖ω))) + e[n]
QAM estimate
如果沒有誤差，那麼很容易計算系統頻譜。
為了應付誤差，輸出做amplitude modulation。
輸出分別乘上餘弦波和正弦波，反推原始振幅、原始相位，
稱作quadrature amplitude modulation。
cos(ωn)
╷
e ↓× ┌───────────────┐
╷ ┌─→⊕─→│N-point average│──→ Ic(ω)
α cos(ωn) ┌─────┐ ↓+ │ └───────────────┘
x ──────────→│ A │─→⊕─→ y ─┤
└─────┘ │ ┌───────────────┐
└─→⊕─→│N-point average│──→ Is(ω)
↑× └───────────────┘
╵
sin(ωn)
x[n] = α cos(ωn)
y[n] = α |A(exp(𝑖ω))| cos(ωn + ∠A(exp(𝑖ω))) + e[n]
Ic(ω) = sum { y[n] cos(ωn) } = +½ α |A(exp(𝑖ω))| cos(ωn)
n=1⋯N
Is(ω) = sum { y[n] sin(ωn) } = -½ α |A(exp(𝑖ω))| sin(ωn)
n=1⋯N
積化和差公式
cos(a) cos(b) = ½ cos(a-b) + ½ cos(a+b)
推導過程
Ic(ω) = sum { y[n] cos(ωn) }
= sum { (......) cos(ωn) }
= ½ α |A(exp(𝑖ω))| cos(ωn) ①
+ ½ α |A(exp(𝑖ω))| (1/N) sum { cos(2ωn + ∠A(exp(𝑖ω))) } ②
+ (1/N) sum { e[n] cos(ωn) } ③
②→0 as n→∞. since cos() has zero mean.
③→0 as n→∞. since e and x are independent by assumption.
hence only ① remains.
QAM estimate:
|Â(exp(𝑖ω))| = sqrt(Ic²(ω) + Is²(ω)) / (α/2)
∠Â(exp(𝑖ω)) = -tan⁻¹(Is(ω) / Ic(ω))
empirical transfer function estimate
專著《System Identification: Theory for the User》。
empirical transfer function estimate
MAX model:
Y = AX + E skip (z)
empirical transfer function estimate:
 = Y/X
阿就測量一下輸出頻譜、輸入頻譜，兩者相除，即得系統頻譜。
重劍無鋒大巧不工。
系統頻譜＝輸出頻譜／輸入頻譜
 = Y/X
頻譜相除＝振幅相除＆相位相減
複數相除＝長度相除＆角度相減
|Â| = |Y| / |X| amplitude
∠Â = ∠Y - ∠X phase
mean / variance / covariance
empirical transfer function estimate:
 = Y/X = (AX + E)/X = A + E/X
statistics:
(1) mean: E[Â] = E[A + E/X] = E[A] + E[E/X] = E[A]
(2) variance: E[|Â-A|²] = (Rₑₑ + constant) / E[|X|²]
(3) covariance:
E[(Â(exp(𝑖ω₁))-A(exp(𝑖ω₁)))* (Â(exp(𝑖ω₂))-A(exp(𝑖ω₂)))] = 0
assumptions:
(1) zero-mean noise: E[E(exp(𝑖ω))] = 0
(2) white noise: E[E(exp(𝑖ω₁))E(exp(𝑖ω₂))] = 0
(3) independence: E[X(exp(𝑖ω))E(exp(𝑖ω))] = 0
(4) spectrum of system is asymptotically uncorrelated:
E[A(exp(𝑖ω₁))A(exp(𝑖ω₂))] = 0 when M→∞
推導過程省略。
頻譜都是複數，相乘之前記得取共軛複數。
複數乘以共軛複數，恰是絕對值平方。
ETFT的答案絕對是錯的。
畢竟一眼看上去就是亂算一通，根本不考慮誤差項。
ETFT的答案的某些統計學指標至少是對的：
一、平均值是對的。
二、變異數是錯的，但是誤差有上限。
三、共變異數是對的。
（當系統頻譜的頻率種類M趨近無限多的情況下。）
即使統計學指標是對的，
那也只是自我安慰、精神勝利，沒啥屁用。
因此才會需要發明其他演算法，
像是H₁ estimate、H₂ estimate、QAM estimate。
教科書誤植為bias和variance。
bias和variance是指大量實驗的結果呈現哪種分布。
此處是談僅做一次實驗的結果會是如何。
https://people.ee.ethz.ch/~rsmith/idfiles/SysID_lecture05_small.pdf
spectral smoothing
spectral coherence:
A(exp(𝑖ω₁)) and A(exp(𝑖ω₂)) are asymptotically uncorrelated.
即便ω取樣間距變小，函數曲線也不會變得比較連續平滑。
改善方法是平滑化。做大量實驗，取平均值。
A(exp(𝑖ω)) the curve is spiky
↑ ﹏〰〰〰﹏
│﹏〰 〰﹏﹏
└──────┬─┬───────→ ω
ω₁ ω₂
spectral smoothing:
輸入訊號、輸出訊號事先平滑化，其頻譜也隨之平滑化。
(1) k-fold average 時域訊號，相鄰N窗取平均。
(2) window function 時域訊號，套用窗函數。
bias–variance tradeoff:
窗越窄窗越多，bias越大variance越小。
stochastic LTI IIR system
system model
output error model Y = (A/B)X + E skip (z)
ARX model Y = (A/B)X + (1/B)E
ARMAX model Y = (A/B)X + (C/B)E
Box–Jenkins model Y = (A/B)X + (C/D)E
output error model:
e (zero-mean white noise)
╷
┌─────┐ w ↓+
x ───→│ A/B │───→⊕───→ y
└─────┘
⎧ w[n] = a₀ x[n] + a₁ x[n-1] + ... + aₚ x[n-p]
⎨ + b₁ w[n-1] + ... + b₉ w[n-q]
⎩ y[n] = w[n] + e[n]
ARX model:
e
╷
┌─────┐ ↓+ ┌─────┐
x ───→│ A │──→⊕──→│ 1/B │───→ y
└─────┘ └─────┘
y[n] = a₀ x[n] + a₁ x[n-1] + ... + aₚ x[n-p]
+ b₁ y[n-1] + ... + b₉ y[n-q] + e[n]
ARMAX model:
┌─────┐
e ───→│ C │───┐
└─────┘ │
┌─────┐ ↓+ ┌─────┐
x ───→│ A │──→⊕──→│ 1/B │───→ y
└─────┘ └─────┘
y[n] = a₀ x[n] + a₁ x[n-1] + ... + aₚ x[n-p]
+ b₁ y[n-1] + ... + b₉ y[n-q]
+ c₀ e[n] + c₁ e[n-1] + ... + cᵣ e[n-r]
Box–Jenkins model:
┌─────┐ v
e ───→│ C/D │───┐
└─────┘ │
┌─────┐ w ↓+
x ───→│ A/B │──→⊕──→ y
└─────┘
⎧ w[n] = a₀ x[n] + a₁ x[n-1] + ... + aₚ x[n-p]
⎪ + b₁ w[n-1] + ... + b₉ w[n-q]
⎨ v[n] = c₀ e[n] + c₁ e[n-1] + ... + cᵣ e[n-r]
⎪ + d₁ v[n-1] + ... + dₛ v[n-s]
⎩ y[n] = w[n] + v[n]
least squares estimation
linear least squares estimation
linear regression:
(1) linear regression
(2) linear regression with zero-mean Gaussian white error
兩者公式解恰巧相同。
公式解都是normal equation。
因此，系統模型可以省略誤差項。【尚待確認】
請見本站文件「regression」、「estimation」。
least squares estimation with ARMAX model:
θ̂ = argmin ε(θ)
ε(θ) = sum ‖y[n] - ŷ(n;θ)‖²
n=k⋯k+N
ŷ(n;θ) = a₀ x[n] + a₁ x[n-1] + ... + aₚ x[n-p]
+ b₁ y[n-1] + ... + b₉ y[n-q] + e[n]
θ = (a₀, ..., aₚ, b₁, ..., b₉)
k ≥ max(p,q)
linear regression:
⎡ y[n] ⎤ ⎡ x[n] ... x[n-p] y[n-1] ... y[n-q] ⎤ ⎡ a₀ ⎤
⎢ : ⎥ ⎢ : : : : ⎥ ⎢ : ⎥
⎢ : ⎥ = ⎢ : : : : ⎥ ⎢ aₚ ⎥
⎢ : ⎥ ⎢ ⎥ ⎢ b₁ ⎥
⎣ y[n+N] ⎦ ⎣ ⎦ ⎢ : ⎥
⎣ b₉ ⎦
y thin matrix A θ
where n ≥ max(p,q) and N+1 ≥ p+1+q
由於不知道初始條件，千萬不能採用n = 0、在矩陣右上角補零。
solution:
θ = (Aᵀ A)⁻¹ Aᵀ y
model validation:
得到正確答案之後，驗證系統模型是否正確。
檢查各個時刻的誤差值，照理來說，整體呈現常態分布。
如果不是常態分布，那麼系統模型不正確。
nonlinear least squares estimation
optimization:
也可以將迴歸問題化作最佳化問題。
演算法有gradient descent、Newton's method，隨便你用。
請見本站文件「optimization」、「multivariate optimization」。
least squares estimation with Box–Jenkins model:
θ̂ = argmin ε(θ)
ε(θ) = sum ‖y[n] - ŷ(n;θ)‖²
n=k⋯k+N
⎧ ŷ(n;θ) = w[n] + v[n]
⎪ w[n] = a₀ x[n] + a₁ x[n-1] + ... + aₚ x[n-p]
⎨ + b₁ w[n-1] + ... + b₉ w[n-q]
⎪ v[n] = c₀ e[n] + c₁ e[n-1] + ... + cᵣ e[n-r]
⎩ + d₁ v[n-1] + ... + dₛ v[n-s]
θ = (a₀, ..., aₚ, b₁, ..., b₉, c₀, ..., cᵣ, d₁, ..., dₛ)
k ≥ max(p,q,r,s)
model validation:
同前。
frequency domain method:
理論上，最佳化目標可以是時域訊號，也可以是頻域頻譜。
我不知道效果有沒有比較好。
least squares estimation with Box–Jenkins model:
θ̂ = argmin ε(θ)
ε(θ) = sum ‖Yₙ(exp(𝑖ω)) - Ŷₙ(exp(𝑖ω);θ)‖²
n=k⋯k+N
ω=0⋯M
Ŷₙ = (Aₙ/Bₙ)Xₙ + (Cₙ/Dₙ)Eₙ skip (exp(𝑖ω))
model validation:
得到正確答案之後，驗證系統模型是否正確。
檢查特定時刻誤差項Eₙ(exp(𝑖ω))的分布，照理來說，整體呈現均勻分布（白雜訊）。
如果不是均勻分布，那麼系統模型不正確。
