CN105608057A - 一种分时复用硬件资源的信号子空间分解的fpga实现模块及其fpga实现方法 - Google Patents

Abstract

本发明公开了一种分时复用硬件资源的信号子空间分解的FPGA实现模块及其FPGA实现方法，其特征是包括自动排序的过关Jacobi算法单元和空间信号源数目估计计算单元；自动排序的过关Jacobi算法单元用于获得N阶的对角矩阵；空间信号源数目估计计算单元用于获得信号源的数目估计值。本发明能在满足速度要求的前提下，实现分时复用FPGA中的同一组硬件结构，从而节约硬件资源并提高灵活性，进而实现MUSIC算法中信号子空间分解的准确、快速计算。

Description

一种分时复用硬件资源的信号子空间分解的FPGA实现模块及其FPGA实现方法

技术领域

本发明属于阵列信号处理技术领域，具体的说是一种分时复用硬件资源的信号子空间分解的FPGA实现模块及其FPGA实现方法。

背景技术

20世纪70年代末开始，在空间谱估计方面涌现出了大量的研究成果和文献，波达方向估计(DOA)是空间谱估计研究的重要课题。对于同一参考信号源来说，各天线阵元所接受到的信号之间存在着相位差，从而形成谐波，不同方向的信号源对应着不同的谐波频率，只要估计出各谐波频率，即可求出各对应信号源的波达方向。最初的波达方向估计方法是基于傅里叶变化的线性谱估计方法，主要包括BT法和周期法。由于受到瑞利极限的限制，无法获得超高分辨率性能，且抗噪声能力差，未能取得满意的效果。后出现基于统计分析的极大似然谱估计方法，因其具有很高分辨性能而受到人们关注，然而，最大似然估计法要对高维参量空间进行搜索，运算量极大，难于在实践中得到应用。1967年，Burg提出了最大熵谱估计方法，开始了现代谱估计的研究。八十年代以后，学术界提出了一类基于矩阵特征值分解的谱估计方法，其中由美国Schmidt等人提出的多重信号分类(MUSIC)算法最为突出。其核心思想是对阵列输出的协方差矩阵做特征值分解，把整个特征空间分为信号子空间和噪声子空间，利用两者互相正交的特性，定义一个伪谱，搜索到的谱峰对应的角度即为入射信号的波达方向。

在MUSIC算法中，信号和噪声子空间分解是硬件实现的一个难题，矩阵特征值分解是算法硬件实现中最为复杂的问题。针对这一问题，目前，矩阵特征值分解主要有Jacobi算法和QR算法，由于Jacobi比QR算法更加准确，因此设计中常常采用Jacobi算法进行矩阵的特征值分解，再根据估计信号源数目计算，将特征向量空间分解为信号子空间和噪声子空间，这其中包含了大量的乘法和非线性的运算，如何合理的实现这些运算是在硬件实现时亟需考虑的问题。

子空间分解需要大量的计算，在矩阵特征值分解和信号源估计中又涉及不同的运算方式，需要采用不同的硬件处理结构，目前所提出的几种硬件实现方法上，存在以下的缺陷和不足：

第一，对于高速实现子空间分解问题，国内目前大多使用数字信号处理器(DSP)实现，对于八元阵而言，能够把MUSIC整个算法的计算时间降低至少一半，加速比达到2.37，在ms量级，但在很多实际场合仍难达到要求。

第二，在矩阵特征值分解步骤采用FPGA实现，信号源数目估计计算采用DSP实现。对于计算八元阵的矩阵特征值而言，目前可在us量级完成，时间明显少于DSP，但信号源数目估计计算使用DSP，一方面使得整个计算时间不得不停留在ms量级，另一方面导致整个系统涉及到FPGA和DSP的互联互接问题，涉及到FPGA和DSP运算之间的相互转换，增加了设计难度。

第三，采用FPGA来实现整个子空间分解，能够完成的子空间分解更加快捷的实现技术，具有灵活性强、并行度高、便于调试和扩展接口的特点。通过Jacobi算法实现特征值分解，对于其中的非线性运算，在保证一定精确度的条件下，可以仅用加法和移位就能实现，计算速度大大增加。目前，完成四元阵的MUSIC算法需要30.59us，使用逻辑资源是12007个；完成八元阵的MUSIC算法需要373.99us，使用逻辑资源是29472个。而对基于信息论方法的AIC准则信号源数目估计计算，涉及大量对数运算，需要硬件系统中存在对数运算的协处理器，所以基于FPGA的子空间分解实现，需要两个不同的结构来实现Jacobi算法和信号源数目估计计算，消耗了大量硬件资源，同时也降低了其分解的速度，因此在资源消耗和速度上有可以改进的空间。

发明内容

本发明是针对现有FPGA实现MUSIC中子空间分解算法的不足，并结合自动排序的过关Jacobi算法的特点，利用等价变换的信号源数估计计算，提出一种分时复用硬件资源的信号子空间分解的FPGA实现模块及其FPGA实现方法，以期能在满足速度要求的前提下实现分时复用FPGA的硬件资源，从而节约硬件资源并提高灵活性，进而实现信号子空间分解的准确、快速计算。

本发明为解决技术问题采用如下技术方案：

本发明一种分时复用硬件资源的信号子空间分解的FPGA实现模块的特点是包括：自动排序的过关Jacobi算法单元和空间信号源数目估计计算单元；

所述自动排序的过关Jacobi算法单元读取空间信号源对应的N阶协方差矩阵Α＝[a_pq]_n×n并进行运算，获得N阶的对角矩阵Α_k；所述对角矩阵Α_k的主对角线上的元素为所述N阶协方差矩阵Α的特征值λ₁,λ₂,…,λ_p,…,λ_q,…,λ_n；且λ₁,λ₂,…,λ_p,…,λ_q,…,λ_n按降序排列；1≤p≤n；1≤q≤n；p≠q；

所述空间信号源数目估计计算单元将所述特征值λ₁,λ₂,…,λ_p,…,λ_q,…,λ_n作为源操作数并进行运算，获得2n个运算结果，包括n个特征值和的对数ln_sum₁,ln_sum₂,…，ln_sum_n及n个特征值积的对数ln_acc₁,ln_acc₂,…,ln_acc_n；再利用等价变换的AIC准则法对所述2n个运算结果进行运算，获得一个先递减后递增的AIC序列AIC₁,AIC₂,…,AIC_p,…,AIC_q,…,AIC_n；从所述AIC序列AIC₁,AIC₂,…,AIC_p,…,AIC_q,…,AIC_n中选出一个最小值，以最小值所对应的下标并减去1即为信号源的数目估计值num。

本发明一种分时复用硬件资源的信号子空间分解的FPGA实现方法的特点是按如下步骤进行：

步骤1、对接收到的空间信号源所对应的N阶协方差矩阵进行自动排序的过关Jacobi算法处理，获得N阶的对角矩阵A_k；

步骤1.1、选取一个单调递减且逐渐趋于零的数列a₁,a₂,…,a_i,…,a_k作为关值，并满足a₁＞a₂＞…a_i…＞a_k且l $\underset{k\→\mathrm{\∞}}{\lim }{a}_k=0,$ 并有 $a_i=a_{i-1}/n=\sqrt{2\underset{p=1}{\overset{n-1}{\Σ}}\underset{q=p+1}{\overset{n}{\Σ}}{a}_{pq}}/n;1\≤i\≤k;$

步骤1.2、初始化i＝1；

步骤1.3、初始化p＝1；

步骤1.4、初始化q＝1；

步骤1.5、判断p＝q是否成立，若成立，则将q+1赋值给q；并执行步骤1.6，若不成立，则直接执行步骤1.6；

步骤1.6、判断|a_pq|＞a_i是否成立，若成立，则对N阶协方差矩阵Α利用自动排序的过关Jacobi算法进行旋转变换，获得第i次旋转变换矩阵A_i；若不成立，则将q+1赋值给q；并返回步骤1.6执行，直到q＝n为止；

步骤1.7、将p+1赋值给p；并返回步骤1.5执行；直到p＝n为止；

步骤1.8、将i+1赋值给i；并返回步骤1.3执行；直到i＝k为止；从而获得N阶的对角矩阵Α_k；且所述对角矩阵Α_k的主对角线上的元素为所述N阶协方差矩阵Α的特征值λ₁，λ₂，…λ_p，…λ_q，…λ_n；且λ₁，λ₂，…λ_p，…λ_q，…λ_n按降序排列；

步骤2、利用等价变换的AIC准则估计信号源数目

步骤2.1、对式(1)所示的AIC表达式进行等价变换，获得等价变换的AIC表达式如式(2)所示：

AIC(num+1)＝2K(N-num)lnΛ(num)+2num(2N-num)(1)

$AIC\left(num+1\right)=-2\left(M-num\right)K\log 10\left(\right(\left(\underset{n=num+1}{\overset{N}{\Π}}{\mathrm{\λ}}_j^{\left(L/\left(M-num\right)\right)}\right)\left(\frac{\underset{n=num+1}{\overset{N}{\Σ}}{\mathrm{\λ}}_j}{M-num}\right))+2num(2M-num)---\left(2\right)$

式(1)和式(2)中，K表示信号采样的快拍数，num表示待估计的信号源的个数，N为阵元数，Λ(num)为似然函数；λ_j表示第j个特征值，1≤j≤n；

步骤2.2、利用式(3)所示的泰勒级数展开式获得第j个特征值的对数ln(λ_j)：

$\ln \left({\mathrm{\λ}}_j\right)=ln\left(\frac{x+1}{1-x}\right)=2x(1+\frac{1}{3}{x}^2+\frac{1}{5}{x}^4+\frac{1}{7}{x}^6+...+\frac{1}{2m+1}{x}^{2m}---\left(3\right)$

式(3)中，x表示操作数，且x∈[-1,1]；m表示泰勒级数展开式的项数；

步骤2.3、利用式(4)获得第j个特征值λ_j所对应的N个累加和对数ln_sum_j，从而获得n个特征值和的对数ln_sum₁,ln_sum₂,…,ln_sum_j,…,ln_sum_n：

$\ln \_{\mathrm{sum}}_j=ln\left(\underset{m=j}{\overset{N}{\Σ}}{\mathrm{\λ}}_j\right)---\left(4\right)$

步骤2.4、利用式(5)获得第j个特征值λ_j所对应的N个乘积的对数ln_acc_j，从而获得n个特征值和的对数ln_acc₁,ln_acc₂,…,ln_acc_j,…,ln_acc_n：

$\ln \_{\mathrm{acc}}_j=ln\left(\underset{m=j}{\overset{N}{\Σ}}{\mathrm{\λ}}_j\right)---\left(5\right)$

步骤2.5、将n个特征值和的对数ln_sum₁,ln_sum₂,…,ln_sum_j,…,ln_sum_n以及n个特征值和的对数ln_acc₁,ln_acc₂,…,ln_acc_j,…,ln_acc_n代入式(2)中，从而获得一个先递减后递增的AIC序列AIC₁,AIC₂…,AIC_p…，AIC_q…,AIC_n；

步骤2.6、从所述AIC序列AIC₁,AIC₂…,AIC_p…,AIC_q…,AIC_n中选出一个最小值，以最小值所对应的下标并减去1即为信号源的数目估计值num。

本发明一种分时复用硬件资源的信号子空间分解的FPGA实现电路的特点是包括：第一处理模块PE1、第二处理模块PE2、第三处理模块PE3、寄存器组0register_group0、寄存器组1register_group1、寄存器组2register_group2、寄存器组3register_group3、寄存器组4register_group4、比较单元Comp_PE、存储器D_RAM；

所述存储器D_RAM存入空间信号源对应的N阶协方差矩阵Α＝[a_pq]_n×n；

所述寄存器组0register_group0从所述存储器D_RAM中读取N阶协方差矩阵Α的p行元素；

所述寄存器组1register_group1从所述存储器D_RAM中读取N阶协方差矩阵Α的q行元素；

所述寄存器组4register_group4从所述存储器D_RAM中读取N阶协方差矩阵Α的p列和q列元素；

所述比较单元Comp_PE从所述存储器D_RAM中依次读取N阶协方差矩阵Α中的非主对角线元素a_pq并与所设定的关值进行比较，若非主对角线元素a_pq的绝对值大于设定的关值，则利用第一处理模块PE1对N阶协方差矩阵中主对角线上的元素进行计算，获得自动排序的过关Jacobi算法中的cosθ和sinθ；再根据cosθ和sinθ、所述寄存器组0register_group0、寄存器组1register_group1、寄存器组4register_group4中存储的元素，当非主对角线元素a_pq的绝对值大于设定的关值时，利用第二处理模块PE2和所述第三处理模块PE3进行旋转变换；从而获得所述N阶协方差矩阵Α变换后的对角矩阵Α_k并传递给所述存储器D_RAM进行存储；所述变换后的对角矩阵Α_k的主对角线上的元素为按降序排列的特征值λ₁，λ₂，…λ_p，…λ_q，…λ_n；

所述第一处理模块PE1和第二处理模块PE2分别从所述存储器D_RAM中读取特征值λ₁，λ₂，…λ_p，…λ_q，…λ_n并进行计算，获得特征值对数lnλ₁，lnλ₂，…lnλ_p，…lnλ_q，…lnλ_n并存入所述寄存器组1register_group1中；

所述寄存器组2register_group2从所述存储器D_RAM中读取特征值λ₁，λ₂，…λ_p，…λ_q，…λ_n并传递给所述第三处理模块PE3；所述第三处理模块PE3对所述特征值λ₁，λ₂，…λ_p，…λ_q，…λ_n进行计算获得运算结果ln_acc₁,ln_acc₂,…ln_acc_n存入寄存器组0register_group0中；

所述第二处理模块PE2从所述寄存器组1register_group1所读取特征值对数lnλ₁，lnλ₂，…lnλ_p，…lnλ_q，…lnλ_n并进行计算，获得运算结果ln_sum₁,ln_sum₂,…ln_sum_n并存入寄存器组3register_group3中；

所述第一处理模块PE1和所述第二处理模块PE2分别从所述寄存器组0register_group0和寄存器组3register_group3中读取特征值积的对数ln_acc₁,ln_acc₂,…ln_acc_n和特征值和的对数ln_sum₁,ln_sum₂,…ln_sum_n并进行计算，获得AIC序列AIC₁,AIC₂…,AIC_p…,AIC_q…,AIC_n后传递给所述比较单元Comp_PE；

所述比较单元Comp_PE对所述AIC序列AIC₁,AIC₂…,AIC_p…,AIC_q…,AIC_n中的各个元素进行比较，获得信号源的数目估计值num。

与已有技术相比，本发明有益效果体现在：

1本发明选择了自动排序的过关Jacobi算法，对传统的Jacobi算法进行改进，通过对传统的过关Jacobi算法中旋转变换cosθ和sinθ取值，得到最终的运算结果对角矩阵，其对角线上元素自动按照降序排列，即协方差矩阵的特征值按降序排列，使用改进的Jacobi算法实现了特征值的自动排序，节省了特征值排序的所需的计算时间和硬件资源。

2、本发明对估计信号源数的信息论方法中的一些运算作了等价代换，根据对数运算规律，从设计高效硬件电路的角度出发，将计算AIC函数值的理论公式进行了等价变形，主要是将似然函数的lnΛ(num)的运算中乘方运算和除法运算转化为乘积运算和减法运算，最大程度地降低了运算复杂度，减少了资源消耗并且提高了运算速度。

3、本发明通过自动排序的过关Jacobi算法实现特征值分解，对于其中的非线性运算，在保证一定精确度的条件下，可以仅用加法和移位就能实现，使得Jacobi算法能够很好的运用在FPGA中，计算速度大大增加。

4、本发明通过泰勒级数展开的方法求解信号源估计计算中所需要的对数运算，通过运算得出满足运算精度要求的对数运算结果，不需要使用协处理器进行对数运算，使得AIC准则可以成功地运用到FPGA硬件电路中，实现了信号源个数估计，缩短了实现信号源个数估计所需要的总时间。

5、本发明基于FPGA实现两个计算单元的运算，可以不受硬件运算单元数量(如乘法累加器数量)的影响，将自动排序过关Jacobi计算和空间信号源数目估计计算这两个步骤按一定结构映射到FPGA上，不需要考虑不同硬件结构实现时，每个硬件之间的接口通讯问题，且FPGA提供了更多的带宽，运算速度也有较大提高。

6、本发明实现自动排序过关Jacobi计算单元和空间信号源数目估计计算单元，在改进算法和等价变换的基础上，使得两者有大量相似的运算过程，从而不需要使用两个不同的硬件实现，分时复用FPGA上的三个处理模块PE1、PE2、PE3、比较单元Comp_PE和寄存器组，与其他的设计相比，在满足相同精度的情况下，计算速度大大提高，消耗的硬件资源面积很小。

7、本发明分时复用FPGA中的同一组硬件结构实现自动排序过关Jacobi算法和空间信号源数目估计计算，采用串行计算结构依次进行整个过程，灵活度强、并行度高，能够快速、准确的执行Jacobi算法和信号源数目估计的计算要求。

附图说明

图1为本发明雅克比算法的执行流程图；

图2为本发明FPGA实现总体结构图；

图3为本发明硬件中处理单元PE1内部结构图；

图4为本发明硬件中处理单元PE2、PE3内部结构图；

图5为本发明自动排序的过关Jacobi算法硬件结构图；

图6为本发明空间信号源数目估计计算硬件结构图。

具体实施方式

本实施例中，一种分时复用硬件资源的信号子空间分解的FPGA实现模块，包括自动排序的过关Jacobi算法单元和空间信号源数目估计计算单元；

本实施例中，自动排序的过关Jacobi算法总体执行流程图如图1所示，自动排序的过关Jacobi单元读取空间信号源对应的N阶协方差矩阵Α＝[a_pq]_n×n并进行运算，获得N阶的对角矩阵Α_k；对角矩阵Α_k的主对角线上的元素为N阶协方差矩阵Α的特征值λ₁,λ₂,…,λ_p,…,λ_q,…,λ_n；且λ₁,λ₂,…,λ_p,…,λ_q,…,λ_n按降序排列；1≤p≤n；1≤q≤n；p≠q；

具体实施中，设定阵元数为8，则协方差矩阵Α为一个8阶矩阵，Α中有64个元素，输入Α中64个元素，经过自动排序的过关Jacobi算法单元，获得一个8阶的对角矩阵Α_k，Α_k的主对角线上8个元素即为8个特征值λ₁,λ₂,…,λ₈；

空间信号源数目估计计算单元将特征值λ₁,λ₂,…,λ_p,…,λ_q,…,λ_n作为源操作数并进行运算，获得2n个运算结果，包括n个特征值和的对数ln_sum₁,ln_sum₂,…,ln_sum_n及n个特征值积的对数ln_acc₁,ln_acc₂,…,ln_acc_n；再利用等价变换的AIC准则法对2n个运算结果进行运算，获得一个先递减后递增的AIC序列AIC₁,AIC₂,…,AIC_p,…,AIC_q,…,AIC_n；从AIC序列AIC₁,AIC₂,…,AIC_p,…,AIC_q,…,AIC_n中选出一个最小值，以最小值所对应的下标并减去1即为信号源的数目估计值num；

具体实施中，将特征值λ₁,λ₂,…λ₈作为源操作数并进行运算，获得16个运算结果，包括8个特征值和的对数ln_sum₁,ln_sum₂,…,ln_sum₈及8个特征值积的对数ln_acc₁,ln_acc₂,…,ln_acc₈；利用等价变换的AIC准则法对16个运算结果进行运算，获得一个先递减后递增的AIC序列AIC₁,AIC₂,…,AIC₈；从AIC序列AIC₁,AIC₂,…,AIC₈中选出一个最小值，以最小值所对应的下标并减去1即为信号源的数目估计值num。

本实施例中，一种分时复用硬件资源的信号子空间分解的FPGA实现方法，按如下步骤进行：

步骤1、对接收到的空间信号源所对应的N阶协方差矩阵进行自动排序的过关Jacobi算法处理，获得N阶的对角矩阵Α_k；

步骤1.1、选取一个单调递减且逐渐趋于零的数列a₁,a₂,…,a_i,…,a_k作为关值，并满足a₁＞a₂＞…a_i…＞a_k且l $\underset{k\→\mathrm{\∞}}{\lim }{a}_k=0,$ 并有 $a_i=a_{i-1}/n=\sqrt{2\underset{p=1}{\overset{n-1}{\Σ}}\underset{q=p+1}{\overset{n}{\Σ}}{a}_{pq}}/n;1\≤i\≤k;$

步骤1.2、初始化i＝1；

步骤1.3、初始化p＝1；

步骤1.4、初始化q＝1；

步骤1.5、判断p＝q是否成立，若成立，则将q+1赋值给q；并执行步骤1.6，若不成立，则直接执行步骤1.6；

步骤1.6、判断|a_pq|＞a_i是否成立，若成立，则对N阶协方差矩阵Α利用自动排序的过关Jacobi算法进行旋转变换，获得第i次旋转变换矩阵A_i；若不成立，则将q+1赋值给q；并返回步骤1.6执行，直到q＝n为止；

步骤1.7、将p+1赋值给p；并返回步骤1.5执行；直到p＝n为止；

步骤1.8、将i+1赋值给i；并返回步骤1.3执行；直到i＝k为止；从而获得N阶的对角矩阵Α_k；且所述对角矩阵Α_k的主对角线上的元素为所述N阶协方差矩阵Α的特征值λ₁，λ₂，…λ_p，…λ_q，…λ_n；且λ₁，λ₂，…λ_p，…λ_q，…λ_n按降序排列；

具体实施中，利用自动排序的过关Jacobi算法，过关Jacobi算法的基本思想是：对N阶实对称矩阵Α，一定存在一个正交矩阵Q，使Q^-1AQ＝Q^TAQ＝D,其中对角矩阵D＝diag(λ₁,λ₂,…,λ_n)，λ₁,λ₂,…,λ_n为N阶实对称矩阵Α的特征值，而正交矩阵Q的第j列(1≤j≤n)就是对应于λ_j的特征向量；对N阶实对称协方差矩阵Α，记Α₀＝Α，首先对选择N阶实对称协方差矩阵计算第一个关值N阶实对称协方差矩阵Α中的元素，按a₁₂,a₁₃,…，a_1n,a₂₃,a₂₄,…,a_2n,a₁₂,…,a_n-1,n进行扫描，若|a_pq|＞a_i，则选取旋转矩阵G对A做旋转变换 $\left\{\begin{array}{c}{A}_1={\mathrm{GA}}_0{G}^T\\ A_0=A\end{array}\right.,$ 是个n×n正交的平面旋转方阵，其中g_pp＝cosθ，g_pq＝sinθ，g_qp＝-sinθ，g_qq＝cosθ，sinθ、cosθ的取值由条件(a_qq-a_pp)sin2θ+2a_pq(cos²θ-sin²θ)＝0定出；由于某些绝对值小于a₁的元素的绝对值可能在后面的旋转变换中增长，所有要进行多次扫描，直到Α₁的所有非零非对角线元素的绝对值都小于a₁为止；再分别设置第二个关值a₂＝a₁/n、第三个关值a₃＝a₂/n、……、第k个关值a_k＝a_k-1/n，重复上述过程，直至所有的非对角线元素都小于a_k时(此时的a_k应满足系统所给定的误差限ε)，迭代停止，最终化为对角矩阵此时对角矩阵Α_k的主对角线元素就是N阶协方差矩阵Α的特征值的λ₁，λ₂，…λ_p，…λ_q，…λ_n；

具体实施中，进一步采用自动排序Jacobi方法进行矩阵特征值自动排序，对于多次迭代得到的对角矩阵Α_k，由Jacobi方法推导得到 $a_{pp}^k-a_{qq}^k=a_{pp}^{k-1}({\cos }^2\mathrm{\&theta;}-{\sin }^2\mathrm{\&theta;})+a_{qq}^{k-1}({\sin }^2\mathrm{\&theta;}-{\cos }^2\mathrm{\&theta;})+4a_{pq}^{k-1}cos\mathrm{\&theta;}sin\mathrm{\&theta;},$ 令 $c=2a_{pq}^{k-1},$ $d=abs(a_{pp}^{k-1}-a_{qq}^{k-1}),$ 得到 $a_{pp}^k-a_{qq}^k=2cos\mathrm{\&theta;}sin\mathrm{\&theta;}x(d^2/c^2+1),$ 其中只要使 $a_{pp}^k-a_{qq}^k>0,$ 即cosθsinθ与c符号相同，就能够使对角矩阵Α_k中的对角线元素从大到小排列；取 $\cos \mathrm{\&theta;}=sqrt\left(\frac{sqrt(c^2+d^2)+d}{2sqrt(c^2+d^2)}\right),\sin \mathrm{\&theta;}=sign\left(c\right)sqrt\left(\frac{sqrt(c^2+d^2)-d}{2sqrt(c^2+d^2)}\right),$ 其中 $sign\left(c\right)=\left\{\begin{array}{c}1,c>0\\ -1,c<0\end{array}\right.,$ 实现矩阵特征值λ₁，λ₂，…λ_p，…λ_q，…λ_n按降序排列。

步骤2、利用等价变换的AIC准则估计信号源数目

步骤2.1、对式(1)所示的AIC表达式进行等价变换，获得等价变换的AIC表达式如式(2)所示：

AIC(num+1)＝2K(N-num)lnΛ(num)+2num(2N-num)(1)

$AIC\left(num+1\right)=-2\left(M-num\right)K\log 10\left(\right(\left(\underset{n=num+1}{\overset{N}{\&Pi;}}{\mathrm{\&lambda;}}_j^{\left(L/\left(M-num\right)\right)}\right)\left(\frac{\underset{n=num+1}{\overset{N}{\&Sigma;}}{\mathrm{\&lambda;}}_j}{M-num}\right))+2num(2M-num)---\left(2\right)$

式(1)和式(2)中，K表示信号采样的快拍数，num表示待估计的信号源的个数，N为阵元数，Λ(num)为似然函数；λ_j表示第j个特征值，1≤j≤n；

步骤2.2、利用式(3)所示的泰勒级数展开式获得第j个特征值的对数ln(λ_j)：

$ln\left({\mathrm{\&lambda;}}_j\right)=ln\left(\frac{x+1}{1-x}\right)=2x(1+\frac{1}{3}{x}^2+\frac{1}{5}{x}^4+\frac{1}{7}{x}^6+...+\frac{1}{2m+1}{x}^{2m}---\left(3\right)$

式(3)中，x表示操作数，且x∈[-1,1]；m表示泰勒级数展开式的项数；

步骤2.3、利用式(4)获得第j个特征值λ_j所对应的N个累加和对数ln_sum_j，从而获得n个特征值和的对数ln_sum₁,ln_sum₂,…,ln_sum_j,…,ln_sum_n：

$\ln \_{\mathrm{sum}}_j=ln\left(\underset{m=j}{\overset{N}{\&Sigma;}}{\mathrm{\&lambda;}}_j\right)---\left(4\right)$

步骤2.4、利用式(5)获得第j个特征值λ_j所对应的N个乘积的对数ln_acc_j，从而获得n个特征值和的对数ln_acc₁,ln_acc₂,…，ln_acc_j,…,ln_acc_n：

$\ln \_{\mathrm{acc}}_j=ln\left(\underset{m=j}{\overset{N}{\&Sigma;}}{\mathrm{\&lambda;}}_j\right)---\left(5\right)$

步骤2.5、将n个特征值和的对数ln_sum₁,ln_sum₂,…,ln_sum_j,…,ln_sum_n以及n个特征值和的对数ln_acc₁,ln_acc₂,…,ln_acc_j,…,ln_acc_n代入式(2)中，从而获得一个先递减后递增的AIC序列AIC₁,AIC₂…,AIC_p…,AIC_q…,AIC_n；

步骤2.6、从所述AIC序列AIC₁,AIC₂…,AIC_p…,AIC_q…,AIC_n中选出一个最小值，以最小值所对应的下标并减去1即为信号源的数目估计值num；

本实施例中，利用泰勒级数展开来求解对数运算，对数函数的泰勒级数展开中，操作数x需满足x∈[-1,1]，这样级数才能有效收敛，所以对泰勒级数展开作下面的变换，

当x∈[-1,1]时，根据泰勒级数展开式可得 $\left\{\begin{array}{c}ln(x+1)=x-\frac{1}{2}{x}^2+\frac{1}{3}{x}^3-\frac{1}{4}{x}^4+\frac{1}{5}{x}^5+...\\ ln(1-x)=-x-\frac{1}{2}{x}^2-\frac{1}{3}{x}^3-\frac{1}{4}{x}^4-\frac{1}{5}{x}^5+...\end{array}\right.,$ 对任意的x∈R,(x+1)/(1-x)∈[-1,1]；当需要计算某个任意数y∈R的对应得对数值ln(y)时，令y＝(x+1)/(1-x)，解得x＝(y-1)/(y+1)，任意数y∈R,|x|＜1，操作数x满足要泰勒级数展开式要求，可得 $\ln \left(y\right)=ln\left(\frac{x+1}{1-x}\right)=2x(1+\frac{1}{3}{x}^2+\frac{1}{5}{x}^4+\frac{1}{7}{x}^6+...+\frac{1}{2m+1}{x}^{2m},$ 其中取m＝30，即把对数展开到第30项，以此满足系统给定的精度要求。

本实施例中，一种分时复用硬件资源的信号子空间分解的FPGA实现电路包括第一处理模块PE1、第二处理模块PE2、第三处理模块PE3、寄存器组0register_group0、寄存器组1register_group1、寄存器组2register_group2、寄存器组3register_group3、寄存器组4register_group4、比较单元Comp_PE、存储器D_RAM；

具体实施中，在FPGA上实现的总体结构如图2所示，设定阵元数N为8，其中D_RAM为单端口RAM，存放64个64位的系统初始值，高32位是协方差矩阵元素，低32为存放8×8的对角阵元素。四个寄存器组寄存器组0register_group0、寄存器组1register_group1、寄存器组2register_group2、寄存器组3register_group3分别由8个64位寄存器组成，寄存器组4register_group4由16个32位寄存器组成，用于存放计算的中间数据。选择控制MUX_generate模块用来对各个多路选择器MUX进行不同的配置，得到不同的数据路径，从而实现不同的运算功能。

本实施例中，第一处理模块PE1的内部结构如图3，由1个加法、1个乘法、1个除法、1个根号运算单元及6个选择器构成；第二处理模块PE2和第三处理模块PE3内部结构相同，结构图如图4，由1个加法、2个乘法运算单元及4个选择器构成；

本实施例中，存储器D_RAM存入空间信号源对应的N阶协方差矩阵Α＝[a_pq]_n×n；

寄存器组0register_group0从存储器D_RAM中读取N阶协方差矩阵Α的p行元素；

寄存器组1register_group1从存储器D_RAM中读取N阶协方差矩阵Α的q行元素；

寄存器组4register_group4从存储器D_RAM中读取N阶协方差矩阵Α的p列和q列元素；

比较单元Comp_PE从存储器D_RAM中依次读取N阶协方差矩阵Α中的非主对角线元素a_pq并与所设定的关值进行比较，若非主对角线元素a_pq的绝对值大于设定的关值，则利用第一处理模块PE1对N阶协方差矩阵中主对角线上的元素进行计算，获得自动排序的过关Jacobi算法中的cosθ和sinθ；再根据cosθ和sinθ和寄存器组0register_group0、寄存器组1register_group1、寄存器组4register_group4中存储的元素，当非主对角线元素a_pq的绝对值大于设定的关值时，利用第二处理模块PE2和第三处理模块PE3进行旋转变换；若非主对角线元素a_pq的绝对值小于等于设定的关值，则所述第二处理模块PE2和所述第三处理模块PE3不做任何变换，将结果反馈到D_RAM中；重复上述过程，多次旋转变换后，获得N阶协方差矩阵Α的变换矩阵Α_k并传递给存储器D_RAM进行存储；变换矩阵Α_k的主对角线上的元素为按降序排列的特征值λ₁,λ₂,…,λ_p,…,λ_q,…,λ_n；

具体实施中，自动排序的过关Jacobi算法实现结构如图5，步骤如下：

步骤a：读入空间信号源对应的8阶协方差矩阵Α存入存储器D_RAM；

步骤b：寄存器组0register_group0、寄存器组1register_group1、寄存器组4register_group4分别用来存储旋转变换中p行、q行和p列、q列的数据；

步骤c：对8阶协方差矩阵Α中的第一个非主对角线元素a₁₂进行扫描，进入比较单元Comp_PE与过关Jacobi算法所设定的第一关值a₁进行比较；

步骤d：8协方差矩阵Α中的对角线元素a₁₁、a₂₂输入第一处理模块PE1，计算自动排序的过关Jacobi算法中cosθ、sinθ值，输出到第二处理模块PE2和第三处理模块PE3；

步骤e：若比较单元Comp_PE得到第一个非主对角线元素a₁₂的绝对值大于a₁，根据cosθ和sinθ及寄存器组0register_group0、寄存器组1register_group1、寄存器组4register_group4中存储的元素，第二处理模块PE2和第三处理模块PE3进行一次旋转变换得到变换后的8阶矩阵Α′，一次选择变换后的8阶矩阵Α′中的两个非主对角线元素a₁₂、a₂₁为0；

步骤f：若比较单元Comp_PE得到第一个非主对角线元素a₁₂的绝对值小于等于a₁，第二处理模块PE2和第三处理模块PE3不做任何变换，直接将8阶协方差矩阵Α反馈到存储器D_RAM中；

步骤g：对8阶协方差矩阵Α从行到列、从左到右进行多次扫描，扫描下一个非主对角线元素，重复步骤c、d、e、f，直到第一关值a₁对应的8阶矩阵Α₁中，所有非零非对角线元素的绝对值都小于设定的关值a₁为止；

步骤h：步骤c中取设定的其他关值a₂、a₃…做类似处理，直至所有的非对角线元素都小于a_k时，得到一个8阶的对角矩阵Α_k，8阶的对角矩阵Α_k中主对角线上的元素即为计算后自动降序排列的特征值λ₁，λ₂，…λ₈，计算结果返回存储器D_RAM中。

本实施例中，第一处理模块PE1和第二处理模块PE2分别从存储器D_RAM中读取特征值λ₁,λ₂,…,λ_p,…,λ_q,…,λ_n并进行计算，获得特征值对数lnλ₁,lnλ₂,…,lnλ_p,…,lnλ_q,…,lnλ_n并存入寄存器组1register_group1；

寄存器组2register_group2从存储器D_RAM中读取特征值λ₁,λ₂,…,λ_p,…,λ_q,…,λ_n并传递给第三处理模块PE3；第三处理模块PE3对特征值λ₁,λ₂,…,λ_p,…,λ_q,…,λ_n进行计算获得运算结果ln_acc₁,ln_acc₂,…,ln_acc_n存入寄存器组0register_group0；

第二处理模块PE2从寄存器组1所读取特征值对数lnλ₁,lnλ₂,…,lnλ_p,…,lnλ_q,…,lnλ_n并进行计算，获得运算结果ln_sum₁,ln_sum₂,…,ln_sum_n并存入寄存器组3register_group3；

第一处理模块PE1和第二处理模块PE2分别从寄存器组0register_group0和寄存器组3register_group3中读取特征值积的对数ln_acc₁,ln_acc₂,…,ln_acc_n和特征值和的对数ln_sum₁,ln_sum₂,…,ln_sum_n并进行计算，获得AIC序列AIC₁,AIC₂,…,AIC_p,…,AIC_q,…,AIC_n后传递给比较单元Comp_PE；

比较单元Comp_PE对AIC序列AIC₁,AIC₂,…,AIC_p,…,AIC_q,…,AIC_n中的各个元素进行比较，获得信号源的数目估计值num；

本实施例中，基于AIC准则等价变换的信号源数目估计计算实现结构如图6，设定阵元数N＝8,步骤如下：

步骤a：寄存器组2register_group2从存储器D_RAM中读入8阶协方差矩阵Α对应的8个降序特征值λ₁,λ₂,…,λ₈；

步骤b：第一处理模块PE1和第二处理模块PE2从存储器D_RAM依次读入8阶协方差矩阵Α对应的8个降序特征值λ₁,λ₂,…,λ₈，计算得到特征值对数lnλ₁,lnλ₂,…,lnλ₈，存入寄存器组1register_group1中；

步骤c：第三处理模块PE3从寄存器组2register_group2读入8阶协方差矩阵Α对应的8个降序特征值λ₁,λ₂,…,λ₈，计算得到特征值积的对数ln_acc₁,ln_acc₂,…,ln_acc₈，存入寄存器组0register_group0中；

步骤d：第二处理模块PE2从寄存器组1register_group1中读入特征值对数lnλ₁,lnλ₂,…lnλ₈，计算得到特征值和的对数ln_sum₁,ln_sum₂,…,ln_sum₈并存入寄存器组3register_group3中；

步骤e：第一处理模块PE1和第二处理模块PE2分别从寄存器组0register_group0和寄存器组3register_group3中读取特征值积的对数ln_acc₁,ln_acc₂,…,ln_acc₈和特征值和的对数ln_sum₁,ln_sum₂,…,ln_sum₈并进行计算，获得AIC序列AIC₁,AIC₂,…,AIC₈后传递给比较单元Comp_PE；

步骤f：比较单元Comp_PE对AIC序列AIC₁,AIC₂,…,AIC₈中的各个元素进行比较，获得信号源的数目估计值num。

本实施例中，设定阵元数N为8的均匀圆阵，阵元间距为半波长，系统采用Verilog_HDL进行描述，运算规则符合单精度ieee-754标准，仿真环境为Modelsim10.1a，并在Xilinx公司的xc5vlx110t-1ff1136上实现。另一方面，基于相同的算法和输入，在Matlab中采用双精度浮点进行软件模拟。本发明中信号子空间分解算法的FPGA设计的资源消耗及系统工作频率如表1所示：

表1FPGA设计的资源消耗及系统工作频率

设定三个非相干信号的方位角分别为-160.5°、10.10°、50.80°，俯仰角分别为50.20°、60.10°、70.70°，信噪比为10dB。为了体现硬件计算的正确性，本发明做了如下实验：对同等源数据进行上位机软件算法处理，软硬件实现的特征值分解计算结果数据对比如表2所示，信号源数目估计计算结果数据对比如表3所示，根据实验的结果可知，本发明可以有效地估计出信号源数。

表2特征值计算结果

序号	软件	硬件	相对误差
				1	5.314883E3	5.314881E3	4.618492E-7
2	3.641558E3	3.641557E3	5.901926E-8
				3	3.044003E3	3.044003E3	8.348386E-811 -->
4	-2.230000E-10	-2.240000E-10	4.484304E-3
				5	-3.666770E-7	3.665450E-7	3.599898E-4
6	-4.594620E-7	4.597480E-7	6.224671E-4
				7	7.745600E-7	7.739630E-7	7.707601E-4
8	7.987330E-7	7.990120E-7	3.493032E-4

表3信号源数目估计计算结果

序号	软件	硬件	相对误差
				1	4.754456E4	4.754457E4	3.812213E-7
2	4.733234E4	4.733235E4	2.442832E-7
				3	4.686936E4	4.686937E4	2.366951E-7
4	3.901296E1	3.900424E1	2.234656E-4
				5	4.800518E1	4.801527E1	2.101905E-4
6	5.500416E1	5.500636E1	4.004997E-5
				7	6.000396E1	6.001696E1	2.167183E-4
8	6.300000E1	6.300000E1	0

本发明与目前用不同硬件实现信号子空间分解算法的性能进行比较，结果如表4所示。郑洪、肖先赐的《MUSIC算法在高速并行处理机上的实现》中，利用多片DSP实现8X8矩阵的特征分解所需的时间约为20ms；徐德琛、刘志文、徐友根的《某侧向系统中MUSIC算法的FPGA实现》中，使用P-BLV阵列计算12X12矩阵的特征值分解需要0.187ms，所消耗的Slice数量为22661；本发明实现8X8矩阵的特征值分解只需要86.83us，实现信号源数估计仅需要24.66us，且仅需要消耗5670个Slice。

表4不同硬件实现方法性能比较

Claims (3)

1.一种分时复用硬件资源的信号子空间分解的FPGA实现模块，其特征是包括：自动排序的过关Jacobi算法单元和空间信号源数目估计计算单元；

所述自动排序的过关Jacobi算法单元读取空间信号源对应的N阶协方差矩阵Α＝[a_pq]_n×n并进行运算，获得N阶的对角矩阵Α_k；所述对角矩阵Α_k的主对角线上的元素为所述N阶协方差矩阵Α的特征值λ₁,λ₂,…,λ_p,…,λ_q,…,λ_n；且λ₁,λ₂,…,λ_p,…,λ_q,…,λ_n按降序排列；1≤p≤n；1≤q≤n；p≠q；

所述空间信号源数目估计计算单元将所述特征值λ₁,λ₂,…,λ_p,…,λ_q,…,λ_n作为源操作数并进行运算，获得2n个运算结果，包括n个特征值和的对数ln_sum₁,ln_sum₂,…,ln_sum_n及n个特征值积的对数ln_acc₁,ln_acc₂,…,ln_acc_n；再利用等价变换的AIC准则法对所述2n个运算结果进行运算，获得一个先递减后递增的AIC序列AIC₁,AIC₂,…,AIC_p,…,AIC_q,…,AIC_n；从所述AIC序列AIC₁,AIC₂,…,AIC_p,…,AIC_q,…,AIC_n中选出一个最小值，以最小值所对应的下标并减去1即为信号源的数目估计值num。

2.一种分时复用硬件资源的信号子空间分解的FPGA实现方法，其特征是按如下步骤进行：

步骤1、对接收到的空间信号源所对应的N阶协方差矩阵进行自动排序的过关Jacobi算法处理，获得N阶的对角矩阵A_k；

步骤1.1、选取一个单调递减且逐渐趋于零的数列a₁,a₂,…,a_i,…,a_k作为关值，并满足a₁＞a₂＞…a_i…＞a_k且 $\underset{k\&RightArrow;\mathrm{\&infin;}}{\lim }{a}_k=0,$ 并有 $a_i=a_{i-1}/n=\sqrt{2\underset{p=1}{\overset{n-1}{\&Sigma;}}\underset{q=p+1}{\overset{n}{\&Sigma;}}{a}_{pq}}/n;$ 1≤i≤k；

步骤1.2、初始化i＝1；

步骤1.3、初始化p＝1；

步骤1.4、初始化q＝1；

步骤1.5、判断p＝q是否成立，若成立，则将q+1赋值给q；并执行步骤1.6，若不成立，则直接执行步骤1.6；

步骤1.6、判断|a_pq|＞a_i是否成立，若成立，则对N阶协方差矩阵Α利用自动排序的过关Jacobi算法进行旋转变换，获得第i次旋转变换矩阵A_i；若不成立，则将q+1赋值给q；并返回步骤1.6执行，直到q＝n为止；

步骤1.7、将p+1赋值给p；并返回步骤1.5执行；直到p＝n为止；

步骤1.8、将i+1赋值给i；并返回步骤1.3执行；直到i＝k为止；从而获得N阶的对角矩阵Α_k；且所述对角矩阵Α_k的主对角线上的元素为所述N阶协方差矩阵Α的特征值λ₁，λ₂，…λ_p，…λ_q，…λ_n；且λ₁，λ₂，…λ_p，…λ_q，…λ_n按降序排列；

步骤2、利用等价变换的AIC准则估计信号源数目

步骤2.1、对式(1)所示的AIC表达式进行等价变换，获得等价变换的AIC表达式如式(2)所示：

AIC(num+1)＝2K(N-num)lnΛ(num)+2num(2N-num)(1)

$AIC(num+1)=-2(M-num)K\log 10\left(\right(\left(\underset{n=num+1}{\overset{N}{\&Pi;}}{\mathrm{\&lambda;}}_j^{(L/(M-num\left)\right)}\right)\left(\frac{\underset{n=num+1}{\overset{N}{\&Sigma;}}{\mathrm{\&lambda;}}_j}{M-num}\right))+2num(2M-unm)---\left(2\right)$

式(1)和式(2)中，K表示信号采样的快拍数，num表示待估计的信号源的个数，N为阵元数，Λ(num)为似然函数；λ_j表示第j个特征值，1≤j≤n；

步骤2.2、利用式(3)所示的泰勒级数展开式获得第j个特征值的对数ln(λ_j)：

$ln\left({\mathrm{\&lambda;}}_j\right)=ln\left(\frac{x+1}{1-x}\right)=2x(1+\frac{1}{3}{x}^2+\frac{1}{5}{x}^4+\frac{1}{7}{x}^6+\mathrm{...}+\frac{1}{2m+1}{x}^{2m}---\left(3\right)$

式(3)中，x表示操作数，且x∈[-1,1]；m表示泰勒级数展开式的项数；

步骤2.3、利用式(4)获得第j个特征值λ_j所对应的N个累加和对数ln_sum_j，从而获得n个特征值和的对数ln_sum₁,ln_sum₂,…,ln_sum_j,…,ln_sum_n：

$\ln \_{\mathrm{sum}}_j=\ln \left(\underset{m=j}{\overset{N}{\&Sigma;}}{\mathrm{\&lambda;}}_j\right)---\left(4\right)$

步骤2.4、利用式(5)获得第j个特征值λ_j所对应的N个乘积的对数ln_acc_j，从而获得n个特征值和的对数ln_acc₁,ln_acc₂,…,ln_acc_j,…,ln_acc_n：

$\ln \_{\mathrm{acc}}_j=\ln \left(\underset{m=j}{\overset{N}{\&Pi;}}{\mathrm{\&lambda;}}_j\right)---\left(5\right)$

步骤2.5、将n个特征值和的对数ln_sum₁,ln_sum₂,…,ln_sum_j,…,ln_sum_n以及n个特征值和的对数ln_acc₁,ln_acc₂,…,ln_acc_j,…,ln_acc_n代入式(2)中，从而获得一个先递减后递增的AIC序列AIC₁,AIC₂…,AIC_p…,AIC_q…,AIC_n；

步骤2.6、从所述AIC序列AIC₁,AIC₂…,AIC_p…,AIC_q…,AIC_n中选出一个最小值，以最小值所对应的下标并减去1即为信号源的数目估计值num。

3.一种分时复用硬件资源的信号子空间分解的FPGA实现电路，其特征是包括：第一处理模块PE1、第二处理模块PE2、第三处理模块PE3、寄存器组0register_group0、寄存器组1register_group1、寄存器组2register_group2、寄存器组3register_group3、寄存器组4register_group4、比较单元Comp_PE、存储器D_RAM；

所述存储器D_RAM存入空间信号源对应的N阶协方差矩阵Α＝[a_pq]_n×n；

所述寄存器组0register_group0从所述存储器D_RAM中读取N阶协方差矩阵Α的p行元素；

所述寄存器组1register_group1从所述存储器D_RAM中读取N阶协方差矩阵Α的q行元素；

所述寄存器组4register_group4从所述存储器D_RAM中读取N阶协方差矩阵Α的p列和q列元素；

所述比较单元Comp_PE从所述存储器D_RAM中依次读取N阶协方差矩阵Α中的非主对角线元素a_pq并与所设定的关值进行比较，若非主对角线元素a_pq的绝对值大于设定的关值，则利用第一处理模块PE1对N阶协方差矩阵中主对角线上的元素进行计算，获得自动排序的过关Jacobi算法中的cosθ和sinθ；再根据cosθ和sinθ、所述寄存器组0register_group0、寄存器组1register_group1、寄存器组4register_group4中存储的元素，当非主对角线元素a_pq的绝对值大于设定的关值时，利用第二处理模块PE2和所述第三处理模块PE3进行旋转变换；从而获得所述N阶协方差矩阵Α变换后的对角矩阵Α_k并传递给所述存储器D_RAM进行存储；所述变换后的对角矩阵Α_k的主对角线上的元素为按降序排列的特征值λ₁，λ₂，…λ_p，…λ_q，…λ_n；

所述第一处理模块PE1和第二处理模块PE2分别从所述存储器D_RAM中读取特征值λ₁，λ₂，…λ_p，…λ_q，…λ_n并进行计算，获得特征值对数lnλ₁，lnλ₂，…lnλ_p，…lnλ_q，…lnλ_n并存入所述寄存器组1register_group1中；

所述寄存器组2register_group2从所述存储器D_RAM中读取特征值λ₁，λ₂，…λ_p，…λ_q，…λ_n并传递给所述第三处理模块PE3；所述第三处理模块PE3对所述特征值λ₁，λ₂，…λ_p，…λ_q，…λ_n进行计算获得运算结果ln_acc₁,ln_acc₂,…ln_acc_n存入寄存器组0register_group0中；

所述第二处理模块PE2从所述寄存器组1register_group1所读取特征值对数lnλ₁，lnλ₂，…lnλ_p，…lnλ_q，…lnλ_n并进行计算，获得运算结果ln_sum₁,ln_sum₂,…ln_sum_n并存入寄存器组3register_group3中；

所述第一处理模块PE1和所述第二处理模块PE2分别从所述寄存器组0register_group0和寄存器组3register_group3中读取特征值积的对数ln_acc₁,ln_acc₂,…ln_acc_n和特征值和的对数ln_sum₁,ln_sum₂,…ln_sum_n并进行计算，获得AIC序列AIC₁,AIC₂…,AIC_p…,AIC_q…,AIC_n后传递给所述比较单元Comp_PE；

所述比较单元Comp_PE对所述AIC序列AIC₁,AIC₂…,AIC_p…,AIC_q…,AIC_n中的各个元素进行比较，获得信号源的数目估计值num。