基于软件无线电的机会信号实时盲分离系统

2020-06-08叶万洋郎荣玲

导航定位与授时 2020年3期

叶万洋，郎荣玲

(北京航空航天大学电子信息工程学院，北京 100191)

0 引言

在复杂的电磁环境中，存在多种异质随机信号共存的情况，无法直接对机会信号进行识别、评估以及定位。为了解决这些问题，需要对接收端的多源信号进行分离。一般情况下，存在对机会信号的先验知识未知的情况。因此，需要借助盲源分离(Blind Source Separation，BSS)技术解决这个问题。

BSS起源于著名的鸡尾酒会问题[1]，其试图在未知信号源特性与通道传输特性的前提下，仅根据接收到的信号来恢复原始信号。盲代表我们对信号的形式并不关心，而仅利用预先假设信号所具有的某一特性来建立目标函数，从而实现混合信号的分离处理。独立成分分析(Independent Component Analysis, ICA)是一种常用的BSS技术[2]，ICA的核心是寻找一个可以最大化混合信号各分量相互独立性的线性变换，从混合信号中分离出各个独立的源信号。基于ICA技术的BSS算法包括ICA算法[3]、快速独立成分分析(Fast-ICA)算法[4]和等变自适应盲分离(Equivariant Adaptive Source Separation via Independence，EASI)算法[5]。目前，基于ICA技术的BSS算法在语音信号、生物医学信号、水声信号、移动通信、地震信号处理等领域中均得到了广泛应用[6-7]。但是大多数BSS算法都使用计算机或服务器执行离线分析，并没有应用于信号的实时处理中。由于BSS算法的计算中涉及大量矩阵代数运算，因此非常适合在大规模并行处理的现场可编程门阵列(Field-Programmable Gate Array，FPGA)和专用集成电路(Application Specific Integrated Circuit，ASIC)芯片上实现，并行的计算方式和大量的逻辑单元能够保证BSS算法的实时实现。

目前已经提出了几种BSS算法的ASIC或FPGA实现[8-9]，以加速BSS处理。但是这些研究成果大多都是针对实数域的声音信号或胎儿心电信号，并不能够直接用于复数域混合信号的BSS。Van等提出了一种面向ASIC芯片内存不足的在线Fast-ICA算法和用于八通道脑电图(Electroencephalography，EEG)信号分离的硬件体系结构[10]。Charoensak等提出了基于FPGA的ICA算法[11]，该算法能够实时地分离两种声音的混合信号，但是只能够处理低频的声音信号。Kim等提出了用于BSS和自适应噪声消除(Adaptive Noise Cancelling，ANC)的ICA算法的FPGA实现[12]，在嘈杂环境中能够分离出噪声信号从而增强语言信号，但是该算法在实际应用中会存在一定的延时。Shyu等提出了用于实时信号BSS处理的Fast-ICA算法的FPGA实现方案[13]，设计采用了流水线的架构和浮点运算处理来提高工作性能。ODOM在文献[14]中对比了FPGA实现ICA、Fast-ICA和EASI三种BSS算法的性能和所需资源的数量，其中EASI算法由于具有收敛速度快、算法简单、有效性高等优点，被认为是最佳算法。

在本文中，为了能够实时地分离复数域的多源混合机会信号，设计了一套基于软件无线电的机会信号实时盲分离系统。系统软件部分利用C++和QT编写的信号采集程序，能够在上位机中设置采集信号的频带范围、下变频过程中的本振频率大小和信号上传的通道选择，并且能够实时存储和显示分离信号的时域波形和频域波形。硬件部分的主要功能是在FPGA中实时实现两通道的复数域EASI算法。算法采用流水线架构和浮点运算处理设计，其中流水线架构能够提高算法的实时处理能力；相比于数字信号处理中的定点设计，浮点算术设计能够提供更高的精度和动态性能。在实际实验中，该系统能够实时地分离两路复数域混合机会信号。

1 盲源分离

1.1 系统模型

盲源分离的系统模型如图1所示。

图1 盲源分离系统模型Fig.1 Model of blind source separation system

假设天线阵列由M个天线单元组成，则瞬时线性混合模型可以表示为

x=As+N

(1)

其中，x表示接收到的混合信号，s是源信号，A是满秩的混合矩阵，N是加性高斯白噪声。BSS的核心思想就是在s和A未知的情况下，仅根据x来估计一个分离矩阵W，使得x与W相乘后的分离信号y尽可能地接近s，可以表示为

y=Wx=WAs+WN

(2)

在盲源分离过程中，由于仅利用源信号之间相互统计独立这一条性质，并不知道s和A的先验信息。因此根据文献[2]可知，分离信号y和源信号s之间存在着排列顺序和幅度不确定性，即有

C=WA=PΛ

(3)

其中，C为广义置换矩阵(Generalized Permutation Matrix)；P为置换矩阵(即初等矩阵)，表示了分离信号各分量排列顺序的不确定性；Λ为满秩对角矩阵，表示了分离信号各分量幅度的不确定性。由于分离信号各分量的排列顺序和幅度的不确定性并不影响信号的波形，因此可以认为源信号能够被完全恢复。

1.2 EASI算法

EASI算法是将白化过程和高阶去相关过程同时进行的一种等变化的自适应算法[5]。其核心思想是通过自适应地更新W，使得y的各分量尽可能地相互独立，即各个分量的互信息代价函数最小。EASI算法中W的迭代更新方向为自然梯度方向，即代价函数的最陡下降方向。

典型的代价函数的形式为

φ(W)=E(f(y))

(4)

但其必须在去相关约束下进行优化

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算法不能恢复源信号的真实幅度，因此可以假设源信号具有单位方差和零协方差,即

Rs=E[ssT]=I

(9)

因为U是正交矩阵，y=Uz,为了满足去相关约束条件：z需要具有单位方差和零协方差，即

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的互信息代价函数进行优化

φ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表示相对梯度矩阵，非线性函数φ2(y)的计算公式通常选择为φ2(y)=y3。

2 EASI算法的FPGA实现

本文使用Xilinxx公司提供的EDA工具ISE手动编写Verilog代码，并在FPGA中实时实现了两通道的复数域EASI算法。根据式(21)将EASI算法分为3个模块，y计算模块，H计算模块，W更新模块。其中y、W和H都是复数域矩阵，表示形式如下

(22)

其中，i表示实部，q表示虚部，j表示虚数单位。

2.1 FPGA架构设计

EASI算法采用自适应的更新方法，即每一时刻的输入x都会更新一次W，并且每一个更新过程中都利用了上一个时刻x更新得到的W。然而在FPGA中更新一次W往往需要多个时钟周期，并不能在1个时钟周期内完成。为了在FPGA设计中满足自适应算法的要求，采用流水线架构来保证1个时钟周期更新一次W。流水线的原理如图2所示，将EASI算法分成若干个时间上均衡的子模块，从流水线的起点连续地输入数据，流水线的各子模块以重叠方式执行。这种处理方式能够使得W的更新速度只与x的输入速度有关，而与处理所需的时间无关。采用流水线的架构设计虽然能够保证每一时刻更新一次W，但是W的更新过程仍然需要多个时钟周期。为了解决这个问题，将迭代步长μ的值设置为0.0004，这样在一次W的更新过程内，W的值基本不变，因此可以用当前时刻的W来代替多个时钟周期后更新的W，从而保证了算法中整个W的更新过程的正确性和稳定性。

图2 EASI算法流水线处理Fig.2 EASI algorithm pipeline processing

由于将迭代步长μ的值设置的很小，对每一次更新过程中的数据的精度要求很高，因此在本文中整个算法的运算处理都采用浮点运算方式来保证数据精度。浮点运算采用ISE中的浮点IP核进行处理，包括浮点加法器、浮点减法器和浮点乘法器。除此以外，由于AD量化后的数字信号是定点数据，因此需要使用定点转浮点IP核将输入的定点数转化为浮点数；由EASI模块输出的到FIFO中缓存的数据是定点数，而实际运算的结果是浮点数，因此还需要使用浮点转定点IP核将输出信号的浮点数转化为定点数上传。EASI算法的硬件整体结构如图3所示。

图3 EASI算法硬件实现结构图Fig.3 EASI algorithm hardware implementation structure diagram

2.2 EASI算法模块设计

为了节省FPGA的逻辑元素，整个算法采用了时序硬件设计，其余计算将按顺序共享浮点运算单元。下面将详细描述EASI算法的3个主要模块

(1)y计算模块

此模块使用W和x计算复数y1和y2，运算流程如图3所示。W的初始值是单位矩阵，将式(2)改写成复数域形式如下

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实现中，此模块需要3个IP核时钟周期，其中y1的运算流程如图4所示。输出结果y分别是H计算模块的输入和EASI算法的输出。

图4 复数域矩阵运算硬件流程图Fig.4 Hardware flowchart of complex matrix operation

(2)H计算模块

根据式(22)可知，H函数的计算公式如下

H(t)=I+y(t)φ2(y(t))T-y(t)y(t)T-

φ2(y(t))y(t)T

(24)

由于式(24)中包含非线性函数φ2(y)=y3，因此在硬件中直接求解相对梯度矩阵H会消耗大量的逻辑资源。为了降低H的计算复杂度，从矩阵运算的角度将式(22)中y的复数域表示形式代入式(24)中，推导并简化相对梯度矩阵H各分量的求解公式，结果如下

(25)

其中，由于H11q和H22q的计算结果始终为0，因此在实际的硬件处理过程中不需要考虑这2个参数以节省硬件资源。此模块的输出H是W更新模块的输入。

(3)W更新模块

根据式(21)利用上一个模块得到H和设定的迭代步长μ来更新W。首先计算中间矩阵M=I+μH，由于H11q和H22q的计算结果始终为0，因此M11q和M22q在整个算法的运算过程中也始终为0，这样在硬件计算中同样可以不用考虑这4个参数，以节省硬件资源。则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的计算。

2.3 系统总体结构

机会信号实时盲分离系统采用软件无线电的设计结构，系统总体结构如图5所示。软件部分在上位机中利用C++和QT编写完成，硬件部分由接收天线、AD9361、FPGA组成。其中FPGA芯片型号为Xc7k325t，AD9361芯片的采样率为40MHz， AD9361芯片前端的混频模块能够对信号做下变频处理，并将中频信号量化为两路I、Q信号，量化位数为12bit。

图5 系统总体结构Fig.5 Overall system structure

系统的工作过程如图6所示：阵列天线接收两路混合的机会信号，FPGA配置本振频率大小，将AD前端的射频信号下变频为中频信号，通过AD转换后得到两路I、Q信号，并传送到FPGA芯片内，接下来在FPGA中启动EASI算法模块对信号做BSS处理，最后通过USB3.0数据线将FIFO中缓存的数据上传到上位机中进行识别和实时显示。上位机通过RS485接口实现与FPGA的通信，以配置信号采集过程中的各种参数信息，包括本振频率大小、信号采集频段范围和上传数据的通道选择等。

图6 系统工作流程Fig.6 System workflow

3 实验验证

为了能够验证本文提出的基于FPGA的EASI算法的实际分离效果，在实际环境中测试了基于软件无线电的机会信号实时盲分离系统，实际环境中测试的实物流程如图7所示。首先利用本实验室开发的机会信号仿真系统发射2个机会信号，然后利用阵列天线采集射频信号，数字中频板中FPGA板卡通过FMC(FPGA连接通用模块)与AD9361芯片相连，FPGA硬件板卡负责对AD9361进行编程控制以产生相应的本振频率用于射频信号的下变频，中频信号经过模数转换后传送到FPGA中进行BSS处理，最后将分离信号上传到上位机中。信号上传的通道选择和存储方式选择可以在上位机中设置，同时在上位机界面中能够显示信号的时域波形和频域波形，系统的分离效果通过上位机中的识别结果和载频估计结果来评价。

图7 系统实际测试流程Fig.7 System actual test process

1)实验1

机会信号仿真系统产生1275.52MHz的单频信号和载频为1276.52MHz、带宽为5MHz的BPSK信号，经射频天线发出。在上位机界面中设置接收频率为1268.52MHz，中频为1MHz，即在AD前端产生1267.52MHz的本振频率用于下变频。同时在QT界面中设置接收信号通道，点击参数发送将参数信息下发至FPGA中。之后系统开始对接收的混合信号做实时的BSS处理，并将分离后的信号通过USB3.0数据线上传，上位机接收FPGA上传的信号数据。最后点击单路信号识别按键能够实时识别信号、计算信号载频，并显示信号的时域波形和频域波形。图8所示为通道设置选择通道一时的接收信号波形，图9所示为通道设置选择通道二时的接收信号波形。

图8 通道一上传的CW信号Fig.8 CW signal uploaded on channel one

图9 通道二上传的BPSK信号Fig.9 BPSK signal uploaded on channel two

从图8和图9中可以发现，通道一上传的信号是分离后的CW信号，通道二上传的信号是分离后的BPSK信号。CW信号的载频估计误差为0.1%，BPSK信号的载频估计误差为2.5%，识别结果和载频估计结果表明了系统能够实时地分离CW和BPSK的混合信号。

2)实验2

机会信号仿真系统产生载频为1276.52MHz、带宽为5MHz的FM信号，和载频为1275.52MHz的CW信号，其他实验配置和实验1相同。

从图10和图11中可以发现，通道一上传的信号是分离后的FM信号，通道二上传的信号是分离后的CW信号。FM信号的载频估计误差为2.5%，CW信号的载频估计误差为0.2%，识别结果和载频估计结果表明了系统能够实时地分离CW和FM的混合信号。

图10 通道一上传的FM信号Fig.10 FM signal uploaded on channel one

3)实验3

机会信号仿真系统产生载频为1274.52MHz、带宽为5MHz的BPSK信号，和载频为1276.52MHz、带宽为5MHz的BPSK信号，其他实验配置和实验1相同。

从图12和图13中可以发现，通道一上传的信号是分离后的载频为1276.52MHz的BPSK信号，通道二上传的信号是分离后的载频为1274.52MHz的BPSK信号。通道一上传的BPSK信号的载频估计误差为3.3%，通道二上传的BPSK信号的载频估计误差为1.9%，识别结果和载频估计结果表明了系统能够实时地分离两种相同调制体制的机会信号。