自適應定階的快速Burg 算法設計與FPGA 實現*

郭鳴晗 ，陳立平 ，張 浩 ，趙 坤 ，柏 偉

(1.中國科學院大學，北京 100049；2.中國科學院微電子研究所，北京 100029)

0 引言

現代功率譜估計的AR 模型法使用有限長的數據序列來估計假設模型的參數，再將參數帶入功率譜密度模型中，可獲得較好的功率譜估計結果[1-3]。

Burg 算法是一種常見的AR 模型求功率譜的方法，其主導思想是利用前后向預測誤差功率之和最小的方法來計算反射系數k，然后帶入Levinson 遞推，求解AR模型參數[4-6]。此方法在處理短數據時具有較高的頻率分辨率[7-8]，但求解反射系數計算量較大。為了改進這一問題，Vos 提出一種快速Burg 算法[9]，通過一系列矩陣變換降低了反射系數求解時的計算量，但是不能確定AR模型的階數，并且串行算法的執行耗時較長。針對上述問題，本文將快速Burg 算法與FPE 準則[10]結合，對短序列的功率譜估計實現自適應定階的功能，達到較高頻率分辨率，并使用Verilog 硬件描述語言設計電路，達到硬件加速功能[11]。電路結構在二級流水線的基礎上[12]，結合自適應定階方案，提出一種新的流水線結構，并設置狀態機靈活控制。本文對計算單元進行并行化處理加速計算。考慮速度與面積的折中，針對算法特點設計內存讀寫方案，減少數據存儲長度，從而減小了存儲單元的面積。

1 自適應定階的快速Burg 算法

快速Burg 算法過一系列矩陣變換將反射系數計算過程化簡[9]，避免了Burg 算法[4]中計算反射系數前需計算出前、后向預測誤差ef、eb而帶來的較大計算量。在本算法中，通過自適應定階的方式，對長度為N 的輸入序列xn的遞歸計算步驟如下：

(1)計算自相關ci。(本實驗中自適應定階在10 階以內，在此只計算10 個自相關值)。

(2)初始化。初始為第一階，此時迭代次數i=0；用于計算反射系數的列向量g0=[2c0-|x0|2-|xN-1|2；2c1]T；用于計算列向量gi的列向量r1=c1；預測系數a0=1；反射系數。

(3)計算初始預測誤差功率Pmin(0)，預測誤差FPE 值和初始閾值TH：

(4)判斷迭代條件(TH＜FPE(i))，即若滿足閾值TH 小于預測誤差FPE 值，則將當前閾值TH 更新為預測誤差值TH=FPE(i)，之后i+1，繼續執行步驟(5)；若不滿足，則計算階數，輸出預測系數ai。

(5)更新列向量ri+1：

(6)更新用于計算列向量gi的矩陣：

(7)更新列向量gi：

其中，J 是i+1 行的方陣：

(8)計算初始反射系數：

(9)計算最小預測誤差功率Pmin(i)與最終預測誤差準則中的FPE(i)值：

(10)更新AR 模型參數ai，跳轉到步驟(3)。

2 硬件電路設計與實現

2.1 并行二級流水線結構

本算法完成硬件加速的關鍵思路是使用流水線的結構。自相關ci的計算只和輸入序列xn有關，與AR 參數的計算無關，而AR 參數的計算需要利用自相關的值，故電路結構將劃分為二級流水線，第一級為自相關ci的計算，第二級為AR 參數和閾值的計算。為了提高計算速度，本設計將采用并行的二級流水線結構，時序關系如圖1 所示。

圖1 并行二級流水線時序圖

本電路結構的頂層電路圖如圖2 所示，計算單元由兩個并行連接的自相關計算單元和AR 參數、閾值計算單元構成，控制單元為狀態機來控制迭代次數以及計算單元。存儲單元包括兩個存放xn的RAM，寬度為16，深度為序列長度N(在本設計中N=128)，一個用于自相關計算，一個用于AR 參數計算；第二部分用來存儲AR 參數、閾值計算單元產生的中間結果。

圖2 快速Burg 算法頂層電路圖

2.2 自適應定階與狀態機控制

為了達到自適應定階目的和提高流水線使用效率，采用狀態機來控制計算單元，根據閾值FPE 是否達到最小值來確定是否停止運算。在本算法中，AR 參數的階數自適應得確定在2～10 階內。

本狀態機的輸入信號為外部輸入的算法使能信號和計算單元的結束狀態標識信號，輸出為計算單元的使能控制信號。

根據并行二級流水線結構和自適應定階算法，控制單元狀態機如圖3 所示，共分為6 個狀態。

圖3 快速Burg 算法控制單元狀態機

算法啟動后，初始階數為1，迭代次數i=0，進入S1狀態進行初始自相關計算；之后進入S2 啟動后續自相關計算及AR 參數和閾值計算。計算完成后若FPE 達到最小值，跳至S5 結束運算，最終定階為i+1 階；否則跳轉至S4，啟動下一階AR 參數模型和閾值計算。計算完成后若FPE 達到最小值，則跳轉至S5 結束運算，最終定階為i+2 階；否則跳轉至S4，等待當前的自相關計算結束后再次進入S2 進行下一輪的迭代，直至FPE 達到最小值，迭代停止，RAM 中的AR 參數即為所求。

2.3 計算單元的設計

2.3.1自相關的計算單元

自相關的計算使用流水線結構。產生地址、讀取數據、乘累加運算，三拍完成一次乘累加輸出，兩個自相關計算單元分別進行N-i 次和N-i-1 次乘累加后，輸出自相關值ci、ci+1給下一級計算單元，當前階數下消耗的總計算時長為N-i+3 個時鐘周期。

2.3.2 AR 參數和閾值計算單元

本計算單元電路包括4 個模塊：預備參數計算模塊、反射系數k 值和AR 參數a 計算模塊、預測誤差功率Pmin計算閾值計算模塊，以及判斷閾值是否達到最小值的定階模塊。本單元頂層結構如圖4 所示。

圖4 AR 參數和閾值計算單元頂層電路圖

快速Burg 算法的關鍵即為計算預備變量，包括r、ΔR和g 3 個變量。

r 列向量的計算電路為三級流水線結構，其電路如圖5 所示。第一拍鎖存當前階數下輸入序列的第一個乘數，并將讀入的數據寫入Xf_buffer 和Xb_buffer；第二拍并行完成兩次乘加運算；第三拍做選通輸出c1或乘加值，并倒序輸出r 列向量的行值。當前階數下消耗的總計算時長為i+3 個時鐘周期。

圖5 列向量r 的運算電路

中間變量矩陣ΔR 和列向量g 的計算電路采用并行化處理，矩陣ΔR 和列向量g 的四級流水線電路如圖6所示。第一拍從兩個buffer 中各依次讀出一個x 值送入乘法器；第二拍做乘加運算，得到矩陣ΔR 的一行值；第三拍ΔR 的行值分別與上一階求得的AR 參數依次相乘并求和；第四拍選通輸出列向量g 的行值。當前階數下消耗的總計算時長為i+5 個時鐘周期。

圖6 矩陣ΔR 和列向量g 的運算電路

2.4 存儲單元設計

由于參數計算過程中部分變量需要在每階下遞歸運算，因此需為其分配內存空間。為了起到加速的作用，使用雙口RAM 對中間結果進行存儲，提高讀寫數據的速度。

3 結果分析

3.1 功率譜估計的分辨率

在本節中，當輸入數據采樣率為512 Hz，輸入數據點數為512 時，通過ModelSim 實驗平臺對基于定點數[13-14]的算法進行仿真，并將FFT 算法和本算法的頻率分辨率進行比較。仿真實驗中，AR 參數計算結果a 和最小預測誤差Pmin仿真結果如圖7 所示。

圖7 ModelSim 仿真結果

為了更直觀地觀察頻譜結果，將硬件電路的結果傳入MATLAB 中，經過處理打印功率譜圖，并將硬件電路計算出的定點數模型的功率譜結果與MATLAB 直接計算的浮點數功率譜結果進行比較，如圖8 所示，譜峰一致說明量化良好。

圖8 定浮點算法功率譜結果比較

實驗設置輸入數據為15 Hz 和20 Hz 正弦波的疊加信號，采樣率為512 Hz。

當輸入數據長度為512 時，FFT 算法和自適應定階的快速Burg 算法的譜分析結果分別如圖9(a)和9(b)所示，兩種方式都可分辨出15 Hz 和20 Hz 兩個譜峰，快速Burg 算法的譜分析結果更為平滑，階數為7 階；當輸入數據長度為128 時，上述兩種算法的譜分析結果分別如圖9(c)和9(d)所示，FFT 算法無法得到兩個獨立譜峰，而自適應定階的快速Burg 算法依然可以產生15 Hz 和20 Hz 兩個獨立譜峰，自適應定階為5 階。

圖9 FFT 算法與自適應定階快速Burg 算法譜分析結果

綜上，自適應定階的快速Burg 算法在短序列譜分析中具有更高的頻率分辨率，在實際應用中，對實時性要求較高時，具有較為重要的意義。

3.2 計算時間

快速Burg 算法與Burg 算法相比，從硬件加速設計上分析，Burg 算法的前后向預測誤差ef、eb，反射系數k，AR 參數a，最小預測誤差功率Pmin的計算4 組變量完全是串行關系，無法并行化[15]。Burg 算法計算單元在每階下所消耗的時間如表1 所示。

表1 Burg 算法計算單元消耗時鐘數

對于快速Burg 算法，自相關的計算與其他變量無關，只與輸入序列和當前階數相關，它作為耗時最長的運算單元，可以提前一個流水線周期計算。在設計中采用并行的二級流水線結構，并對流水線進行靈活的啟動終止控制，使該算法與加速前相比可以將整體計算速度提高約50%。其各個計算單元所消耗的時鐘數如表2 所示。

表2 快速Burg 算法計算單元消耗時鐘數

最后，本文對兩種算法在不同階數下的硬件實現時計算單元所消耗的時鐘數進行比較。以輸入數據長度N=128 為例，如表3 所示，快速Burg 算法的硬件加速實現方案確實可以起到更為有效的硬件加速作用，在AR 階數在10 階以內時，與Burg算法相比，可將計算時間降低75%。

表3 各算法在不同最終階數時計算單元消耗時鐘數

3.3 內存空間分析

快速Burg 算法的硬件實現方案主要從反射系數k的計算中減少RAM 的使用。其中自相關的計算由于流水線的設計無需存儲器進行存儲，列向量r 和g 需要進行存儲，其存儲深度等于階數，在本設計中階數在2～10之間，故深度為10，寬度分別為4 B 和6 B。而Burg 算法在計算反射系數k 時，需要為前后向預測誤差ef、eb這兩個參量分配存儲空間，所需RAM 的寬度為2 B，深度和輸入數據長度N 有關。因此，快速Burg 算法的硬件實現明顯更加節省面積。

3.4 硬件實現資源使用情況

采用如圖4 所示的自適應定階快速Burg 算法的實現方案在Xilinx Virtex7 實驗平臺上做了仿真測試，綜合結果如表4 所示。

表4 資源使用情況

4 結論

本文以Xilinx 公司的Virtex7 開發板作為硬件平臺，完成了自適應定階的快速Burg 算法設計與實現。提出靈活控制的并行二級流水線結構和并行化計算單元結構，達到了硬件加速的目標。同Burg 算法相比，其對短序列的處理在達到較高分辨率的同時，計算時間可降低75%以上，同時減小了內存空間，并自適應定階。實驗發現，AR 模型階數是否達到最優決定了功率譜質量的好壞，但是此最優階數的確定也與輸入序列長度、采樣率等相關，這一關系有待進一步研究。