1. 從“黑盒”到“白盒”為什么我們需要對隨機信號建模在信號處理的世界里我們每天都要和各種信號打交道。有些信號是確定性的比如一個正弦波它的頻率、幅度、相位都是清清楚楚的我們可以用一個精確的數(shù)學(xué)公式來描述它。但更多時候尤其是在現(xiàn)實世界的工程應(yīng)用中我們面對的是隨機信號。比如一段語音、一段腦電圖、一段股票價格波動或者是你手機接收到的無線信號。這些信號在任何一個具體時刻的取值我們無法用一個固定的公式來預(yù)測它充滿了不確定性像一個“黑盒”。那么面對這樣一個“黑盒”我們是不是就束手無策了當(dāng)然不是。雖然我們無法預(yù)測其每一個具體的瞬時值但我們可以研究它的統(tǒng)計特性比如它的平均值、方差、以及不同時刻取值之間的關(guān)聯(lián)性也就是自相關(guān)函數(shù)。更進一步一個非常強大的思路是我們嘗試用一個相對簡單的、參數(shù)化的數(shù)學(xué)模型來“模仿”或“逼近”這個復(fù)雜的隨機信號。這個過程就是隨機信號的參數(shù)建模。參數(shù)建模法的核心思想是把一個復(fù)雜的隨機過程看作是由一個簡單的、確定性的系統(tǒng)我們的模型在受到一個簡單的、隨機的輸入通常是白噪聲激勵后產(chǎn)生的輸出。一旦我們找到了這個模型就等于把這個“黑盒”打開了一個口子。我們不再需要存儲或處理冗長的原始信號數(shù)據(jù)只需要記住模型的幾個關(guān)鍵參數(shù)就能在很大程度上“復(fù)現(xiàn)”或“理解”這個信號。這帶來的好處是巨大的數(shù)據(jù)壓縮、信號預(yù)測、特征提取、系統(tǒng)辨識、故障診斷……幾乎所有高級信號處理應(yīng)用都建立在有效的參數(shù)模型之上。而在眾多參數(shù)模型中自回歸模型也就是我們常說的AR 模型無疑是應(yīng)用最廣泛、理論最成熟、也最直觀的一種。它背后的邏輯非常符合人的直覺一個信號當(dāng)前時刻的值很大程度上可以由它過去若干個時刻的值線性組合來預(yù)測再加上一點無法預(yù)測的隨機“新息”。這就像預(yù)測明天的天氣我們會參考今天、昨天甚至前幾天的天氣情況再考慮一些突發(fā)的、不可控的因素。AR模型正是將這種思想數(shù)學(xué)化了。接下來我們就深入這個“白盒”內(nèi)部看看AR模型是如何構(gòu)建、如何工作以及在實際中我們?nèi)绾斡盟鼇怼榜{馭”隨機信號。2. AR模型的核心原理用過去預(yù)測現(xiàn)在理解了建模的必要性我們現(xiàn)在聚焦于AR模型本身。它的全稱是AutoRegressive Model中文譯為自回歸模型。這個名字本身就揭示了它的核心“Auto”指自身“Regressive”指回歸合起來就是用信號自身的歷史值來回歸預(yù)測當(dāng)前值。2.1 數(shù)學(xué)定義與直觀理解一個p階的AR模型其數(shù)學(xué)表達式非常簡潔x[n] -Σ_{i1}^{p} a_i * x[n-i] w[n]讓我們來拆解這個公式里的每一個符號x[n]這是我們觀測到的隨機信號在時刻n的取值也就是我們想要建模的對象。p模型的階數(shù)。它決定了我們用過去多少個時刻的數(shù)據(jù)來預(yù)測現(xiàn)在。p的選擇至關(guān)重要太小了模型太粗糙太大了又會引入過擬合和計算復(fù)雜度。a_i(i1, 2, ..., p)這就是AR模型的參數(shù)也稱為自回歸系數(shù)。它們是整個建模過程要求解的核心。a_i前面的負號是習(xí)慣寫法有時也省略但含義不變。這些系數(shù)本質(zhì)上是一組權(quán)重告訴我們過去的每一個值x[n-i]對當(dāng)前值x[n]的“影響力”有多大。w[n]這是驅(qū)動整個模型的輸入通常被假設(shè)為一個均值為0、方差為σ2的白噪聲序列。你可以把它理解為我們模型無法解釋的、完全隨機的“創(chuàng)新”或“擾動”。正是這個w[n]為整個輸出信號x[n]注入了隨機性。這個公式的直觀理解非常強當(dāng)前信號值 ≈ 過去p個信號值的加權(quán)和 一個隨機噪聲。模型的任務(wù)就是找到那一組最優(yōu)的權(quán)重{a_i}使得這個線性預(yù)測的誤差——也就是那個隨機噪聲w[n]——的功率盡可能小通常是方差最小。當(dāng)這組權(quán)重找得好時w[n]就真的像一個白噪聲不包含任何可預(yù)測的結(jié)構(gòu)信息所有可預(yù)測的部分都已經(jīng)被a_i和過去的數(shù)據(jù)x[n-i]捕捉到了。2.2 模型背后的系統(tǒng)視角一個全極點濾波器如果我們把上述公式稍微變個形從系統(tǒng)輸入輸出的角度來看會得到更深刻的見解。將公式改寫為w[n] x[n] Σ_{i1}^{p} a_i * x[n-i]這可以看作白噪聲w[n]作為輸入通過一個線性時不變系統(tǒng)后得到了輸出信號x[n]。這個系統(tǒng)的傳遞函數(shù)H(z)是什么對等式兩邊進行Z變換假設(shè)初始條件為0W(z) X(z) * (1 a_1*z^{-1} a_2*z^{-2} ... a_p*z^{-p})因此系統(tǒng)的傳遞函數(shù)為H(z) X(z) / W(z) 1 / (1 a_1*z^{-1} a_2*z^{-2} ... a_p*z^{-p})這是一個典型的全極點濾波器。它的極點完全由AR模型的參數(shù){a_i}決定。這個視角極其重要因為它將AR模型與信號的頻譜特性直接聯(lián)系了起來。為什么這一點很關(guān)鍵因為一個隨機信號的功率譜密度描述了信號功率在不同頻率上的分布。而對于上述系統(tǒng)當(dāng)輸入是白噪聲其功率譜是平坦的時輸出信號x[n]的功率譜P_x(ω)就等于輸入白噪聲的功率譜σ2乘以系統(tǒng)頻率響應(yīng)H(e^{jω})的模平方。即P_x(ω) σ2 / |1 Σ_{i1}^{p} a_i e^{-jωi}|2這意味著一旦我們通過建模估計出了AR參數(shù){a_i}和噪聲方差σ2我們就直接得到了該隨機信號的一個功率譜估計。這種譜估計方法被稱為AR譜估計或最大熵譜估計。與傳統(tǒng)的基于傅里葉變換的周期圖法相比AR譜估計在數(shù)據(jù)記錄短、分辨率要求高的場合如雷達、聲納、生物醫(yī)學(xué)信號處理有著顯著優(yōu)勢因為它隱含著對數(shù)據(jù)范圍外的外推假設(shè)能提供更高的頻率分辨率。3. 如何為你的信號“量身定制”AR模型參數(shù)估計實戰(zhàn)理論很優(yōu)美但落到實操上我們面對一段具體的信號數(shù)據(jù)x[0], x[1], ..., x[N-1]如何找到那組最優(yōu)的AR參數(shù){a_i}和噪聲方差σ2呢這就是AR模型參數(shù)估計要解決的問題。主要有三種經(jīng)典方法它們基于不同的優(yōu)化準(zhǔn)則但核心思想相通。3.1 尤爾-沃克方程法從自相關(guān)函數(shù)出發(fā)這是最經(jīng)典、最直接的方法它建立在“使前向預(yù)測誤差功率最小”的準(zhǔn)則上。推導(dǎo)過程涉及一些線性代數(shù)但其最終形式非常規(guī)整——尤爾-沃克方程[ r[0] r[1] ... r[p-1] ] [ a_1 ] [ -r[1] ] [ r[1] r[0] ... r[p-2] ] [ a_2 ] [ -r[2] ] [ ... ... ... ... ] * [ ... ] [ ... ] [ r[p-1] r[p-2] ... r[0] ] [ a_p ] [ -r[p] ]其中r[m] E{ x[n] * x[nm] }是信號的理論自相關(guān)函數(shù)。在實際中我們用樣本數(shù)據(jù)估計自相關(guān)函數(shù)例如\hat{r}[m] (1/N) * Σ_{n0}^{N-1-m} x[n] * x[nm], 對于 m 0并且r[-m] r[m]?？梢钥吹椒匠探M的系數(shù)矩陣是一個托普利茨矩陣沿對角線元素相同并且是正定的。這使得我們可以用高效的萊文森-德賓遞推算法來求解該算法復(fù)雜度僅為O(p2)避免了直接求逆矩陣的O(p3)復(fù)雜度。求解步驟計算自相關(guān)根據(jù)觀測數(shù)據(jù)估計出自相關(guān)序列\(zhòng)hat{r}[0], \hat{r}[1], ..., \hat{r}[p]。萊文森-德賓遞推初始化a_1(1) -r[1]/r[0],σ_12 (1 - |a_1(1)|2) * r[0]對于 k2 到 p:κ_k - ( r[k] Σ_{i1}^{k-1} a_i(k-1) * r[k-i] ) / σ_{k-1}2a_k(k) κ_ka_i(k) a_i(k-1) κ_k * a_{k-i}(k-1), for i1,..., k-1σ_k2 (1 - |κ_k|2) * σ_{k-1}2最終得到的a_i(p)(i1..p) 就是AR(p)模型的參數(shù)σ_p2就是白噪聲方差估計。注意尤爾-沃克法在數(shù)據(jù)量較大時表現(xiàn)穩(wěn)健但它有一個隱含的假設(shè)在計算自相關(guān)時對觀測窗口外的數(shù)據(jù)做了補零假設(shè)。這可能導(dǎo)致在短數(shù)據(jù)情況下譜估計出現(xiàn)偏差。3.2 協(xié)方差法與修正協(xié)方差法更精確的數(shù)據(jù)匹配為了克服尤爾-沃克法在數(shù)據(jù)邊界處的假設(shè)問題協(xié)方差法直接基于原始數(shù)據(jù)最小化前向預(yù)測誤差的平方和。其正則方程中的矩陣元素計算如下c_{ij} (1/(N-p)) * Σ_{np}^{N-1} x[n-i] * x[n-j], 其中 i, j 0, 1, ..., p (這里定義a_01)。這個矩陣C不再是托普利茨矩陣但仍然是正定的。求解這個方程通常使用喬里斯基分解或奇異值分解等線性代數(shù)方法無法使用萊文森-德賓遞推。協(xié)方差法通常能給出比尤爾-沃克法更高的頻率分辨率。修正協(xié)方差法則同時最小化前向預(yù)測誤差和后向預(yù)測誤差的平方和進一步提升了數(shù)據(jù)的利用率和平穩(wěn)性特別適用于短數(shù)據(jù)序列其譜估計特性往往更優(yōu)。3.3 伯格法在保證模型穩(wěn)定的前提下遞推伯格方法非常巧妙它通過遞推的方式在每一步都同時滿足前向和后向預(yù)測誤差最小化并且強制保證最終得到的AR模型是穩(wěn)定的即其對應(yīng)的系統(tǒng)極點都在單位圓內(nèi)。這是伯格法最大的優(yōu)點。伯格遞推的核心是計算反射系數(shù)或稱偏相關(guān)系數(shù)κ_k。其步驟與萊文森-德賓類似但計算κ_k的公式不同它直接基于前向和后向預(yù)測誤差能量來計算。由于保證了穩(wěn)定性伯格法在實際中應(yīng)用非常廣泛許多軟件工具如MATLAB的arburg函數(shù)默認(rèn)采用的就是伯格算法。方法選擇經(jīng)驗談追求穩(wěn)健和快速數(shù)據(jù)較長且信噪比較高時尤爾-沃克法萊文森-德賓遞推是首選。追求高分辨率數(shù)據(jù)較短時協(xié)方差法或修正協(xié)方差法通常能給出更尖銳的譜峰。保證模型穩(wěn)定當(dāng)模型階數(shù)p較高或者你需要確保生成的合成信號不發(fā)散時伯格法是最安全的選擇。我在處理語音信號合成時就曾因為使用其他方法在高階時得到不穩(wěn)定模型導(dǎo)致合成語音爆炸幅度無限增大改用伯格法后問題迎刃而解。4. 模型階數(shù)p一個至關(guān)重要的超參數(shù)無論采用哪種估計方法你都必須事先指定一個階數(shù)p。p選得太小模型過于簡單無法捕捉信號中復(fù)雜的相關(guān)性這稱為“欠擬合”會導(dǎo)致譜估計平滑、細節(jié)丟失。p選得太大模型會開始擬合信號中的隨機噪聲成分這稱為“過擬合”會導(dǎo)致譜估計出現(xiàn)虛假的峰值模型參數(shù)方差增大。那么如何確定這個“恰到好處”的階數(shù)p呢沒有絕對正確的答案但有以下幾種實用的準(zhǔn)則和方法4.1 信息論準(zhǔn)則AIC與MDL這類準(zhǔn)則在擬合優(yōu)度和模型復(fù)雜度之間進行折衷。它們會計算不同階數(shù)p下的一個準(zhǔn)則函數(shù)值選擇使該函數(shù)值最小的p作為最佳階數(shù)。赤池信息量準(zhǔn)則AIC(p) N * ln(σ_p2) 2p其中σ_p2是p階模型下的白噪聲方差估計。AIC傾向于選擇稍高階的模型。最小描述長度準(zhǔn)則MDL(p) N * ln(σ_p2) p * ln(N)MDL的懲罰項比AIC更重因此傾向于選擇比AIC更低的階數(shù)在樣本量N較大時MDL準(zhǔn)則具有一致性即當(dāng)N趨于無窮時能選出真實階數(shù)。在實際操作中我會計算從1到一個預(yù)設(shè)最大階數(shù)比如N/3或N/2范圍內(nèi)所有p對應(yīng)的AIC和MDL值然后畫出曲線尋找明顯的拐點或最小值點。4.2 最終預(yù)測誤差準(zhǔn)則FPE準(zhǔn)則FPE(p) σ_p2 * (Np1)/(N-p-1)FPE準(zhǔn)則的目標(biāo)是最小化一步預(yù)測的均方誤差也是一個常用的參考。4.3 觀察預(yù)測誤差方差或反射系數(shù)的變化這是一個更直觀的方法隨著階數(shù)p增加白噪聲方差估計σ_p2通常會單調(diào)遞減。當(dāng)p達到或超過真實階數(shù)后σ_p2的下降會變得非常緩慢出現(xiàn)一個“肘部”。同樣在伯格算法中反射系數(shù)|κ_k|的絕對值會隨著k增大而減小。當(dāng)k超過真實階數(shù)后|κ_k|通常會趨近于0。你可以將|κ_k|首次低于某個閾值比如0.05時的k作為階數(shù)估計。我的實戰(zhàn)經(jīng)驗不要迷信單一準(zhǔn)則。最好的做法是多方法交叉驗證。例如同時觀察AIC、MDL的曲線并結(jié)合σ_p2下降的“肘部”位置。然后用選出的幾個候選p值分別進行AR譜估計觀察其功率譜圖。一個“好”的譜圖應(yīng)該具有清晰的物理可解釋的譜峰而沒有太多雜亂無章的小毛刺。例如在分析一個包含50Hz和120Hz工頻干擾的腦電信號時如果AR譜在50Hz和120Hz處出現(xiàn)了尖銳且合理的峰值而在其他頻率很平坦那這個階數(shù)p可能就是合適的。如果譜圖上出現(xiàn)了很多密集的、無法解釋的小峰那很可能就是過擬合了。5. AR模型的力量從譜估計到預(yù)測與合成當(dāng)我們成功估計出AR模型的參數(shù){a_i}和σ2后這個模型就成為了我們理解和操作該隨機信號的有力工具。它的應(yīng)用遠不止于“理解”更在于“創(chuàng)造”和“預(yù)測”。5.1 高分辨率功率譜估計如前所述將估計出的參數(shù)代入公式P_x(ω) σ2 / |1 Σ_{i1}^{p} a_i e^{-jωi}|2即可得到信號的AR譜估計。與傳統(tǒng)的周期圖法相比AR譜估計尤其適用于短數(shù)據(jù)記錄傳統(tǒng)方法分辨率受限于數(shù)據(jù)長度1/T而AR譜估計可以突破這個限制。銳峰頻譜對于由多個正弦波疊加而成的信號AR譜能呈現(xiàn)出非常尖銳的譜線便于頻率檢測。平滑背景上的譜峰能有效區(qū)分寬頻帶背景噪聲上的窄帶信號。在雷達目標(biāo)速度估計、語音共振峰分析、腦電節(jié)律提取等領(lǐng)域AR譜估計是標(biāo)準(zhǔn)工具之一。5.2 線性預(yù)測與信號濾波AR模型本身就是一個線性預(yù)測器。給定過去p個樣本x[n-1], ..., x[n-p]我們對當(dāng)前值的最優(yōu)線性預(yù)測在最小均方誤差意義下就是\hat{x}[n] -Σ_{i1}^{p} a_i * x[n-i]預(yù)測誤差e[n] x[n] - \hat{x}[n]理論上應(yīng)該接近于白噪聲w[n]。這個性質(zhì)被廣泛應(yīng)用于語音編碼如線性預(yù)測編碼傳輸預(yù)測誤差殘差和模型參數(shù)而非原始語音樣本實現(xiàn)高效壓縮。信號去噪如果信號符合AR模型而噪聲是加性的可以通過預(yù)測和相減來增強信號。異常檢測在平穩(wěn)運行的系統(tǒng)如旋轉(zhuǎn)機械中其振動信號可以用AR模型描述。一旦模型建立實時計算預(yù)測誤差。當(dāng)系統(tǒng)出現(xiàn)故障時信號特性改變預(yù)測誤差e[n]的功率會突然增大從而觸發(fā)報警。5.3 隨機信號合成這是AR模型一個非?！翱帷钡膽?yīng)用。既然我們認(rèn)為信號是由白噪聲w[n]通過一個傳遞函數(shù)為H(z)的系統(tǒng)產(chǎn)生的那么反過來我們也可以用計算機生成一段白噪聲序列然后讓它通過我們估計出的AR模型系統(tǒng)H(z)來合成一段與原始信號統(tǒng)計特性相似的新信號。具體步驟估計原始信號x[n]的AR(p)模型參數(shù){a_i}和噪聲方差σ2。生成一個方差為σ2、均值為0的白噪聲序列w_synth[n]。用差分方程進行濾波x_synth[n] -Σ_{i1}^{p} a_i * x_synth[n-i] w_synth[n]。忽略前若干點的瞬態(tài)響應(yīng)得到平穩(wěn)的合成信號x_synth[n]。合成的信號x_synth[n]與原始信號x[n]具有相同的自相關(guān)函數(shù)和功率譜密度二階統(tǒng)計特性相同。這在需要大量具有特定統(tǒng)計特性的仿真數(shù)據(jù)時非常有用例如通信系統(tǒng)仿真、金融風(fēng)險蒙特卡洛模擬、以及音頻合成中的背景噪聲生成。一個我踩過的坑在合成信號時務(wù)必確保模型的穩(wěn)定性。不穩(wěn)定的AR模型其系統(tǒng)極點有在單位圓外的會導(dǎo)致合成信號幅度指數(shù)增長迅速溢出。這就是為什么在合成應(yīng)用中我強烈推薦使用伯格法來估計參數(shù)因為它能保證穩(wěn)定性。如果用了其他方法在合成前一定要檢查系統(tǒng)極點即多項式1 a_1 z^{-1} ... a_p z^{-p} 0的根是否全部在單位圓內(nèi)。6. 超越基礎(chǔ)AR模型的局限與擴展AR模型雖然強大但并非萬能鑰匙。理解它的局限性才能知道何時該用它何時該尋求其他工具。6.1 主要局限性對信號特性的假設(shè)AR模型最適合建模全極點譜的信號。也就是說它的功率譜密度可以通過全極點濾波器很好地匹配。對于在頻譜上有深谷即“零點”的信號AR模型需要很高的階數(shù)才能近似效率低下。對噪聲敏感估計過程特別是基于自相關(guān)的方法對觀測噪聲比較敏感。如果信號被加性白噪聲污染估計出的AR參數(shù)和譜峰會有所偏差。階數(shù)選擇的主觀性如前所述最佳階數(shù)p的選擇沒有黃金標(biāo)準(zhǔn)需要經(jīng)驗和多種準(zhǔn)則輔助判斷。6.2 模型家族的擴展為了克服AR模型的局限更一般的參數(shù)模型被提出滑動平均模型x[n] Σ_{i0}^{q} b_i * w[n-i] 其系統(tǒng)函數(shù)只有零點適合表征具有凹槽頻譜的信號。自回歸滑動平均模型x[n] -Σ_{i1}^{p} a_i * x[n-i] Σ_{i0}^{q} b_i * w[n-i] 這是AR模型和MA模型的結(jié)合系統(tǒng)函數(shù)既有極點也有零點理論上可以用更低的階數(shù)(p, q)擬合更廣泛的信號。但ARMA模型的參數(shù)估計比AR模型復(fù)雜得多。自回歸積分滑動平均模型在ARMA基礎(chǔ)上引入了差分運算專門用于處理非平穩(wěn)時間序列如具有趨勢或季節(jié)性的經(jīng)濟數(shù)據(jù)。在實際工作中AR模型因其概念簡單、計算高效、算法成熟仍然是首選的“第一模型”。當(dāng)AR模型效果不佳時我們才會考慮更復(fù)雜的ARMA等模型。通常一個實用的策略是先嘗試用AR模型如果發(fā)現(xiàn)需要的階數(shù)p異常高或者殘差檢驗顯示預(yù)測誤差還不是白噪聲再考慮引入MA部分。從我多年的工程實踐來看AR模型及其譜估計是每個信號處理工程師工具箱里的必備品。它的價值在于將看似不可捉摸的隨機信號轉(zhuǎn)化為幾個具有物理或數(shù)學(xué)意義的參數(shù)從而打開了分析、預(yù)測、合成信號的大門。掌握它不僅僅是學(xué)會幾個算法函數(shù)更是建立起一種“建?！钡乃季S方式——用簡單的數(shù)學(xué)結(jié)構(gòu)去理解和駕馭復(fù)雜的世界這正是工程技術(shù)的魅力所在。當(dāng)你下次再面對一段嘈雜的、看似無規(guī)律的信號時不妨試著用AR模型去“問詢”它你很可能會得到一幅清晰得多的頻率“肖像”。