CN1746697A - 一种芯片可实现的多信号分类算法 - Google Patents

Abstract

一种芯片可实现的多信号分类算法，属于通信技术领域，具体应用于阵列信号处理中的波达方向估计。包括以下步骤：采集数据，构造数据矩阵；计算协方差矩阵；特征值分解；构造噪声子空间；谱峰搜索。其本质是使用了本发明提出的改进坐标旋转算法来实现顺序Jacobi算法中的反正切运算和坐标旋转，而通常使用CORDIC算法需分别实现顺序Jacobi算法中的反正切运算和坐标旋转；同时改进了顺序Jacobi算法的迭代结构，把原来需要进行两次坐标旋转的计算复杂度降低为只需要进行一次坐标旋转。本发明可以提高算法并行性，大大提高运算速度。在芯片中进行EVD运算的时间是原来运行时间的1/3左右，而进行整个多信号分类算法时间降低为1/2左右。

Description

一种芯片可实现的多信号分类算法

技术领域

一种芯片可实现的多信号分类算法，属于通信技术领域，具体应用于阵列信号处理中的波达方向估计(Direction Of Arrival，以下简称DOA)。

背景技术

随着移动通信的蓬勃发展，用户数量迅速增加，频谱资源越来越紧张，如何利用现有频谱资源进一步扩展容量已成为移动通信发展的关键问题。智能天线技术可以产生多个空间定向波束，动态改变覆盖区域形状，使天线主波束对准用户信号到达方向，旁瓣或零陷对准干扰信号到达方向，从而能够降低系统干扰，提高系统容量和频谱效率。

实现智能天线的前提是首先要确定用户信号到达方向，DOA是其中重要的一种方法。现有的DOA方法有很多，如：多信号分类(MUltiple SIgnal Classification，以下简称MUSIC)算法、ESPRIT算法、最大似然估计算法等。MUSIC算法是一种综合考虑性能和复杂度后得到的一种比较优化的算法。现在将MUSIC算法简洁描述如下：

1.采集数据，构造数据矩阵：

当信号从远场(一般认为信号源与阵列天线距离大于信号源信号的波长10倍以上)入射到阵列天线时候，可以近似认为该信号到各天线阵元的方向角一致。常用的阵列天线是均匀的线阵或圆阵，如图1、2所示。

首先在阵列天线的N个阵元并行采集M个数据，形成数据矩阵A：

$A=\left[\begin{array}{cccccc}{x}_1\left(1\right)& x_1\left(2\right)& \·\·\·& x_1\left(j\right)& \·\·\·& x_1\left(M\right)\\ \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·\\ x_i\left(1\right)& x_i\left(2\right)& \·\·\·& x_i\left(j\right)& \·\·\·& x_i\left(M\right)\\ \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·\\ x_N\left(1\right)& x_N\left(2\right)& \·\·\·& x_N\left(j\right)& \·\·\·& x_N\left(M\right)\end{array}\right]---\left(1\right)$

其中x_i(j)中的i代表第i(i≤N)个阵元，j代表第j(j≤N)个采集数据。

2.计算协方差矩阵：

在得到数据矩阵A后，则可以计算协方差矩阵R：

R＝AA^T                                            (2)

从(2)可确定矩阵R是N×N对称实矩阵。

3.特征值分解(EVD)：

对协方差矩阵R进行特征值分解，对R的特征值和特征向量矩阵按照R的特征值大小排序，得到：

{λ₁，λ₂…λ_j…λ_N}其中λ₁≥λ₂…λ_j…≥λ_N≥0，

和对应特征向量矩阵

$V=\left[\begin{array}{cccccc}{v}_{11}& v_{12}& \·\·\·& v_{1j}& \·\·\·& v_{1N}\\ \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·\\ V_{i1}& v_{i2}& \·\·\·& v_{\mathrm{ij}}& \·\·\·& v_{\mathrm{iN}}\\ \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·\\ v_{N1}& v_{N2}& \·\·\·& v_{\mathrm{Nj}}& \·\·\·& v_{\mathrm{NN}}\end{array}\right]---\left(3\right)$

记矩阵V第j列矢量v_j为λ_j对应的特征向量，则V＝[V₁ V₂…V_j…V_N]。

4.构造噪声子空间：

对于D个信号入射的情况，R的特征值有如下性质：

${\mathrm{\λ}}_k=\left\{\begin{array}{cc}{\mathrm{\μ}}_k+{\mathrm{\σ}}_k& k=\mathrm{1,2}\·\·\·D\\ {\mathrm{\σ}}_k& k=D+1,\·\·\·N\end{array};\right.$

上式中的μ_k代表入射信号的特征值(μ₁≥μ₂…≥μ_D≥0)，σ_k代表噪声方差。一般说来，λ_D和λ_D+1在信噪比不恶劣的情况下，都有数量级上的差距，从而能够很好区分出信号个数D；

估计出D大小，从V中构造噪声子空间E_N ^H，其中 $E_N^H=[V_{D+1}\·\·\·V_N].$

5.谱峰搜索：

阵列天线的方向矩阵由信号入射方向θ_k(k＝1，2…D)、天线阵列几何参数(均匀圆阵、均匀线阵)等决定。设方向矩阵的列矢量是α(θ)(当θ一定时候，该值为常数)，算法证明，当θ＝θ_k时 $E_N^{H}a\left({\mathrm{\θ}}_k\right)=0.$

于是构造空间谱函数P(θ)，在θ∈[0，180]的范围内进行谱峰搜索：

$P\left(\mathrm{\θ}\right)=\frac{1}{{\left|\right|E_N^{H}a\left(\mathrm{\θ}\right)\left|\right|}^2}$

‖·‖表示矩阵的模2范数，P(θ)中D个峰值对应的D个θ就是待求的信号入射方向θ_k。

MUSIC算法第2步中的特征值分解的计算量占到了整个MUSIC算法的60％。其中，涉及特征值分解的算法有很多，如Jacobi算法、顺序Jacobi算法、QR算法等等。这些算法的特点无一例外的是运算繁杂，都存在非线性运算，不利于芯片实现。相对来说，顺序Jacobi算法是计算MUSIC算法中第2步的特征值分解最理想算法，并行计算时收敛速度最快，非线性计算相对较少，便于芯片实现。

顺序Jacobi算法求解特征值和特征向量的基本步骤如下：

设待特征值分解矩阵是R，特征向量矩阵V(矩阵V的初值是N×N单位矩阵)，

1、从R选择消去对象r_pq：

按照r₁₂，r₁₃，…r_1N；r₂₃，r₂₄…r_2N；……；r_(N-1)N矩阵元素所在的位置的先后顺序从R选择一个消去对象r_pq；当选择r_(N-1)N后，下次选择r₁₂并重复上述顺序重新进行选择；

2、计算反正切：

${\mathrm{\θ}}_{\mathrm{opt}}=\frac{1}{2}{\tan }^{-1}\left[\frac{{2r}_{\mathrm{pq}}}{r_{\mathrm{pp}}-r_{\mathrm{qq}}}\right]---\left(4\right)$

在确定第1步的r_pq后，按照r_pq从R中提取r_pp和r_qq，计算反正切得到旋转角θ_opt；

3、计算坐标旋转：

$R^{\′}=W_{\mathrm{pq}}^T{\mathrm{RW}}_{\mathrm{pq}}---\left(6\right)$

         V′＝W_pq ^TV                             (7)

在第2步计算得到旋转角θ_opt后，按照公式(5)构造W_pq ^T，然后按照公式(6)，对R做2次坐标旋转后得到矩阵R′；在计算R′同时，按照公式(7)，对V做1次坐标旋转就得到V′，具体计算时候可以用如下表达式

$\left[\begin{array}{cc}{r_{\mathrm{pp}}}^{\′}& {r_{\mathrm{pq}}}^{\′}\\ {r_{\mathrm{qp}}}^{\′}& {r_{\mathrm{qq}}}^{\′}\end{array}\right]=\left[\begin{array}{cc}\cos {\mathrm{\θ}}_{\mathrm{opt}}& {-\sin \mathrm{\θ}}_{\mathrm{opt}}\\ {\sin \mathrm{\θ}}_{\mathrm{opt}}& {\cos \mathrm{\θ}}_{\mathrm{opt}}\end{array}\right]\left[\begin{array}{cc}{r}_{\mathrm{pp}}& r_{\mathrm{pq}}\\ r_{\mathrm{qp}}& r_{\mathrm{qq}}\end{array}\right]\left[\begin{array}{cc}{\cos \mathrm{\θ}}_{\mathrm{opt}}& {\sin \mathrm{\θ}}_{\mathrm{opt}}\\ {-\sin \mathrm{\θ}}_{\mathrm{opt}}& {\cos \mathrm{\θ}}_{\mathrm{opt}}\end{array}\right]$

$\left[\begin{array}{c}{r_{\mathrm{pj}}}^{\′}\\ {r_{\mathrm{qj}}}^{\′}\end{array}\right]=\left[\begin{array}{cc}{\cos \mathrm{\θ}}_{\mathrm{opt}}& {-\sin \mathrm{\θ}}_{\mathrm{opt}}\\ {\sin \mathrm{\θ}}_{\mathrm{opt}}& {\cos \mathrm{\θ}}_{\mathrm{opt}}\end{array}\right]\left[\begin{array}{c}{r}_{\mathrm{pj}}\\ r_{\mathrm{qj}}\end{array}\right]j\≠p,q;$

$\left[\begin{array}{c}{v_{\mathrm{jp}}}^{\′}\\ {v_{\mathrm{jq}}}^{\′}\end{array}\right]=\left[\begin{array}{cc}{\cos \mathrm{\θ}}_{\mathrm{opt}}& {-\sin \mathrm{\θ}}_{\mathrm{opt}}\\ {\sin \mathrm{\θ}}_{\mathrm{opt}}& {\cos \mathrm{\θ}}_{\mathrm{opt}}\end{array}\right]\left[\begin{array}{c}{v}_{\mathrm{jp}}\\ v_{\mathrm{jq}}\end{array}\right]$ j任意值

4.判断R′

如果R′不是三角型矩阵，则令R＝R′、V＝V′后，返回第1步选择下一个消去对象，然后重复步骤1-4；如果R′已经转化为三角型，此时R′的主对角线元素就是特征值，V′就是特征向量矩阵。

5.排序

对R′中的特征值按从大到小顺序，对R′，V′进行对应调整，使得{λ₁，λ₂…λ_j…λ_N}其中λ₁≥λ₂…λ_j…≥λ_N≥0

V＝{v₁ v₂…v_j…v_N}(v_j是λ_j对应的特征向量)

使用上述步骤实现特征值分解时候，运算量主要集中在第2步的反正切和第3步的坐标旋转，因为采用了串行方式来实现反正切和2次坐标旋转，导致整个实现算法并行性很低，运算时间较长，无法满足智能天线在工程应用中的要求。

工程应用中实现反正切和坐标旋转常用的一种算法是可芯片实现的CORDIC算法，它可以在仅使用加法器和移位器的情况下在一个芯片上就实现这两种非线性运算。在使用CORDIC算法实现这两种运算的时候，其运算时间和运算原理基本一致，下面仅就坐标旋转的原理和步骤进行描述，反正切不再赘述。

CORDIC算法实现坐标旋转的原理可以解释如下：

已知

和旋转角θ，需要计算：

$\left[\begin{array}{c}{y}^{\′}\\ x^{\′}\end{array}\right]=\left[\begin{array}{cc}\cos \mathrm{\θ}& -\sin \mathrm{\θ}\\ \sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right]\left[\begin{array}{c}y\\ x\end{array}\right]---\left(8\right)$

首先定义CORDIC基本的角度集：

{α₁，α₂…α_k}＝{arctg2^-1，arctg2^-2…arctg2^-k}(k根据结果精度选择)和CORDIC符号集：

{ξ₁，ξ₂…ξ_k，ξ_i＝1或-1，i＝1，2…k}(k根据结果精度选择)，

CORDIC算法实现过程就是确定符号集{ξ₁，ξ₂…ξ_k，ξ_i＝1或-1，i＝1，2…k}，使之满足： $\mathrm{\θ}\≈\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\ξ}}_i{\mathrm{\α}}_i,$ 从而把坐标

旋转角θ近似等效为坐标

旋转1组CORDIC基本的角{ξ₁α₁，ξ₂α₂…ξ_kα_k}，所以式(8)可以近似等效

$\left[\begin{array}{c}{y}^{\′}\\ x^{\′}\end{array}\right]=\left[\begin{array}{cc}\cos {\mathrm{\ξ}}_k{\mathrm{\α}}_k& {-\sin \mathrm{\ξ}}_k{\mathrm{\α}}_k\\ \sin {\mathrm{\ξ}}_k{\mathrm{\α}}_k& \cos {\mathrm{\ξ}}_k{\mathrm{\α}}_k\end{array}\right]\·\·\·\left[\begin{array}{cc}\cos {\mathrm{\ξ}}_1{\mathrm{\α}}_1& -\sin {\mathrm{\ξ}}_1{\mathrm{\α}}_1\\ \sin {\mathrm{\ξ}}_1{\mathrm{\α}}_1& \cos {\mathrm{\ξ}}_1{\mathrm{\α}}_1\end{array}\right]\left[\begin{array}{c}y\\ x\end{array}\right]$

下面介绍CORDIC实现坐标旋转具体步骤：

把原坐标 和旋转角θ作为计算步骤初值，令i＝1，则步骤如下：

1.确定ξ_i的值

根据θ的正负，判断ξ_i的值。具体判定标准如下：如果θ大于0，令ξ_i＝1；反之θ小于0，令ξ_i＝-1。

2.计算

旋转i个CORDIC基本角之后的坐标

$\left[\begin{array}{c}y\\ x\end{array}\right]=\left[\begin{array}{cc}1& -{\mathrm{\ξ}}_i{2}^{-i}\\ {\mathrm{\ξ}}_i{2}^{-i}& 1\end{array}\right]\left[\begin{array}{c}y\\ x\end{array}\right]---\left(9\right)$

式(9)使用加法和移位就可以实现。

3.判定是否旋转完毕

如果i＝k，说明旋转完毕，转到第4步进行，否则计算θ＝θ-ξ_iα_i(α_i为固定常数arctg2^-i)然后i＝i+1，后回到第1步重新计算；

4.模值修正：

把步骤3计算出的坐标

乘以固定常数const即可 $(\mathrm{const}=\underset{i=1}{\overset{k}{\mathrm{\Π}}}\cos {\mathrm{\α}}_i).$ 此时输出值就是所求。

发明内容：

本发明的主要任务就是提供一种基于快速、芯片可实现的特征值算法的MUSIC算法来实现波达方向估计。

一种芯片可实现的多信号分类算法，其步骤如下：

步骤一：采集数据，构造数据矩阵

首先在阵列天线的N个阵元并行采集M个数据，形成数据矩阵A：

$A=\left[\begin{array}{cccccc}{x}_1\left(1\right)& x_1\left(2\right)& \·\·\·& x_1\left(j\right)& \·\·\·& x_1\left(M\right)\\ \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·\\ x_i\left(1\right)& x_i\left(2\right)& \·\·\·& x_i\left(j\right)& \·\·\·& x_i\left(M\right)\\ \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·\\ x_N\left(1\right)& x_N\left(2\right)& \·\·\·& x_N\left(j\right)& \·\·\·& x_N\left(M\right)\end{array}\right]$

其中x_i(j)中的i代表第i(i≤N)个阵元，j代表第j(j≤N)个采集数据；数据矩阵A可用存储器RAM存储。

步骤二：计算协方差矩阵：

在得到数据矩阵A后，则可以计算协方差矩阵R：

              R＝AA^T

协方差矩阵R是一个对称实矩阵，具体表达形式如下：

其中(p＜q≤N)

该步骤可用

个乘法累加器实现。

步骤三：特征值分解(EVD)：

已知协方差矩阵是R，设特征向量矩阵V的初值是N×N的单位矩阵。定义CORDIC基本的角度集：

{α₁，α₂…α_ORDER}＝{arctg2^-1，arctg2^-2…arctg2^-ORDER}和CORDIC符号集：

{ξ₁，ξ₂…ξ_ORDER，ξ_i＝1或-1，i＝1，2…ORDER}

其中，ORDER是CORDIC阶数，根据精度选择，是固定常数。

具体步骤如下：

1.从R选择消去对象r_pq：

按照r₁₂，r₁₃，…r_1N；r₂₃，r₂₄…r_2N；……；r_(N-1)N矩阵元素所在的位置的先后顺序从R选择一个消去对象r_pq；当选择r_(N-1)N后，下次选择r₁₂并重复上述顺序重新进行选择；本步骤可用基于状态机的数据选择器实现。

2.改进的坐标旋转：

(1)计算γ＝2*r_pq、β＝r_pp-r_qq；本步骤可用移位器和减法器实现；

(2)根据γ和β确定ξ_i。如果γ*β≥0，ξ_i＝1，否则ξ_i＝-1；

(3)计算γ′和β′，

$\left[\begin{array}{c}{\mathrm{\γ}}^{\′}\\ {\mathrm{\β}}^{\′}\end{array}\right]=\left[\begin{array}{cc}1-2^{-2i}& -{\mathrm{\ξ}}_i{2}^{-(i-1)}\\ {\mathrm{\ξ}}_i{2}^{-(i-1)}& 1-2^{-2i}\end{array}\right]\left[\begin{array}{c}\mathrm{\γ}\\ \mathrm{\β}\end{array}\right];$

本步骤可用移位器和加法器实现；

(4)旋转第i个CORDIC基本旋转角，

$\left[\begin{array}{c}{r_{\mathrm{pj}}}^{\′}\\ {r_{\mathrm{qj}}}^{\′}\end{array}\right]=\left[\begin{array}{cc}1& -{\mathrm{\ξ}}_i{2}^{-i}\\ {\mathrm{\ξ}}_i{2}^{-i}& 1\end{array}\right]\left[\begin{array}{c}{r}_{\mathrm{pj}}\\ r_{\mathrm{qj}}\end{array}\right]j\≠p,q;$

$\left[\begin{array}{c}{v_{\mathrm{jp}}}^{\′}\\ {v_{\mathrm{jq}}}^{\′}\end{array}\right]=\left[\begin{array}{cc}1& -{\mathrm{\ξ}}_i{2}^{-i}\\ {\mathrm{\ξ}}_i{2}^{-i}& 1\end{array}\right]\left[\begin{array}{c}{v}_{\mathrm{jp}}\\ v_{\mathrm{jq}}\end{array}\right]$ j任意值

本步骤可用移位器和加法器实现；

(5)判定是否旋转完所有基本旋转角：

如果i＝ORDER，则执行第(6)步；否则i+1，令

$\mathrm{\&gamma;}={\mathrm{\&gamma;}}^{\&prime;};\mathrm{\&beta;}={\mathrm{\&beta;}}^{\&prime;};\left[\begin{array}{c}{r}_{\mathrm{pj}}\\ r_{\mathrm{qj}}\end{array}\right]=\left[\begin{array}{c}{r_{\mathrm{pj}}}^{\&prime;}\\ {r_{\mathrm{qj}}}^{\&prime;}\end{array}\right](j\&NotEqual;p,q<mrow>\; <mtable>\; </mtable>\; )</mrow>;\left[\begin{array}{c}{v}_{\mathrm{jp}}\\ v_{\mathrm{jq}}\end{array}\right]=\left[\begin{array}{c}{v_{\mathrm{jp}}}^{\&prime;}\\ {v_{\mathrm{jq}}}^{\&prime;}\end{array}\right]$ (j是任意值)后，回到第(2)

步重新执行；本步骤可用比较器实现；

(6)模值修正

$\left[\begin{array}{c}{r_{\mathrm{pj}}}^{\&prime;}\\ {r_{\mathrm{qj}}}^{\&prime;}\end{array}\right]=\mathrm{const}*\left[\begin{array}{c}{r_{\mathrm{pj}}}^{\&prime;}\\ {r_{\mathrm{qj}}}^{\&prime;}\end{array}\right](j\&NotEqual;p,q),\left[\begin{array}{c}{v}_{\mathrm{jp}}\\ v_{\mathrm{jq}}\end{array}\right]=\mathrm{const}*\left[\begin{array}{c}{v_{\mathrm{jp}}}^{\&prime;}\\ {v_{\mathrm{jq}}}^{\&prime;}\end{array}\right]$ (j任意值)

(其中 $\mathrm{const}=\underset{i=1}{\overset{\mathrm{order}}{\mathrm{\&Pi;}}}\cos {\mathrm{\&alpha;}}_i,$ 在ORDER确定条件下，const近似等于常数)；本步骤可用移位器和加法器实现；

(7)在(6)计算的同时，通过查表(由CORDIC符号集{ξ₁，ξ₂…ξ_ORDER}映射到值tgθ_opt)和乘法器后求出r_pp′、r_qq′，具体公式如下：

r_pp′＝r_pp-r_pqtgθ_opt

r_qq′＝r_qq+r_pqtgθ_opt

而对r_pq′直接置0：

r_pq′＝0；

3.判断R′

如果R′不是三角型矩阵，则令R＝R′、V＝V′后，返回第1步选择下一个消去对象，然后重复步骤1-2；如果R′已经转化为三角型，此时R′的主对角线元素就是特征值，V′就是特征向量矩阵；本步骤可用比较器实现；

4.排序

对R′中的特征值按从大到小顺序，对R′，V′进行对应调整，使得{λ₁，λ₂…λ_j…λ_N}(λ₁≥λ₂…λ_j…≥λ_N)和V＝{v₁ v₂…v_j…v_N}，其中V_j是λ_j对应的特征向量；本步骤可用比较器和数据选择器实现；

步骤四.构造噪声子空间：

根据特征值分解后得到的特征值估计出D大小，然后从V中去构造噪声子空间E_N ^H，其中 $E_N^H=[v_{D+1}\&CenterDot;\&CenterDot;\&CenterDot;v_N];$ 本步骤可用比较器和数据选择器实现；

步骤五.谱峰搜索：

阵列天线的方向矩阵由信号入射方向θ_k(k＝1，2…D)、天线阵列几何位置(均匀圆阵、均匀线阵)等决定。设方向矩阵的列矢量是α(θ)(当θ一定时候，该值为常数)，算法证明，当θ＝θ_k时 $E_N^{H}a\left({\mathrm{\&theta;}}_k\right)=0.$

于是构造空间谱函数P(θ)，在θ∈[0，180]的范围内进行谱峰搜索：

$P\left(\mathrm{\&theta;}\right)=\frac{1}{{\left|\right|E_N^{H}a\left(\mathrm{\&theta;}\right)\left|\right|}^2}$

P(θ)中D个峰值对应的D个θ就是待求的信号入射方向θ_k。本步骤可用加法器、乘法器、移位器、比较器和数据选择器实现。

本发明的本质在于：本发明使用了由本发明提出的改进的坐标旋转算法来实现顺序Jacobi算法中的反正切运算和坐标旋转，而通常使用CORDIC算法需要分别实现顺序Jacobi算法中的反正切运算和坐标旋转(在本发明中不需要使用CORDIC反正切求出坐标旋转的旋转角)；同时改进了顺序Jacobi算法的迭代结构，把原来需要进行2次坐标旋转的计算复杂度降低为只需要进行1次坐标旋转。以上两点改进，把顺序Jacobi算法单次迭代中的3次CORDIC运算(2次CORDIC坐标旋转和1次CORDIC反正切运算)大约降低到1.2次CORDIC运算(1次改进的坐标旋转、1次乘法和1次查表)。

本发明正是基于上述两点改进，提出了本发明所述的一种芯片可实现的多信号分类算法，在不增加芯片复杂程度和保证同等精度要求的情况下，可以提高算法并行性，大大提高运算速度，减小运算时间。在芯片中实现EVD的运算时间降低到原来运算时间的1/3左右，而完成整个芯片可实现的多信号分类算法的运算时间降低到1/2左右。比如：对于4个阵元的阵列天线来说，在芯片中实现整个EVD运算时间约为5.76微秒，实现整个MUSIC算法的时间在20微秒以下；而现有的多信号分类算法在芯片实现中整个EVD运算时间约为16.7微秒以上，整个算法运算时间在30微秒以上。

附图说明：

图1：8阶均匀线阵示意图。

图2：8阶均匀圆阵示意图。

图3：MUSIC算法步骤框图。

图4：基于CORDIC算法的顺序Jacobi特征值分解流程图。

图5：本发明所述的基于改进坐标旋转算法的顺序Jacobi特征值分解流程图。

Claims (4)

1、一种芯片可实现的多信号分类算法，其步骤如下：

步骤一：采集数据，构造数据矩阵

首先在阵列天线的N个阵元并行采集M个数据，形成数据矩阵A：

$A=\left[\begin{array}{cccccc}{x}_1\left(1\right)& x_1\left(2\right)& \&CenterDot;\&CenterDot;\&CenterDot;& x_1\left(j\right)& \&CenterDot;\&CenterDot;\&CenterDot;& x_1\left(M\right)\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ x_i(11)& x_i\left(2\right)& \&CenterDot;\&CenterDot;\&CenterDot;& x_i\left(j\right)& \&CenterDot;\&CenterDot;\&CenterDot;& x_i\left(M\right)\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ x_N\left(2\right)& x_N\left(2\right)& \&CenterDot;\&CenterDot;\&CenterDot;& x_N\left(j\right)& \&CenterDot;\&CenterDot;\&CenterDot;& x_N\left(M\right)\end{array}\right]$

其中x_i(j)中的i代表第i(i≤N)个阵元，j代表第j(j≤N)个采集数据；数据矩阵A可用存储器RAM存储；

步骤二：计算协方差矩阵：

在得到数据矩阵A后，则可以计算协方差矩阵R：

                   R＝AA^T

协方差矩阵R是一个对称实矩阵，具体表达形式如下：

其中(p＜q≤N)

该步骤可用

个乘法累加器实现；

步骤三：特征值分解(EVD)：

已知协方差矩阵是R，设特征向量矩阵V初值是N×N的单位矩阵，定义CORDIC基本的角度集：

{α₁，α₂…α_ORDER}＝{arctg2^-1，arctg2^-2…arctg2^-ORDER}和CORDIC符号集：

{ξ₁，ξ₂…ξ_ORDER，ξ_i＝1或-1，i＝1，2…ORDER}

其中，ORDER是CORDIC阶数，根据精度选择，是固定常数；

具体步骤如下：

1.从R选择消去对象r_pq：

按照r₁₂，r₁₃，…r_1N；r₂₃，r₂₄…r_2N；……；r_(N-1)N矩阵元素所在的位置的先后顺序从R选择消去对象r_pq；当选择r_(N-1)N后，下次选择r₁₂并重复上述顺序重新进行选择；本步骤可用基于状态机的数据选择器实现；

2.改进的坐标旋转：

(1)计算γ＝2*r_pq、β＝r_pq-r_qq，令i＝1；本步骤可用移位器和减法器实现；

(2)根据γ和β确定ξ_i。如果γ*β≥0，ξ_i＝1，否则ξ_i＝-1；

(3)计算γ′和β′， $\left[\begin{array}{c}{\mathrm{\&gamma;}}^{\&prime;}\\ {\mathrm{\&beta;}}^{\&prime;}\end{array}\right]=\left[\begin{array}{cc}1-2^{-2i}& -{\mathrm{\&xi;}}_i{2}^{-(i-1)}\\ {{\mathrm{\&xi;}}_{i}2}^{-(i-1)}& 1-2^{-2i}\end{array}\right]\left[\begin{array}{c}\mathrm{\&gamma;}\\ \mathrm{\&beta;}\end{array}\right];$ 本步骤可用移位器和加法器实现；

(4)旋转第i个CORDIC基本旋转角，

$\left[\begin{array}{c}{r_{\mathrm{pj}}}^{\&prime;}\\ {r_{\mathrm{qj}}}^{\&prime;}\end{array}\right]=\left[\begin{array}{cc}1& -{\mathrm{\&xi;}}_i{2}^{-i}\\ {{\mathrm{\&xi;}}_{i}2}^{-i}& 1\end{array}\right]\left[\begin{array}{c}{r}_{\mathrm{pj}}\\ r_{\mathrm{qj}}\end{array}\right]---j\&NotEqual;p,q;$

$\left[\begin{array}{c}{v_{\mathrm{jp}}}^{\&prime;}\\ {r_{\mathrm{jq}}}^{\&prime;}\end{array}\right]=\left[\begin{array}{cc}1& -{\mathrm{\&xi;}}_i{2}^{-i}\\ {{\mathrm{\&xi;}}_{i}2}^{-i}& 1\end{array}\right]\left[\begin{array}{c}{r}_{\mathrm{jp}}\\ r_{\mathrm{jq}}\end{array}\right]$ j任意值；

本步骤可用移位器和加法器实现；

(5)判定是否旋转完所有基本旋转角：

如果i＝ORDER，则执行第(6)步；否则i+1，令

γ＝γ′；β＝β′； $\left[\begin{array}{c}{r}_{\mathrm{pj}}\\ r_{\mathrm{qj}}\end{array}\right]=\left[\begin{array}{c}{r_{\mathrm{pj}}}^{\&prime;}\\ {r_{\mathrm{qj}}}^{\&prime;}\end{array}\right](j\&NotEqual;p,q);\left[\begin{array}{c}{v}_{\mathrm{jp}}\\ v_{\mathrm{jq}}\end{array}\right]=\left[\begin{array}{c}{v_{\mathrm{jp}}}^{\&prime;}\\ {v_{\mathrm{jq}}}^{\&prime;}\end{array}\right]$ (j是任意值)后，回到第(2)

步重新执行；本步骤可用比较器实现；

(6)模值修正

$\left[\begin{array}{c}{r_{\mathrm{pj}}}^{\&prime;}\\ {r_{\mathrm{qj}}}^{\&prime;}\end{array}\right]=\mathrm{const}*\left[\begin{array}{c}{r_{\mathrm{pj}}}^{\&prime;}\\ {r_{\mathrm{qj}}}^{\&prime;}\end{array}\right](j\&NotEqual;p,q),\left[\begin{array}{c}{v}_{\mathrm{jp}}\\ v_{\mathrm{jq}}\end{array}\right]=\mathrm{const}*\left[\begin{array}{c}{v_{\mathrm{jp}}}^{\&prime;}\\ {v_{\mathrm{jq}}}^{\&prime;}\end{array}\right]$ (j任意值)

(其中 $\mathrm{const}=\underset{i=1}{\overset{\mathrm{order}}{\mathrm{\&Pi;}}}\cos {\mathrm{\&alpha;}}_i,$ 在ORDER确定条件下，const近似等于常数)；本步骤可用移位器和加法器实现；

(7)在(6)计算的同时，通过查表(由确定后的CORDIC符号集{ξ₁，ξ₂…ξ_ORDER，}映射到tgθ_opt)和乘法器后求出r_pp′、r_qq′，具体公式如下：

r_pp′＝r_pp-r_pqtgθ_opt

r_qq′＝r_qq+r_pqtgθ_opt

而对r_pq′直接置0：

r_pq′＝0；

3.判断R′

如果R′不是三角型矩阵，则令R＝R′、V＝V′后，返回第1步选择下一个消去对象，然后重复步骤1-2；如果R′已经转化为三角型，此时R′的主对角线元素就是特征值，V′就是特征向量矩阵；本步骤可用比较器实现；

4.排序

对R′中的特征值按从大到小顺序，对R′，V′进行对应调整，使得{λ₁，λ₂…λ_j…λ_N}(λ₁≥λ₂…λ_j…≥λ_N)和V＝{v₁ v₂…v_j…v_N}，其中v_j是λ_j对应的特征向量；本步骤可用比较器和数据选择器实现；

步骤四.构造噪声子空间：

根据特征值分解后得到的特征值估计出D大小，然后从V中构造噪声子空间E_N ^H，其中 $E_N^H=[v_{D+1}\&CenterDot;\&CenterDot;\&CenterDot;v_N];$ 本步骤可用比较器和数据选择器实现；

步骤五.谱峰搜索：

阵列天线的方向矩阵由信号入射方向θ_k(k＝1，2…D)、天线阵列几何位置(均匀圆阵、均匀线阵)等决定。设方向矩阵的列矢量是a(θ)(当θ一定时候，该值为常数)，算法证明，当θ＝θ_k时 $E_N^{H}a\left({\mathrm{\&theta;}}_k\right)=0;$

于是构造空间谱函数P(θ)，在θ∈[0，180]的范围内进行谱峰搜索：

$P\left(\mathrm{\&theta;}\right)=\frac{1}{{\left|\right|E_N^{H}a\left(\mathrm{\&theta;}\right)\left|\right|}^2}$

P(θ)中D个峰值对应的D个θ就是待求的信号入射方向θ_k；本步骤可用加法器、乘法器、移位器、比较器和数据选择器实现。

Images