本发明公开了一种基于角度‑空间频率域快速傅里叶变换的互质阵列波达方向估计方法，主要解决现有技术计算复杂度高的问题。其实现步骤是：接收端架构互质阵列；利用互质阵列接收入射信号并建模；根据互质阵列接收信号推导对应于增广虚拟均匀线性阵列的二阶等价虚拟信号；定义角度‑空间频率并用其描述虚拟均匀线阵的二阶等价虚拟信号；对用角度‑空间频率描述的二阶等价虚拟信号进行预加窗操作；对加窗预处理后的二阶等价虚拟信号进行补零操作；对补零后的二阶等价虚拟信号进行快速傅里叶变换操作，并构建空间谱；根据所构建空间谱进行波达方向估计。本发明在提升信号波达方向估计自由度性能的同时降低了计算复杂度，易于在实际系统上实现。
1.一种基于角度-空间频率域快速傅里叶变换的互质阵列波达方向估计方法，其特征在于，包含以下步骤：(1)接收端使用2M+N-1个阵元，并按照互质阵列结构进行架构；其中M与N为互质整数；
(2)有L个来自θ＝[θ1,θ2,…,θL]T方向的远场窄带非相干信号源，[·]T表示转置操作，利用步骤(1)中构建的互质阵列对入射信号进行接收，则在t时刻的互质阵列接收信号x(t)可建模为：其中，x(t)为(2M+N-1)×1维向量，sl(t)为第l个入射信号的波形，n(t)为与各信号源相互独立的噪声分量，a(θl)为对应于θl方向信号源的互质阵列导引向量，可表示为：其中，μi,i＝1,2,3,…,2M+N-1表示互质阵列中第i个物理阵元的实际位置，且首个物理阵元的位置为μ1＝0，λ为入射窄带信号的波长，j为虚数单位；
根据互质阵列接收信号x(t)构建协方差矩阵：使用T个采样快拍得到的采样协方差矩阵 对理论协方差矩阵Rx进行近似替代；
(3)根据互质阵列接收信号推导对应于增广虚拟均匀线性阵列的二阶等价虚拟信号：通过向量化采样协方差矩阵 得到虚拟阵列等价虚拟信号y：
1)2×L维导引矩阵， 为包含L个入射信号源功率的向量， 为噪声功率， I2M+N-1表示(2M+N-1)×(2M+N-1)维单位矩阵，vec(·)表示向量化操作，即把矩阵中的各列依次堆叠以形成一个新的向量，(·)*表示共轭操作， 表示克罗内克积；
SD＝{±(Mn-Nm)d,0≤n≤N-1,0≤m≤2M-1}，其中d为入射窄带信号波长的一半，即
选取非均匀虚拟阵列SD中最大连续部分的虚拟阵元，组成一个包含2V+1个虚拟阵元的虚拟均匀线阵SV＝{-Vd,-(V-1)d,…,0,…,(V-1)d,Vd}，V＝MN+M-1；
从向量y中选择对应于SV中各虚拟阵元位置的等价虚拟信号，组成虚拟均匀线阵对应的二阶等价虚拟信号zθ，可表示为：其中，B(θ)＝[b(θ1),b(θ2),…b(θL)]，其第l列为对应于
θl方向信号源的虚拟均匀线阵SV的导引向量，e是由从 中选择对应于虚拟均匀线阵SV的元素组成的向量；
(4)定义角度-空间频率并用其描述虚拟均匀线阵的二阶等价虚拟信号：定义角度-空间频率为空间中来自方向θ的窄带信号在相邻阵元间的传播距离差范围内所传播的信号周期数；步骤(3)中虚拟均匀线阵的二阶等价虚拟信号zθ可在角度-空间频率域中等价表示为：其中，B(ξ)＝[b(ξ1),b(ξ2),…b(ξL)]，其第l列ξ＝[ξ1,ξ2,…,ξL]T对
(5)对用角度-空间频率描述的二阶等价虚拟信号zξ进行预加窗操作，得到预加窗后的二阶等价虚拟信号 可表示为：其中，w表示窗函数组成的与zξ维度相同的向量，表示哈达玛(Hadamard)积；
(6)对预加窗处理后的二阶等价虚拟信号 进行补零操作，使补零后的向量中元素的为K，其中K满足2的整数次幂，从而得到补零后的等价虚拟信号(7)对补零后的二阶等价虚拟信号 进行快速傅里叶变换操作，并构建空间谱：通过对补零后的二阶等价虚拟信号 进行快速傅里叶变换得到K×1维空间响应构建一个空间谱，该谱的横轴表示角度θ，其与空间响应 的第k个元素的关系可表示为：其中，k＝0,1,…,K-1，arccos(·)为反余弦函数，h为保证 满足反余弦函数的定义域的系数，当 时，h＝1，当 时，h＝0；该谱的纵轴表示空间响应 中第k个元素的模P(k)；
(8)根据所构建的空间谱进行波达方向估计：对步骤(7)中所构建的空间谱进行谱峰搜索操作，将幅度最大的前L个峰值对应的角度，作为L个入射信号的波达方向估计。
2.根据权利要求1所述的基于角度-空间频率域快速傅里叶变换的互质阵列波达方向估计方法，其特征在于：步骤(1)所述的互质阵列结构可具体描述为：首先，选取一组互质的整数M、N，构造一对稀疏均匀线性子阵列；第一个子阵列包含2M个间距为Nd的阵元，其位置为0,Nd,…,(2M-1)Nd；第二个子阵列包含N个间距为Md的阵元，其位置为0,Md,…,(N-1)Md；
之后将两个子阵列按照首个阵元重合的方式进行组合，得到包含2M+N-1个物理阵元的非均匀互质阵列。
3.根据权利要求1所述的基于角度-空间频率域快速傅里叶变换的互质阵列波达方向估计方法，其特征在于：步骤(3)中，若y中存在多个不同元素与SV中同一虚拟阵元位置对应的情况，则选择这些元素中的任意一个作为组成向量zθ的元素。
4.根据权利要求1所述的基于角度-空间频率域快速傅里叶变换的互质阵列波达方向估计方法，其特征在于：步骤(4)中的角度-空间频率是一种与信号入射角度有关的频率，通过如下方式定义：空间中来自方向θ的窄带信号到达间距为d的相邻阵元会产生一段传播距离差u，可表示为：u＝dcosθ
在阵元间距固定的情况下，该传播距离差u仅随着入射信号角度θ变化，因此得到步骤(4)中角度-空间频率的定义；角度-空间频率ξ与入射信号角度θ的关系可表示为：
5.根据权利要求1所述的基于角度-空间频率域快速傅里叶变换的互质阵列波达方向估计方法，其特征在于：步骤(5)中所采用的预加窗处理可以采用凯瑟(Kaiser)窗，其窗函数可表示为：其中，I0(·)为零阶修正贝塞尔函数，β为一个非负系数用于决定窗函数的形状，v＝0,
1,2,…,2V，则窗函数组成的向量w与向量zξ具有相同的维度，可表示为：w＝[w(0),w(1),…,w(2V)]T。
6.根据权利要求1所述的基于角度-空间频率域快速傅里叶变换的互质阵列波达方向估计方法，其特征在于：所述步骤(7)中，通过快速傅里叶变换得到补零后的二阶等价虚拟信号 的空间响应 可表示为：其中， 表示快速傅里叶变换运算，FK可表示为：
7.根据权利要求1所述的基于角度-空间频率域快速傅里叶变换的互质阵列波达方向估计方法，其特征在于：所述步骤(7)中，所构建的空间谱反映了空间中各角度上的响应幅度，其中存在对应于L个入射信号的L个峰值。
8.根据权利要求6所述的基于角度-空间频率域快速傅里叶变换的互质阵列波达方向估计方法，其特征在于：步骤(7)中，空间谱中存在对应于L个入射信号的L个峰值，该结论是通过如下方式得到的：用于构建空间谱的关系P(k)具体可表示为：其中， 为第l个信号的功率，δ(·)表示冲激函数，r为整数，用于表示该冲激函数为周期性冲激串序列；根据冲激函数的性质可知，仅当 时，P(k)为峰值；由于L个入射信号具有不同的角度-空间频率，每个角度-空间频率能且仅能在空间谱中引起一个峰值，因此空间谱中存在对应于L个入射信号的L个峰值。
基于角度-空间频率域快速傅里叶变换的互质阵列波达方向
估计方法
技术领域
[0001] 本发明属于信号处理技术领域，尤其涉及对雷达信号、声学信号及电磁信号的波达方向估计，具体是一种基于角度-空间频率域快速傅里叶变换的互质阵列波达方向估计方法，可用于无源定位和目标探测。
背景技术
[0002] 波达方向(Direction-of-Arrival,DOA)估计是阵列信号处理领域的一个基本问题，它是指利用传感器阵列接收信号，并通过一系列信号处理方法对接收信号进行统计处理，从而获得信号中所包含的波达方向信息，在雷达、声呐、语音、无线通信等领域均有着广泛的应用。
[0003] DOA估计方法的自由度是指其能够估计的入射信号源个数。传统的DOA估计方法通常采用均匀线性阵列进行信号的接收和建模，但基于均匀线性阵列的DOA估计方法的自由度受限于物理阵元个数。当入射的信号源个数大于阵列中物理阵元个数时，基于均匀线性阵列的DOA估计方法将无法得到有效的估计结果。互质阵列是一种具有系统化结构的非均匀稀疏阵列，可以突破传统均匀线性阵列自由度的限制，实现在物理阵元数一定的情况下自由度性能的提升，因此近年来得到学术界的广泛关注。基于互质阵列的DOA估计方法的原理为利用质数的性质将互质阵列推导至虚拟域，并构建一个对应于虚拟均匀线性阵列的二阶等价虚拟信号用于DOA估计。该虚拟均匀线性阵列中虚拟阵元数量大于物理阵元数量，因此实现了自由度的提升。
[0004] 现有基于互质阵列的DOA估计方法大多在互质阵列二阶等价虚拟信号的基础上进行较为复杂的运算，包括求逆、特征值分解等复杂的矩阵运算，以及凸优化问题的设计与求解等过程，这些运算过程导致了较高的计算复杂度，在实时性要求较高的应用场景下面临一定的挑战，且在实际系统中的硬件实现较为困难。
发明内容
[0005] 本发明的目的在于针对上述现有技术存在的不足，提出一种基于角度-空间频率域快速傅里叶变换的互质阵列波达方向估计方法，基于采用角度-空间频率描述的互质阵列二阶等价虚拟信号进行快速傅里叶变换操作，实现DOA估计自由度性能提升，且本发明所提方法计算复杂度低，易于在实际系统中硬件实现。本发明的目的是通过以下技术方案来实现的：一种基于角度-空间频率域快速傅里叶变换的互质阵列波达方向估计方法，包含以下步骤：
[0006] (1)接收端使用2M+N-1个阵元，并按照互质阵列结构进行架构；其中M与N为互质整数；
[0007] (2)假设有L个来自θ＝[θ1,θ2,…,θL]T方向的远场窄带非相干信号源，[·]T表示转置操作，利用步骤(1)中构建的互质阵列对入射信号进行接收，则在t时刻的互质阵列接收信号x(t)可建模为：
[0008]
[0009] 其中，x(t)为(2M+N-1)×1维向量，sl(t)为第l个入射信号的波形，n(t)为与各信号源相互独立的噪声分量，a(θl)为对应于θl方向信号源的互质阵列导引向量，可表示为[0010]
[0011] 其中，μi,i＝1,2,3,…,2M+N-1表示互质阵列中第i个物理阵元的实际位置，且首个物理阵元的位置为μ1＝0，λ为入射窄带信号的波长，j为虚数单位；
[0012] 根据互质阵列接收信号x(t)构建协方差矩阵：使用T个采样快拍得到的采样协方差矩阵 对理论协方差矩阵Rx进行近似替代；
[0013] (3)根据互质阵列接收信号推导对应于增广虚拟均匀线性阵列的二阶等价虚拟信号：通过向量化采样协方差矩阵 得到虚拟阵列等价虚拟信号y：
[0014]
[0015] 其中， 为(2M+N-1)2×L维导引矩阵， 为包含L个入射信号源功率的向量， 为
噪声功率，i＝vec(I2M+N-1)，I2M+N-1表示(2M+N-1)×(2M+N-1)维单位矩阵，vec(·)表示向量化操作，即把矩阵中的各列依次堆叠以形成一个新的向量，(·)*表示共轭操作， 表示克罗内克积；
[0016] 向量y对应的非均匀虚拟阵列SD表示为：
[0017] SD＝{±(Mn-Nm)d,0≤n≤N-1,0≤m≤2M-1}，
[0018] 其中d为入射窄带信号波长的一半，即
[0019] 选取非均匀虚拟阵列SD中最大连续部分的虚拟阵元，组成一个包含2V+1个虚拟阵元的虚拟均匀线阵SV＝{-Vd,-(V-1)d,…,0,…,(V-1)d,Vd}，V＝MN+M-1；
[0020] 从向量y中选择对应于SV中各虚拟阵元位置的等价虚拟信号，组成虚拟均匀线阵对应的二阶等价虚拟信号zθ，可表示为：
[0021]
[0022] 其中，B(θ)＝[b(θ1),b(θ2),…b(θL)]，其第l列
[0023]为对应于θl方向信号源的虚拟均匀线阵SV的导引向量，e是由从i中选择对应于虚拟均匀线阵SV的元素组成的向量；
[0024] (4)定义角度-空间频率并用其描述虚拟均匀线阵的二阶等价虚拟信号：定义角度-空间频率为空间中来自方向θ的窄带信号在相邻阵元间的传播距离差范围内所传播的信号周期数；步骤(3)中虚拟均匀线阵的二阶等价虚拟信号zθ可在角度-空间频率域中等价表示为：
[0025]
[0026] 其中，B(ξ)＝[b(ξ1),b(ξ2),…b(ξL)]，其第l列
[0027]对应于θ中所包含的L个角度的角度-空间频率；
[0028] (5)对用角度-空间频率描述的二阶等价虚拟信号zξ进行预加窗操作，得到预加窗后的二阶等价虚拟信号 可表示为：
[0029]
[0030] 其中，w表示窗函数组成的与zξ维度相同的向量， 表示哈达玛(Hadamard)积；
[0031] (6)对预加窗处理后的二阶等价虚拟信号 进行补零操作，使补零后的向量中元素的为K，其中K满足2的整数次幂，从而得到补零后的等价虚拟信号
[0032] (7)对补零后的二阶等价虚拟信号 进行快速傅里叶变换操作，并构建空间谱：通过对补零后的二阶等价虚拟信号 进行快速傅里叶变换得到K×1维空间响应
[0033] 构建一个空间谱，该谱的横轴表示角度θ，其与空间响应 的第k个元素的关系可表示为：
[0034]
[0035] 其中，k＝0,1,…,K-1，arccos(·)为反余弦函数，h为保证 满足反余弦函数的定义域的系数，当 时，h＝1，当 时，h＝0；该谱的纵轴表示空间响应 中第k个元素的模P(k)；
[0036] (8)根据所构建的空间谱进行波达方向估计：对步骤(7)中所构建的空间谱进行谱峰搜索操作，将幅度最大的前L个峰值对应的角度，作为L个入射信号的波达方向估计。
[0037] 进一步地，步骤(1)所述的互质阵列结构可具体描述为：首先，选取一组互质的整数M、N，构造一对稀疏均匀线性子阵列。第一个子阵列包含2M个间距为Nd的阵元，其位置为0,Nd,…,(2M-1)Nd；第二个子阵列包含N个间距为Md的阵元，其位置为0,Md,…,(N-1)Md。之后将两个子阵列按照首个阵元重合的方式进行组合，得到包含2M+N-1个物理阵元的非均匀互质阵列。
[0038] 进一步地，步骤(3)中，若y中存在多个不同元素与SV中同一虚拟阵元位置对应的情况，则选择这些元素中的任意一个作为组成向量zθ的元素。
[0039] 进一步地，步骤(4)中的角度-空间频率是一种与信号入射角度有关的频率，通过如下方式定义：空间中来自方向θ的窄带信号到达间距为d的相邻阵元会产生一段传播距离差u，可表示为：
[0040] u＝dcosθ。
[0041] 在阵元间距固定的情况下，该传播距离差u仅随着入射信号角度θ变化，因此得到步骤(4)中角度-空间频率的定义；角度-空间频率ξ与入射信号角度θ的关系可表示为：
[0042]
[0043] 进一步地，步骤(5)中所采用的预加窗处理可以采用凯瑟(Kaiser)窗，其窗函数可表示为：
[0044]
[0045] 其中，I0(·)为零阶修正贝塞尔函数，β为一个非负系数用于决定窗函数的形状，v＝0,1,2,…,2V，则窗函数组成的向量w与向量zξ具有相同的维度，可表示为：
[0046] w＝[w(0),w(1),…,w(2V)]T。
[0047] 进一步地，所述步骤(7)中，通过快速傅里叶变换得到补零后的二阶等价虚拟信号的空间响应 可表示为：
[0048]
[0049] 其中， 表示快速傅里叶变换运算，FK可表示为：
[0050]
[0051] 所得的空间响应 是K×1维向量。
[0052] 进一步地，所述步骤(7)中，所构建的空间谱反映了空间中各角度上的响应幅度，其中存在对应于L个入射信号的L个峰值。
[0053] 进一步地，所述步骤(7)中，空间谱中存在对应于L个入射信号的L个峰值，该结论是通过如下方式得到的：用于构建空间谱的关系P(k)具体可表示为：
[0054]
[0055] 其中， 为第l个信号的功率，δ(·)表示冲激函数，r为整数，用于表示该冲激函数为周期性冲激串序列；根据冲激函数的性质可知，仅当 时，P(k)为峰值；由于L个入射信号具有不同的角度-空间频率，每个角度-空间频率能且仅能在空间谱中引起一个峰值，因此空间谱中存在对应于L个入射信号的L个峰值。
[0056] 本发明与现有技术相比具有以下优点：
[0057] (1)本发明所提方法基于采用角度-空间频率描述的互质阵列二阶等价虚拟信号进行快速傅里叶变换操作，获得空间响应并基于此构建一个空间谱，通过对所构建空间谱的谱峰搜索过程得到波达方向估计结果，避免了传统波达方向估计方法中常用的凸优化问题设计求解以及矩阵求逆、矩阵特征值分解等复杂计算，在保证互质阵列波达方向估计性能的基础上，有效地降低波达方向估计的计算复杂度，使其更适用于对波达方向估计结果实时性要求较高的应用场景。
[0058] (2)本发明中快速傅里叶变换过程可以仅通过加法及乘法完全实现，有效地降低了计算复杂度，更好地满足了实时性要求较高的应用需求，同时易于在实际系统中硬件实现。
附图说明
[0059] 图1是本发明的方法总体流程框图；
[0060] 图2是本发明中组成互质阵列的一对稀疏均匀子阵列结构示意图；
[0061] 图3是本发明中互质阵列的结构示意图；
[0062] 图4是用于体现本发明所提方法自由度性能的空间谱示意图。
具体实施方式
[0063] 以下参照附图，对本发明的技术方案和效果作进一步的详细说明。
[0064] 互质阵列由于可以对其二阶等价虚拟阵列信号进行统计处理，实现在物理阵元数一定的情况下自由度性能的提升，近年来得到学术界的广泛关注。在实际系统应用中，现有的DOA估计方法计算复杂度较高，难以满足实时性要求较高的应用场景，且复杂的运算过程在实际系统中的硬件实现存在一定的困难。针对以上问题，本发明提供了一种基于角度-空间频率域快速傅里叶变换的互质阵列波达方向估计方法，参照图1，本发明的实现步骤如下：
[0065] 步骤一：在接收端使用2M+N-1个实际阵元构建互质阵列。首先，选取一组互质的整数M、N，构造一对稀疏均匀线性子阵列。第一个子阵列包含2M个间距为Nd的阵元，其位置为0,Nd,…,(2M-1)Nd；第二个子阵列包含N个间距为Md的阵元，其位置为0,Md,…,(N-1)Md，其中单位间距d为入射窄带信号的半波长即d＝λ/2。之后将两个子阵列按照首个阵元重合的方式进行组合，得到包含2M+N-1个实际阵元的非均匀互质阵列。
[0066] 步骤二：采用互质阵列接收信号并对信号建模。假设有L个来自θ＝[θ1,θ2,…,θL]T方向的远场窄带非相干信号源，[·]T表示转置操作，利用步骤一中构建的互质阵列对入射信号进行接收，则在t时刻的互质阵列接收信号x(t)可建模为：
[0067]
[0068] 其中，x(t)为(2M+N-1)×1维向量，sl(t)为第l个入射信号的波形，n(t)为与各信号源相互独立的噪声分量，a(θl)为对应于θl方向信号源的互质阵列导引向量，可表示为[0069]
[0070] 其中，μi,i＝1,2,3,…,2M+N-1表示互质阵列中第i个物理阵元的实际位置，且首个物理阵元的位置为μ1＝0，λ为入射窄带信号的波长，j为虚数单位。实际中，使用T个采样快拍得到的采样协方差矩阵 对理论协方差矩阵Rx进行近似替代， 和Rx可分别表示为[0071]
[0072] Rx＝E[x(t)xH(t)]，
[0073] 其中，(·)H表示共轭转置，E[·]表示数学期望。
[0074] 步骤三：根据互质阵列接收信号推导对应于增广虚拟均匀线性阵列的二阶等价虚拟信号。虚拟阵列等价虚拟信号y可以通过向量化步骤二中的采样协方差矩阵 得到：
[0075]
[0076] 其中， 为(2M+N-1)2×L维虚拟阵列导引矩阵， 为包含L个入射信号源功率
的向量， 为噪声功率，i＝vec(I2M+N-1)，I2M+N-1表示(2M+N-1)×(2M+N-1)维单位矩阵，vec(·)表示向量化操作，即把矩阵中的各列依次堆叠以形成一个新的向量，(·)*表示共轭操作， 表示克罗内克积。向量y对应的非均匀虚拟阵列SD可表示为：
[0077] SD＝{±(Mn-Nm)d,0≤n≤N-1,0≤m≤2M-1}，
[0078] 其中d为入射窄带信号波长的一半，即 选取非均匀虚拟阵列SD中最大连续部分的虚拟阵元，组成一个包含2V+1个虚拟阵元的虚拟均匀线阵SV＝{-Vd,-(V-1)d,…,0,…,(V-1)d,Vd}，其中V＝MN+M-1。从向量y中选择对应于SV中各虚拟阵元位置的等价虚拟信号，组成虚拟均匀线阵对应的二阶等价虚拟信号zθ，可表示为：
[0079]
[0080] 其中，B(θ)＝[b(θ1),b(θ2),…b(θL)]，其第l列为对应于
θl方向信号源的虚拟均匀线阵导引向量，e是由从i中选择对应于SV中阵元的元素组成的向量。若y中存在多个不同元素与SV中同一虚拟阵元位置对应的情况，则选择这些元素中的任意一个作为组成向量zθ的元素。
[0081] 步骤四：定义角度-空间频率并用其描述虚拟均匀线阵的二阶等价虚拟信号。定义角度-空间频率为空间中来自方向θ的窄带信号在相邻阵元间的传播距离差范围内所传播的信号周期数。具体地，空间中来自方向θ的窄带信号到达间距为d的相邻阵元会产生一段传播距离差，可表示为：
[0082] u＝dcosθ,
[0083] 在阵元间距一定的情况下，该传播距离差随着入射角度θ变化，角度-空间频率定义为在该传播距离差上窄带信号传播的周期数，即：
[0084]
[0085] 相应地，步骤三中虚拟均匀线阵的二阶等价虚拟信号zθ可在角度-空间频率域中等价表示为：
[0086]
[0087] 其中ξ＝[ξ1,ξ2,…,ξL]T对应于θ中所包含的L个角度的角度-空间频率。
[0088] 步骤五：对用角度-空间频率描述的二阶等价虚拟信号进行预加窗操作。
[0089] 对步骤四中采用角度-空间频率表示的二阶等价虚拟信号zξ进行预加窗操作。以采用凯瑟(Kaiser)窗为例，其窗函数w(k)可表示为：
[0090]
[0091] 其中，I0(·)为零阶修正贝塞尔函数，β为一个非负系数用于决定窗函数的形状，v＝0,1,2,…,2V。在该窗函数作用下，得到加窗后的二阶等价虚拟信号 可表示为：
[0092]
[0093] 其中，w＝[w(0),w(1),…,w(2v)]T表示窗函数组成的向量， 表示哈达玛(Hadamard)积。
[0094] 步骤六：对预加窗处理后的二阶等价虚拟信号进行补零操作。根据快速傅里叶变换的性质，对待变换序列的末尾补充一定数量的零，可以在不改变谱分辨率的情况下增大谱密度，因此为保证谱的足够多细节，在加窗后的等价虚拟信号向量末尾进行补零操作，使得其向量中有K个元素，其中K满足2的整数次幂。得到补零后的等价虚拟信号 可表示为：
[0095]
[0096] 其中，0表示满足要求的适当长度的零向量。
[0097] 步骤七：对补零后的二阶等价虚拟信号进行快速傅里叶变换操作，构建空间谱。通过对补零后的二阶等价虚拟信号 进行快速傅里叶变换得到K×1维空间响应 可表示为：
[0098]
[0099] 其中， 表示快速傅里叶变换运算，FK可表示为：
[0100]
[0101] 构建一个空间谱，该谱的横轴表示角度θ，其与空间响应第k个元素的关系可表示为：
[0102]
[0103] 其中，k＝0,1,…,K-1，arccos(·)为反余弦函数，h为保证 满足反余弦函数的定义域的系数，当 时，h＝1，当 时，h＝0；该谱的纵轴表示空间响应向量中第k个元素的模P(k)，可表示为：
[0104]
[0105] 其中，[·]k表示向量中第k个元素，|·|表示复数的模。具体地，P(k)可表示为：
[0106]
[0107] 其中， 为第l个信号的功率，δ(·)表示冲激函数，r为整数，用于表示该冲激函数为周期性冲激串序列。根据冲激函数的性质可知，仅当 时，P(k)为峰值。由于L个入射信号具有不同的角度-空间频率，每个角度-空间频率能且仅能在空间谱中引起一个峰值，因此空间谱中存在对应于L个入射信号的L个峰值。
[0108] 步骤八：根据所得空间谱进行波达方向估计。对步骤七中所构建的空间谱进行谱峰搜索操作，将其峰值按照由高到低排列，则最大的前L个峰值对应的角度即为L个入射信号的波达方向估计。
[0109] 本发明所提波达方向估计方法基于采用角度-空间频率描述的互质阵列二阶等价虚拟信号进行快速傅里叶变换，获得空间响应并基于此构建一个空间谱，通过对所构建空间谱的谱峰搜索过程得到波达方向估计结果。相比于传统的基于均匀线性阵列的波达方向估计方法，本发明所提方法在保证了波达方向估计自由度性能提升的同时，计算复杂度仅为 更好地满足了对估计实时性有较高要求的应用场景。进一步地，由于快速傅里叶变换的计算仅由复数加法和乘法组成，本发明所提波达方向估计方法相比于传统基于互质阵列的波达方向估计方法而言更容易在实际硬件系统中实现。
[0110] 下面结合仿真实例对本发明的效果做进一步的描述。
[0111] 采用互质阵列接收入射信号，其参数选取为M＝9，N＝10，即架构的互质阵列共包含2M+N-1＝27个物理阵元。信噪比为10dB，采样快拍数为T＝500，凯瑟窗参数β＝2。假定入射窄带信号数量为30，入射方向均匀分布于45°至135°这一空间角度域范围内，补零后快速傅里叶变换点数K＝2048。本发明所提方法得到的空间谱如图4所示，其中垂直虚线代表入射信号源的真实方向。可以看出，本发明所提方法能够有效分辨这30个入射信号源。对于传统采用均匀线性阵列的波达方向估计方法，利用27个物理阵元最多只能分辨26个入射信号，以上结果体现了本发明所提方法实现了自由度的增加。
[0112] 综上所述，本发明所提方法利用快速傅里叶变换提高了计算效率，在保证互质阵列波达方向估计自由度性能提升的同时，避免了现有方法优化问题的设计求解以及矩阵求逆、特征值分解等高复杂度的计算过程，更有利于在实际硬件系统上实现，也更好地满足了对估计实时性要求较高的应用场景。
[0113] 上述实施例用来解释说明本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明作出的任何修改和改变，都落入本发明的保护范围。
