獨立分量分析中NLPCA－RLS算法IP核的設計

徐長輝，李 宏，常 帥

（西北工業大學 電子信息學院，陜西 西安710129）

0 引 言

盲信號分離是指在信號混迭參數未知的情況下，從多個混合在一起的觀測信號中恢復出源信號。通常，觀測信號來自一組傳感器，每個傳感器的接收信號是源信號的不同組合。獨立分量分析是在盲源分離研究中發展起來的一種多通道信號處理技術，近年來，該技術在語音降噪、通信、生物醫學信號處理、圖像處理中的巨大優勢引起了廣泛的關注。

隨著FPGA技術的快速發展，有關獨立分量分析算法硬件實現的方案也不斷涌現。2006年Shyu和Li首先通過硬件描述語言（HDL）實現了FastICA算法，之后對該方案進行改進，實現了算法的流水線操作［1－2］。獨立分量分析的信息極大化算法和 EASI（equivariant adaptive separation via independence）等其他算法硬件方案的研究也已完成［3－6］。隨著相關算法和技術的不斷完善，獨立分量分析在工程實際中的應用研究也已展開，Huang等人實現了四通道腦電信號的實時分離［7］，而Thuong在FPGA中將獨立分量分析算法應用于音頻水印技術［8］。在ICA眾多算法中，NLPCA－RLS算法是一種隨機在線算法，更適用于對實時性處理要求較高的場合。本文研究了非線性主分量分析算法的原理、主要功能和分離過程，依照算法的流程進行了功能模塊的劃分，利用Simulink/DSP builder與QuartusII組合開發工具及Altera CycloneIII硬件仿真平臺，首次完成了NLPCA－RLS算法硬件IP的實現方案，并將其應用于語音信號的分離與降噪中。

1 獨立分量分析的基本模型及NLPCA－RLS算法

1.1 獨立分量分析的基本模型

獨立分量分析是一種基于統計量的信號處理方法，它將多通道觀測信號按照統計獨立的原則進行優化運算，分解成若干相互獨立的成分。其模型可描述為［9］

式中：x（t）——觀測信號矢量，s（t）——統計獨立的源信號矢量，A——混合矩陣。獨立分量分析的目的是在s（t）和A均未知的條件下，尋求解混矩陣W，使得

如圖1所示。

圖1 獨立分量分析模型

根據Darmois定理［10］，在如下兩個假設成立的條件下：

（1）源信號各分量之間統計獨立；

（2）源信號分量中最多只有一個高斯分量；

如果我們能找到矩陣W，使得（2）中的y（t）各分量統計獨立，那么y（t）即為s（t）的一個拷貝，即

式中：∧——一個對角矩陣，P——一個置換矩陣。即通過獨立分量分析得到的分離信號y（t），與源信號s（t）之間存在幅值和排列的不確定性。

1.2 NLPCA－RLS算法

NLPCA算法是獨立分量分析的一種典型算法，它是在主分量分析（principal component analysis，PCA）中引入非線性函數而得。PCA的一個基本概念是最小均方誤差壓縮，即對于隨機輸入向量x，尋找標準正交基向量W，使得x和它在該標準正交基向量張成的子空間上的投影之間的均方誤差最小，其目標函數為

當我們對于白化后的x用非線性函數g（Wx）代替（Wx），即用x在子空間上的非線性投影代替線性投影，上述方程的解在一定程度上能夠使y＝W＊x的各分量獨立［11］，這種方法叫做NLPCA算法。

NLPCA－RLS算法是在NLPCA算法中采用近似遞歸最小二乘。使用下標t代表迭代步數，對輸入數據x（t），迭代算法如下

式中：向量變量q（t），h（t），m（t），r（t）及矩陣變量P（t）——算法內部的輔助變量，g——NLPCA指標中的非線性函數，此處取反正切函數。參數β是一值接近1的“遺忘常量”。符號Tri的數學意義是只計算矩陣的上三角部分，然后將上三角部分求轉置并復制到下三角部分。經過迭代之后的矩陣W 即為式（2）中的解混矩陣。

上述NLPCA－RLS算法的完整實現步驟如圖2所示。

圖2 NLPCA－RLS算法流程

由于使用近似遞歸最小二乘方法，算法可以根據輸入數據決定學習速率參數。相比于其他算法，NLPCA－RLS具有更快的收斂速度［12］。

2 基于DSP builder的硬件設計

2.1 設計工具與流程

目前出現的基于FPGA的DSP開發工具，使得設計者能夠遵循類似于軟件設計流程的開發方法進行FPGA的設計。尤其當設計包含諸多數學算法時，可使效率大為提高。

本文利用 Matlab、DSP Builder、Modelsim、QuartusII等組合工具完成整個設計。①在Matlab/Simulink中進行設計輸入，即在Matlab的Simulink環境中建立模型文件，用圖形方式調用DSP Builder和其他Simulink庫中的圖形模塊，構成系統級和算法級設計框圖，或稱Simulink設計模型。②利用Simulink的圖形化仿真、分析功能，分析此設計模型的正確性，完成模型仿真。③由于EDA工具軟件不能直接處理Matlab的模型文件，需要通過DSP builder中的SignalCompiler模塊將模型文件自動轉換成通用的硬件描述語言文件，即VHDL文件。第四步，對以上頂層設計產生的VHDL文件在Modelsim中進行RTL級仿真。④對功能仿真正確的硬件代碼在Quartus中進行綜合、編譯、時序仿真并下載到FPGA中進行驗證。

2.2 整體設計方案及技術指標

本方案基于2x2的數據模型搭建，即源信號向量s（t）維數為2，混合后的信號向量x（t）維數為2。為減少復雜度，所處理數據采用24位固定點類型，其中整數部分占12比特，小數部分占12比特。數據處理采用批處理模式。根據MATLAB計算機仿真可知，算法在200次迭代之內即可收斂，每個數據RAM的大小由算法的總迭代步數決定， 即為200x24bit。設計框圖及關鍵信號連接如圖3所示。

圖3 NLPCA－RLS算法硬件實現框架

其中，SCU（System Contrl Unit）模塊控制系統的復位、啟動、運行及算法迭代流程的狀態機。程序運行時，待分離的兩路混合信號x（t）數據首先分別載入輸入端兩個雙端口RAM（RAM＿IN＿1和RAM＿IN＿2）中，接著NLPCA＿RLS算法模塊計算出分離矩陣W。混合信號與分離矩陣W 相乘后得到分離后的信號y（t），該數據存入輸出端雙端口RAM（RAM＿OUT＿1和RAM＿OUT＿2）中。在down＿o＝1，即運算完成后，外部控制器可讀取輸出端RAM中的數據。

2.3 遞歸最小二乘迭代算法的實現

由于系統級開發工具DSP Builder有專用模塊供調用，可大大提高開發速度，因此算法模塊采用該工具進行建模。圖3給出的迭代算法主要涉及三個方面：迭代流程的控制，矩陣乘法運算和反正切函數運算。其中矩陣乘法采用FPGA中的專用乘法模塊和加法模塊實現。各模塊的在DSP builder中的具體建模方法如下：

（1）算法流程控制

本文中采用有限狀態機的方式來控制算法頂層的迭代與循環，如圖4所示。

圖4 算法流程控制狀態機

程序啟動后首先進入Idel狀態；在外部輸入信號start＿i置1時跳轉到Clear狀態，此時模塊中的乘法器和加法器進行一次復位，完成后進入Load狀態；在Load狀態下，迭代式（5）中的W 和P矩陣將上一次迭代后的值載入，作為本次迭代的起始值（若為第一次迭代，則載入初值）；Calc狀態下程序進行一次迭代運算，當運算結束后狀態機進入Save狀態；Save狀態將本次計算得到的W 和P保存起來，之后程序判斷迭代次數是否達到，如果沒有達到，狀態機跳轉至Clear狀態繼續進行下一次迭代運算；否則跳轉到Down狀態，并將down＿o信號置1，告知外部模塊運算已完成。外部模塊將運算結果讀取之后復位start＿i信號，狀態機返回Idel狀態。

（2）矩陣乘法的實現

式（5）中涉及多個2位矩陣乘法運算，其一般格式為

圖5即為式（6）在DSP builder中的具體實現方案，圖中的乘法器單元可選用Altera FPGA的專用乘法器單元LPM－MULT，它是一個參數化的宏函數［13］。圖中各運算單元均使用統一的使能信號（ena）與異步清零信號（aclr），便于頂層狀態機控制整個流程。使用該矩陣乘法單元實現的整個算法流程，每次獨立迭代運算共需22個系統時鐘周期。

圖5 Simulink中二維矩陣乘法的實現

（3）反正切函數的實現

反正切函數q＝atan（y）相對比較復雜，如果通過計算的方式實現必然會占用較多的資源，同時延緩算法的速度。本方案采用查表法實現，查找表地址位寬為8，輸出數據位寬為14，故可存放128個14位固定點類型數據。變量y在之前的運算中已將范圍限制在 ［－4，4］，因此我們在8位地址線的查找表中線性存放 ［－4，4］的反正切函數值。根據計算可得，查找表的地址add與輸入變量y的關系為

在實際設計中，我們依然采用LPM參數化宏模塊例化LUT（lookup－table），即將生成的反正切函數數據存放在由嵌入式RAM構成的ROM中，這樣可以節省大量邏輯單元。圖6給出了反正切函數模塊的實現框圖。

圖6 Simulink中反正切函數的實現

3 系統仿真與分析

3.1 系統性能指標

按照DSP Builder的設計流程，首先在Simulink中對設計模塊進行仿真。然后進行編譯綜合，選用Altera CycloneIII EP3C25芯片開發板進行硬件仿真。該IP核在QuartusII中綜合后的結果見表1。系統采用批處理模式，每次處理200個數據，輸入輸出數據寬度為24位。由于每次迭代運算需22個時鐘周期，因此在系統輸入時鐘頻率為50MHz時，對外部數據的處理速率約為1MHz，完全可以實現常規采樣率下語音信號的實時處理。

表1 FPGA硬件資源占用

3.2 系統盲信號分離功能驗證

實驗中分別采用周期信號與語音信號驗證IP核的有效性。在周期信號的驗證中，我們將一段正弦信號與方波信號進行混合，混合矩陣為混合后的信號如圖7（a）所示。為了便于仿真，我們將上述待處理的混合信號數據寫入hex文件，作為輸入端RAM的初始化數據。圖7（b）顯示了該IP核對混合信號的處理結果，可以看出，程序能夠正確的提取出原始正弦信號與方波信號。

圖7 周期信號分離試驗

接下來，我們采用一路真實的語音信號與一路白噪聲信號的混合進行驗證。語音信號來源于SPIB（Signal Processing Information Base）數據庫［14］，采樣率為44kHz，采樣點數為200個。混合矩陣仍為分離結果如圖8所示。其中，圖8（a）所示為語音信號與周期信號的混合信號，圖8（b）所示為混合信號經IP核處理后復原出的原語音信號和白噪聲信號。

從以上仿真結果中我們可以看出，無論對于人工合成的周期信號還是真實語音信號，該IP核都能很好的完成NLPCA－RLS算法，對滿足盲信號模型數據進行處理，達到分離的目的。

3.3 系統準確性驗證

為了評價算法的性能，通常需要相應的性能指標。在獨立分量分析算法中，誤差指數E［15］是衡量算法收斂精度的常用指標，它用來表征算法對信道矩陣的估計值與真實值之間的相似性，相似性越大，E越小，信號分離性能越好

式中：cij——矩陣C＝WA 的第i行、第j列的元素，A——信道混合矩陣，W——信道解混矩陣。

圖8 語音信號分離試驗

在語音信號分離實驗中，我們記錄了誤差指數E隨算法收斂次數的變化，并將其與算法的matlab仿真結果進行了對比。為避免數據隨機性帶來的誤差，結果取1000次非重復實驗平均值，如圖9所示。

圖9 IP核與matlab結果比較

從圖9中可以看出，隨著收斂次數的增加，E逐漸收斂，并與matlab結果趨近。但由于在硬件設計中為提高運算速度，所有數據均采用24bit固定點類型，存在舍入誤差，因此最終結果與matlab理論值之間亦存在一定誤差。

4 結束語

本文提出了一種獨立分量分析中NLPCA－RLS算法的硬件方案，以IP核的形式在FPGA中實現。相比于傳統基于FastICA算法的模型，該方案不需要在每次迭代過程中使用全部樣本數據，大大減少了對硬件資源中存儲空間的占用。系統建模時采用了Altera優化的IP核和庫函數，可以更好的發揮FPGA的數據處理能力。所有數據均采用24位固定點類型，既滿足常規數據處理精度要求，又提高了運算速度。使用VHDL語言編寫的頂層系統狀態機，使各運算模塊更加簡潔高效，數據處理速率能達到1MHz。系統能夠準確完成盲信號分離功能，滿足工程設計中對實時性的要求，具有較高的實用價值。

［1］Shyu K K，LI M H.FPGA implementation of FastICA based on floating point arithmetic design for real－time blind source separation ［C］.International Joint Conference on the Neural Networks，2006：2785－2792.

［2］Shyu K K，LI M H，WU Y T，et al.Implementation of Pipelined FastICA on FPGA for Real－Time Blind Source Separation［J］.IEEE Transactions on Neural Networks，2008，19（6）：958－970.

［3］DU H，QI H，WANG X.Comparative study of VLSI solutions to independent component analysis ［J］.IEEE Trans Ind Electron，2007，54（1）：548－558.

［4］Charoensak C，Sattar F.Design of low－cost FPGA hardware for real－time ICA－based blind source separation algorithm ［J］.EURASIP Journal on Applied Signal Processing，2005，doi＞10.1155/ASP.2005.3076.

［5］Kim S J，Ken Umeno，Ryo Takahashi.FPGA implementation of EASI algorithm ［J］.IEICE Electron.Express，2007，22：707－711.

［6］Moreno L N O，Arce M A A，Lamont J G.Implementation of infomax ICA algorithm for blind source separation ［C］.Electronics，Robotics and Automotive Mechanics Conference，2008：447－451.

［7］HUANG W H，HUNG S H，CHUNG J F，et al.FPGA implementation of 4－channel ICA for on－line EEG signal separation ［C］.Biomedical Circuits and Systems Conference，2008：65－68.

［8］Thuong L T，Dien V N，LAN N H，et al.Independent component analysis applied to watermark extraction and its implemented model on FPGAs ［C］.Fifth IEEE International Symposium on Electronic Design，Test＆Applications，2010：71－76.

［9］Hyvarinen A，Karhunen J，Erkki O.Independent component analysis［D］.周宗譚，董國華，徐昕，等譯.北京：電子工業出版社，2007.

［10］Pierre Comon，Christian Jutten.Handbook of blind source separation：Independent component analysis and applications［D］.Academic Press，2010.

［11］Aleksandr Nikolaevich Gorban，Donald C Wunsch.Principal manifolds for data visualization and dimension reduction ［D］.Berlin Heidelberg：Springer－Verlag，2007.

［12］ZHANG Zhilin，PI Yiming，SUN Zhijian.Independent component analysis based denoising technology ［J］.Journal of UEST of China，2005，34（3）：296－299（in Chinese）. ［張智林，皮亦鳴，孫志堅.基于獨立分量分析的降噪技術 ［J］.電子科技大學學報，2005，34（3）：296－299.］

［13］PAN Song，HUANG Jiye.EDA technology and VHDL ［D］.Beijing：Tsinghua University Press，2007（in Chinese）.［潘松，黃繼業.EDA技術與VHDL［D］.北京：清華大學出版社，2007.］

［14］Signal processing society and the national science foundation.signal processing information base（SPIB）［EB/OL］.http：//www.spib.rice.edu/spib/select＿noise.html，1995.

［15］MA Jiancang，NIU Yilong，CHEN Haiyang.Blind singal processing.beijing：National defence industial press，2006：88－89（in Chinese）. ［馬建倉，牛奕龍，陳海洋.盲信號處理［D］.北京：國防工業出版社，2006：88－89.］