面向DSP平臺(tái)的CiSSA-CSP特征提取算法的移植與優(yōu)化

劉哲賢，趙金庫(kù)，趙玉峰，王 鵬

(1.清華大學(xué) 精密儀器系，北京 100084；2.黑龍江北方工具有限公司，黑龍江 牡丹江 157000)

0 引言

腦機(jī)接口(BCI，brain computer interfaces)是一種人與外部設(shè)備的交互系統(tǒng)，通過(guò)信號(hào)解析算法，解析人的大腦活動(dòng)，從而根據(jù)人的大腦活動(dòng)生成相應(yīng)的控制信號(hào)，幫助人類在不通過(guò)神經(jīng)與肌肉的情況下與外界實(shí)現(xiàn)交互[1-2]。BCI系統(tǒng)主要由四個(gè)部分組成，分別是信號(hào)監(jiān)測(cè)、特征提取、控制信號(hào)輸出和外部設(shè)備響應(yīng)與交互。

從信號(hào)監(jiān)測(cè)的手段來(lái)看，BCI系統(tǒng)主要分為兩類，分別是侵入式與非侵入式。侵入式腦機(jī)接口通過(guò)外科手術(shù)將微電極陣列植入大腦中，而非侵入式則將電極布置在頭皮上已檢測(cè)神經(jīng)活動(dòng)。相比于非侵入式，侵入式BCI檢測(cè)到的信號(hào)具有更高的信噪比，可用于病理分析、運(yùn)動(dòng)控制等領(lǐng)域；但侵入式BCI也面臨成本高、有醫(yī)療風(fēng)險(xiǎn)、使用不便等問(wèn)題，同時(shí)侵入式BCI的應(yīng)用場(chǎng)景主要局限于癱瘓病人的協(xié)助，而非侵入式則具有更廣闊的應(yīng)用前景，不僅能用于運(yùn)動(dòng)障礙人士的康復(fù)治療也能為健康人帶來(lái)不一樣的生活體驗(yàn)。非侵入式BCI中最常用的檢測(cè)手段是腦電圖(EEG，electroencephalography)[3]。

EEG信號(hào)為多通道數(shù)據(jù)，由緊貼在頭皮上的若干個(gè)電極采集得到，是記錄大腦皮層神經(jīng)元活動(dòng)的電生理信號(hào)[4]。EEG相比其他檢測(cè)手段成本更低，便攜性更好，但EEG的信號(hào)質(zhì)量是困擾研究人員的問(wèn)題[5]。由于神經(jīng)元的電生理信號(hào)需穿過(guò)顱骨和頭皮才能被位于頭皮表面的電極檢測(cè)到，因此EEG信號(hào)具有幅值低(微伏量級(jí))、信噪比低、極易受外界干擾、信號(hào)存在偽影[6]等特點(diǎn)。

從信號(hào)類型來(lái)看，腦電信號(hào)主要分為自發(fā)腦電和誘發(fā)腦電。其中誘發(fā)腦電指人接收到外界刺激(視覺(jué)、聽(tīng)覺(jué)等)的情況下大腦特定部位產(chǎn)生的具有一定特征的腦電信號(hào)。自發(fā)腦電則由人類自身心理活動(dòng)產(chǎn)生。

運(yùn)動(dòng)想象(MI，motor imagine)是一種典型的自發(fā)腦電信號(hào)。用戶想象肢體運(yùn)動(dòng)時(shí)，大腦皮層中會(huì)產(chǎn)生感覺(jué)運(yùn)動(dòng)節(jié)律(SMR，sensorimotor rhythms)。SMR包含mu節(jié)律和beta節(jié)律，分別位于8～13 Hz和13～30 Hz。感覺(jué)刺激、運(yùn)動(dòng)行為和心理意象可以改變皮層內(nèi)的功能連接，并導(dǎo)致mu節(jié)律和beta節(jié)律出現(xiàn)事件相關(guān)同步(ERS，event-related synchronization)和事件相關(guān)去同步(ERD，event-related desynchronization)現(xiàn)象[7]。運(yùn)動(dòng)想象任務(wù)主要包括左手運(yùn)動(dòng)想象、右手運(yùn)動(dòng)想象、以及想象腳和舌頭的運(yùn)動(dòng)。因此常用的運(yùn)動(dòng)想象解析算法為二分類或三分類算法，一般為對(duì)左右手想象的分類，對(duì)手和腳想象分類，以及對(duì)左右手和舌頭想象的分類。

典型的運(yùn)動(dòng)想象腦電(MI EEG)信號(hào)解析算法流程包括特征提取和分類。其中，由于MI EEG信號(hào)具有非平穩(wěn)性和低信噪比的性質(zhì)，特征提取具有一定挑戰(zhàn)性，同時(shí)也是研究熱點(diǎn)。最常見(jiàn)的算法是共空間模式算法(CSP，common spatial pattern)及其改進(jìn)方法。CSP是一種空間濾波方法，用于將EEG數(shù)據(jù)轉(zhuǎn)換到一個(gè)新空間中，其中一個(gè)類的方差最大化，而另一類的方差最小化。EEG信號(hào)的不同頻帶包含不同的信息，而CSP 能夠從特定頻帶中提取這些信息。研究人員嘗試了許多基于CSP算法的改進(jìn)方法。共空間頻譜模式(CSSP，common spatio-spectral patterns)將有限脈沖響應(yīng)濾波器(FIR，finite impulse response)集成到 CSP 算法中[8]。共稀疏頻譜空間模式(CSSSP，common sparse spatio-spectral patterns)是對(duì)CSSP的進(jìn)一步改進(jìn)，旨在找到不同通道間共同的頻譜模式，而非各通道獨(dú)立的頻譜模式[9]。子頻帶共空間模式(SBCSP，sub-band common spatial pattern)將腦電信號(hào)花費(fèi)為若干個(gè)子頻帶，再對(duì)不同的子帶進(jìn)行濾波并對(duì)每個(gè)子帶用CSP提取特征，然后使用線性判別分析方法來(lái)降低子帶的維數(shù)[10]。這些改進(jìn)方法相比于原始CSP算法，都能適當(dāng)提高分類準(zhǔn)確度。

循環(huán)奇異譜分析-共空間模式(CiSSA-CSP)方法可以提取非平穩(wěn)信號(hào)的最優(yōu)時(shí)頻空間特征，相比于其他特征提取方法，該方法可以顯著提高M(jìn)I EEG的分類精度[14]。CiSSA將信號(hào)分解為一系列頻率成分已知的重構(gòu)分量，然后對(duì)每一組分量分別用CSP提取特征。CiSSA可以去除相同頻率中的噪聲和偽影成分，因此CiSSA-CSP方法提取到的特征具有更好的分類效果。

當(dāng)前大多數(shù)腦機(jī)接口系統(tǒng)都是基于計(jì)算機(jī)實(shí)現(xiàn)的，相比較于計(jì)算機(jī)，嵌入式系統(tǒng)具有體積小、功耗低、成本低、響應(yīng)速度快和便攜式等特點(diǎn)，因此腦機(jī)接口系統(tǒng)的嵌入式實(shí)現(xiàn)也是腦機(jī)接口系統(tǒng)應(yīng)用化的必經(jīng)之路。Bueno等人利用16位數(shù)字信號(hào)控制器實(shí)現(xiàn)了運(yùn)動(dòng)想象三分類(左右手想象和心里想一句話)，分類準(zhǔn)確率約為60%[11]。使用SOMANN神經(jīng)網(wǎng)絡(luò)實(shí)現(xiàn)分類算法，從空間濾波到分類結(jié)果輸出的完整算法耗時(shí)2 ms。Ma等人出了一種基于FPGA的EEG-MI分類的純硬件系統(tǒng)[12]。EEG信號(hào)被處理為連續(xù)時(shí)域中的一系列多通道圖像，顯示受試者 MI 期間大腦皮層的能量變化。分類準(zhǔn)確度達(dá)到80.5%，系統(tǒng)耗時(shí)為計(jì)算機(jī)的，功耗為計(jì)算機(jī)的。Belwafi等人利用Stratix-IV(FPGA)實(shí)現(xiàn)了用于運(yùn)動(dòng)想象分類的嵌入式腦機(jī)接口系統(tǒng)[13]。使用動(dòng)態(tài)濾波算法對(duì)MI EEG信號(hào)進(jìn)行預(yù)處理，動(dòng)態(tài)濾波算法相比靜態(tài)濾波算法復(fù)雜度更低，更有利于FPGA嵌入式實(shí)現(xiàn)，能有效降低延時(shí)與功耗。動(dòng)態(tài)濾波算法通過(guò)選擇包含最活躍 ERS/ERS 模式的時(shí)間窗口，為每個(gè)用戶動(dòng)態(tài)選擇 EEG 試驗(yàn)的時(shí)期。其利用CSP進(jìn)行特征提取，LDA完成分類。相比于SVM、KNN等方法，LDA具有更低的算法復(fù)雜度。公共數(shù)據(jù)庫(kù)數(shù)據(jù)測(cè)試，離線訓(xùn)練時(shí)分類準(zhǔn)確率為77%，在線訓(xùn)練分類準(zhǔn)確度達(dá)到80.25%，單次實(shí)驗(yàn)計(jì)算耗時(shí)0.43 s(22通道，采樣時(shí)間5 s，采樣頻率100 Hz)。

嵌入式腦機(jī)接口系統(tǒng)面臨的挑戰(zhàn)包括高分類精度、穩(wěn)健性和自適應(yīng)性、功耗、時(shí)間響應(yīng)和成本效益。在運(yùn)動(dòng)想象分類任務(wù)中，分類準(zhǔn)確度應(yīng)當(dāng)達(dá)到80%以上，才算是可靠的腦機(jī)接口系統(tǒng)。系統(tǒng)應(yīng)當(dāng)具有較好的穩(wěn)健性和自適應(yīng)性，需要有良好的預(yù)處理算法以消除偽跡、眼電等噪聲，同時(shí)不同用戶的腦電信號(hào)特征具有明顯差異，系統(tǒng)應(yīng)當(dāng)能適用于不同用戶的信號(hào)特征。

嵌入式DSP處理器專門用于信號(hào)處理，通過(guò)對(duì)系統(tǒng)架構(gòu)和指令算法的特殊設(shè)計(jì)，實(shí)現(xiàn)高效率的編譯效率和數(shù)字運(yùn)算。相比于FPGA、ARM等其他嵌入式平臺(tái)，DSP具有更強(qiáng)的信號(hào)運(yùn)算能力，同時(shí)可使用C語(yǔ)言進(jìn)行編譯，易于實(shí)現(xiàn)。

本文通過(guò)優(yōu)化CiSSA-CSP算法流程，合理使用DSP庫(kù)函數(shù)和關(guān)鍵字等優(yōu)化工具，成功將CiSSA-DSP特征提取算法移植實(shí)現(xiàn)在DSP平臺(tái)中，并利用公共數(shù)據(jù)庫(kù)數(shù)據(jù)驗(yàn)證了其分類有效性和低延時(shí)的特性。

1 CiSSA-CSP特征提取算法原理

1.1 CiSSA原理

循環(huán)奇異譜分析(CiSSA，circulant singular spectrum analysis)是對(duì)奇異譜分析算法(SSA，singular spectrum analysis)的改進(jìn)，最初由Juan Bógalo等人[15]提出，其可以有效去除與有效信號(hào)頻率重疊的部分噪聲和偽影。CiSSA可以將信號(hào)分解為一系列頻率成分已知的重構(gòu)分量，其包含四個(gè)步驟：序列嵌入、分解、分組和對(duì)角線平均。

步驟1：序列嵌入

輸入為單通道腦電時(shí)序數(shù)據(jù)(x1，x2，…，…，xN)，使用窗口長(zhǎng)度為L(zhǎng)的滑動(dòng)窗口將輸入數(shù)據(jù)映射到軌跡矩陣X上。其中K=N-L+1，一般有L?K。

(1)

步驟2：分解

與SSA中的SVD分解所不同的是，CiSSA引入了循環(huán)矩陣CL，利用循環(huán)矩陣的特征向量將軌跡矩陣分解為一系列代表不同頻段的秩為1的基本矩陣。

(2)

相比于一般矩陣需要進(jìn)行SVD分解以計(jì)算特征向量和特征值，循環(huán)矩陣的特征向量和特征值計(jì)算相對(duì)簡(jiǎn)單，其第k個(gè)特征值及其對(duì)應(yīng)的特征向量表達(dá)式為：

(3)

(4)

(5)

循環(huán)矩陣CL的值由序列樣本協(xié)方差決定。

(6)

(7)

與SSA相同，改進(jìn)SSA算法中軌跡矩陣X被分解為L(zhǎng)個(gè)子矩陣相加。

X=X1+…+XL

(8)

Xi=wiwiHX，i=1，…，L

(9)

步驟3：分組

由于循環(huán)矩陣的特征向量為復(fù)向量，因此需要通過(guò)分組配對(duì)，將復(fù)矩陣Xi轉(zhuǎn)化為實(shí)矩陣。配對(duì)方法如下：

(10)

(11)

由此將原有軌跡矩陣X分解得到了一組新的子矩陣XBk。

步驟4：對(duì)角線平均

(12)

(13)

(14)

依上所述，CiSSA的算法流程如圖1所示。

圖1 CiSSA算法流程圖

1.2 CSP原理

共空間模式算法(CSP，common spatial pattern)算法是一種專用于二分類的多通道信號(hào)特征提取方法，CSP通過(guò)最大化兩類想象任務(wù)腦電信號(hào)方差的差異性，以提取出差異性較大的特征用于分類[16]。

首先計(jì)算各樣本的協(xié)方差陣并繼續(xù)歸一化，然后按照所屬任務(wù)類別，分別對(duì)兩類任務(wù)信號(hào)的歸一化協(xié)方差陣求均值，然后相加得到混合協(xié)方差陣。其中Xi為第i類樣本數(shù)據(jù)，大小為NC*NS的矩陣。

(15)

(16)

通過(guò)白化矩陣P，將不同方向的特征值完成歸一化。

PRPT=I

(17)

P=∑-1/2UT

(18)

R=U∑UT

(19)

UUT=I

(20)

此時(shí)兩類任務(wù)的協(xié)方差陣特征值對(duì)應(yīng)相加為1，因此某一類任務(wù)的方差大時(shí)，另一類任務(wù)的方差就小，接下來(lái)需要調(diào)整投影方向使得兩類方差的差異盡可能大[17]。方法是將每類任務(wù)的平均協(xié)方差陣白化，并沿特征向量方向進(jìn)行投影。其中特征向量的排序按照對(duì)應(yīng)特征值降序排布。

(21)

Si=BΛiBT，i=1，2

(22)

BBT=I

(23)

白化矩陣沿特征向量方向投影即得到空間濾波矩陣W?？臻g濾波矩陣W的每一行對(duì)應(yīng)一種空間投影模式，將W與樣本X相乘得到空間濾波后的信號(hào)，再對(duì)該信號(hào)的每一行求方差并取對(duì)數(shù)[18]，即得到此種投影方式下的特征。

W=BTP

(24)

Z=WX

(25)

Z=[Z1，…，ZNC]T

(26)

Fi=ln(var(Zi))

(27)

其中：Λi應(yīng)當(dāng)按照特征值降序排布，因此Z1，…，ZNC的順序與不同投影方向方差的差異性順序一致，Z1代表的投影方向最大化了第一類樣本投影后方差減第二類樣本投影后方差，ZNC代表的投影方向則最大化了投影后第二類樣本方差減第一類樣本方差。因此一般選取Z的前m行和后m行來(lái)計(jì)算用于分類的特征，這2m行是將類間方差差異最大化的2m個(gè)投影方向。

(28)

依上所述，CSP算法流程圖見(jiàn)圖2。

圖2 CSP算法流程圖

2 CiSSA-CSP特征提取算法的優(yōu)化移植

2.1 移植平臺(tái)和優(yōu)化方法介紹

C66x系列DSP是德州儀器公司近年來(lái)推出的高端DSP芯片，具有強(qiáng)大的定點(diǎn)浮點(diǎn)運(yùn)算能力，能用于圖像處理、音頻分析、信號(hào)處理等領(lǐng)域[20]。定點(diǎn)使用時(shí)，C66x核的乘加運(yùn)算能力(MAC，multiply accumulate)為C64x核的4倍。同時(shí)集成了浮點(diǎn)運(yùn)算能力和各核的原始運(yùn)算能力，達(dá)到了行業(yè)領(lǐng)先的44.8 GMACS和22.4 GFLOPS(在1.4 GHz工作頻率下)。單個(gè)周期內(nèi)可以執(zhí)行8次單精度浮點(diǎn)MAC操作，也可以執(zhí)行雙精度和混合精度操作。本文采用的是C66x系列的TMS320C6678DSP，C6678集成了8個(gè)C66x內(nèi)核，其主頻率最高可到1.4 GHz[19]。為方便實(shí)現(xiàn)，沒(méi)有自己設(shè)計(jì)接口電路，直接選用現(xiàn)有的開(kāi)發(fā)板以初步實(shí)現(xiàn)功能。本文選用的是廣州創(chuàng)龍公司的TL6678-EasyEVM開(kāi)發(fā)板。開(kāi)發(fā)板上已預(yù)先燒入了初始化程序，并引出GPIO、SPI等外設(shè)接口，使用開(kāi)發(fā)板可以大大方便我們的調(diào)試工作。

軟件開(kāi)發(fā)平臺(tái)采用軟件編譯平臺(tái)CCS(code composer studio)，版本為v5.5。CCS中集成了多個(gè)開(kāi)發(fā)套件，同時(shí)支持C語(yǔ)言、C++等語(yǔ)言，且提供了方便的調(diào)試環(huán)境，利于用戶調(diào)試代碼程序。同時(shí)還需要仿真器用于連接電腦和DSP開(kāi)發(fā)板。仿真器通過(guò)USB接口與電腦相連，通過(guò)JTAG接口與開(kāi)發(fā)板相連。通過(guò)仿真器才可以將CCS中程序?qū)懭隓SP并利用CCS控制調(diào)試DSP。選用的仿真器是創(chuàng)龍公司的TL-XDS200仿真器。

DSP軟件優(yōu)化主要包含兩個(gè)方面，算法和編譯。算法優(yōu)化指重新設(shè)計(jì)算法流程，通過(guò)數(shù)學(xué)手段降低計(jì)算復(fù)雜度。編譯優(yōu)化指優(yōu)化代碼實(shí)現(xiàn)方式，通過(guò)合理編譯，減少計(jì)算耗時(shí)。CCS中有編譯器優(yōu)化選項(xiàng)，分別為off，0，1，2，3五檔，off代表無(wú)優(yōu)化，0代表寄存器級(jí)別優(yōu)化，1代表局部?jī)?yōu)化，2代表全局優(yōu)化，3代表工程文件整體優(yōu)化[21]。我們選用～2級(jí)別優(yōu)化，其可以實(shí)現(xiàn)高效率的編譯同時(shí)也方便進(jìn)行程序調(diào)試。C語(yǔ)言優(yōu)化也可以提高編譯效率，其主要從兩個(gè)角度著手，其一是使用各類關(guān)鍵字或編譯指令告知編譯器更多信息，讓編譯器能更好的對(duì)C語(yǔ)言進(jìn)行編譯；其二是使用專用于相應(yīng)型號(hào)DSP的C語(yǔ)言函數(shù)，這些函數(shù)由德州儀器公司編寫提供，已專門針對(duì)相應(yīng)DSP做了特定優(yōu)化[22]。

2.2 CiSSA算法優(yōu)化

MI EEG信號(hào)的特征頻段集中在mu節(jié)律和beta節(jié)律，同時(shí)我們并不需要各分量的功率譜信息，因此，我們無(wú)需計(jì)算循環(huán)矩陣的特征值。由式(3)～(5)可知，循環(huán)矩陣的特征向量只與L有關(guān)，而與循環(huán)矩陣中元素?cái)?shù)值無(wú)關(guān)，因此在DSP上實(shí)現(xiàn)改進(jìn)SSA實(shí)現(xiàn)時(shí)，無(wú)需計(jì)算cm，求出特征向量后，利用式(9)，即可求出子矩陣。但由于特征向量為復(fù)向量，因此子矩陣Xi為復(fù)矩陣，這在DSP硬件上實(shí)現(xiàn)并不劃算。因此實(shí)際移植到DSP上時(shí)，希望跳過(guò)復(fù)矩陣的計(jì)算步驟，直接得到配對(duì)后的實(shí)矩陣。利用式(5)、(11)和(12)，計(jì)算得到XBk。XBk的公式如下：

(29)

(30)

(31)

(32)

由式(29)、(31)計(jì)算得到XBk后，對(duì)角平均化即可得到原時(shí)間序列的相應(yīng)分量。由上述公式推得的改進(jìn)SSA算法流程見(jiàn)圖3。相比于圖1，此處算法流程省去了計(jì)算循環(huán)矩陣和復(fù)特征向量的步驟，得到了顯著簡(jiǎn)化，更適于移植到DSP上。

圖3 CiSSA優(yōu)化后算法流程圖

2.3 CiSSA算法移植與復(fù)雜度分析

是MI EEG信號(hào)特征集中在mu節(jié)律(8～13 Hz)和beta節(jié)律(13～30 Hz)上，因此在利用改進(jìn)SSA算法處理運(yùn)動(dòng)想象腦電信號(hào)時(shí)，只需要提取出8～30 Hz范圍內(nèi)信號(hào)分量，其他頻段的信號(hào)都不需要。在DSP實(shí)現(xiàn)時(shí)，應(yīng)當(dāng)首先篩選出需要的分量，然后只用計(jì)算相應(yīng)的分量即可，不用計(jì)算所有的分量，以降低運(yùn)算量。設(shè)所需分量k的取值集合記為Range，由式(13)推得Range計(jì)算公式如下：

(33)

由式(31)注意到VBk的取值只與k和L有關(guān)，與待處理信號(hào)的值無(wú)關(guān)。因此在DSP上改進(jìn)SSA算法時(shí)，首先根據(jù)上式計(jì)算出Range，再計(jì)算需要用到的VBk并保存。這樣在處理不同信號(hào)樣本時(shí)，VBk可以直接調(diào)用已有的數(shù)據(jù)，不用重復(fù)計(jì)算，以節(jié)省計(jì)算時(shí)間。

在實(shí)際應(yīng)用中，有兩個(gè)階段：通過(guò)大量樣本訓(xùn)練模型階段和單次實(shí)驗(yàn)信號(hào)解析階段。訓(xùn)練模型階段指利用大量已有信號(hào)樣本訓(xùn)練得到算法中的各類模型參數(shù)。

圖4描述的是移植到DSP上的用于腦電信號(hào)預(yù)處理的改進(jìn)SSA算法流程，用于對(duì)單個(gè)樣本信號(hào)進(jìn)行預(yù)處理。不過(guò)不論是單次實(shí)驗(yàn)信號(hào)解析還是大量樣本訓(xùn)練分類模型，預(yù)處理都是必不可少的一步。因此在單次實(shí)驗(yàn)信號(hào)解析與訓(xùn)練分類模型過(guò)程中，預(yù)處理都按照?qǐng)D4的流程進(jìn)行。

圖4 CiSSA移植后算法流程圖

假設(shè)運(yùn)動(dòng)想象腦電信號(hào)中包含Nc個(gè)通道數(shù)據(jù)，采樣頻率為fsHz，單次實(shí)驗(yàn)采樣時(shí)間為T秒，則單次實(shí)驗(yàn)信號(hào)為Nc×Ns的數(shù)據(jù)矩陣。假設(shè)Range的元素個(gè)數(shù)為numK。

Ns=fs×T

(34)

Range的算法復(fù)雜度可忽略不計(jì)，計(jì)算numK個(gè)Vbk的復(fù)雜度為O(numK×L2)；嵌入軌跡矩陣的復(fù)雜度為O(L(Ns-L+1))=O(LNs)；計(jì)算numK個(gè)XBk的復(fù)雜度為O(numK×L2(Ns-L+1))=O(numK×L2Ns)，對(duì)角平均化的復(fù)雜度為O(numK×LNs)。

通過(guò)上述分析可知，移植到DSP上的改進(jìn)SSA算法中，計(jì)算量最大的部分就是XBk的計(jì)算，這一步是單純的矩陣乘法，共需要L2(Ns-L+1)次浮點(diǎn)乘加運(yùn)算，而在單次實(shí)驗(yàn)信號(hào)處理中，對(duì)單通道下Range中每一個(gè)k要計(jì)算一次，因此其一共需要L2(Ns-L+1)Nc×numK次浮點(diǎn)乘加運(yùn)算。

2.4 CSP算法移植與復(fù)雜度分析

為提高運(yùn)行效率，在移植過(guò)程中應(yīng)當(dāng)多使用DSP的庫(kù)函數(shù)，這些函數(shù)都專門針對(duì)DSP結(jié)構(gòu)進(jìn)行了優(yōu)化，執(zhí)行效率很高。CSP算法中有兩步特征分解，而德州儀器公司所編寫的dsplib庫(kù)(3.4.0.0版本)中提供了單精度數(shù)值矩陣SVD分解函數(shù)，函數(shù)名為DSPF_sp_svd。而對(duì)于方陣而言，SVD分解就是特征分解，因此移植時(shí)使用該函數(shù)來(lái)完成特征分解的步驟。

如2.3節(jié)所述，實(shí)際應(yīng)用中分為模型訓(xùn)練和單次實(shí)驗(yàn)信號(hào)解析兩個(gè)階段，模型訓(xùn)練階段算法與2.2節(jié)所述一致。移植后算法流程見(jiàn)圖5。而在單次實(shí)驗(yàn)信號(hào)階段，已有空間濾波矩陣W，此時(shí)只需將W的前m行和后m行與樣本相乘，并對(duì)結(jié)果的每一行方差取對(duì)數(shù)即可完成特征提取，移植后算法流程見(jiàn)圖6。

圖5 CSP移植后模型訓(xùn)練階段流程圖

圖6 CSP移植后單次實(shí)驗(yàn)信號(hào)解析階段流程圖

模型訓(xùn)練階段中，對(duì)nT個(gè)樣本計(jì)算歸一化協(xié)方差矩陣的復(fù)雜度為O(nT×Nc2Ns)，特征分解和矩陣乘法的復(fù)雜度為O(Nc3)，其余步驟算法復(fù)雜度可忽略不計(jì)。在MI EEG信號(hào)中，一般有Nc?Ns，因此由上述分析可知，CSP整體算法復(fù)雜度為O(nT×Nc2Ns)。而在歸一化協(xié)方差矩陣計(jì)算中，矩陣乘法XXT需要Nc2Ns次浮點(diǎn)乘加運(yùn)算，trace(XXT)需要Nc次浮點(diǎn)加法運(yùn)算，歸一化需要Nc2次浮點(diǎn)除法運(yùn)算。由于浮點(diǎn)除法運(yùn)算相比于乘法更消耗時(shí)間，因此移植到DSP時(shí)，為節(jié)省時(shí)間，可以將除法運(yùn)算轉(zhuǎn)化為乘法運(yùn)算，即先求出trace(XXT)的倒數(shù)，之后每次與該倒數(shù)相乘即可。這樣只用計(jì)算1次除法和Nc2次乘法。單次實(shí)驗(yàn)信號(hào)解析階段的算法復(fù)雜度為O(2m×Nc×Ns)，相比于2.3節(jié)所述CiSSA復(fù)雜度可忽略。

2.5 CiSSA-CSP整體算法移植與內(nèi)存分配

模型訓(xùn)練階段的CiSSA-CSP移植后整體算法流程如圖7所示。注意到樣本數(shù)據(jù)是存儲(chǔ)在MSMCRAM中的，經(jīng)過(guò)預(yù)處理后，再存儲(chǔ)到DDR3上，這是因?yàn)镃SP需要nT個(gè)樣本的數(shù)據(jù)，而MSMCRAM可能無(wú)法提供足夠的存儲(chǔ)空間。以公開(kāi)數(shù)據(jù)庫(kù)數(shù)據(jù)為例，nT=280，Nc=17，Ns=350，假設(shè)Range={3，4，5，6，7，8}，即numK=6，于是280個(gè)樣本經(jīng)過(guò)CiSSA預(yù)處理后，有numK×nT×Nc×Ns=9 996 000個(gè)數(shù)需要存儲(chǔ)。假設(shè)在DSP中把數(shù)據(jù)存為單精度格式，一個(gè)單精度數(shù)需要4B存儲(chǔ)空間，則結(jié)果預(yù)處理得到的所有數(shù)據(jù)至少需要39 047 MB的存儲(chǔ)空間，而在C6678上MSMCRAM只有4 096 MB的空間，而DDR3有1 GB空間，可以滿足要求。因此需要將經(jīng)過(guò)CiSSA預(yù)處理的數(shù)據(jù)存儲(chǔ)到DDR3上。

圖7 模型訓(xùn)練階段CiSSA-CSP移植后算法流程圖

單次實(shí)驗(yàn)信號(hào)解析階段中已有分類參數(shù)，CSP模塊算法得到極大簡(jiǎn)化，而CiSSA模塊的算法不變。此時(shí)整體算法流程圖見(jiàn)圖8。CiSSA模塊中無(wú)需計(jì)算Range，Vbk，其已經(jīng)計(jì)算好并存儲(chǔ)在內(nèi)存中了，同樣的CSP模塊中各分量的空間濾波矩陣Wk也已儲(chǔ)存在內(nèi)存中，需要時(shí)直接調(diào)用即可。單次實(shí)驗(yàn)數(shù)據(jù)處理中則無(wú)需使用DDR3存儲(chǔ)空間，所有內(nèi)存均存放在MSMCRAM上。

圖8 單次實(shí)驗(yàn)信號(hào)解析階段CiSSA-CSP移植后算法流程圖

3 移植測(cè)試實(shí)驗(yàn)結(jié)果

C語(yǔ)言算法通過(guò)CCS平臺(tái)編譯并移植到C6678DSP開(kāi)發(fā)板上，本節(jié)將對(duì)移植到DSP上的算法進(jìn)行測(cè)試，一方面與在電腦上使用Matlab計(jì)算的結(jié)果對(duì)比以驗(yàn)證計(jì)算準(zhǔn)確性，另一方面測(cè)試使用DSP完成計(jì)算所消耗的時(shí)間。其中PC端處理器為Intel i5-8250 CPU，機(jī)帶RAM為8 GB，Matlab版本為2021b；C6678開(kāi)發(fā)板上僅使用核0進(jìn)行運(yùn)算，另外7核未啟動(dòng)，芯片主頻設(shè)定為1.0 GHz，信號(hào)存儲(chǔ)及運(yùn)算均在32位浮點(diǎn)數(shù)精度下實(shí)現(xiàn)。

3.1 數(shù)據(jù)來(lái)源

用于測(cè)試的數(shù)據(jù)來(lái)源于BCI Competition Ⅲ，選用的數(shù)據(jù)集是data set ⅠVa，其中包含5個(gè)健康受試者的數(shù)據(jù)[23]。受試者進(jìn)行右手和右腳的運(yùn)動(dòng)想象實(shí)驗(yàn)，使用國(guó)際通用的10～20系統(tǒng)電極排布方式，共有118個(gè)通道[15]。每個(gè)受試者進(jìn)行280次實(shí)驗(yàn)，每個(gè)想象任務(wù)140次，單次實(shí)驗(yàn)為持續(xù)3.5秒的運(yùn)動(dòng)想象，兩次想象間隔為1.75～2.25秒，采樣頻率為100 Hz。處理數(shù)據(jù)時(shí)選擇所有通道中與運(yùn)動(dòng)想象相關(guān)的17個(gè)通道(FC3，F(xiàn)C1，F(xiàn)Cz，F(xiàn)C2，F(xiàn)C4，C5，C3，C1，Cz，C2，C4，C6，CP3，CP1，CPz，CP2，CP4，如圖9所示)[24]，因此單次實(shí)驗(yàn)樣本數(shù)據(jù)為17*350的矩陣。

圖9 公共數(shù)據(jù)庫(kù)信號(hào)通道分布圖

3.2 CiSSA-CSP特征提取算法移植誤差測(cè)試結(jié)果

取L=24，可得Range=(3，4，5，6，7，8)，numK=6。于是單個(gè)信號(hào)樣本經(jīng)過(guò)CiSSA算法預(yù)處理后，得到6個(gè)與原信號(hào)大小一致的矩陣。取m=2，即一個(gè)信號(hào)樣本矩陣通過(guò)CSP特征提取后，得到2m=4個(gè)特征值。于是原信號(hào)一共提取出24個(gè)特征。對(duì)比DSP上運(yùn)行CiSSA-CSP算法得到的特征與PC機(jī)上MATLAB計(jì)算的特征，結(jié)果見(jiàn)表1，平均相對(duì)誤差在10-4量級(jí)。其中絕對(duì)誤差為DSP運(yùn)算得到特征減去PC機(jī)上MATLAB計(jì)算得到特征的絕對(duì)值，相對(duì)誤差為DSP運(yùn)算得到特征減去PC機(jī)上MATLAB計(jì)算得到特征的絕對(duì)值再除以PC機(jī)上MATLAB計(jì)算得到特征。

表1 CiSSA-CSP誤差測(cè)試

注意到CiSSA-CSP算法步驟中，就計(jì)算類型而言主要分為兩種：矩陣乘法和特征分解。矩陣乘法采用DSP庫(kù)函數(shù)實(shí)現(xiàn)，誤差相對(duì)較小，我們用生成的隨機(jī)方陣檢驗(yàn)了特征分解的誤差，結(jié)果見(jiàn)表2，其中方陣大小為N×N。單次特征分解的相對(duì)誤差在10-5量級(jí)，而CiSSA-CSP運(yùn)算中有兩步特征分解，兩步誤差累積可達(dá)10-4量級(jí)，這與表1的結(jié)果相吻合，說(shuō)明移植后CiSSA-CSP特征提取算法的主要誤差來(lái)源為特征分解。

表2 特征分解誤差測(cè)試

3.3 CiSSA-CSP特征提取算法移植耗時(shí)測(cè)試結(jié)果

表3測(cè)試了CiSSA運(yùn)算耗時(shí)與式(29)單步運(yùn)算耗時(shí)。在單次實(shí)驗(yàn)信號(hào)解析中，CiSSA運(yùn)算中需要計(jì)算17×6=102次式(29)，由表3結(jié)果可得式(29)的計(jì)算耗時(shí)占比達(dá)到預(yù)處理模塊計(jì)算耗時(shí)的94.49%。這與我們?cè)?.3節(jié)中復(fù)雜度分析的結(jié)果相吻合。

表3 CiSSA耗時(shí)測(cè)試 ms

矩陣乘法函數(shù)是CiSSA-CSP運(yùn)算過(guò)程中的核心步驟，我們使用關(guān)鍵字、編譯器優(yōu)化等方式，對(duì)乘法函數(shù)文件進(jìn)行了優(yōu)化。表3中也給出了優(yōu)化后CiSSA的耗時(shí)測(cè)試結(jié)果，優(yōu)化后耗時(shí)僅為優(yōu)化前耗時(shí)的14%。

訓(xùn)練模型階段，CSP運(yùn)算的耗時(shí)測(cè)試結(jié)果見(jiàn)表4，其中單次協(xié)方差計(jì)算為式(15)。此時(shí)特征提取模塊算法中，共包含280次協(xié)方差計(jì)算(共280個(gè)樣本)，所有樣本的協(xié)方差計(jì)算耗時(shí)在特征提取模塊中耗時(shí)占比達(dá)94.03%。在訓(xùn)練模型階段，特征提取模塊對(duì)時(shí)間消耗并不敏感，對(duì)耗時(shí)敏感的是單次實(shí)驗(yàn)信號(hào)解析。而在單次實(shí)驗(yàn)信號(hào)解析階段中，其耗時(shí)測(cè)試結(jié)果見(jiàn)表5。與CiSSA耗時(shí)相比，單次實(shí)驗(yàn)信號(hào)解析中CSP的耗時(shí)可以忽略，這與2.4節(jié)分析的結(jié)果一致。

表4 CSP耗時(shí)測(cè)試(模型訓(xùn)練階段) ms

表5 CSP耗時(shí)測(cè)試(單次實(shí)驗(yàn)信號(hào)解析階段) ms

3.4 移植后分類測(cè)試結(jié)果

CiSSA-CSP特征提取算法移植到DSP后的有效性需要依賴MI EEG信號(hào)分類精度來(lái)確認(rèn)，因此我們將在DSP上還移植實(shí)現(xiàn)了LDA分類算法，并對(duì)完整算法流程進(jìn)行了測(cè)試，比較相同樣本在DSP上運(yùn)算的分類精度和MATLAB上運(yùn)算的分類精度。分類精度測(cè)試采用10折交叉驗(yàn)證，結(jié)果見(jiàn)表6。DSP運(yùn)算的分類精度相比于PC機(jī)上Matlab計(jì)算結(jié)果下降小于0.5%，且平均分類精度達(dá)到93.2%。同時(shí)LDA的耗時(shí)相比CiSSA-CSP可以忽略不計(jì)，而由3.3節(jié)結(jié)果可知，單次實(shí)驗(yàn)信號(hào)解析階段中CiSSA-CSP在DSP運(yùn)算的耗時(shí)小于150 ms，滿足MI EEG分類實(shí)時(shí)性的要求。

表6 DSP運(yùn)算與PC機(jī)上Matlab運(yùn)算分類精度對(duì)比

4 結(jié)束語(yǔ)

本文通過(guò)優(yōu)化CiSSA算法、合理規(guī)劃整體算法流程、使用CCS編譯器優(yōu)化工具等方法成功在C6678DSP開(kāi)發(fā)板中移植實(shí)現(xiàn)了CiSSA-CSP特征提取算法，并使用MI EEG公共數(shù)據(jù)庫(kù)數(shù)據(jù)驗(yàn)證了算法的準(zhǔn)確性，分類準(zhǔn)確度相比于PC機(jī)上Matlab運(yùn)算結(jié)果下降小于0.5%。同時(shí)單次實(shí)驗(yàn)信號(hào)解析耗時(shí)少于0.15 s，滿足了解析實(shí)時(shí)性的要求。