正交匹配追蹤算法的優(yōu)化設計與FPGA實現(xiàn)*

莫禹鈞 ，柏正堯 ，黃 振 ，董 亮 ，2，周 燕

(1.云南大學 信息學院，云南 昆明 650091；2.中國科學院云南天文臺，云南 昆明 650011)

0 引言

2006年，CANDES D E等人提出了壓縮感知（Compressed Sensing，CS）理論[1]，CS 理論利 用與表達基 不相干的觀測矩陣，以低于奈奎斯特的采樣速率非自適應地采樣可稀疏表示的信號，得到低維的離散信息矢量，該信息矢量包含了原始信號的全部信息，然后通過非線性重建算法完美地重建信號。

壓縮感知理論主要包含了三大核心部分：信號的稀疏表示、測量矩陣的構造和信號重構算法的設計。在壓縮感知理論的三個核心問題中，如何設計并用硬件實現(xiàn)根據(jù)離散信息樣點準確重構原始信號的行之有效的算法是該理論中較為重要的一環(huán)。目前，壓縮感知信號重構算法主要分為兩類：基于凸松弛的優(yōu)化算法，如基追蹤（Basis Pursuit，BP）算法；基于貪婪迭代的匹配追蹤算法，如OMP算法[2]。這兩類算法各有優(yōu)缺點：凸松弛算法具有很好的魯棒性，然而由于需要將求解問題轉化為線性規(guī)劃問題，計算量大，信號重構效率低；貪婪算法雖然不具有強保證性，但實現(xiàn)簡單，重構效率高，在工程應用中得到廣泛使用[3]。

首次對壓縮感知恢復算法進行VLSI設計是在參考文獻[4]中，而之后，有文獻進行優(yōu)化設計。參考文獻[5]根據(jù)OMP算法必須按照特定順序執(zhí)行這一特征，采用資源復用技術，提高了資源利用率。參考文獻[6]設計了一個快速求逆平方根算法，在矩陣分解部分采用QR算法。參考文獻[7]對OMP算法進行優(yōu)化，減少了計算延時。參考文獻[8]同時進行了OMP算法和AMP算法的VLSI設計。本文先對OMP算法進行理論分析，然后對OMP算法進行改進，通過增加一個閾值來減少乘法運算次數(shù)，使運算速度更快。在矩陣分解部分采用ACD方法避免開方運算，同時在硬件實現(xiàn)上也進行了相應的優(yōu)化。仿真結果驗證了設計的可行性。

1 OMP算法

1.1 基本OMP算法

在壓縮感知中，原始信號x的稀疏度為k，觀測矢量y是所采集的數(shù)據(jù)，y可通過測量矩陣Φ與x相乘而得。本設計的目的是在已知y和Φ的前提下恢復出x。OMP算法主要分為兩部分，即尋找稀疏矢量中非零元素的位置和計算非零元素的值。

在OMP算法中殘差r是一個很關鍵的參數(shù)，殘差是通過當前選取的列向量和原始信號的線性組合不能對壓縮測量值進行表示的部分。

1.2 改進OMP算法

令原始信號x的稀疏度為k，測量矩陣Φ大小為M×N，那么y為M維的離散信息矢量。本文提出一種新的方法，即加閾值法，通過添加一個閾值來減少乘法運算次數(shù)，閾值定為內積和的平均值的α倍，內積小于閾值的那些列在下一次迭代中不再求內積。每次迭代計算后都要對閾值進行更新。信號估計的均方誤差隨著α的增大而增大，當α為0時均方誤差最小。改進的OMP算法步驟如下：

循環(huán)步驟：

(1)尋找測量矩陣Φ中與殘差最匹配的原子，即與殘差的內積最大的列：

(2)更新索引集和列的集：

(3)更新殘差：

(4)計算內積的平均值：

(5)計算閾值，并去掉內積小于閾值的列：

(6)t=t+1，如果 t＜k，則返回第(1)步。

(7)解如下最小二乘問題，得到新的估計值：

2 計算步驟

本文利用硬件實現(xiàn)重構長度N=256、稀疏度k=8的原始信號，觀測矢量長度M=64。

改進后的OMP算法可分為4個模塊。第1個模塊對應重建過程的第(1)和第(2)步，也就是在剩余列的集Φ^中尋找對殘差貢獻最大的列為最匹配原子。

第2個模塊對應重建過程的第(3)步，即計算新殘差，為下次迭代做準備。

第 3個模塊對應重建過程的第(4)和第(5)步，即計算新的閾值并除去剩余列的集Φ^中和殘差求內積小于閾值的列。求閾值前要先求內積的平均值。第t次迭代的內積平均值可用以下公式計算：

通過以上公式可以看出，只要將Csumt向右移8位，就可以得到Cavgt。

為解決對Φ的列的定位問題，用一個256位的標志位來追蹤Φ的列，標志位的第i位對應Φ的第i列。在第i列和殘差求內積后，下一個時鐘和殘差求內積的就是下一個標志位為非零所對應的列，跳過標志位為零對應的列。開始前先把標志位的每一位全部初始化成1，在每一次迭代之后對標志位進行更新。

第4個模塊對應重構過程的第(7)步，求解非零元素的值，即解決最小二乘問題。對于這類運算一般用Moore-Penrose偽逆的方法求解：

矩陣C的求逆可以通過如下求得：

求出C的逆矩陣后，就可以求得原始信號的估計：

由于OMP算法的迭代性質，4個模塊是不能并行執(zhí)行的，只能每個模塊依次執(zhí)行。

3 硬件設計

硬件電路主要由以上4個模塊組成，分為兩個部分。整體硬件電路如圖1所示。

圖1 整體電路結構圖

首先用觀測矢量y對殘差r進行初始化。y用寄存器組存儲，而觀測矩陣Φ用多個RAM存儲，這樣就能在一個時鐘內讀出y的所有值和Φ的一列值。數(shù)據(jù)用24位定點數(shù)表示，10位整數(shù)，14位小數(shù)。設計 64個 24位乘法器并行工作來求內積，然后找到內積最大值來更新。 矩陣的大小變化從 N×1～N×8。

每次迭代后會把Φ中和殘差內積小于閾值的列過濾掉，根據(jù)式(9)、(10)和(11)，剩余列的集中的每一列和殘差的內積都送到累加器進行求和，然后通過求內積平均值求得閾值。閾值參數(shù)α設置為一個常數(shù)。

256位標志位作為Φ的地址尋址，標志位每一位對應Φ每一列，初始化為所有位為1。每次迭代后對標志位進行更新，把Φ中和殘差內積小于閾值的列所對應的標志位賦為零，否則保持為1。然后在下一次迭代時跳過標志位為零所對應的Φ的列，也就是直接用下一個非零標志位所對應的列與殘差進行求內積。通過把標志位的前32位送到一個32位前導零計算器可以找出下一個非零位。

在尋找非零元素位置的部分迭代8次后，就開始計算非零元素的值。首先要計算矩陣C=，可通過以下等式計算：

此處復用之前的64個乘法器。C是一個對稱矩陣，所以只需要計算C的對角線上8個元素和對角線下半部（或上半部）的28個元素。

然后要對C進行交替的柯列斯基分解，矩陣分解要求出下三角矩陣L和對角矩陣D。從式(13)和(14)可以看出，L和D是相互依存的，必須以特定的順序計算。本設計中稀疏度k=8，L和D可以按照圖2箭頭所指順序計算。設計7個乘法器并行計算D中的元素，那么每計算一個元素需要一個時鐘周期。計算D-1時采用參考文獻[9]的方法進行除法運算。由于L的同一列的各個元素并不是相互依存的，所以求L的每一列值都設計為并行計算各個元素，那么每一列的計算只需要一個時鐘周期。

圖2 矩陣L和矩陣D的計算順序

矩陣L的求逆需要迭代進行，如式(18)：

由于L的逆矩陣的各列的各個元素是相互依存的，所以列和列可以并行運算，每一列要按照特定的順序運算，那么計算L-1需要7個時鐘周期。

求 C-1=(L-1)T×D-1×L-1時可以先求 A=(L-1)T×D-1，然后再計算 C-1=A×L-1。

4 仿真及結果分析

考慮到兩個模塊的最大運行頻率不一樣，本設計在尋找非零元素部分采用85 MHz的時鐘，在求解非零元素值部分采用65 MHz的時鐘。為了進行更好的對比，在MATLAB上用相同的算法、測量矩陣和觀測矢量來重構原始估計值。當α=0.25時，軟件和硬件的重構結果進行歸一化后的對比如圖3所示。

圖3 軟硬件仿真結果對比（歸一化后）

當α取值為零時，尋找非零元素部分共需要2 100個時鐘周期，而僅僅是計算內積就需要256×8=2 048個時鐘周期，計算非零元素部分共需要110個時鐘周期，總的重構時間為26.40μs。當α取值為0.25時，計算內積所需減少到約1 300個時鐘周期，總的重構時間減少到約 16.99μs。在相同條件下，參考文獻[7]重構時間為17.61μs。而在參考文獻[4]中，測量矩陣維數(shù)為 32×128，觀測向量維數(shù)為 32×1，原始信號的稀疏度為 5，總的重構時間就需要 24μs。

但是改進OMP算法歸一化誤差會隨著α的增大而增大，當α取值為零時，歸一化均方誤差為 0.001 5，取α=0.25時，歸一化均方誤差增加到0.007 1。

5 結論

本文采用一種閾值法，使得OMP恢復算法的求內積次數(shù)大大減少，從而縮短了信號重構所需要的時間，提高了恢復速率。同時，本文在硬件結構設計上也進行了一些優(yōu)化，較好地平衡了占用資源和運算時間。本設計采用VHDL對改進的OMP算法進行了RTL級描述，在Quartus II上針對Altera公司的 Cyclone II EP2C70F672C6進行設計和仿真，結果表明信號能夠以更少的重構時間較好地恢復。

[1]DONOHO D L.Compressed sensing[J].Information Theory，IEEE Trans.on，2006，52(4)：1289-1306.

[2]TROPP J A，GILBERT A C.Signal recovery from random measurements via orthogonal matching pursuit[J].Information Theory，IEEE Trans.on，2007，53(12)：4655-4666.

[3]趙貽玖.稀疏模擬信號壓縮采樣與重構算法研究[D].成都：電子科技大學，2012.

[4]SEPTINUS A，STEINBERG R.Compressive sampling hardware reconstruction[C].Circuits and Systems(ISCAS)，Proc.of 2010 IEEE International Symposium on.IEEE，2010：3316-3319.

[5]BLACHE P，RABAH H，AMIRA A.High level prototyping and FPGA implementation of the orthogonal matching pursuit algorithm[C].Information Science，Signal Processing and their Applications(ISSPA)，2012 11th International Conference on.IEEE，2012：1336-1340.

[6]STANISLAUS J L V M，MOHSENIN T.High performance compressive sensing reconstruction hardware with QRD process[C].Circuits and Systems(ISCAS)，2012 IEEE International Symposium on.IEEE，2012：29-32.

[7]STANISLAUS J，MOHSENIN T.Low-complexity fpga implementation of compressive sensing reconstruction[C].International Conference on Computing，Networking and Communications.2013.

[8]BAI L，MAECHLER P，MUEHLBERGHUBER M，et al.High-speed compressed sensing reconstruction on FPGA using OMP and AMP[C].Proc.19th Int.Conf.Electronics，Circuits and Systems(ICECS)，Dec.2012：53-56.

[9]周殿鳳，王俊華.基于 FPGA的 32位除法器設計[J].信息化研究，2010(3)：26-28.