葉萬洋,郎榮玲
(北京航空航天大學(xué)電子信息工程學(xué)院,北京 100191)
在復(fù)雜的電磁環(huán)境中,存在多種異質(zhì)隨機信號共存的情況,無法直接對機會信號進(jìn)行識別、評估以及定位。為了解決這些問題,需要對接收端的多源信號進(jìn)行分離。一般情況下,存在對機會信號的先驗知識未知的情況。因此,需要借助盲源分離(Blind Source Separation,BSS)技術(shù)解決這個問題。
BSS起源于著名的雞尾酒會問題[1],其試圖在未知信號源特性與通道傳輸特性的前提下,僅根據(jù)接收到的信號來恢復(fù)原始信號。盲代表我們對信號的形式并不關(guān)心,而僅利用預(yù)先假設(shè)信號所具有的某一特性來建立目標(biāo)函數(shù),從而實現(xiàn)混合信號的分離處理。獨立成分分析(Independent Component Analysis, ICA)是一種常用的BSS技術(shù)[2],ICA的核心是尋找一個可以最大化混合信號各分量相互獨立性的線性變換,從混合信號中分離出各個獨立的源信號。基于ICA技術(shù)的BSS算法包括ICA算法[3]、快速獨立成分分析(Fast-ICA)算法[4]和等變自適應(yīng)盲分離(Equivariant Adaptive Source Separation via Independence,EASI)算法[5]。目前,基于ICA技術(shù)的BSS算法在語音信號、生物醫(yī)學(xué)信號、水聲信號、移動通信、地震信號處理等領(lǐng)域中均得到了廣泛應(yīng)用[6-7]。但是大多數(shù)BSS算法都使用計算機或服務(wù)器執(zhí)行離線分析,并沒有應(yīng)用于信號的實時處理中。由于BSS算法的計算中涉及大量矩陣代數(shù)運算,因此非常適合在大規(guī)模并行處理的現(xiàn)場可編程門陣列(Field-Programmable Gate Array,F(xiàn)PGA)和專用集成電路(Application Specific Integrated Circuit,ASIC)芯片上實現(xiàn),并行的計算方式和大量的邏輯單元能夠保證BSS算法的實時實現(xiàn)。
目前已經(jīng)提出了幾種BSS算法的ASIC或FPGA實現(xiàn)[8-9],以加速BSS處理。但是這些研究成果大多都是針對實數(shù)域的聲音信號或胎兒心電信號,并不能夠直接用于復(fù)數(shù)域混合信號的BSS。Van等提出了一種面向ASIC芯片內(nèi)存不足的在線Fast-ICA算法和用于八通道腦電圖(Electroencephalography,EEG)信號分離的硬件體系結(jié)構(gòu)[10]。Charoensak等提出了基于FPGA的ICA算法[11],該算法能夠?qū)崟r地分離兩種聲音的混合信號,但是只能夠處理低頻的聲音信號。Kim等提出了用于BSS和自適應(yīng)噪聲消除(Adaptive Noise Cancelling,ANC)的ICA算法的FPGA實現(xiàn)[12],在嘈雜環(huán)境中能夠分離出噪聲信號從而增強語言信號,但是該算法在實際應(yīng)用中會存在一定的延時。Shyu等提出了用于實時信號BSS處理的Fast-ICA算法的FPGA實現(xiàn)方案[13],設(shè)計采用了流水線的架構(gòu)和浮點運算處理來提高工作性能。ODOM在文獻(xiàn)[14]中對比了FPGA實現(xiàn)ICA、Fast-ICA和EASI三種BSS算法的性能和所需資源的數(shù)量,其中EASI算法由于具有收斂速度快、算法簡單、有效性高等優(yōu)點,被認(rèn)為是最佳算法。
在本文中,為了能夠?qū)崟r地分離復(fù)數(shù)域的多源混合機會信號,設(shè)計了一套基于軟件無線電的機會信號實時盲分離系統(tǒng)。系統(tǒng)軟件部分利用C++和QT編寫的信號采集程序,能夠在上位機中設(shè)置采集信號的頻帶范圍、下變頻過程中的本振頻率大小和信號上傳的通道選擇,并且能夠?qū)崟r存儲和顯示分離信號的時域波形和頻域波形。硬件部分的主要功能是在FPGA中實時實現(xiàn)兩通道的復(fù)數(shù)域EASI算法。算法采用流水線架構(gòu)和浮點運算處理設(shè)計,其中流水線架構(gòu)能夠提高算法的實時處理能力;相比于數(shù)字信號處理中的定點設(shè)計,浮點算術(shù)設(shè)計能夠提供更高的精度和動態(tài)性能。在實際實驗中,該系統(tǒng)能夠?qū)崟r地分離兩路復(fù)數(shù)域混合機會信號。
盲源分離的系統(tǒng)模型如圖1所示。
圖1 盲源分離系統(tǒng)模型Fig.1 Model of blind source separation system
假設(shè)天線陣列由M個天線單元組成,則瞬時線性混合模型可以表示為
x=As+N
(1)
其中,x表示接收到的混合信號,s是源信號,A是滿秩的混合矩陣,N是加性高斯白噪聲。BSS的核心思想就是在s和A未知的情況下,僅根據(jù)x來估計一個分離矩陣W,使得x與W相乘后的分離信號y盡可能地接近s,可以表示為
y=Wx=WAs+WN
(2)
在盲源分離過程中,由于僅利用源信號之間相互統(tǒng)計獨立這一條性質(zhì),并不知道s和A的先驗信息。因此根據(jù)文獻(xiàn)[2]可知,分離信號y和源信號s之間存在著排列順序和幅度不確定性,即有
C=WA=PΛ
(3)
其中,C為廣義置換矩陣(Generalized Permutation Matrix);P為置換矩陣(即初等矩陣),表示了分離信號各分量排列順序的不確定性;Λ為滿秩對角矩陣,表示了分離信號各分量幅度的不確定性。由于分離信號各分量的排列順序和幅度的不確定性并不影響信號的波形,因此可以認(rèn)為源信號能夠被完全恢復(fù)。
EASI算法是將白化過程和高階去相關(guān)過程同時進(jìn)行的一種等變化的自適應(yīng)算法[5]。其核心思想是通過自適應(yīng)地更新W,使得y的各分量盡可能地相互獨立,即各個分量的互信息代價函數(shù)最小。EASI算法中W的迭代更新方向為自然梯度方向,即代價函數(shù)的最陡下降方向。
典型的代價函數(shù)的形式為
φ(W)=E(f(y))
(4)
但其必須在去相關(guān)約束下進(jìn)行優(yōu)化
Ry=E(yyT)=I
(5)
其中,I表示單位矩陣。
基于此可將W分解為
W=UV
(6)
其中,V是n×m的白化矩陣,U是n×n的正交矩陣。
白化信號的形式為
z=Vx
(7)
分離信號的形式為
y=Uz
(8)
因此,將W的更新過程分為兩步,首先獲取V和U的串行更新,然后將其組合成W的串行更新。
1)白化矩陣的串行更新
由于BSS算法不能恢復(fù)源信號的真實幅度,因此可以假設(shè)源信號具有單位方差和零協(xié)方差,即
Rs=E[ssT]=I
(9)
因為U是正交矩陣,y=Uz,為了滿足去相關(guān)約束條件:z需要具有單位方差和零協(xié)方差,即
Rz=E[zzT]=E[VAssTATVT]
=(VA)E[ssT](VA)T=I
(10)
其中,VA是正交矩陣。V可以通過最小化Rz和I之間的距離獲得,即
(11)
φ1(V)的相對梯度為
(12)
因此,白化矩陣串行迭代更新算法為
V(t+1)=V(t)+μ(I-zzT)V(t)
(13)
其中,μ是迭代步長。
2)正交矩陣的串行更新
U的更新是利用自然梯度法對U的互信息代價函數(shù)進(jìn)行優(yōu)化
φ2(U)=-lg|det(U)|-E{lg[f(y)]}
(14)
得到φ2(U)的自然梯度,即
(15)
U是正交矩陣,則有U-1=UT,即
(16)
因此,正交矩陣的串行迭代更新算法為
U(t+1)=U(t)+μ{I-E[φ2(y)yT]}U(t)
(17)
將式(17)改寫成U(t+1)=U(t)+ΔUU(t)的形式,由于U(t+1)也需要保留正交性,即
(18)
(U(t)+ΔUU(t))(U(t)+ΔUU(t))T=I+o(ΔU)
(19)
U(t+1)=U(t)+u{E[y(t)φ2(y(t))T-
φ2(y(t))y(t)T]}U(t)
(20)
W的全局更新算法由式(17)和式(20)組合得到,即
W(t+1)=U(t+1)V(t+1)
=(I+μ[I+y(t)φ2(y(t))T-
y(t)y(t)T-φ2(y(t))y(t)T])W(t)
=W(t)+μH(t)W(t)
(21)
其中,H表示相對梯度矩陣,非線性函數(shù)φ2(y)的計算公式通常選擇為φ2(y)=y3。
本文使用Xilinxx公司提供的EDA工具ISE手動編寫Verilog代碼,并在FPGA中實時實現(xiàn)了兩通道的復(fù)數(shù)域EASI算法。根據(jù)式(21)將EASI算法分為3個模塊,y計算模塊,H計算模塊,W更新模塊。其中y、W和H都是復(fù)數(shù)域矩陣,表示形式如下
(22)
其中,i表示實部,q表示虛部,j表示虛數(shù)單位。
EASI算法采用自適應(yīng)的更新方法,即每一時刻的輸入x都會更新一次W,并且每一個更新過程中都利用了上一個時刻x更新得到的W。然而在FPGA中更新一次W往往需要多個時鐘周期,并不能在1個時鐘周期內(nèi)完成。為了在FPGA設(shè)計中滿足自適應(yīng)算法的要求,采用流水線架構(gòu)來保證1個時鐘周期更新一次W。流水線的原理如圖2所示,將EASI算法分成若干個時間上均衡的子模塊,從流水線的起點連續(xù)地輸入數(shù)據(jù),流水線的各子模塊以重疊方式執(zhí)行。這種處理方式能夠使得W的更新速度只與x的輸入速度有關(guān),而與處理所需的時間無關(guān)。采用流水線的架構(gòu)設(shè)計雖然能夠保證每一時刻更新一次W,但是W的更新過程仍然需要多個時鐘周期。為了解決這個問題,將迭代步長μ的值設(shè)置為0.0004,這樣在一次W的更新過程內(nèi),W的值基本不變,因此可以用當(dāng)前時刻的W來代替多個時鐘周期后更新的W,從而保證了算法中整個W的更新過程的正確性和穩(wěn)定性。
圖2 EASI算法流水線處理Fig.2 EASI algorithm pipeline processing
由于將迭代步長μ的值設(shè)置的很小,對每一次更新過程中的數(shù)據(jù)的精度要求很高,因此在本文中整個算法的運算處理都采用浮點運算方式來保證數(shù)據(jù)精度。浮點運算采用ISE中的浮點IP核進(jìn)行處理,包括浮點加法器、浮點減法器和浮點乘法器。除此以外,由于AD量化后的數(shù)字信號是定點數(shù)據(jù),因此需要使用定點轉(zhuǎn)浮點IP核將輸入的定點數(shù)轉(zhuǎn)化為浮點數(shù);由EASI模塊輸出的到FIFO中緩存的數(shù)據(jù)是定點數(shù),而實際運算的結(jié)果是浮點數(shù),因此還需要使用浮點轉(zhuǎn)定點IP核將輸出信號的浮點數(shù)轉(zhuǎn)化為定點數(shù)上傳。EASI算法的硬件整體結(jié)構(gòu)如圖3所示。
圖3 EASI算法硬件實現(xiàn)結(jié)構(gòu)圖Fig.3 EASI algorithm hardware implementation structure diagram
為了節(jié)省FPGA的邏輯元素,整個算法采用了時序硬件設(shè)計,其余計算將按順序共享浮點運算單元。下面將詳細(xì)描述EASI算法的3個主要模塊
(1)y計算模塊
此模塊使用W和x計算復(fù)數(shù)y1和y2,運算流程如圖3所示。W的初始值是單位矩陣,將式(2)改寫成復(fù)數(shù)域形式如下
y1i=W11i*x1i-W11q*x1q+W12i*x2i-
W12q*x2q
y1q=W11i*x1q+W11q*x1i+W12i*x2q+
W12q*x2i
y2i=W21i*x1i-W21q*x1q+W22i*x2i-
W22q*x2q
y2q=W21i*x1q+W21q*x1i+W22i*x2q+
W22q*x2i
(23)
在FPGA實現(xiàn)中,此模塊需要3個IP核時鐘周期,其中y1的運算流程如圖4所示。輸出結(jié)果y分別是H計算模塊的輸入和EASI算法的輸出。
圖4 復(fù)數(shù)域矩陣運算硬件流程圖Fig.4 Hardware flowchart of complex matrix operation
(2)H計算模塊
根據(jù)式(22)可知,H函數(shù)的計算公式如下
H(t)=I+y(t)φ2(y(t))T-y(t)y(t)T-
φ2(y(t))y(t)T
(24)
由于式(24)中包含非線性函數(shù)φ2(y)=y3,因此在硬件中直接求解相對梯度矩陣H會消耗大量的邏輯資源。為了降低H的計算復(fù)雜度,從矩陣運算的角度將式(22)中y的復(fù)數(shù)域表示形式代入式(24)中,推導(dǎo)并簡化相對梯度矩陣H各分量的求解公式,結(jié)果如下
(25)
其中,由于H11q和H22q的計算結(jié)果始終為0,因此在實際的硬件處理過程中不需要考慮這2個參數(shù)以節(jié)省硬件資源。此模塊的輸出H是W更新模塊的輸入。
(3)W更新模塊
根據(jù)式(21)利用上一個模塊得到H和設(shè)定的迭代步長μ來更新W。首先計算中間矩陣M=I+μH,由于H11q和H22q的計算結(jié)果始終為0,因此M11q和M22q在整個算法的運算過程中也始終為0,這樣在硬件計算中同樣可以不用考慮這4個參數(shù),以節(jié)省硬件資源。則W的更新公式有
W11i←M11i*W11i+M12i*W21i-M12q*W21q
W11q←M11i*W11q+M12i*W21q+M12q*W21i
W12i←M11i*W12i+M12i*W22i-M12q*W22q
W12q←M11i*W12q+M12i*W22q+M12q*W22i
(26)
W21i、W21q、W22i、W22q的計算公式和式(26)類似,這里不再贅述。更新后的W使用寄存器寄存后用于下一時刻y的計算。
機會信號實時盲分離系統(tǒng)采用軟件無線電的設(shè)計結(jié)構(gòu),系統(tǒng)總體結(jié)構(gòu)如圖5所示。軟件部分在上位機中利用C++和QT編寫完成,硬件部分由接收天線、AD9361、FPGA組成。其中FPGA芯片型號為Xc7k325t,AD9361芯片的采樣率為40MHz, AD9361芯片前端的混頻模塊能夠?qū)π盘栕鱿伦冾l處理,并將中頻信號量化為兩路I、Q信號,量化位數(shù)為12bit。
圖5 系統(tǒng)總體結(jié)構(gòu)Fig.5 Overall system structure
系統(tǒng)的工作過程如圖6所示:陣列天線接收兩路混合的機會信號,F(xiàn)PGA配置本振頻率大小,將AD前端的射頻信號下變頻為中頻信號,通過AD轉(zhuǎn)換后得到兩路I、Q信號,并傳送到FPGA芯片內(nèi),接下來在FPGA中啟動EASI算法模塊對信號做BSS處理,最后通過USB3.0數(shù)據(jù)線將FIFO中緩存的數(shù)據(jù)上傳到上位機中進(jìn)行識別和實時顯示。上位機通過RS485接口實現(xiàn)與FPGA的通信,以配置信號采集過程中的各種參數(shù)信息,包括本振頻率大小、信號采集頻段范圍和上傳數(shù)據(jù)的通道選擇等。
圖6 系統(tǒng)工作流程Fig.6 System workflow
為了能夠驗證本文提出的基于FPGA的EASI算法的實際分離效果,在實際環(huán)境中測試了基于軟件無線電的機會信號實時盲分離系統(tǒng),實際環(huán)境中測試的實物流程如圖7所示。首先利用本實驗室開發(fā)的機會信號仿真系統(tǒng)發(fā)射2個機會信號,然后利用陣列天線采集射頻信號,數(shù)字中頻板中FPGA板卡通過FMC(FPGA連接通用模塊)與AD9361芯片相連,F(xiàn)PGA硬件板卡負(fù)責(zé)對AD9361進(jìn)行編程控制以產(chǎn)生相應(yīng)的本振頻率用于射頻信號的下變頻,中頻信號經(jīng)過模數(shù)轉(zhuǎn)換后傳送到FPGA中進(jìn)行BSS處理,最后將分離信號上傳到上位機中。信號上傳的通道選擇和存儲方式選擇可以在上位機中設(shè)置,同時在上位機界面中能夠顯示信號的時域波形和頻域波形,系統(tǒng)的分離效果通過上位機中的識別結(jié)果和載頻估計結(jié)果來評價。
圖7 系統(tǒng)實際測試流程Fig.7 System actual test process
1)實驗1
機會信號仿真系統(tǒng)產(chǎn)生1275.52MHz的單頻信號和載頻為1276.52MHz、帶寬為5MHz的BPSK信號,經(jīng)射頻天線發(fā)出。在上位機界面中設(shè)置接收頻率為1268.52MHz,中頻為1MHz,即在AD前端產(chǎn)生1267.52MHz的本振頻率用于下變頻。同時在QT界面中設(shè)置接收信號通道,點擊參數(shù)發(fā)送將參數(shù)信息下發(fā)至FPGA中。之后系統(tǒng)開始對接收的混合信號做實時的BSS處理,并將分離后的信號通過USB3.0數(shù)據(jù)線上傳,上位機接收FPGA上傳的信號數(shù)據(jù)。最后點擊單路信號識別按鍵能夠?qū)崟r識別信號、計算信號載頻,并顯示信號的時域波形和頻域波形。圖8所示為通道設(shè)置選擇通道一時的接收信號波形,圖9所示為通道設(shè)置選擇通道二時的接收信號波形。
圖8 通道一上傳的CW信號Fig.8 CW signal uploaded on channel one
圖9 通道二上傳的BPSK信號Fig.9 BPSK signal uploaded on channel two
從圖8和圖9中可以發(fā)現(xiàn),通道一上傳的信號是分離后的CW信號,通道二上傳的信號是分離后的BPSK信號。CW信號的載頻估計誤差為0.1%,BPSK信號的載頻估計誤差為2.5%,識別結(jié)果和載頻估計結(jié)果表明了系統(tǒng)能夠?qū)崟r地分離CW和BPSK的混合信號。
2)實驗2
機會信號仿真系統(tǒng)產(chǎn)生載頻為1276.52MHz、帶寬為5MHz的FM信號,和載頻為1275.52MHz的CW信號,其他實驗配置和實驗1相同。
從圖10和圖11中可以發(fā)現(xiàn),通道一上傳的信號是分離后的FM信號,通道二上傳的信號是分離后的CW信號。FM信號的載頻估計誤差為2.5%,CW信號的載頻估計誤差為0.2%,識別結(jié)果和載頻估計結(jié)果表明了系統(tǒng)能夠?qū)崟r地分離CW和FM的混合信號。
圖10 通道一上傳的FM信號Fig.10 FM signal uploaded on channel one
3)實驗3
機會信號仿真系統(tǒng)產(chǎn)生載頻為1274.52MHz、帶寬為5MHz的BPSK信號,和載頻為1276.52MHz、帶寬為5MHz的BPSK信號,其他實驗配置和實驗1相同。
從圖12和圖13中可以發(fā)現(xiàn),通道一上傳的信號是分離后的載頻為1276.52MHz的BPSK信號,通道二上傳的信號是分離后的載頻為1274.52MHz的BPSK信號。通道一上傳的BPSK信號的載頻估計誤差為3.3%,通道二上傳的BPSK信號的載頻估計誤差為1.9%,識別結(jié)果和載頻估計結(jié)果表明了系統(tǒng)能夠?qū)崟r地分離兩種相同調(diào)制體制的機會信號。
圖12 通道一上傳的BPSK信號Fig.12 BPSK signal uploaded on channel one
圖13 通道二上傳的BPSK信號Fig.13 BPSK signal uploaded on channel two
本文針對多機會源信號識別難度較大的問題,設(shè)計了一套基于軟件無線電的機會信號實時盲分離系統(tǒng)。實際實驗結(jié)果表明:
1)該系統(tǒng)能夠?qū)崟r有效地分離兩通道的混合機會信號,相比于其他硬件實現(xiàn)BSS算法的技術(shù),該系統(tǒng)能夠?qū)﹃嚵行盘柼幚碇谐R姷膹?fù)數(shù)域混合信號進(jìn)行盲源分離。
2)基于軟件無線電的系統(tǒng)結(jié)構(gòu)能夠靈活地配置接收信號的頻段范圍和中頻大小,使得系統(tǒng)能夠適用于各種頻段范圍內(nèi)信號的盲源分離處理。除此以外,在上位機界面中對分離后的信號進(jìn)行了實時的識別、載頻估計和波形顯示,以評價系統(tǒng)的分離效果。
3)EASI算法模塊的設(shè)計采用了流水線的架構(gòu)和浮點運算單元,提高了整個算法的實時處理能力和數(shù)據(jù)計算精度。
4)目前,由于AD采集通道數(shù)目的限制,系統(tǒng)只能夠分離兩路混合機會信號,而實際環(huán)境中可能會存在更多的機會信號,從而影響系統(tǒng)的分離效果。今后會將兩通道的復(fù)數(shù)域EASI算法推廣到四通道,進(jìn)一步提高系統(tǒng)的工作性能。