基于FPGA的改進的排序QR分解實現

陳健，莊耀宇，楊丹，張俊杰

［1.上海大學通信與信息工程學院，上海 200444；2.特種光纖與光接入網重點實驗室（上海大學），上海200444］

隨著人們對移動通信需求的急劇增加，復雜的數據通信和處理將成為移動通信中一個重要挑戰［1-4］.多輸入多輸出（Multiple-Input Multiple-Output，MIMO）技術的高帶寬、低延遲、強魯棒性的特性使其成為下一代移動通信的關鍵技術［5］.

雖然MIMO 技術在不增加額外頻譜和功率消耗的情況下能顯著提高通信系統的容量和質量，但是與傳統的單輸入單輸出（Single-Input Single-Output，SISO）通信系統相比，MIMO檢測需要包含空間維度在內的多維度信號處理，復雜度更高.與此同時，MIMO系統發射天線所引起的同信道干擾為信號的檢測增加了難度［6］.因此，近些年來，人們針對不同的MIMO傳輸方案提出了不同的檢測算法.

MIMO 檢測算法主要分為線性檢測算法和非線性檢測算法兩大類.雖然線性檢測算法復雜度較低，但檢測性能不如非線性檢測算法.迫零（Zero-Force，ZF）檢測算法和最小均方誤差（Minimum Mean Square Error，MMSE）檢測算法［7］是兩種典型的線性檢測算法.它們利用信道矩陣得到加權矩陣，對接收信號進行判決.為了提高線性檢測的精度，排序的連續干擾消除（Ordered Successive Interference Cancellation，OSIC）方法［8］采用一組線性接收機，每個接收機檢測數據流中的一個符號，其次，在每個階段都能夠成功地從接收信號中刪除已檢測出的信號成分，從而減少對后續待檢測信號的干擾.

非線性檢測算法中，最大似然（Maximum Likelihood，ML）檢測是一種最優的檢測算法［9-10］.傳統的ML 檢測通過計算接收信號向量和所有可能的后處理向量之間的歐氏距離，找到一個最小距離來進行判決.由于ML 要遍歷搜索所有可能的發送信號，其復雜度隨信號調制階數和發射天線數呈指數上升.文獻［1］提出了一種基于Cholesky分解和k-best檢測相結合的檢測方案.其中k-best［11］是一種次優檢測算法，對矩陣進行Cholesky［12］分解是為k-best檢測做預處理.通過選擇合適的k值使得檢測性能逐漸逼近ML.然而對矩陣進行Cholesky 分解，該方法要求該矩陣正定對稱，適用范圍受限.QR 分解（QR Decomposition，QRD）算法［13-14］不但可以為諸如k-best 檢測算法做預處理，其本身也可作為一種非線性檢測手段進行信號檢測.QR 分解通過對接收到的信號進行分層檢測從而消除符號間的干擾.同時QR分解避免了復雜的矩陣求逆運算且適用性最廣.但是在分層檢測的過程中，容易出現誤差累積效應.因此，基于排序的QR分解算法被提出［15-16］.

研究者們針對QR 分解的實現提出了不同方案.矩陣進行QR 分解主要包括格拉姆-施密特（Gram-Schmidt，GS）正交化［17］，豪斯霍爾德（Householder）變換［18］和吉文斯（Givens）變換［19］三種方式.由于Givens變換只包含酉矩陣的乘法且運算邏輯簡單，更適合硬件實現.文獻［14］和文獻［20］提出了一種基于Givens 變換的QR 分解的FPGA 實現方案，但是都沒有引入排序過程.本文主要工作如下：

1）通過對各類QR 分解實現方案的對比，選用Givens 變換作為QR 分解的實現方案并基于脈動陣列設計相應的電路結構.

2）由于L2 范數的計算復雜度較高，本文改進了排序QR 分解中的基于L2 范數的排序方案，采用L1范數進行排序.Matlab 性能仿真顯示，改進后的方案對系統性能的影響幾乎不變.同時L1 范數在單個列范數計算中組合邏輯資源比L2 范數至少節省29.2%，觸發器資源至少節省32.4%.

3）針對排序QR分解的實現特點，本文設計了基于串行迭代的排序單元.利用該排序單元對N個數據進行排序，只需要N個比較器，且迭代延遲為N+1個時鐘周期.排序的QR分解中只需要進行最小值交換，因此最終只會產生1個時鐘的排序延遲.

1 MIMO系統檢測

1.1 系統模型

圖1給出了MIMO系統中一個典型的信號傳輸模型.發射天線數為NT，接收天線數為NR.發射信號矢量，接收信號矢量

圖1 MIMO系統模型Fig.1 The model of MIMO system

信道傳輸矩陣H表示為：

式中：hji表示第i根發射天線和第j根接收天線之間的信道增益，它們統計獨立且服從高斯隨機分布［7，21］.在本文中，我們將一直假定信道信息已知.則該NR×NTMIMO系統可以表示為：

式中：z=是均值為0，方差為的加性高斯白噪聲.

1.2 QR分解檢測

相較于ZF 檢測和MMSE 檢測，QR 分解檢測避免了復雜的求解偽逆矩陣的過程，受到了人們的廣泛關注.NR×NT（NR≥NT）的MIMO 信道傳輸矩陣H，進行QR分解后可以得到：

式中：Q是NR×NR的酉矩陣且滿足QHQ=I；R是NT×NT的上三角矩陣.接收信號y左乘QH得到：

式中：=QHz.為方便闡述，我們假設NR=NT=N，并將式（4）改寫為：

由于Q是酉矩陣，因此QHz不會放大噪聲.考慮到第N層信號的檢測不受其他層的干擾，先從第N層開始檢測，得到第N層信號的估計值：

式中：Q[·]表示按照信號調制方式作判決.利用第N層得到的估計值消除對第N-1 層的干擾，從而得到第N-1 層的估計值依此類推，那么第k層的估計值可以表示為：

顯然，這種方式容易存在持續的誤差累積.為了降低誤差傳播，理想的檢測是使每一層的SNR 盡可能大，第k層的信噪比可以表示為：

從式（9）可以看出，每層的SNR 正比于R矩陣的主對角線元素，因此應使QR 分解后得到的R矩陣的主對角線元素盡可能從左往右單調不減，即rNN≥rN-1N-1≥… ≥r22≥r11.這便引出了基于排序的QR分解.

1.3 基于Givens變換的排序的QR分解

圖2 以4 × 4 復矩陣的分解為例，展示了不同分解方法的計算復雜度［22］.整體來看，Householder 變換復雜度最高，GS 變換復雜度較低.由于Givens 旋轉只包含酉矩陣的乘法，被認定為是一個數值穩定的算法，且可以通過旋轉的方式來代替除法和開方運算.因此，在本文中，我們采用Givens 旋轉的方式來實現QR分解.

圖2 分解方法的復雜度Fig.2 The complexity of decomposition methods

以2 × 2 復矩陣為例，闡述Givens 旋轉過程.給定復矩陣H：

為了使左下角的元素H10變為0，復數Givens 旋轉運算如下：

式中：θ0==-∠H00；θ2=-∠H10；R10=0.排序的QR 分解檢測是在分解的過程中引入排序，使得R矩陣的主對角線元素盡可能按照升序排列.對矩陣進行排序的QR 分解時，首先通過對比原始矩陣列向量范數的大小，對矩陣的列向量進行交換.利用交換后的矩陣進行Givens 旋轉.在每輪旋轉完成后，需要再次排序并進行列向量交換.一個N列的矩陣H一共需要N-1 次迭代排序.在本文中，我們以列向量的L1 范數的大小作為參考進行排序.Givens 旋轉排序的QR 分解的偽代碼邏輯如表1 所示.此外，當矩陣H列欠秩時，對矩陣H進行QR 分解后，會導致分解得到的R矩陣主對角線上的元素有0的情況.結合QR 分解的檢測流程，當其中某一層主對角線為0 時，會產生嚴重的誤判.因此，QR 分解不太適用于欠秩的MIMO信道.

表1 Givens旋轉的排序QR分解Tab.1 Sorted QR decomposition based on Givens rotation

2 FPGA實現方案

2.1 整體結構設計

為了硬件實現排序的QR 分解，我們提出的整體結構如圖3 所示.該結構主要由三個部分組成：信道矩陣的存儲RAM、范數排序單元以及Givens 迫零旋轉單元.為了方便范數計算以及Givens旋轉，信道矩陣以列向量的形式存儲到RAM 中.范數排序單元對RAM 地址中的列向量重新進行映射.映射的結果同樣按列向量的形式保存到由RAM 構成的Buffer 中.Givens 迫零旋轉單元一次只對一列進行迫零.對于NR×NT的信道矩陣H，則需要NT次迭代且在每一次迭代之前重新對列向量排序.

圖3 排序QR分解結構Fig.3 The architecture of sorted QR decomposition

2.2 排序方案設計

相較于傳統的QR 分解算法，排序的QR 分解在實現過程中不可避免地會引入排序延遲.目前研究者們基于FPGA 的并行排序已經設計出了不同的方案［23-26］.考慮在QR 分解的實現方案中，以矩陣列向量串行輸入，通過構建脈動陣列的方式，達到節省資源的目的.因此，本文構建了基于迭代的串行排序單元.圖4 展示了該單元的硬件結構.它由2 個數據寄存器、2個多路選擇器和1個比較器組成.

圖4 串行迭代排序單元Fig.4 Sorting unit based on serial iteration

我們假設當前按照升序排序.寄存器A 在初始化時，存儲當前位寬下所能表示的最大無符號數.當待排序的數據被送入該排序單元，立刻與寄存器A中的數據進行比較.較小的數據被寫入寄存器A，較大的數據被寫入寄存器B 用于后續的迭代比較.在排序的QR 分解中，我們需要找到L1 范數最小的列向量.因此我們做了進一步的簡化.簡化后的排序單元只需要1個寄存器、1個比較器和1個多路選擇器.但多了一個索引計數器.索引計數器的計數值范圍為0～3，初始化為0.待排序的數據以每個時鐘1 個數據的方式串行送入該排序單元.信號input_valid 高電平表明當前輸入數據有效.同時索引計數器和寄存器A 利用該信號完成數據的更新.當數據D 被送入該單元后，立刻與寄存器A 中的數進行比較.如果D 小于A，則將數據D 連同當前的索引計數值一同寫入寄存器A，否則寄存器A中的數保持不變.檢測input_valid 的下降沿作為排序結束標志，寄存器A 和索引計數器同時重新初始化.

2.3 列向量映射交換設計

如2.1 小節所述，整體結構中的范數排序單元，除了完成對信道矩陣H列L1 范數的排序，還需要依據排序的結果對矩陣列向量進行重新映射.這里以4 × 4 的H矩陣第一輪迭代為例說明映射的過程.為了方便描述，我們定義矩陣的列向量從左至右依次為第一列、第二列，并以此類推.排序QR分解的目的在于使得分解后的R矩陣的主對角線元素從左至右盡可能按照升序排序.因此我們選擇L1 范數最小的列與第一列進行交換.在計算L1 范數時，我們從左至右依次讀取每一列數據，然后按照讀順序，通過修改RAM 寫地址的方式，完成矩陣的列交換.圖5 給出了對應的讀寫流程.如果排序得到的最小范數由第一列產生，那么列向量的寫地址依次為｛0，1，2，3｝；如果由第二列產生，列向量的寫地址依次為｛1，0，2，3｝；如果由第三列產生，寫地址依次為｛2，1，0，3｝；如果由第四列產生，寫地址依次為｛3，1，2，0｝.

圖5 排序映射邏輯Fig.5 Sorting and mapping logic

2.4 Givens迫零旋轉設計

任何n階復非奇異矩陣H=（hj，i）可通過左乘有限個復初等旋轉矩陣轉化為上三角矩陣，且對角線元素除最后一個外都是正數［27］.參考式（11），在對第i列第j行元素hj，i迫零時，利用hj，i和hj-1，i得到旋轉角度，且迫零運算過程中，矩陣的乘法只對第j行元素和第j-1 行元素產生影響.圖6 為4 × 4 矩陣的Givens 旋轉運算的過程，首先利用X31和X41對X41進行迫零，然后利用第一次迫零的結果R31和X21對R31進行迫零，最后利用第二次迫零的結果R21和X11對R21進行迫零.至此第一列屬于下三角的元素全部變為0.同樣的流程，第二列需要2 次迫零運算，第三列需要1次迫零運算.

圖6 4 × 4矩陣Givens旋轉過程Fig.6 The process of Givens rotation on 4 × 4 matrix

為了得到旋轉所需要的角度以及完成相應的矩陣乘法運算，設計了圖7 所示的硬件單元，包括邊界單元和迫零運算單元.在邊界單元中，通過調用Xilinx Cordic IP 核得到旋轉矩陣所需要的角度，迫零運算單元主要利用邊界單元得到的角度，完成復數乘法運算.

圖7 QR分解硬件組成單元Fig.7 The hardware units of QR decomposition

結合式（11），圖8 展示了邊界單元具體的角度運算流程.由于θ0的計算需要先得到H00和H10的模，然后再求解反三角得到角度，最后才能得到正弦和余弦的值.相較于θ1和θ2的計算，多出了求解反三角的延遲.我們對（11）進行改寫，得到：

圖8 旋轉角度計算流程Fig.8 The calculation process of rotation angle

我們通過在迫零計算單元中先計算(cosθ1+jsinθ1)H00和(cosθ2+jsinθ2)H10，從而可以將求解反三角函數的延遲隱藏在迫零運算中，達到減少等待延遲的目的.圖9 展示了基于邊界單元和迫零計算單元的脈動陣列結構.該陣列結構對矩陣第一列除主對角線元素進行迫零運算.通過逐列讀取RAM 的方式，完成一輪迫零迭代.

圖9 Givens旋轉脈動陣列Fig.9 The systolic array of Givens rotation

3 仿真與結果分析

3.1 Matlab仿真分析

本文利用Matlab 對不同算法的性能進行對比分析.實驗仿真環境選擇數據幀長度為10 000，調制方式采用16-QAM（Quadrature Amplitude Modulation，QAM）調制，信道為平坦瑞利衰落信道，噪聲為均值為0 的加性高斯白噪聲.圖10 和圖11 以4 × 4 MIMO模型給出了檢測算法在不同的SNR 下的誤符號率.從圖10 中可以看出，系統誤符號率從低到高依次為：ML，Givens 旋轉的排序QR 分解（Givens SQR），QR分解，MMSE，ZF.

圖10 各類算法性能仿真Fig.10 The performance simulation of different algorithms

圖11 L1范數和L2范數性能對比Fig.11 The performance comparison between L1-norm and L2-norm

圖11 給出了分別利用L1 范數和L2 范數進行排序對系統性能的影響，其中Givens SQR1 表示基于L1 范數排序，Givens SQR2 表示基于L2 范數排序.雖然在Givens 旋轉時引入了排序，但是并不能使得分解后的R矩陣的主對角線元素從左至右嚴格按照遞增的要求排序.從圖11 的仿真結果來看，兩種排序的方案所起到的效果幾乎相同.

表2 給 出 了 在4 × 4、32 × 32、64 × 64 三 種MIMO 模型下的不同SNR 的誤符號率.從仿真結果來看，所提出的算法在多種場景下，其算法性能具有一致性.當收發天線規模擴大后，算法的檢測性能也進一步提高.

表2 算法性能對比Tab.2 The comparison of algorithm performance

由于浮點數運算會占用較多的硬件資源，因此在本文中采用定點數的方式實現整個邏輯的設計.圖12 以相同的仿真環境展示了不同精度的數據表示對系統性能的影響.8 位定點小數性能較差.10 位定點小數性能略低于浮點表示，但基本滿足檢測需要.為了保證運算的精度以及防止數據溢出，我們采用16 bit量化的方式.其中1 bit符號位，5 bit整數位，10 bit小數位.數值的表示范圍為-32≤n≤31.999 023 4，精度為1/210.

圖12 數據精度對比Fig.12 The comparison on different precisions

3.2 硬件邏輯仿真

為了驗證本文所提結構的正確性.利用Mentor Graphics 的Modelsim 仿真平臺搭建仿真環境并與Maltab 的結果進行對比.在Matlab 端產生了4 × 4 的平坦瑞利衰落信道并按照3.1小節所述進行量化，量化后的結果作為激勵導入到所提設計中.

圖13 展示了對讀取的矩陣列向量求解L1 范數并以其為參考進行排序以完成矩陣列向量交換的過程.Norm_1_result_valid 表明當前求解的范數有效，Norm_1_result_temp 表示某一列向量的L1 范數值.按照RAM 讀取順序，分別為第一列、第二列、第三列和第四列.Recent_min_norm 表明了當前排序中最小的列向量范數，其中低2 bit 為該范數所對應的列向量標號.Norm_sorting_finish 標志著列向量排序結束.Ram_data_o 表示按照讀順序需要重新寫入RAM 的數據.Ram_data_addr 為當前Ram_data_o 所要寫入的RAM 地址.由于當前迫零迭代下，第一列的范數最小，因此按照2.3 節的地址映射邏輯，寫地址依次為0，1，2，3.參考Matlab，L1 范數的求解和排序交換映射邏輯正確.

圖13 排序并交換列向量Fig.13 Sorting and exchanging matrix columns

圖14 展示了QR 分解的Modelsim 仿真和Matlab仿真的對比結果.此時Matlab 仿真得到的R矩陣主對角線元素恰好是單調不減的.在Modelsim 仿真圖中，Rotate_wr_ram_valid_temp 表明當前寫RAM 數據和地址有效.Rotate_wr_ram_addr_temp 表示當前時鐘周期下寫數據所對應的寫地址.Rotate_wr_ram_dout_temp 表示分解完成后矩陣列向量中的4 個元素.RAM地址0對應矩陣的第一列，地址1對應矩陣的第二列，地址2對應第三列，地址3對應第四列.仿真誤差范圍控制在10-3內.基本滿足作為原理模型的需要.

圖14 QR分解的結果Fig.14 The results of QR decomposition

3.3 FPGA布局實現

本文主要是利用矩陣L1 范數來替換L2 范數.對于4 × 4 復矩陣，每列中包含4 個復數.L1 范數可以直接利用Xilinx Cordic IP 核進行求解.L2 范數在忽略開方的情況下，可以利用復數乘法器乘共軛復數或者在模值的基礎上做平方運算獲得.本文基于Xilinx FPGA VC709 開發套件，以默認綜合選項對相關資源進行了評估.表3給出了Cordic 核進行8 次旋轉迭代，復數乘法器以及無符號數乘法器所需要的硬件資源，其中數據位寬都是16 bit.FF 表示觸發器資源，LUT表示查找表資源.因此，對于一個4 × 4 矩陣求某一列向量L1 范數時，需要2 424 個LUT，2 504個FF.而求L2 范數時，至少需要3 428 個LUT，3 708個FF.表4 中給出了整體資源使用情況及電路性能.為了便于流水線實現，我們對每一級運算結果都做了緩存，因此占用了較多的存儲資源.

表3 計算單元資源利用情況Tab.3 The resource utilization of computing units

表4 整體資源消耗情況Tab.4 The resource utilization of the entire design

圖15 展示了基于VC709 FPGA 設計的QR 分解硬件電路板級測試結果.利用Matlab 將經過量化后的兩個4 × 4 復矩陣導出成Coefficient 文件，作為板級測試的初始化文件，并保存在只讀存儲器ROM中.當FPGA 上電后，QR 分解處理單元循環地從ROM 中讀取數據進行運算.通過調用Xilinx 提供的集成邏輯分析儀IP 核對運算結果實時抓取.OUT_VALID 高電平表明當前時刻運算結果有效.Aij_I 表示矩陣第i行第j列元素的實部，Aij_Q 表示矩陣第i行第j列元素的虛部.圖中以矩陣第四列的結果進行了展示.板級測試結果與仿真預期結果一致，結果正確.

圖15 FPGA板級測試結果Fig.15 The test results on FPGA

4 結論

本文基于矩陣列向量的L1 范數設計了排序的QR 分解硬件電路結構.利用L1 范數進行排序，單列計算中LUT 資源至少節省29.2%，FF 資源至少節省32.4%.通過Matlab 仿真，所提出的優化算法在大規模MIMO 場景下仍有適用空間.針對范數排序問題，我們設計了基于串行迭代排序單元，只需要利用一個比較器便可完成一次迭代中的排序需要.在每次Givens 旋轉角度計算中，我們依靠6 個Cordic 核完成相應的運算.并基于脈動陣列構建了矩陣乘法運算單元.對4 × 4 復矩陣進行分解時，本文所提出的結構主頻性能得到明顯改善并通過仿真以及板級測試分析，電路功能正確.