張金鵬，錢 慧，王仁平

(福州大學(xué) 物理與信息工程學(xué)院，福州 350116)

1 引 言

壓縮感知(Compressed Sensing，CS)[1，2]是近十多年來信號處理領(lǐng)域最著名的研究進(jìn)展之一.CS主要由壓縮投影和重構(gòu)兩個步驟組成.在壓縮投影階段，利用稀疏基和觀測矩陣構(gòu)建CS矩陣，將稀疏信號在CS矩陣中進(jìn)行投影，獲得信號的觀測向量.在重構(gòu)階段，CS通過貪婪追蹤、凸優(yōu)化等方法從觀測向量中重構(gòu)出原始信號.這些CS重構(gòu)算法以非線性計算為基礎(chǔ)，是典型的計算密集型算法.在早期CS系統(tǒng)中，一般通過離線計算重構(gòu)原始樣本.但是隨著物聯(lián)網(wǎng)、人工智能技術(shù)的發(fā)展，便攜式應(yīng)用逐漸普及，實現(xiàn)CS快速重構(gòu)成為研究領(lǐng)域的關(guān)注重點[3].

正交匹配追蹤算法(Orthogonal Matching Pursuit，OMP)[4,5]是CS系統(tǒng)中討論最廣泛的重構(gòu)算法之一.近幾年，研究學(xué)者對面向OMP算法的硬件加速計算進(jìn)行了大量的研究.文獻(xiàn)[6-10]使用改進(jìn)的Cholesky分解算法優(yōu)化最小二乘(Least Squares，LS)問題，避免了平方根運(yùn)算，減少了計算延時.其中在文獻(xiàn)[6]中提出了第一個OMP算法的硬件實現(xiàn)設(shè)計，獲得了比CPU和GPU[11]更快的重構(gòu)速度.文獻(xiàn)[7]提出了一個新的高度流水線化的OMP硬件實現(xiàn)架構(gòu)，通過適當(dāng)水平的展開最大化并行度，通過共享硬件降低資源使用量.文獻(xiàn)[12，13]使用的是QR分解算法實現(xiàn)最小二乘求解.其中文獻(xiàn)[12]中使用了一種基于改進(jìn)的Gram-Schmidt算法的增量式QR分解算法優(yōu)化LS問題，通過每次迭代只更新Q矩陣和R矩陣的一列減少了存儲資源消耗，并設(shè)計了一種高并行度的硬件架構(gòu)，該架構(gòu)的硬件復(fù)雜度幾乎不受稀疏度的影響.文獻(xiàn)[13]提出了一種無方根增量式QR分解算法，進(jìn)一步消除了平方根操作，轉(zhuǎn)換為一些更簡單的操作，降低了計算的復(fù)雜度.另外文獻(xiàn)[14-16]使用了一種矩陣求逆旁路(Matrix inversion bypass，MIB)技術(shù)，通過分解矩陣求逆的過程有效降低了算法的計算復(fù)雜度.上述方法都是針對LS問題優(yōu)化算法，然后用硬件實現(xiàn)加速重構(gòu)，獲得了一定的加速效果，但是重構(gòu)速度會隨著迭代次數(shù)的增加而迅速降低.

從減少矩陣求逆次數(shù)的角度出發(fā)，本文提出了一種近似OMP重構(gòu)算法.在前K-1次迭代中用最大值計算方法近似替代LS計算去索引最大列，只在第K次迭代時進(jìn)行LS計算重構(gòu)原始信號估計值，僅需要做一次矩陣求逆計算.設(shè)計了一個近似OMP算法的硬件實現(xiàn)架構(gòu)，主要包含計算模塊、存儲模塊和控制模塊3個部分.該架構(gòu)用Verilog HDL語言進(jìn)行描述，在Xilinx公司的Vivado軟件上進(jìn)行綜合仿真實驗.實驗結(jié)果表明，本文設(shè)計的硬件架構(gòu)重構(gòu)速度快同時具備一定的可擴(kuò)展性.

2 正交匹配追蹤算法及近似實現(xiàn)

2.1 正交匹配追蹤算法

OMP算法基本思想是在每次迭代中，通過計算當(dāng)前殘差和CS矩陣的內(nèi)積值找到CS矩陣?yán)锱c當(dāng)前殘差最匹配的一列，然后通過去掉這一列的貢獻(xiàn)計算一個新的殘差，通過迭代計算最后輸出原始信號的估計值.OMP算法流程如下：

輸入：CS矩陣A=ΦΨ，隨機(jī)解調(diào)架構(gòu)矩陣Φ∈RM×N，傅里葉變換基Ψ∈RN×N，觀測向量y∈RM，稀疏度為K；

初始化：迭代次數(shù)t=1，支撐索引集Λ為空集，r0=y；

迭代過程：

2.更新索引集Λt=Λt-1∪λt，更新原子集At=At∪φλt；

5.判斷，若迭代次數(shù)t

2.2 近似正交匹配追蹤算法

輸入：CS矩陣A=ΦΨ，隨機(jī)解調(diào)架構(gòu)矩陣Φ∈RM×N，傅里葉變換基Ψ∈RN×N，觀測向量y∈RM，稀疏度為K；

初始化：迭代次數(shù)t=1，支撐索引集Λ為空集，r0=y；

迭代過程：

5.判斷，若迭代次數(shù)t

C=LDLT

(1)

其中L是對角線元素為1的下三角矩陣，D是對角矩陣，按以下公式計算求出C-1：

(2)

(3)

(4)

C-1=(L-1)TD-1L-1

(5)

矩陣D-1可以通過將矩陣D對角線上的元素求倒數(shù)獲得.

本文在MATLAB中進(jìn)行近似OMP算法的仿真驗證.原始信號是由兩個不同頻率的正弦波相加得到的，一個正弦波的頻率從20Hz增加到119Hz，另一個從150Hz增加到159Hz，步進(jìn)都是1，兩個信號可以組合出1000個樣本，這1000個樣本有4000個最大列，經(jīng)過仿真驗證，近似OMP算法找對了全部的最大列并成功重構(gòu)原始信號，驗證了近似OMP算法的可行性.

3 硬件系統(tǒng)設(shè)計

3.1 整體架構(gòu)

算法整體硬件架構(gòu)如圖1所示，主要分為計算模塊、存儲模塊和控制模塊.計算模塊完成算法迭代過程中的所有運(yùn)算，比如矩陣內(nèi)積、最大值計算、殘差計算、矩陣求逆的運(yùn)算等等.存儲模塊用來存儲所有的矩陣、向量和數(shù)據(jù).控制模塊負(fù)責(zé)控制算法的各個狀態(tài)的轉(zhuǎn)換以及從存儲模塊中讀取相應(yīng)的數(shù)據(jù)到計算模塊中和將結(jié)果寫回到存儲模塊中.

圖1 近似OMP算法硬件架構(gòu)Fig.1 Structure of approximate OMP algorithm

3.2 計算模塊

計算模塊又分為矩陣內(nèi)積單元、最大值計算單元、殘差計算單元、C矩陣求逆單元和最小二乘計算單元.

3.2.1 矩陣內(nèi)積單元

矩陣內(nèi)積單元由256個復(fù)數(shù)乘法器構(gòu)成的并行乘法器和255個復(fù)數(shù)加法器構(gòu)成的加法器樹和一個比較器組成.一個復(fù)數(shù)乘法器通過調(diào)用4個DSP資源實現(xiàn)，總共調(diào)用了1024個DSP資源.加法器通過調(diào)用slices資源實現(xiàn).每個周期從存儲模塊中讀取A的一列和y或(r)進(jìn)行內(nèi)積操作，結(jié)果送入到比較器進(jìn)行比較，并記錄較大列的列序號.最終在比較完所有列的內(nèi)積后，輸出內(nèi)積值最大的列的列序號.硬件架構(gòu)如圖2所示.最大值計算單元、殘差計算單元和最小二乘計算單元的大部分運(yùn)算都可以復(fù)用這里的并行乘法器和加法器樹結(jié)構(gòu)實現(xiàn)，從而降低硬件資源使用量.

圖2 矩陣內(nèi)積單元架構(gòu)Fig.2 Structure of matrix inner product unit

3.2.2C矩陣求逆單元

矩陣C和C-1是正定對稱矩陣，只需要求出C-1對角線上的元素和下三角部分的元素.首先需要根據(jù)公式(2)和公式(3)求出L矩陣和D矩陣的元素.其中L是對角線元素為1的下三角矩陣，D是對角矩陣.然后L矩陣和D矩陣分別求逆，最后利用公式(5)求出C-1.硬件實現(xiàn)架構(gòu)如圖3所示.

圖3 C矩陣求逆單元架構(gòu)Fig.3 Structure of C matrix inversion unit

在求矩陣D-1的元素時涉及到除法，本文選用的是比較經(jīng)典的牛頓-拉夫森方法求倒數(shù).該方法是找到一個在零點有解X=1/M的函數(shù)f(X)，一般這個函數(shù)如式(6)所示：

f(X)=1/X-M

(6)

這個函數(shù)在零點的解就是M的倒數(shù).牛頓-拉夫森方法給出了迭代公式：

(7)

根據(jù)文獻(xiàn)[7]中確定的初始值X0=3-2M，利用公式(7)進(jìn)行迭代計算，Xi的值就會收斂于M的倒數(shù).用這種方法將除法運(yùn)算轉(zhuǎn)化成了乘法和減法運(yùn)算.硬件實現(xiàn)架構(gòu)如圖4.

圖4 牛頓-拉夫森求倒數(shù)硬件架構(gòu)Fig.4 Structure of Newton-Raphson inversion

然后利用公式(4)求L-1矩陣的元素.最后可以根據(jù)公式(5)計算矩陣C-1.如果想要設(shè)計能重構(gòu)稀疏度更高的信號的硬件架構(gòu)，只需要修改最小二乘單元的硬件架構(gòu)，而不需要改變其它部分的硬件架構(gòu)，具有一定的可擴(kuò)展性.

3.3 存儲模塊

存儲模塊主要用來存儲算法計算所需要的數(shù)據(jù)，包括CS矩陣、觀測向量、稀疏基、殘差等和從計算模塊寫回的各種數(shù)據(jù).各個數(shù)據(jù)的存儲具體情況如表1所示.為了能夠在一個周期讀出矩陣一列的數(shù)據(jù)并行的送入到并行乘法器中，本文將矩陣一列的數(shù)據(jù)存在BRAM同一個地址中.其中CS矩陣A用單獨(dú)的一塊BRAM存儲，地址是0到1023，傅里葉變換基Ψ用單獨(dú)的一塊BRAM存儲，地址是0到1023，觀測向量y和殘差r存儲到相同的一塊BRAM中，觀測向量存到0地址，殘差r按迭代順序存在y向量后面.本文設(shè)計的架構(gòu)比其它相似設(shè)計多存儲了一個傅里葉變換基Ψ，所以消耗的BRAM資源會多一些.

表1 內(nèi)存資源使用情況Table 1 Memory resources usage

3.4 控制模塊

控制模塊主要通過狀態(tài)機(jī)完成對各個運(yùn)算子模塊工作條件的控制.狀態(tài)機(jī)的狀態(tài)轉(zhuǎn)換圖見圖5，各狀態(tài)說明如表2所示.每個子模塊都有一個控制信號，在子模塊需要工作的時候狀態(tài)機(jī)將該子模塊的控制信號置1，當(dāng)任務(wù)完成時，子模塊反饋給控制模塊一個信號，然后狀態(tài)機(jī)將該子模塊的控制信號置0，進(jìn)入下一個狀態(tài).當(dāng)?shù)螖?shù)t

圖5 狀態(tài)轉(zhuǎn)換圖Fig.5 Diagram of FSM

表2 狀態(tài)機(jī)的各狀態(tài)說明Table 2 Description of FSM

4 實驗結(jié)果與分析

4.1 實驗結(jié)果

本文實驗步驟如下：1)輸入原始信號f由兩個頻率分別為50Hz和100Hz，幅值分別0.3和0.6的正弦波相加得到，采樣序列1：1024，采樣頻率為1024，然后對原始信號f進(jìn)行歸一化處理；2)在MALTAB中運(yùn)行近似OMP算法程序，獲得CS矩陣A，傅里葉變換基Ψ，觀測向量y的原始數(shù)據(jù)，并對數(shù)據(jù)進(jìn)行定點量化處理；3)將得到的數(shù)據(jù)導(dǎo)入coe格式文件，在Xilinx公司的開發(fā)工具Vivado 2016.2中進(jìn)行仿真測試；4)將FPGA中仿真得到的結(jié)果導(dǎo)入MATLAB中，在MATLAB中進(jìn)行反傅里葉變換得到重構(gòu)信號，并與原始信號進(jìn)行對比，驗證是否重構(gòu)成功.

當(dāng)K=4時，近似OMP算法在FPGA中輸出原始信號在頻域的值，是放大了2^5倍的定點數(shù)據(jù)，分別為235-305i、214+308i、57+176i、58-168i，對應(yīng)的在MATLAB中仿真的數(shù)據(jù)是209-300i、210+301i、55+173i、54-169i.圖6是K=4時的仿真結(jié)果，依次為：原始信號、在MATLAB中用經(jīng)典的OMP算法重構(gòu)的信號、近似OMP算法在FPGA中的輸出結(jié)果重構(gòu)的信號.

圖6 重構(gòu)結(jié)果圖Fig.6 Diagram of reconstruction results

本文用重構(gòu)信噪比(recovery signal-to-noise ratio，RSNR)對重構(gòu)誤差進(jìn)行分析，公式如下：

(8)

圖7 RSNR曲線圖Fig.7 Curve of RSNR

4.2 分析比較

在表3中列出了與其它相似設(shè)計比較的具體情況.本文的稀疏基是傅里葉變換基，因此設(shè)計的是復(fù)數(shù)硬件架構(gòu)，在矩陣內(nèi)積模塊本文調(diào)用了1024個DSP資源實現(xiàn)復(fù)數(shù)并行乘法器，降低了slices資源的使用量.由于需要存儲額外的傅里葉變換基，增加了BRAM資源的使用量.對比稀疏度K=36時的重構(gòu)時間，本文的重構(gòu)速度快了1.25-1.72倍.

表3 實現(xiàn)結(jié)果與其它相似設(shè)計的比較Table 3 Comparison of implementation results with other similar existing designs

5 結(jié) 論

本文提出了一個新的近似OMP算法，從優(yōu)化算法的最小二乘問題角度出發(fā)，使算法只需在第K次迭代時進(jìn)行最小二乘計算重構(gòu)原始信號估計值.在其它K-1次迭代中用最大值計算方法替代最小二乘計算去索引最大列，從而減少了矩陣求逆的次數(shù).設(shè)計了近似OMP算法的硬件實現(xiàn)架構(gòu)主要包含計算模塊，存儲模塊和控制模塊3個部分.該架構(gòu)用Verilog HDL語言進(jìn)行描述，在Xilinx公司的Vivado軟件上進(jìn)行綜合仿真實驗.對比現(xiàn)有相似實現(xiàn)設(shè)計，本文設(shè)計的架構(gòu)通過調(diào)用更多的DSP資源減少了slices資源的使用量，在208MHz頻率下重構(gòu)速度提升了1.25-1.72倍，同時具備一定的可擴(kuò)展性.