引言
訊號學家自創一個詞彙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 system
開迴路->輸入的加權總和->輸入只取幾項->有限脈衝響應
2. LTI IIR system = ARMA system
閉迴路->輸入與輸出的加權總和->輸出強行展開->輸入取所有項->無限脈衝響應
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
│ └ Y = AX + E
└ stochastic LTI IIR system
├ Y = (A/B)X + E
├ Y = (A/B)X + (1/B)E
├ Y = (A/B)X + (C/B)E
└ Y = (A/B)X + (C/D)E skip (z)
有些訊號學書籍替這些基礎模型命名。
名稱直接引用統計學術語，意義卻完全對不上，牛頭不對馬嘴。
下述名稱取自專著《System Identification: Theory for the User》。
很不幸的，MATLAB也採用這些名稱。
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 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
system model:
e (zero-mean white noise)
╷
┌─────┐ ↓
x ───→│ A │──→⊕──→ y
└─────┘
y[n] = sum { aₖ x[n-k] + e[n] }
k=0⋯∞
Y(z) = A(z)X(z) + E(z)
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 system，
可以直接估計系統參數，
也可以間接估計系統頻譜。
LTI FIR system = MA system。
卷積核恰是系統參數。
系統頻譜做逆向傅立葉轉換得到卷積核。
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
system model:
e (zero-mean white noise)
╷
┌─────┐ ↓+
x ───→│ A │──→⊕──→ y
└─────┘
y[n] = sum { aₖ x[n-k] + e[n] }
k=0⋯∞
Y(z) = A(z)X(z) + E(z)
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
system 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
(1) Y = (A/B)X + E
(2) Y = (A/B)X + (1/B)E
(3) Y = (A/B)X + (C/B)E
(4) Y = (A/B)X + (C/D)E skip (z)
(1)
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]
(2)
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]
(3)
┌─────┐
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]
(4)
┌─────┐ 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 Y = (A/B)X + (C/D)E:
θ̂ = 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 Y = (A/B)X + (C/D)E:
θ̂ = argmin ε(θ)
ε(θ) = sum ‖Yₙ(exp(𝑖ω)) - Ŷₙ(exp(𝑖ω);θ)‖²
n=k⋯k+N
ω=0⋯M
Ŷₙ = (Aₙ/Bₙ)Xₙ + (Cₙ/Dₙ)Eₙ skip (exp(𝑖ω))
model validation:
得到正確答案之後，驗證系統模型是否正確。
檢查特定時刻誤差項Eₙ(exp(𝑖ω))的分布，照理來說，整體呈現均勻分布（白雜訊）。
如果不是均勻分布，那麼系統模型不正確。
